Be Careful in Multimessenger Inference of the Hubble Constant: A Path Forward for Robust Inference

The bright sirens method for measuring the Hubble flow in the local Universe accessible to gravitational-wave (GW) detectors was proposed by B. F. Schutz (1986) based on the robust prediction by general relativity (GR) of the luminosity and time evolution of a binary neutron star (BNS) coalescence. It has been validated at the single-event level by the multimessenger observation of GW170817, which provided luminosity distance and redshift information independently (e.g., B. P. Abbott et al. 2017a) and produced a ∼15% measurement of the Hubble parameter H₀; see B. P. Abbott et al. (2017b). This measurement of the Hubble parameter is further improved to a ∼7% measurement by including the observation of a radio jet and improved peculiar velocity (see B. P. Abbott et al. 2017b; K. Hotokezaka et al. 2019; S. Mukherjee et al. 2021). Though such measurements are theoretically robust, in practice, multiple astrophysical uncertainties such as inclination angle (S. Nissanke et al. 2010), peculiar velocity (T. M. Davis et al. 2019; C. Nicolaou et al. 2020; S. Mukherjee et al. 2021; H. Nimonkar & S. Mukherjee 2023), waveform uncertainty (N. Kunert et al. 2024), nonstationary noise (S. Mozzon et al. 2022), etc. can obscure the inference of H₀ if not taken into account correctly. To accurately infer the value of the Hubble constant from GW data and shed light on the ongoing tension on the true value of the Hubble constant (see E. Abdalla et al. 2022 for a review), it is important to mitigate the influence of any possible systematic due to astrophysical and instrument uncertainties.

In order to mitigate these additional uncertainties, combining GW and electromagnetic (EM) observations can be beneficial. For a BNS or black hole–neutron star (BHNS) merger, the remnant, if matter is still present, is expected to source a highly relativistic, collimated matter outflow (jet) from its poles, which produces a short gamma-ray burst (GRB; ) and a long-lived transient spanning all wavelengths () called the afterglow (see E. Berger 2014; B. P. Abbott et al. 2017a; E. Nakar & T. Piran 2021). Besides allowing for the identification of a host galaxy, and through that a redshift, the EM counterpart can also mitigate a degeneracy between the luminosity distance (d_L ) and the system orientation in the GW analysis, which limits the precision of the d_L measurement (e.g., S. Nissanke et al. 2010; B. F. Schutz 2011).

The angle between the observer line of sight and the orbital angular momentum axis (see Figure 1) is generally referred to as the system inclination, θ_JN. Since the jet comes from a polar outflow of the remnant, θ_JN is expected to be correlated with the jet inclination. ⁵ While (non)observation of the prompt γ-ray emission can provide a weak (lower) upper bound on the viewing angle due to its highly collimated nature, the light curve of the afterglow, which is also detectable further off-axis, is more sensitive to the viewing angle θ_vi (which is defined as the angle between the jet's direction of propagation and the observer’s line of sight ; see Figure 1). In particular, the initial rising slope of the light curve (for bright signals) and the jet break contain orientation information and allow extraction of the viewing angle (e.g., E. Nakar & T. Piran 2021; G. Ryan et al. 2024). The constraint from the light curve is actually on a combination of the viewing angle and opening angle of the jet, and this degeneracy curtails the precision of the measurement (G. Ryan et al. 2020; K. Takahashi & K. Ioka 2020; E. Nakar & T. Piran 2021). A possible way to break the degeneracy in the light-curve data is through astrometry performed on the radio image of the afterglow, which can reveal a lateral apparent superluminal motion of the afterglow centroid that is highly sensitive to the viewing angle ⁶ and thus lead to a 10%-level measurement of the viewing angle; see, e.g., K. P. Mooley et al. (2018). Such a constraint for the viewing angle, combined with assumptions about its correlation with the system inclination, can limit or break the degeneracy between θ_JN and d_L and improve the precision of the d_L –H₀ measurement by a factor of 2 (S. Nissanke et al. 2010; H.-Y. Chen et al. 2019). However, by including EM data, the clear dependence of the inferred d_L on GR alone is lost, and one has to make sure that the additional precision at the single-event level is not accompanied by a reduced accuracy in the inference.

Zoom In Zoom Out Reset image size

In particular, the interpretation of the photometry and astrometry data for the jet afterglow and its connection to the properties of the merging compact binary system requires the assumption of a model that describes the emission from the jet afterglow. However, the emission model for the afterglow that is fit to the EM data to constrain θ_vi is strongly dependent on the assumed jet structure and geometry. Understanding the latter based on the properties of the merging binary, unfortunately, is presently not feasible due to the computational challenges involved in the physical processes at vastly different scales that govern the merger and subsequent remnant evolution with jet formation (L. Combi & D. M. Siegel 2023; E. R. Most & E. Quataert 2023; K. Kiuchi et al. 2024 for BNS and E. R. Most et al. 2021; K. Hayashi et al. 2022 for BHNS) and then jet propagation (R. Fernández et al. 2019; A. Kathirgamaraju et al. 2019) and nonthermal emission (e.g., J. Granot & R. Sari 2002). Thus, one has to assume a particular model for the jet structure to fit the light curve and centroid motion, which cannot be constrained by knowledge about the inspiraling binary from the GW data. We refer to the systematic uncertainty in the inference of H₀, if the inference of θ_vi is incorrect due to a theoretical modeling choice, as model bias.

Furthermore, while the GW and EM emission share the luminosity distance as a common parameter, it is necessary to specify a relation between the viewing angle for the jet and the system inclination for the inspiraling binary. It is often assumed that the jet axis is aligned with the total angular momentum of the postmerger system and that this in turn does not deviate much from the orbital angular momentum of the merging system—see, e.g., B. P. Abbott et al. (2017c), H.-Y. Chen et al. (2019), and A. Farah et al. (2020)—which allows for a direct identification of the two parameters when the posteriors on θ_JN and θ_vi are combined. However, there is presently no theoretical framework for confirming this assumption, and in general, a misalignment of between the three axes might be present due to remnant properties (e.g., N. Stone et al. 2013; K. Hayashi et al. 2022; Y. Li et al. 2023) or collimation in an anisotropic medium (e.g., H. Nagakura et al. 2014). We define the line-of-sight misalignment angle as θ_m = θ_JN − θ_vi ⁷ (see Figure 1) and refer to systematic deviations in H₀ due to an intrinsic mismatch between the jet angle and inclination angle () as misalignment bias.

The effect of model and misalignment biases is similar in the sense that both will shift the EM posterior with respect to the GW posterior and, consequently, the region of joint posterior support, as shown in Figure 2, and a failure to take this into account in the joint inference will bias the combined posterior for d_L . In particular, Figure 2 shows that a positive/negative misalignment angle, which is equivalent to an under-/overestimation of the viewing angle in the model bias scenario, will bias the posterior support toward larger/smaller d_L due to the anticorrelation of the GW posterior in the d_L –θ_JN plane.

Zoom In Zoom Out Reset image size

Previously, the presence of a model bias has been studied in great detail for the EM counterpart of GW170817 (K. P. Mooley et al. 2018; Z. Doctor 2020; E. Nakar & T. Piran 2021; G. Gianfagna et al. 2023, 2024; T. Govreen-Segal & E. Nakar 2023). While these authors study the systematic based on GW170817 and/or future similar events ⁸ with detailed EM modeling, the more general analysis of H.-Y. Chen (2020) uses Gaussian estimates for EM likelihoods. ⁹ Consequently, a systematic exploration of the model bias for systems that differ from GW170817 that employs realistic EM inference models and mock data has been missing so far.

This Letter points out that the misalignment between the jet direction inferred from EM observations and GW observations can be a potential hindrance in achieving accurate H₀ at the single-event level and the ultimate bottleneck when the statistical uncertainty is negligible. By performing joint cosmological inference on a suite of realistic mock data simulations, we show that due to a positive/negative value of the misalignment angle, there can be an under-/overestimation of H₀. We propose a Bayesian technique to mitigate this effect in a multimessenger analysis where light-curve, centroid motion of the jet, and GW data are used to infer the value of H₀.

To study the impact of potential model and misalignment biases on the joint GW+EM inference, we generate a series of mock observations in the GW and EM sectors that are consistent with an equal-mass BNS assuming different system inclinations and jet viewing angles in the range 0°–80°, which allows us to combine posteriors from both sources that either have an intrinsic misalignment between the jet and orbital angular momentum axes () or not.

We restrict the present analysis to BNS systems for simplicity. While BHNS systems may differ in the details of the jet launching mechanism and structure, if they launch a relativistic jet, we expect broadly similar phenomenology to a BNS: a GW inspiral with d_L –θ_vi degeneracy and a nonthermal EM synchrotron afterglow (see also M. Ruiz et al. 2019). While the asymmetries introduced by unequal masses or misaligned component spins have, in principle, the capacity to lead to a misaligned jet, a more complete theoretical understanding of the jet formation is required to leverage such correlations. Consequently, as a first step, we only consider equal-mass binaries in this analysis. We consider a flat Lambda cold dark matter cosmology (see H. Mo et al. 2010) with km s⁻¹ Mpc⁻¹, consistent with B. P. Abbott et al. (2017b), and we fix the redshift to z^(true) = 0.005, which corresponds to a luminosity distance of Mpc, since we will not focus on the peculiar velocity bias. This is an unrealistically small distance (for comparison, GW170817 was observed at around 40 Mpc; see B. P. Abbott et al. 2017b), and it will lead to superideal posteriors; however, we investigate the biases under such optimal conditions to clearly distinguish between statistical fluctuations (which should be strongly suppressed in the considered loud signals) from systematic deviations. We shall use a particular jet structure to generate all of the events; see Section 2.2.2 for details. Below, we outline the statistical framework that is used for the joint GW+EM inference of the Hubble parameter, as well as the methods with which we generated our mock observations in both sectors.

2.1. Bayesian Framework

From Bayes’s theorem (M. Bayes & M. Price 1763), we obtain the joint Hubble posterior through

where is the collection of all relevant data sets from the EM and GW sectors, respectively, and represents the prior/model assumptions under which the inference is performed. Assuming that the measurements in the subsectors are independent and that the data only depend on the event parameters θ, the data likelihood for the joint inference model is

where represents the prior assumptions for the analysis of the combined data. In this work, we are interested in combining posterior knowledge about the redshift z, the luminosity distance d_L , the inclination of the GW source θ_JN, and the viewing angle for the postmerger jet θ_vi to learn about H₀, and both observation channels have access to different but overlapping subsets of this parameter space. As discussed in Sections 2.2.1 and 2.2.2, the respective models for the jet afterglow and the GW emission have a much higher parameter space with complex correlations; however, we shall marginalize the full posteriors to the relevant subspace θ = θ_EM ∪ θ_GW, with and . In the choice of the prior for the joint inference, we assume the general scenario where the jet viewing angle and system inclination differ by a misalignment angle θ_m, whose value is determined by an a priori assumption in . Furthermore, there is a mismatch in the parameter range of θ_vi and θ_JN, since the reflection symmetry of the jet about the spin plane causes the scenarios θ_vi and to be indistinguishable. Thus, the GW inclination parameter will always be converted to a value compatible with θ_vi with (e.g., B. P. Abbott et al. 2017c). In the context of the inference, these are still treated as separate parameters, but we enforce this relation in the Hubble inference prior, which we shall take to be of the following form:

where the priors and are chosen to be Dirac delta distributions to account for the different parameters’ correlations. We have already introduced the misalignment angle θ_m as a free parameter in the prior for θ_vi to account for the possibility that we have a nonzero shift between θ_JN and θ_vi. In this work, we also assume that the redshift can be measured with high precision, ¹⁰ to the point that we can approximate the likelihood/posterior as a Dirac delta distribution. Inserting the prior, expressing the GW and EM data likelihoods in terms of posteriors, and integrating out all of the Dirac delta distributions result in the following reduced form for the combined data likelihood:

We choose to use a uniform prior for the Hubble parameter, km s⁻¹ Mpc⁻¹, where denotes the uniform distribution between two bounds. The priors for the λ_GW and λ_EM models that are used for the posteriors in the disjoint GW and EM analysis are specified in Section 2.2.1 and Appendix A.

2.2. Mock Data Generation

Based on the mock data generated with the procedures outlined in Sections 2.2.1 and 2.2.2 (see also Figure A1 for the mock light curves), we conduct several inference campaigns to probe the model and misalignment bias, respectively. For reference, we compute H₀ posteriors only using the GW data and the knowledge about the redshift, and the results will be labeled with GW+z. On the EM side, we consider four different scenarios, which either use a GJ (as in the injection) or a PLJ model for the fitting, and for each jet structure, we compute posteriors once only from the light-curve fit (denoted GJ/PLJ) and then also with the inclusion of centroid data (denoted GJc/PLJc). Once we combine the EM posteriors with the GW and redshift information, we will label the four different scenarios with GW+GJ+z and GW+PLJ+z for the light-curve-only fits and GW+GJc+z and GW+PLJc+z for the posteriors including the centroid motion. Since we are only working with system inclinations of <90°, the definitions of θ_JN and θ_vi agree with each other, and we shall quote all results in terms of the jet viewing angle θ_vi (and the misalignment angle θ_m if it is nonzero). Below, we describe the mismatch that can arise due to errors in modeling the jet and the intrinsic mismatch between the inclination angle and jet angle in Sections 3.1 and 3.2, respectively.

3.1. Aligned Runs and Afterglow Model Bias

For the analysis of the model bias, we combine GW and EM data without using a misalignment angle, i.e., . The median and 2.5th–97.5th percentiles for these posteriors, as well as relevant statistics, are tabulated in Tables B1 and B2 in Appendix B.

The individual reference runs without EM observations, GW+z, demonstrate a difference in H₀ inference at small viewing angles that decreases to as θ_JN increases past 40°. This can be attributed to the d_L –θ_JN anticorrelation, which is strongest at smaller viewing angles, and the noise in the mock data, which can easily shift the posterior peak along the d_L –θ_JN degeneracy. Since our mock GW data use a single noise realization, all GW+z runs inherit an identical push to larger θ_JN (smaller d_L ) and hence larger H₀. As the viewing angle increases past ≳50°, the two GW polarizations decouple in intensity, which breaks the d_L –θ_JN degeneracy and increases the H₀ precision to (H.-Y. Chen et al. 2019; E. Nakar & T. Piran 2021).

The effect of EM afterglow observations on H₀ inference depends on how well the EM sector constrains θ_vi relative to the GW data (see the results in Tables B1 and B2 in Appendix B). At small inclinations, θ_vi ≤ 20°, the EM data are bright and highly constraining. All runs show a median H₀ within ∼0.5σ of the true H₀ with σ/H₀ ≲ 3%, a substantial improvement over GW+z. The photometric-only EM runs (GW+GJ+z and GW+PLJ+z) perform worse than GW+z for θ_vi ≥ 30°, constraining H₀ to only ∼10% with differences as large as 1.5σ. Runs including centroid position observations (GW+GJc+z and GW+PLJc+z) remain constraining for θ_vi = 30°, 40°, with a similar difference () but improved precision (∼few percent). As θ_vi increases past 50°, even at our very close fiducial luminosity distance, the afterglow emission becomes too faint to observe, and the mock data sets become almost entirely upper limits; see Figure A1. At this point, the only added information from the EM posteriors can be an improved lower bound on the inclination, which translates into a lower bound on H₀.

Figure 3 (the θ_m = 0 row) shows the H₀ posteriors for θ_vi = 30°, 40°, including the less (and more) precise posterior on H₀ for the GJ (and GJc) case. A similar behavior is also identified in the analysis of the GW170817 afterglow in G. Gianfagna et al. (2024), where it is attributed to an attempt of the jet models to accommodate a late-time excess in the light curve. However, since we also find this for the fits of our mock observations that do not include such an effect, the behavior may be more generic, possibly due to parameter degeneracies altering the available prior volume (see G. Ryan et al. 2024 for discussion).

Zoom In Zoom Out Reset image size

To assess the presence of model bias, we compare the H₀ posterior distributions from fits run on the same data set with different assumed jet models. Table B2 contains the difference in median H₀ values from GJ and PLJ runs normalized by the GJ H₀ (ΔH₀) and by the width σ of the GJ H₀ posterior (δ H₀). In the no-centroid runs, GJ and PLJ models show a maximum bias of 2% in median H₀. This occurs in the θ_vi = 20° run and corresponds to a 0.49σ bias. All other runs show less than 0.08σ bias. Including astrometric centroid data increases the constraining power of the EM afterglow, increasing the precision of H₀ inference, which also increases its vulnerability to bias from model mis-specification. Between the GW+GJc+z and GW+PLJ+z runs, we again find a maximum discrepancy in median H₀ of 2%, this time in the θ_vi = 40° run.

3.2. Angular Mismatch Bias

We now turn our attention to systems that have a nonvanishing line-of-sight misalignment angle and compute joint GW+EM posteriors for the Hubble parameter neglecting this fact; i.e., we will still use θ_m = 0 in Equation (3). The mock data in the GW sector were set up with this in mind, and we computed several slightly shifted GW posteriors, scattered around the events that were used in Section 3.1. Consequently, will be fixed, and events with deviating are combined in the joint posterior, resulting in a systematic variation of . As discussed in Section 1, increasing implies that the GW posterior is shifted toward subsequently higher-inclination regions, and the overlap with the EM posterior will occur at larger values of d_L ; thus, we expect the H₀ posterior to shift toward smaller H₀ for greater misalignment . This behavior is demonstrated in Figure 3 for the posteriors from the GW+GJ(c)+z scenarios at and 40°. Two important observations can be made. Clearly, the high-precision posteriors from the GJc analysis () are more sensitive to the misalignment, and already a small shift of ±4° causes a 2σ–4σ bias in the resulting H₀ (see Table B4). Furthermore, we know from Section 3.1 that the posterior for the GJ analysis is biased toward larger H₀, which is evident in Figure 3, and we note that a large and positive will result in a more accurate Hubble posterior, since the two effects cancel. On the other hand, a small negative can amplify the model bias, leading to an even more inaccurate inference. We also confirm this behavior through the change of the 95% credible regions for the cases and 40° (see Figure 4), which seem to suggest the general trend that for an inference that is using light-curve and centroid data, will lead to a bias, which is then quickly increasing for larger magnitudes. ¹¹ For a viewing angle of 20°, we can also state that the inference using only photometry data is less sensitive since will not bias the result beyond ; however, we cannot make a similar result for the higher-inclination systems shown in Figures 3 and 4, since they have an intrinsic bias that can combine with the misalignment to cause more drastic shifts.

Zoom In Zoom Out Reset image size

3.3. Mitigation of Angular Mismatch Bias

Now that we have established that even a small line-of-sight angular mismatch of can lead to a severely biased Hubble posterior, we want to address a possible method of mitigating this bias. Since we have a well-defined mismatch parameter θ_m, we can consider this as a nuisance parameter in the joint inference and co-infer it with the Hubble parameter. Marginalization of the posterior in the H₀–θ_m plane should then lead to an unbiased posterior for the Hubble parameter, and we can also obtain a posterior for the mismatch. The key modification is that we treat θ_m as a variable with an associated prior added to Equation (3). The effect of the misalignment angle is to shift the d_L –θ_vi posterior/likelihood from the EM inference in the θ_vi direction without changing its shape. Consequently, in the joint posterior, we will compute a shifted overlap in θ_JN space in Equation (6). Note that contrary to the behavior of in Section 3.2, θ_m shifts the EM likelihood with respect to the GW+z likelihood; for θ_m > 0, it is shifted to larger θ_vi, leading to a posterior overlap with the GW posterior that occurs at larger inclination than for θ_m = 0 and therefore at smaller d_L , which will increase the inferred H₀. Thus, we expect a direct positive correlation between θ_m and H₀ (see also Figure 2 for a qualitative comparison).

We apply the modified inference to the same cases at , as shown in Figures 3 and 4, using the prior rad for the misalignment angle and considering only the GJ model with or without the centroid. The resulting 2D posteriors in the θ_m–H₀ plane for the cases of , ±4°, −10°, ±16° are shown in Figure 5. Considering the case of first, we observe the expected positive correlation/degeneracy, and both the GW+GJ+z and GW+GJc+z inferences are consistent with the truth to within 1σ, respectively.

Zoom In Zoom Out Reset image size

For the scenarios with intrinsic misalignment, it seems that negative misalignment affects the accuracy worse than positive misalignment, while all posteriors from the GJ data sets lead to a accuracy, and the large negative misalignments are at the edge. This behavior is more pronounced in the inference with the GJc data sets, where all scenarios with are only accurate to within 2σ, while positive misalignment still maintains a accuracy (see Figures B1 and B2 for statistics of the marginal H₀ posteriors). This is related to our choice of using the same noise seed for all GW mock data realizations, ¹² which for the considered data shifts the d_L –θ_vi posterior toward larger inclination angles than the injected value ¹³ ; see the discussion in Section 3.1. Consequently, the combined GW+EM posteriors tend to have higher posterior overlap for positive misalignment angles, thus leading to more posterior support toward θ_m > 0. However, all posteriors have support over a large region in the θ_m–H₀ plane, and once we compute the marginal posteriors for the Hubble parameter, we shall obtain unconstraining posteriors due to the degeneracy of θ_m and H₀. The statistics for the marginal H₀ posteriors in Figure 5 demonstrate that the precision, even for the most constraining EM data sets (GJc), is reduced to within a factor of 1.1–1.2 of the level of the GW+z posteriors (see second column in Figure B2).

Considering the most precise posteriors that are obtained from analyzing light-curve and centroid data, we obtain constraints for H₀ from marginalizing θ_m in the 2D posteriors in Figure 5, which is comparable to the constraints obtained by only fitting the light curve in the cases with zero misalignment in Section 3.1. We also compute 2D posteriors in the θ_m–H₀ plane for the series of events without misalignment that was considered in Section 3.1, where we would expect to obtain an accurate inference of H₀ with the maximum posterior region centered around θ_m = 0. This is indeed what we find; θ_m is constrained to 0 with precision, and we obtain Hubble posteriors that are accurate to within 1.5σ for all viewing angles . The posterior statistics are summarized in Table B5 and the second column of Figure B3 in Appendix B. However, also for these scenarios, the gain in precision from including EM data is reduced due to marginalization of the misalignment angle; see the second column of Figure B3 for more details.

Trusting the theoretical models, we might expect that the physical misalignment will be of at most; consequently, a narrower choice for the θ_m prior could be motivated, e.g., , and this could also reduce the spread in the posterior while still being able to confidently infer the misalignment in the cases with . We have applied the same bias-mitigating procedure as discussed above to data sets with , −42, 0°, 44, 62}, but we applied the narrower prior range for θ_m, and the statistics for the marginal H₀ posterior indeed show that for the posteriors from the GJc data sets, the accuracy behaves similar to the case with a wider θ_m prior, but the precision is again improved to a factor of 2 compared to the GW+z posteriors (Figure B4, as well as the third column in Figure B2). This demonstrates the need for narrow theoretical bounds for the misalignment angle in the joint inference of H₀, for the amount of precision that can be achieved in the single-event-level inference of H₀ in this scenario is directly linked to the prior width in H₀.

In this work, we analyze scenarios for the joint inference of the Hubble parameter from GW and EM data at the single-event level that attempt to mimic a realistic multimessenger campaign in the near future, in terms of available data and models, under ideal conditions. We point out a potential source of bias in the inference of the value of H₀ from multimessenger observation due to incorrect model fitting with EM observations and/or intrinsic misalignment between the inclination angle and the jet angle. For the case of incorrect modeling, we find that over the range θ_vi = 3°–40°, the model bias is a subdominant effect, as both assumptions on the jet structure lead to consistent results to within , extending the findings for a GW170817-like event in G. Gianfagna et al. (2024). GJc/PLJc are precise and accurate for θ_vi = 3°–40°, while GJ/PLJ are precise and accurate for θ_vi = 3°–20° (although already worse by a factor of 3 than the analysis with centroid data) but less precise () and inaccurate at θ_vi = 30°–40°, with an overestimation, which can be prohibitive for cosmological inference.

On the other hand, the physical misalignment between the orbital angular momentum axis of the inspiraling binary and the jet axis was identified as a significant source of bias in the joint inference of GW and EM data, particularly when centroid information is included. Data with line-of-sight misalignment that are analyzed under the assumption of alignment—which is frequently employed (see B. P. Abbott et al. 2017c; H.-Y. Chen et al. 2019; A. Farah et al. 2020), as the misalignment is not constrained by current theoretical models—can lead to biases in H₀, corresponding to 3%–7% deviations for and 5%–12% for . This shows that presently used models for the EM counterpart lead to posteriors that are even more sensitive than what was found by H.-Y. Chen (2020) based on Gaussian estimates. Since we are using superideal events, this is a lower limit, and it is already close to the required accuracy for resolving the Hubble tension at about ≳3σ (e.g., K. Hotokezaka et al. 2019). Ensuring an accurate H₀ inference from bright sirens would therefore require a method to mitigate such systematic biases.

Note that random values of misalignment for a large number of events would also cause random systematic errors in H₀, and if multiple events were combined, their effects might average out to provide an unbiased H₀. However, the prospects for observing many BNS or BHNS mergers that can be analyzed with the bright sirens method are limited, and the expectation is that analyses in the near future will have access to a few events at best; see, e.g., G. Gianfagna et al. (2024). In this case, even a few observations with precise but inaccurate H₀ posteriors from the inclusion of light-curve and centroid data could dominate over other bright siren measurements of H₀ that only have access to the redshift, and the bias of the multievent inference would be driven by the systematic errors in the analysis at the single-event level. Consequently, it is paramount to identify biases at the single-event level and to reduce these if possible.

To mitigate the uncertainty due to the model or misalignment biases, one is required to jointly infer the value of θ_m and H₀ to obtain an unbiased marginal posterior on H₀ and also an inference of θ_m . The multimessenger observations from GW and EM data are capable of providing information on θ_m , and our analysis shows that they should be inferred together, rather than assuming the value of θ_m to be 0. The later assumption can cause a biased H₀ even for a minor deviation that is easily feasible astrophysically (e.g., N. Stone et al. 2013; H. Nagakura et al. 2014; K. Hayashi et al. 2022; Y. Li et al. 2023).

We conclude that co-inferring and marginalizing over the misalignment angle can indeed reduce the bias in the inferred H₀ to within 2σ but at the cost of the precision that is available at the single-event level. If the misalignment angle is given a wide prior range, this can reduce the precision of the GW+EM inference to the level of the GW-only inference, thus removing the benefit of the combined analysis for the considered scenarios. However, this might change at larger redshifts, depending on the quality of the EM data compared to the GW data. Reducing the prior range for θ_m can restore the gain in precision for the combined GW+EM inference up to a factor of 2 with respect to the GW+z posteriors; thus, improvements on the theoretical bracket for this parameter will be paramount for ensuring a highly constraining and unbiased inference.

In the future, it will be important to incorporate this analysis in the bright siren studies for which EM measurements of the viewing angle from light-curve and/or centroid data are available. Along with the inference of H₀, this technique makes it possible to also infer the value of the mismatch angle θ_m . This can be useful in understanding the astrophysical properties of BNS/BHNS mergers better, shedding light on to what extent merger-driven jets are aligned with the orbital angular momentum. With a population of BNS/BHNS events, a population-level inference of θ_m can lead to new insights about the complexity of the postmerger phase.

The authors thank Samsuzzaman Afroz for useful comments on the draft. This research was supported in part by the Perimeter Institute for Theoretical Physics. Research at the Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Colleges and Universities. M.M. acknowledges the support of NSERC, funding reference No. RGPIN-2019-04684. The work of S.M. is part of the 〈data∣theory〉 Universe-Lab, supported by TIFR and the Department of Atomic Energy, Government of India. The authors express gratitude to the computer cluster of 〈data∣theory〉 Universe-Lab. The authors would like to thank the LIGO/Virgo/KAGRA scientific collaboration for providing the data. LIGO is funded by the U.S. National Science Foundation. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN), and the Dutch Nikhef, with contributions by Polish and Hungarian institutes. This material is based upon work supported by NSF’s LIGO Laboratory, which is a major facility fully funded by the National Science Foundation.

GRB afterglows are produced by synchrotron emission from a nonthermal electron population, accelerated at the forward shock driven by the GRB jet. The jet launching mechanism imbues the jet with a distribution of energy E(θ), typically taken to be axisymmetric. The jet sweeps up the ambient medium of density n₀, forms a shock, decelerates, and spreads laterally. As first-principles simulations of BNS mergers have not converged on a jet structure, we instead use simple phenomenological models: a GJ and a PLJ . Both have central energy E₀ and core width θ_c and are truncated at an outer angle θ_w ; the PLJ also includes an index b.

Following standard afterglow theory (e.g., J. Granot & R. Sari 2002), the nonthermal electrons acquire a power-law energy distribution of index −p, made up of a fraction ξ_N of the available electrons. The nonthermal electrons and the magnetic field are each assigned fractions _e and _B , respectively, of the available shock energy. The GRB is placed at redshift z and luminosity distance d_L and oriented θ_vi relative to Earth. For pure photometric data, this is enough to compute a flux; for astrometric centroid data, we also need initial sky coordinates RA₀ and Dec₀ and a PA on the sky.

We use the fiducial GJ model of G. Ryan et al. (2024) for the generation of light curves. This model is consistent with the GW170817 afterglow observations and has the parameters E₀ = 4.8 × 10⁵³ erg, θ_c = 32, θ_w = 224, n₀ = 2.4 × 10⁻³ cm⁻³, p = 2.13, _e = 1.9 × 10⁻³, _B = 5.8 × 10⁻⁴, and ξ_N = 0.95. Emission from this jet model is computed using afterglowpy v0.8.0 with jet spreading enabled; see G. Ryan et al. (2020, 2024) for implementation details.

We generate mock data from the fluxes of the injected GJ by simulating an idealized observation campaign using radio, optical, and X-ray photometric observations as well as radio observations of the centroid position. We designed the data-generation algorithm to roughly model an optimistic observation scenario for a GW EM counterpart.

The procedure goes as follows: for each photometric band, we space observations geometrically between 1 and 10⁴ days postmerger. Starting at t = 1 day, we compute the model flux, add noise, and determine whether our mock observatory would detect the emission. If the source is undetected, an upper limit is recorded, and the next observation is attempted a factor of 10^1/2 ≈ 3.1 later in time, giving these monitoring observations a cadence of two per decade in time. If at any time the source is detected, a flux observation with the relevant uncertainty is recorded, and the next observation is scheduled a factor of 10^1/4 ≈ 1.8 later in time, giving an observing cadence of four per decade in time. Once the afterglow is detected in a band, “observations” continue in that band at the four-per-decade cadence until two successive upper limits are found, at which point observations in that band cease. Radio very long baseline interferometry (VLBI) observations are carried out between 10 and 10⁴ days postmerger at the two-observations-per-decade cadence. This approach generates sufficiently dense observations to measure the initial rising slope of the light curve (if the emission is sufficiently bright) and covers a duration of time over which the jet break is expected to appear. These two features contain the most geometric information in the light curve and, if detectable, are essential for extracting the inclination from GRB afterglow observations.

Our mock radio, optical, and X-ray observatories emulate optimistic observations from the Very Large Array (VLA), the Hubble Space Telescope (HST), and Chandra, respectively, with uncertainties and noise levels similar to what was obtained with GW170817 (e.g., W. Fong et al. 2019; E. Troja et al. 2021). We characterize the radio and optical observatories with a flux limit (representing the noise level of the instrument) and a fractional source uncertainty δ (representing calibration and background-subtraction uncertainties). We assign a total uncertainty to each observation by simply adding these in quadrature: , where F_src is the model afterglow flux. The synthetic observed flux F_obs is F_src plus a random noise contribution drawn from a normal distribution of mean 0 and width σ. If , the observation is declared a nondetection, and we emit an upper limit of . Otherwise, the observation is a detection, and we emit F_obs with the uncertainty σ. In parameter estimation, these observations are assigned a Gaussian likelihood. Our radio photometry is generated at ν = 3 GHz with δ = 0.05 and μJy, representing a long integration by a VLA-like instrument. Our optical photometry is generated at ν = 5.1 × 10¹⁴ Hz with δ = 0.05 and nJy, representing a long integration by an HST-like instrument in the F606W filter.

Since many of these sources are faint, we generate mock X-ray data with Poisson statistics unless we exceed 100 observed counts. For each observation, we compute F_src at 5 keV and convert this to an expected number of counts S by multiplying by a flux-to-counts conversion factor of 1.4 × 10⁸ mJy⁻¹, a typical value for a 100 ks Chandra observation of GW170817 (E. Troja et al. 2021). We then draw a random background count B from a normal distribution of mean 1.0 and width 0.5, truncated to B > 0 (again, typical for GW170817). The observed photon count N is then drawn from a Poisson distribution with mean B + S. If N > 100, we replace this Poisson observation with a normally distributed flux observation of width 0.1 × F_src. In parameter estimation, X-ray observations are assigned a Poisson likelihood following the approach of G. Ryan et al. (2024) unless N > 100, in which case a Gaussian likelihood is used.

Each synthetic VLBI observation produces both a flux and sky coordinates (R.A. and decl.) with relevant uncertainties. The flux observation proceeds identically to the radio photometric routine. The uncertainty σ_c in each direction of the centroid is computed as , where the signal-to-noise ratio (SNR) is F_obs/σ. We compute the true R.A. and decl. from afterglowpy and then add Gaussian noise of width σ_c to produce the observed values. As in the mock photometry, if , the observation is declared a nondetection, a 3σ flux upper limit is emitted, and the R.A. and decl. uncertainties are set to large numbers to make the data points unconstraining. We take the flux uncertainty parameters as above and assign an optimistic σ_beam = 1 mas. In parameter estimation, the flux, R.A., and decl. are given a 3D multivariate normal likelihood with independent uncertainties (a diagonal covariance matrix).

For the fitting of the mock data, we employ the following prior ranges in , θ_c , θ_w , b, , p, , , , RA₀, Dec₀, PA}: , rad, rad, , , , , , , mas, mas, . Furthermore, the EM inference uses the same prior for the viewing angle as is used for the inclination in the λ_GW model, i.e., , but with range θ_vi ∈ [0, π/2], and for the distance prior, we use a Gaussian prior in H₀ that is centered on a fixed value of km s⁻¹ Mpc⁻¹, i.e., (in km s⁻¹ Mpc⁻¹), where denotes a Gaussian distribution with mean μ and variance σ. This prior then induces a prior on the luminosity distance, which can be obtained in closed form if we assume that we are at small redshifts where the luminosity distance is not sensitive to the specific cosmological model and , in which case , for constant redshift. This is chosen to limit the range in d_L beyond the Mpc prior that was specified above, since a large range in d_L reduces the efficiency of the inference. We made sure that the precision of the d_L posteriors from the EM data is always worse than the GW constraints, such that the joint inference is not affected by this choice.

Figure A1 shows the generated mock light-curve observations for this work.

Zoom In Zoom Out Reset image size

In this section, we summarize the results from the different inference campaigns. For every campaign, we quote the median () and the 2.5th–97.5th percentiles of the H₀ posteriors for the different inference scenarios. Furthermore, we provide statistics that help to estimate the accuracy and precision of the respective posteriors, we use to gauge the precision (statistical error) of the measurement, and, since we know the truth in our runs, we define to estimate the accuracy.

To better understand the differences/consistency between different inference models for the EM data, we also define the quantities for the relative deviation from the GJ(c) model and for the deviation from GJ(c) in terms of , in analogy to our accuracy statistic . We check consistency only among light-curve-only and light curve + centroid motion inferences, respectively, and we always use the GJ(c) jet models as a reference, since these should also be consistent with the injection.

We work with unrealistically close and therefore loud sources, which will suppress the statistical error (estimated in terms of ) to be able to identify systematic errors which might be hidden for sources with higher statistical errors (there the bias would only show up on the multievent level, once the statistical error reaches the bias scale). We choose the sources so close that the statistical error is below the bias scale that could still be problematic for the desired precision, any bias that we do not resolve at this point would also not become relevant for the immediate Hubble inference. This implies that and ΔH₀ do not provide representative properties of the posteriors that would generalize to more realistic sources. Only and δ H₀, which are normalized by the statistical error, would provide an idea of how the distributions might change for larger statistical errors.

In Tables B1 and B2, we show the statistics for the inference runs that do not use an intrinsic misalignment parameter and where the only inaccuracy can arise from a biased EM inference.

Table B1. Inferred Hubble Parameter for Different Inference Scenarios with Zero Misalignment Angle in the Datasets

	H0 (Median, 2.5th–97.5th percentiles)
	GW+z	GW+GJ+z	GW+PLJ+z	GW+GJc+z	GW+PLJc+z
3°
10°
20°
30°
40°
50°
60°
70°
80°

Note. This table lists the posterior results for the Hubble parameter from the inference campaign for the different inference scenarios and all viewing angles using the median of the posterior distribution and the 2.5th–97.5th percentiles.

Download table as:  ASCIITypeset image

Table B2. Precision and Accuracy of the Inferred Hubble Parameter for Different Inference Scenarios with Zero Misalignment Angle in the Datasets

	H0 (Median, 2.5th–97.5th percentiles)										Posterior Consistency
	GW+z		GW+GJ+z		GW+PLJ+z		GW+GJc+z		GW+PLJc+z		GJ/PLJ		GJc/PLJc
											ΔH0	δ H0	ΔH0	δ H0
3°	0.1	1.8	0.01	−0.51	0.01	−0.43	0.01	−0.53	0.01	−0.4	−0.0	−0.07	−0.0	−0.12
10°	0.1	1.69	0.01	−0.53	0.02	−0.36	0.01	0.03	0.01	0.03	−0.0	−0.07	−0.0	−0.01
20°	0.11	1.29	0.04	0.24	0.03	−0.26	0.01	0.13	0.01	−0.33	0.02	0.49	0.01	0.46
30°	0.12	0.84	0.09	1.39	0.08	1.45	0.02	−0.32	0.01	−0.65	0.01	0.08	0.0	0.3
40°	0.12	0.33	0.1	1.17	0.1	1.21	0.04	0.48	0.05	0.77	−0.0	−0.03	−0.02	−0.45
50°	0.12	0.15	0.1	0.57	0.1	0.55	0.1	0.5	0.11	0.51	0.0	0.02	−0.0	−0.03
60°	0.09	0.2	0.08	0.49	0.08	0.49	0.08	0.43	0.08	0.45	0.0	0.01	−0.0	−0.03
70°	0.05	0.33	0.04	0.62	0.04	0.63	0.04	0.6	0.05	0.61	−0.0	−0.01	−0.0	−0.02
80°	0.02	0.18	0.02	0.7	0.02	0.72	0.02	0.72	0.02	0.74	−0.0	−0.03	−0.0	−0.02

Note. This table lists statistics for the Hubble posterior from the inference campaign for the different inference scenarios and all viewing angles with two decimal places (0.0 means a number <0.01). We quote (for precision), (for accuracy), (for relative deviation the from GJ(c) model), and (for deviation from GJ(c) in terms of ).

Download table as:  ASCIITypeset image

In Tables B3 and B4, we show the statistics for the H₀ posteriors in Figure 4, which varies the intrinsic misalignment systematically around the aligned case for low-to-intermediate θ_vi.

Table B3. Inferred Hubble Parameter for Different Inference Scenarios with Nonzero Misalignment Angle in the Datasets

		H0 (Median, 2.5th–97.5th percentiles)
		GW+GJ+z	GW+PLJ+z	GW+GJc+z	GW+PLJc+z
20°	−6°
	−3°
	0°
	3°
	6°
	10°
	14°
30°	−16°
	−13°
	−10°
	−7°
	−4°
	0°
	4°
	7°
	10°
	14°
40°	−17°
	−14°
	−10°
	−6°
	−3°
	0°
	4°
	7°
	10°
	13°
	15°

Note. This table lists the median and 2.5th–97.5th percentiles of the posterior results for the Hubble parameter from the data sets at varying , shown in Figure 4.

Download table as:  ASCIITypeset image

Table B4. Precision and Accuracy of the Inferred Hubble Parameter for Different Inference Scenarios with Nonzero Misalignment Angle in the Datasets

		H0 Posterior Statistics								Posterior Consistency
		GW+GJ+z		GW+PLJ+z		GW+GJc+z		GW+PLJc+z		GJ/PLJ		GJc/PLJc
										ΔH0	δ H0	ΔH0	δ H0
20°	−6°	0.04	1.05	0.03	0.62	0.01	3.02	0.01	2.52	0.02	0.49	0.01	0.49
	−3°	0.04	0.69	0.03	0.22	0.01	1.73	0.01	1.25	0.02	0.49	0.01	0.47
	0°	0.04	0.24	0.03	−0.26	0.01	0.13	0.01	−0.33	0.02	0.49	0.01	0.46
	3°	0.04	−0.43	0.03	−0.94	0.01	−2.11	0.01	−2.53	0.02	0.48	0.0	0.43
	6°	0.03	−0.96	0.04	−1.39	0.01	−3.93	0.01	−4.33	0.01	0.42	0.0	0.42
	10°	0.03	−2.03	0.04	−2.26	0.01	−6.8	0.01	−7.16	0.01	0.37	0.0	0.33
	14°	0.03	−3.29	0.04	−3.4	0.01	−11.18	0.01	−11.53	0.01	0.38	0.0	0.34
30°	−16°	0.08	2.53	0.07	2.68	0.02	6.67	0.02	6.97	0.01	0.11	0.01	0.41
	−13°	0.08	2.28	0.07	2.41	0.02	5.95	0.02	6.13	0.01	0.11	0.01	0.4
	−10°	0.08	2.1	0.07	2.21	0.02	5.06	0.02	5.1	0.01	0.12	0.01	0.39
	−7°	0.08	1.89	0.07	2.0	0.02	3.64	0.02	3.51	0.01	0.11	0.01	0.36
	−4°	0.08	1.68	0.08	1.77	0.02	2.4	0.02	2.17	0.01	0.1	0.01	0.35
	0°	0.09	1.39	0.08	1.45	0.02	−0.32	0.01	−0.65	0.01	0.08	0.0	0.3
	4°	0.09	0.97	0.08	0.95	0.01	−4.08	0.01	−4.47	0.01	0.1	0.0	0.29
	7°	0.09	0.85	0.09	0.75	0.01	−6.18	0.01	−6.53	0.01	0.14	0.0	0.22
	10°	0.09	0.52	0.09	0.42	0.01	−9.31	0.01	−9.7	0.01	0.13	0.0	0.16
	14°	0.1	0.2	0.09	0.08	0.01	−14.51	0.01	−14.62	0.01	0.12	0.0	0.01
40°	−17°	0.09	2.23	0.09	2.33	0.04	5.67	0.04	6.26	−0.01	−0.08	−0.01	−0.26
	−14°	0.09	2.05	0.09	2.16	0.05	4.82	0.04	5.4	−0.01	−0.08	−0.01	−0.3
	−10°	0.09	1.82	0.09	1.84	0.05	3.62	0.05	3.92	−0.0	−0.02	−0.02	−0.35
	−6°	0.1	1.38	0.1	1.45	0.05	2.36	0.05	2.58	−0.0	−0.0	−0.02	−0.42
	−3°	0.1	1.36	0.1	1.41	0.05	1.67	0.05	1.86	−0.0	−0.01	−0.02	−0.45
	0°	0.1	1.17	0.1	1.21	0.04	0.48	0.05	0.77	−0.0	−0.03	−0.02	−0.45
	4°	0.1	1.01	0.1	1.01	0.04	−1.36	0.06	−0.57	0.0	0.01	−0.02	−0.53
	7°	0.1	0.86	0.1	0.87	0.04	−2.95	0.07	−1.39	−0.0	−0.02	−0.02	−0.52
	10°	0.09	0.73	0.09	0.75	0.04	−4.09	0.09	−1.6	−0.0	−0.03	−0.02	−0.54
	13°	0.09	0.67	0.09	0.74	0.05	−4.88	0.11	−1.89	−0.01	−0.07	−0.03	−0.56
	15°	0.09	0.57	0.08	0.7	0.06	−4.59	0.12	−1.96	−0.01	−0.11	−0.04	−0.69

Note. This table lists statistics for the Hubble posterior from the data sets at varying , shown in Figure 4, for the different inference campaigns and all viewing angles with two decimal places (0.0 means a number <0.01). We quote (for precision), (for accuracy), (for relative deviation from the GJ(c) model), and (for deviation from GJ(c) in terms of ).

Download table as:  ASCIITypeset images: 1 2

Figure B1 and the second column in Figure B2 show the statistics for the marginal H₀ posteriors from the bias-mitigating inference of H₀ and θ_m with a subsequent marginalization over θ_m (we shall hereafter refer to this as inference). In particular, we consider scenarios where we combine with varying values of , which are also shown in Figure 5. This complements the previous inferences with , where we assumed , and obtained the biased posteriors in Figures 3 and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table B5 and Figure B3 show the statistics corresponding to the marginalized H₀ posteriors from the bias-mitigating inference of all aligned () data sets. Based on these, one can judge if this method can recover zero misalignment angle for the aligned GW and EM data sets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table B5. Inferred Hubble Parameter for Different Inference Scenarios with Zero Misalignment Angle in the Datasets, but using Co-inference of the Misalignment Angle and H₀

	H0 (Median, 2.5th–97.5th Percentiles)
	GW+GJ+z	GW+PLJ+z	GW+GJc+z	GW+PLJc+z
3°
10°
20°
30°
40°
50°
60°
70°
80°

Note. This table lists the median and 2.5th–97.5th percentiles of the marginal posterior results for the Hubble parameter from the data sets in the case of the co-inference of (θ_m, H₀), as discussed in Section 3.3.

Download table as:  ASCIITypeset image

The presentations in Figures B2 and B3 should allow for an easy extraction of and and therefore and . For comparison with the uniformed prior for the misalignment angle, in Figure B4 and the third column of Figure B2, we present statistics for the same analysis as in Figure B1; when we employ a narrower prior on the value of θ_m, in particular, we use .