Investigating the detection rates and inference of gravitational-wave and radio emission from black hole neutron star mergers

Black hole neutron star (BHNS) mergers have recently been detected through their gravitational-wave (GW) emission. BHNS mergers could also produce electromagnetic (EM) emission as a short gamma-ray burst (sGRB), and/or an sGRB afterglow upon interaction with the circummerger medium. Here, we make predictions for the expected detection rates with the Square Kilometre Array Phase 1 (SKA1) of sGRB radio afterglows associated with BHNS mergers. We also investigate the benefits of a multimessenger analysis in inferring the properties of the merging binary. We simulate a population of BHNS mergers and estimate their sGRB afterglow flux to obtain the detection rates with SKA1. We investigate how this rate depends on the GW detector sensitivity, the primary black hole (BH) spin, and the neutron star equation of state. We then perform a multimessenger Bayesian inference study on a fiducial BHNS merger. We simulate its sGRB afterglow and GW emission and take systematic errors into account. The expected rates of a combined GW and radio detection with the current generation GW detectors are likely low. Due to the much increased sensitivity of future GW detectors like the Einstein Telescope, the chances of an sGRB localisation and radio detection increase substantially. The unknown distribution of the BH spin has a big influence on the detection rates, however, and it is a large source of uncertainty. Furthermore, for our fiducial BHNS merger we are able to infer both the binary source parameters as well as the parameters of the sGRB afterglow simultaneously when combining the GW and radio data. The radio data provides useful extra information on the binary parameters such as the mass ratio but this is limited by the systematic errors involved. A better understanding of the systematics will further increase the amount of information on the binary parameters that can be extracted from this radio data.


Introduction
The detection of the first two black hole neutron star (BHNS) mergers, GW200105 and GW200115, by the advanced Laser Interferometer Gravitational-Wave Observatory (aLIGO) and Virgo detector network (Aasi et al. 2015;Acernese et al. 2014) completed the trifecta of compact binary coalescence gravitationalwave (GW) observations . The third GW catalogue (GWTC-3) by the GW detector network contains a total of 90 significant detections from binary black hole (BBH), BHNS, and binary neutron star (BNS) mergers (The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration et al. Of particular interest to astronomers are the GW observations of BNS and BHNS systems because of the possibility of a complementary electromagnetic (EM) signature . This possibility was first confirmed with the BNS source GW170817 (LIGO Scientific Collaboration and Virgo Collaboration et al. 2017), which was accompanied by EM radiation from the dynamical ejecta of a kilonova (see e.g. Chornock et al. 2017;Coulter et al. 2017), the short gamma-ray burst (sGRB) GRB170817A (e.g. Abbott et al. 2017), and the afterglow of the jet interacting with the interstellar environment (e.g. Alexander et al. 2017;Haggard et al. 2017;Hallinan et al. 2017). No such emission was detected from either BHNS mergers GW200105 or GW200115 despite a multitude of follow-up campaigns being performed after their detections (e.g. Antier et al. 2020;Abbasi et al. 2021; Abe et al.  Foucart et al. (2018), showing the dependence on the BH spin ( BH ), mass ratio ( ), and NS tidal deformability (Λ NS ). Best-fit values for the free parameters are assumed. In each of the three panels, BH , , or Λ NS is kept constant, as shown in the title of the panel, and the other two parameters are varied. The coloured regions indicate the amount of remnant mass ( rem ) in solar masses. Paterson et al. 2021;Ridnaia et al. 2020;Kasliwal et al. 2020;Ashkar et al. 2021;Anand et al. 2021) .
The probability of detecting any EM emission of a BHNS merger is tightly connected to the ejecta that might be produced during and after the merger. In BHNS mergers, such ejecta originates from the tidal disruption of the neutron star (NS) if the NS does not directly plunge into the black hole (BH). A small part of the tidally disrupted NS, which is not accreted onto the BH merger remnant, may become bound to a disk around the merger remnant or be ejected as unbounded material (Foucart et al. 2018). The rotational energy of the BH merger remnant in conjunction with the magnetic field of the accretion disk could power an sGRB (e.g. Paschalidis et al. 2015), which would produce a radio afterglow upon interacting with the circum-merger medium (e.g. Metzger & Berger 2012). It has been suggested that some of the sGRBs detected to date originate from BHNS mergers (Gompertz et al. 2020). Other sources of radio emission, such as the afterglow of the unbound dynamical ejecta, have also been proposed (e.g. Nakar & Piran 2011).
The fraction of systems in which tidal disruption occurs is still uncertain as it depends on unknown distributions of the binary mass ratio, NS compactness, and BH spin. Still, this fraction is assumed to be small (Zappa et al. 2019;Fragione 2021). For example, while BHNS mergers with large BH spins aligned with the system angular momentum produce the most ejecta mass, high spins do not seem to be consistent with current GW detections (Abbott et al. 2021b). Furthermore, a hard NS equation of state (EOS) also increases the amount of ejecta, but this is again disfavoured by current GW detections (The LIGO Scientific Collaboration and The Virgo Collaboration et al. 2018;Fragione 2021).
While the rate of BHNS mergers with any EM emission is likely low, this is offset by the high potential science impact of a detection. An EM observation associated with a BHNS merger could, for example, provide information on the NS EOS (Pannarale & Ohme 2014; Ascenzi et al. 2019;Fragione & Loeb 2021), constrain the BH spin (Barbieri et al. 2019), or help determine the progenitor class if it is unclear from the GW detection alone (Hinderer et al. 2019). Such types of multi-messenger analyses have been performed extensively on GW170817 (e.g. Radice See also GCN archive for S200105ae. 2020; GCN archive for S200115j. 2020Coughlin et al. 2018;Radice & Dai 2019;Coughlin et al. 2019;Raaĳmakers et al. 2020;Capano et al. 2020;Dietrich et al. 2020;Breschi et al. 2021b), and various general Bayesian multi-messenger frameworks have been developed with this goal in mind (e.g. Breschi et al. 2021a;Raaĳmakers et al. 2021a;Nicholl et al. 2021).
In this paper we investigate both the expected rates and the parameter inference of radio observations of sGRB afterglows associated with BHNS mergers. We build on recent population synthesis results (Broekgaarden et al. 2021) and connect them to an analytical estimate of the sGRB afterglow flux (Nakar et al. 2002;Duque et al. 2019). This allows us to make predictions for the sGRB afterglow radio detection rates with the Square Kilometre Array Phase 1 (SKA1; Braun et al. 2019). We then perform a comprehensive joint analysis of simulated GW and radio data in a Bayesian framework, taking systematic errors into account. Here, we make use of a recent model of the full afterglow light curve (Ryan et al. 2020) and explore the benefits of a multi-messenger analysis in the inference of the source properties. Notably, we do the parameter inference of the GW data and the radio data simultaneously instead of using a sequential approach (see e.g. Barbieri et al. 2019).
The rest of the paper is structured as follows. In Sect. 2 we give an overview of the connection between the binary properties and the ejecta mass and describe the assumed jet launch mechanism and jet energy. We derive the expected radio afterglow detection rates in Sect. 3. We then describe our Bayesian multi-messenger framework setup in Sect. 4 and show the results for a fiducial BHNS merger in Sect. 5. We discuss our findings in Sect. 6 and end the paper with a summary and conclusion in Sect. 7.

Computation of the ejecta mass
The amount of bound disk and unbound dynamical ejecta in BHNS mergers is heavily dependent on the binary source parameters of the mass ratio = BH / NS , with BH the mass of the BH and NS the mass of the NS, the tidal deformability of the NS Λ NS , and the dimensionless BH spin BH . In general, comparable mass binaries with high BH and a stiff assumed EOS (high Λ NS ) produce the most ejecta (Foucart et al. 2018). The analyses in this paper are limited to non-precessing systems so BH and the dimensionless NS spin NS refer to the component of the spin parallel to the orbital angular momentum.
To fully capture the relevant physics and predict the ejecta properties, numerical relativity (NR) simulations are necessary, but they bring a large computational burden. As an approximation, we used analytical formulae from the literature which are fits to such simulations, covering a range of binary source parameters. These fits incorporate various free parameters, which can change depending on the precise form of the fitting formula and the types of NR simulations used to constrain the free parameters. To avoid any ambiguity on how we implemented these fits, we list the formulae together with the best-fit values of the free parameters below. Foucart et al. (2018) provide a formula, fit to 75 NR simulations, for calculating the remnant mass normalised to the baryonic mass of the NS ( NS ): with best-fit values = 0.406, = 0.139, = 0.255, and = 1.761. Here, = /(1 + ) 2 andˆI SCO = ISCO / BH is the innermost stable circular orbit radius normalised by the BH mass: To compute the compactness of the NS NS , we used the approximately universal C-Love relation for NSs of Eq. 78 in Yagi & Yunes (2017): with best-fit coefficients 0 = 0.360, 1 = −0.0355 and 2 = 0.000705. While this relation is not a perfect substitute for integrating a NS EOS to obtain NS , it performs well across a wide range of EOSs and has much less computational overhead. From NS we could, again approximately, compute NS (Lattimer & Prakash 2001): The dependence of rem on BH , , and Λ NS is visualised in Fig. 1. Given rem , we can calculate the amount of disk mass ( disk ) if we have an estimate of the amount of unbound dynamical ejecta mass ( dyn ): We again turned to fits of NR simulations (Krüger & Foucart 2020) for dyn : with best-fit coefficients 1 = 0.007116, 2 = 0.001436, 4 = −0.02762, 1 = 0.8636, and 2 = 1.6840. Parts of the disk itself can also become unbounded through disk wind outflows that are either thermally or magnetically driven. Raaĳmakers et al. (2021a) derive a simple formula that broadly captures the dependence found in NR simulations of the disk wind ejecta mass ej,disk on (Fernández et al. 2020): Here we assumed average values for the free parameters of 1 = 0.18 and 2 = 0.29.

GRB jet
GRB170817A firmly established BNS mergers as a source of sGRBs (e.g. Abbott et al. 2017). Whether BHNS mergers are also able to produce such ultra-relativistic highly collimated outflows, or jets, is not yet fully understood (see Kyutoku et al. 2021 for a recent review of BHNS mergers and the expected EM emission). If the NS is tidally disrupted, more mass may be accreted in BHNS mergers than in BNS mergers, which could increase the energy budget of the jet (see Gompertz et al. 2020, and references therein). Furthermore, the dynamical part of the ejecta is only present close to the equatorial plane in BHNS mergers (Kyutoku et al. 2013) so it cannot choke a possible jet. Some gas pressure from surrounding ejecta might be necessary, however, to create a highly collimated, and thus jet, outflow (see e.g. Nagakura et al. 2014). Unlike the dynamical ejecta, the disk wind ejecta are present in the polar regions. It is not clear whether an ultra-relativistic outflow can overcome such ejecta with sufficient collimation. Just et al. (2016) still find that BHNS mergers are able to harbour ultra-relativistic jets because of the lack of polar dynamical ejecta. A further topic of debate is the mechanism through which the jet is launched. Two often proposed candidates are the Blandford-Znajek (BZ) mechanism (Blandford & Znajek 1977) and neutrino pair annihilation (Eichler et al. 1989;Mészáros & Rees 1992;Just et al. 2016). Salafia & Giacomazzo (2021, hereafter referred to as SG21) discuss both mechanisms in detail and derive disk accretion-to-jet energy conversion efficiencies. For GW170817, they calculate that both mechanisms have efficiencies that are consistent with GRB170817, and they can therefore not distinguish between the two. Kyutoku et al. (2021) argue that neutrino pair annihilation does not occur on a sufficiently long timescale to explain the observed sGRB duration. Here we thus assumed that a jet gets launched through the BZ mechanism and followed SG21 for the accretion-to-jet energy conversion efficiency (see Barbieri et al. 2019 for a similar derivation). For a total accreted disk mass acc = disk − ej,disk , we took acc ≥ 0.03 as a necessary condition to create an sGRB of ∼ 1s duration (Stone et al. 2013;Pannarale & Ohme 2014).
Taking a fairly typical gamma-ray burst prompt emission efficiency of 10% (Beniamini et al. 2016), the final kinetic energy of the jet responsible for the relativistic shock producing the afterglow becomes: where the factor one-half accounts for an identical counterjet.  where f BH is the spin of the final remnant BH and Ω H is the dimensionless angular frequency at the BH horizon: To capture the dependence of f BH on the binary source parameters, we used the fits to NR simulations in Zappa et al. (2019), which are extensions of fits to BBH mergers (Jiménez-Forteza et al. 2017). Compared to the formulae stated in SG21, in the regime f BH ≤ 0.25 an additional factor of 10 −2 is necessary (Salafia & Giacomazzo 2022).
Analogous to Fig. 1, the dependence of on BH , , and Λ NS is visualised in Fig. 2. Considerable parts of the parameter space do not produce any disk mass or not enough disk mass to satisfy the imposed threshold for sGRB creation. Similar to rem , comparable mass binaries with high Λ NS produce the highest jet energies. There is a strong dependence of on BH as well. This is partly inherited from rem but further strengthened by the dependence of BZ on BH (Zappa et al. 2019). The dependence on BH has a big influence on the expected sGRB afterglow detection rates as we show in Sect. 3.
To compute the afterglow emission of the jet when it shocks the interstellar medium, we took two approaches. For our parameter inference, full light curve models are required, which we detail in Sect. 4.3. In the next section we look at the detection rates for the jet afterglow of BHNS mergers. Here we can limit ourselves to an analytical approximation for the peak flux.

Afterglow detection rates with SKA1
In this section we derive quantitative estimates for the detection rate of sGRB afterglows associated with GW events of BHNS mergers. We rely on a recent population synthesis study done in Broekgaarden et al. (2021, hereafter referred to as B21). A detailed study on the influence of various input parameters such as the star formation rate density is beyond the scope of this work. Instead, we use their fiducial model 'A000' (Sect. 3 of B21) to obtain a baseline estimate of the detection rate given the GW detector type, the BH spin ( BH ), and the NS EOS.

Population synthesis of BHNS mergers
B21 provide both the final properties of the BHNS mergers after star formation and the merger rate per redshift ( ). We followed their methods, which are derived from Neĳssel et al. (2019), and integrated over 250 redshift shells from = 0 to = 0.5 to get the total amount of BHNS mergers per year in this volume of space-time. In each redshift shell, we created BHNS mergers according to the merger rate (which is a function of e.g. the assumed star formation rate density) with various BH and NS given by the mass distribution for model A000. We distributed the mergers uniformly in sky position and orientation in order to calculate the GW signal-to-noise ratio (S/N). For a given GW detector sensitivity and S/N threshold, we then obtained the rate of detected GWs from BHNS mergers.
To calculate the GW S/N, we proceeded similarly to Barrett et al. (2018) and B21. Barrett et al. (2018) computed an interpolated grid of S/Ns for different sets of component masses calculated with GW waveforms IMRPhenomPv2 (Hannam et al. 2014;Khan et al. 2016) and SEOBNRv3 (Pan et al. 2014;Babak et al. 2017). They marginalised over the external parameters (Finn & Chernoff 1993) to compute an averaged detection probability instead of a single S/N value. We omitted the marginalisation and calculated the S/N value instead because we needed the inclination as input for the magnitude of the afterglow flux. We have disregarded the effects of spin and NS tidal disruption on the GW waveform here but take them into account in Sect. 4. We looked at both a network of second-generation (2G) GW detectors at design sensitivity and a future network of third-generation (3G) detectors including the Einstein Telescope (ET) with the ET-D sensitivity curve (Hild et al. 2011). For both networks, we calculated the S/N in a single detector, either aLIGO or ET, and took a detection threshold S/N thresh ≥ 8 as a standard proxy for a detection, and sufficient localisation for follow-up, with the entire network.
For all BHNS mergers we created, we computed the amount, if any, of disk mass and the resulting kinetic energy of the gammaray-burst jet using Eqs. 1 through 12. For the ejecta formulae, we also needed to specify an EOS and a BH spin. We set the EOS by fixing the radius of the NS, in line with B21, to either NS = 11.5 km (consistent with GW observations; see e.g, The GW + SKA1 detections per year aLIGO, RNS = 11.5km aLIGO, RNS = 13.0km ET, RNS = 11.5km ET, RNS = 13.0km Fig. 3: Combined detection rate of BHNS mergers in our simulated population observed through their GW emission and their associated sGRB afterglow emission. The radio emission is observed with SKA1-Mid at nominal frequency and sensitivity. The GW emission is observed with either an aLIGO detector network (blue and green circles) or an ET detector network (orange and red circles). Either a soft NS EOS ( NS = 11.5 km) or a hard NS EOS ( NS = 13.0 km) is assumed. The detection rate is shown as a function of the BH spin ( BH ). The dashed black lines are not computed but connect the points to guide the eye.
use to obtain the peak flux magnitude of the possible sGRB afterglow.

sGRB afterglow
Similar to previous work (Boersma et al. 2021), we followed Duque et al. (2019) for an estimate of the sGRB afterglow peak flux based on the theory in Nakar et al. (2002): where is the opening angle of the jet core, 0 is the circumburst density, e , B and are shock microphysics parameters, is the observing frequency, is the luminosity distance and obs is the observing angle. We assumed the jet is produced orthogonal to the orbital plane so that obs is equal to the inclination angle of the binary merger used in the GW analysis.
As shown in Eq. 13, the observed flux does not directly depend on but on the on-axis isotropic-equivalent jet energy 0 . We used a standard formula for top-hat jets to convert between the two: 0 = /(1 − cos( )), where the opening angle is fixed to a representative value of = 0.1 rad (≈ 5.7 deg) (Beniamini & van der Horst 2017). The literature on gamma-ray-burst afterglow modelling often uses such a simplified top-hat jet approximation instead of a structured jet model (see e.g. Aksulu et al. 2022 for a recent study). The top-hat jet, in contrast to the structured jet models, assumes no angular dependence of the jet energy, which is a good approximation at small inclination angles. This belief breaks down for afterglows associated with GW detections as these will most likely be viewed off-axis. This was first demonstrated by GW170817 (see e.g. Ryan et al. 2020 and references therein) and is expected to largely hold true for future detections as well, as we show in the next section. While the difference in the full afterglow light curve is substantial between top-hat and structured jets, the discrepancy for the peak flux is not as large (Duque et al. 2019). Because we focus on detectability only as a function of the peak flux in this section, we stick to Eq. 13 and employ a structured jet model later in Sect. 4. Similar to Hotokezaka et al. (2019), we fixed = 2.2 and = 0.1 to fiducial values. These values are consistent with the observed gamma-ray-burst population (Beniamini & van der Horst 2017;Aksulu et al. 2022). Both 0 and have a much broader population distributions than or so we did not fix these parameters. We took the same approach as Duque et al. (2019) and considered a log-normal distribution for both parameters with mean = 10 −3 and standard deviation = 0.75. Furthermore, was constrained to the range [10 −4 , 10 −2 ]. We assumed a flat Λ cold dark matter cosmology with parameters given by the WMAP9 dataset (Hinshaw et al. 2013). We observed the afterglows with the SKA1-Mid telescope of SKA, at a nominal frequency of = 1.43 GHz, nominal rms continuum noise of RMS = 2 Jy (Braun et al. 2019) and claimed a detection when ,1.43 passes the 5 RMS threshold. We have used 20 combinations of detector type, NS and BH in total. We ran each combination for 1000 realisations of one year. We show the averaged detection rates in Fig. 3.

Rates
For the aLIGO network, we obtain the same averaged 11 detected GW mergers per year as B21 verifying our methods. We also calculated the percentage of those GW detections with some form of NS tidal disruption. For BH = 0, we again recover similar percentages as B21 of 1.3% and 3.9% with NS = 11.5 km and NS = 13.0 km, respectively. We provide some additional figures of the simulated populations in Appendix A. We focus on Fig. 3 and Fig. 4 below.
For a network of 2G detectors, we find that the SKA1-Mid detection rates of sGRB afterglows associated with BHNS GW events are likely low. Only for a population of BHNSs with prob-Article number, page 5 of 16 A&A proofs manuscript no. 43267corr ably unrealistically high BH ∼ 0.8 (Fragione 2021), the combined detection rate per year nears unity. Transitioning from a soft ( NS = 11.5 km) to a hard ( NS = 13.0 km) NS EOS does increase the rates two-to three-fold but a combined detection remains rare for all but the highest BH .
For the soft EOS, we do not observe any combined events with the 2G network in our 1000 realisations for BH = {0, 0.2}. Similarly for the soft EOS with the 3G network, we do not observe any combined events for BH = 0. In those cases, that would imply an upper limit, assuming a Poisson distribution, at a 95% confidence level of less than one combined detection per ∼ 300 years.
The situation changes considerably for a network of 3G detectors and non-zero BH spin. Now, one combined detection per year is not unlikely for BH ∼ 0.2, assuming a hard NS EOS. Many combined detections are expected per year for BH ∼ 0.4−0.8 but we stress this is most certainly an overestimate of the detection rates.
The stark differences between the rates for the 2G network and the 3G network are a result of the much increased sensitivity of the latter. We find, like Dobie et al. (2021) (see Hotokezaka et al. 2016 for an earlier study on the radio detectability of GW mergers), that SKA1-Mid is able to observe afterglows out to gigaparsec distances, which is beyond the range of aLIGO. Our 3G network, however, is able to observe most of the mergers occurring in our chosen space-time volume, which enables many more combined detections with SKA1-Mid. This is visualised in Fig. 4, where the combined GW and radio detection rates are plotted as a function of the distance and inclination angle. This is done for the ET and BH = {0.4, 0.8} combinations and for the aLIGO and BH = 0.8 combination. The hard NS EOS is used and all rates are normalised to their respective total rate, given in Fig. 3, for ease of comparison. Clearly most of the detected sources with the 3G network and SKA1-Mid are located at much greater distances than those detected with the 2G network. The inclination angle distribution of the observed sources also varies with GW detector type and BH spin. In general, because of the strong dependence of the peak flux on the inclination angle, most sources are detected at relatively small . Nearby sources (i.e. those detected with aLIGO) with comparatively higher fluxes can still be observed at larger 20 deg as well. Similarly, sources with high BH ∼ 0.8 and subsequently larger jet energies are also more readily detected at larger . For the ET and BH = 0.4 combination, we predict about 30% of the afterglows to be observed on-axis ( ≤ ).
A 3G network is expected to detect mergers at significantly larger redshifts than our maximum of = 0.5 (Maggiore et al. 2020). We also find that ∼ 25 − 35% (depending on the chosen BH ) of our detected afterglows in the radio are between 0.4 < < 0.5. This implies SKA1-Mid will be sensitive enough to detect afterglows with > 0.5 as well. While these relatively distant mergers are not taken into account here, some of these will thus be observable both through GWs and their afterglows. Such mergers will increase the rates for the 3G network and SKA1-Mid combination further but we note that a sufficient localisation for follow-up is not guaranteed at larger distances (Maggiore et al. 2020).
The combined detection rates are strongly dependent on the assumed BH and an approximate log-linear relation is visible in Fig. 3. We attribute this partly to the formulation of the BZ accretion-to-jet energy conversion efficiency in Eq. 11. 0 is directly proportional to this efficiency and by extension ,1.43 as well. For our set of BH , we are usually in the regime where f BH > 0.505 (Zappa et al. 2019). As the final BH spin f BH is approximately linear in BH (Zappa et al. 2019), the log-linear dependence on BH follows naturally from Eq. 11 in this regime. The dependence of acc on BH is less intuitively understood, but we observe a roughly log-linear slope there too.
In conclusion, the chances of finding any sGRB afterglow of a BHNS merger in the near future are slim. Fragione (2021) even argue that it is unlikely to detect most forms of EM emission associated with BHNS mergers when looking at the small fraction of mergers with ejecta mass. We remain cautiously optimistic on the basis of our simulations that the low probability of ejecta mass might be offset by the sheer amount of mergers a 3G detector network will be able to observe. Other forthcoming GW detectors built in the more immediate future, such as the LIGO Voyager (Adhikari et al. 2020), also promise significantly increased rates of detection.
We end this section with a few caveats to our results. Any radio source count with a negative slope will find most objects just above the instrument flux density limit (see e.g. Mooley et al. 2013). Any changes in the final sensitivity of SKA1-Mid will thus have a big impact on the detection numbers. If we would set the detection threshold at a more conservative 10 RMS , for example, we would lose ∼ 40% -50% of the previously detected afterglows. A further source of uncertainty is the disk mass required for sGRB creation. We set the threshold at the lower end of the uncertainty range (Stone et al. 2013;Pannarale & Ohme 2014). Increasing the mass required will consequently lower the rates. Furthermore, we would like to iterate that we have used only one population synthesis model of the sizeable 420 model variations in B21. They find a large range in both the predicted GW detections per year as well as the percentage of those detections with any remnant mass over their model variations. We provide point estimates here for the fiducial model but remind the reader that such numbers are inherently surrounded by large population uncertainties given the few BHNS GW detections to date. In the future, even if no EM emission is detected, more GW observations will allow the population of BHNS mergers to be better understood.

Multi-messenger parameter inference
In Sect. 3 we consider single flux measurements of sGRB afterglows associated with BHNS mergers. A sole measurement above the detection threshold could already give some hints on the binary source properties by confirming the presence of ejecta mass. Much more physics, however, is encompassed in a light curve consisting of multiple data points. In the remainder of this paper we investigate, given a GW observation, what the added benefit of a radio light curve could be in inferring the binary source properties. Conversely, we examine the ability of a GW observation to help constrain the parameters of the sGRB afterglow too.

Bayesian framework
We took a Bayesian approach to incorporate both GW and EM radiation in our parameter inference. This method has been used extensively in analysing GW data (see e.g. Romero-Shaw et al. 2020) and we proceeded similarly to Raaĳmakers et al. (2021a) when incorporating the EM data.
To compute the posterior density function of a certain set of signal parameters ì we made use of Bayes' theorem. For a given signal and a model of the signal ℎ( ì ), the posterior density function is proportional to: ( ì | , ℎ( ì )) ∝ L ( | ì , ℎ( ì )) ( ì |ℎ( ì )), where L ( | ì , ℎ( ì )) is the likelihood function and ( ì |ℎ( ì )) is the prior information on ì given the model. For multiple signals, either gravitational or EM in nature, the joint likelihood becomes the product of the individual likelihoods assuming the noise streams in the various detectors are uncorrelated: In practice, we sampled from the joint log-likelihood where the product over the individual likelihoods is replaced with a sum over their natural log counterpart. To sample from log L joint , we employed the MultiNest nested sampler (Feroz et al. 2009) in its Python implementation incorporated into the bilby Python package (Ashton et al. 2019). While bilby is primarily developed for pure GW data analysis, its flexibility and extensive documentation make it well suited for our purposes. In the next sections we discuss the specific forms of the GW and EM likelihoods.

GW likelihood
We used a network of three GW detectors consisting of the two aLIGO detectors (in Hanford and Livingston) and the Virgo detector at their design sensitivities, hereafter referred to as the LHV network. The log-likelihood of this network takes the standard form (see e.g. Cutler & Flanagan 1994): where the inner product for two GW strains ( ) and ( ) is defined as Here the tilde superscript implies a Fourier transform, and ( ) is the power spectral density of the detector's noise. In bilby, the likelihood of Eq. (16)  For our GW waveform model ℎ GW ( ì GW ), we chose a recent model, 'SEOBNRv4_ROM_NRTidalv2_NSBH' (Matas et al. 2020), which is specifically tailored to aligned-spin BHNS mergers. It includes tidal effects and also accounts for NS tidal disruption in the late inspiral, merger and ring-down part of the waveform. The data GW in each detector was generated by projecting the same waveform model onto the detector frame and adding coloured Gaussian noise scaled by the power spectral density. The set of signal parameters ì GW used to generate the waveform consists of various parameters intrinsic and extrinsic to the source. While we have defined most of these parameters in previous sections already, we list them here for clarity: (i) the mass of the BH, BH , and the mass of the NS, NS . As it is relatively inefficient to sample in the two component masses because of strong correlations between them, we transformed them to the chirp mass, M , and the mass ratio, . The chirp mass is defined as M = ( NS · BH ) (3/5) ( NS + BH ) (1/5) , and = NS / BH ; (ii) the spin of the BH, BH , and the spin of the NS, NS ; (iii) the tidal deformability of the NS, Λ NS . We set the tidal deformability of the BH to zero; (iv) the inclination angle, , and the polarisation, ; and (v) the sky position in right ascension, , and declination, , and the luminosity distance, . By assuming an EM counterpart was detected and the source was thus localised to a host galaxy, we did not need to sample in these parameters and we could set them to their true injected values.
To fully specify the waveform, the time of coalescence and the phase of coalescence need to be given as well. Because these parameters are not of interest for our analysis but can drastically increase sampling time, we analytically marginalised over them using built-in functions in bilby.

EM likelihood
For our model ℎ EM ( ì EM ) of the jet synchrotron afterglow at radio frequencies, we turned to the afterglowpy Python package (Ryan et al. 2020  arising from structured jets. It produces reasonably accurate light curves compared to codes using relativistic numerical hydrodynamic simulations such as BoxFit (Eerten et al. 2012) at a small fraction of the compute time (Ryan et al. 2020). This makes it ideal for parameter estimation studies where the model needs to be calculated many times. We refer to Ryan et al. (2020) for further details on afterglowpy. We used the 'power law jet' model light curves from afterglowpy in this work. In general, these light curves are a function of the same parameters as the top-hat jet peak flux in Eq. 13. Some additional parameters are necessary for structured jets, however. We summarise here all the parameters ì EM needed to generate the light curves: (i) the observing frequency, , and observing times, obs ; (ii) the core opening angle of the jet, , which determines the effective width of the jet core; (iii) the on-axis isotropic-equivalent jet energy, 0 ; (iv) the observing angle, obs . As in Sect. 3.2, we assumed obs = ; (v) the truncation angle, , which determines how far the wings of the jet extend; (vi) the power law index, , which determines how energetic the wings of the jet are; (vii) the circumburst density, 0 . We treated 0 as being uniform, though not fixed, in value; (viii) the index of the Lorentz factor distribution, , of the accelerated electrons and the fraction, , of thermal energy in those electrons. As in Sect. 3.2, we set = 2.2 and = 0.1; (ix) the fraction, , of magnetic energy relative to thermal energy; (x) the fraction of accelerated electrons, . We fixed = 1, in line with Ryan et al. (2020); and (xi) the luminosity distance, , which we set to the true injected value, as in Sect. 4.2.
We assumed both Gaussian measurement errors RMS and systematic errors sys on data points EM generated using afterglowpy. This leads to the following likelihood for our EM data: −log(2 ( 2 RMS + 2 sys, )) .

Connecting the EM and GW likelihoods
Similar to Sect. 3, we expressed the kinetic jet energy, , in terms of the binary source parameters. In this way, we connect our EM measurements to our GW measurements. We thus did not sample in , or more precisely 0 , directly but used Eqs. 1 through 12 to calculate from M , , BH and Λ NS . After converting the resulting to 0 using the specific formula for structured jets in the afterglowpy code, the EM data points, EM , could be computed from the rest of the parameters in ì EM . Our complete set of sampling parameters is thus Our joint log likelihood is simply the sum of Eqs. 16 and 19: log L joint = log L GW + log L EM .
We did not assume any other dependence of the parameters in ì EM on ì GW . In reality, these dependences likely do exist to some degree. We return to this point in Sect. 6.

Systematics
The error sys represents the uncertainty that arises when converting the binary source parameters into EM data points using fits to NR data. These fits are not exact and we used their residual errors as input for sys .To incorporate the errors in the remnant mass and the dynamical ejecta, we generated 1000 samples from two Gaussian distributions in rem and dyn each with means given by Eq. 1 and Eq. 8, respectively. The standard deviation for each distribution is assumed to be a relative error of 15% in rem (Foucart et al. 2018) and an absolute error of dyn = 0.0047 in dyn (Krüger & Foucart 2020). For each generated combination of samples, we calculated the data points of our light curve as detailed in Sect. 4.4. We set sys, equal to the standard deviation of the resulting flux distribution for each data point EM . We discuss more potential sources of systematic errors in Sect. 6.

Fiducial BHNS merger
Log-Uniform (10 −5 , 1) 50, 100, 200 Mpc --For our fiducial BHNS merger, we chose binary source parameters comparable to GW200105 ) but adjusted to generate enough ejecta mass for an observable sGRB afterglow. Particularly, we set a high BH = 0.7, which is beneficial for detectability as shown in Fig. 3. We set Λ NS = 400, which implies a hard EOS for the corresponding NS mass. We chose a small NS = 0.02, consistent with the spins observed in galactic BNSs (Burgay et al. 2003). We picked a relatively large inclination angle, = 0.4 rad, consistent with GW170817 (see e.g. Ryan et al. 2020), and a random value for in (0, ). The other parameters, , , , 0 , and , all have fiducial values similar to the inferred values of GW170817 for the power law jet model in Ryan et al. (2020). While we did not sample in , we varied the injected value to generate GW and sGRB afterglow detections with different S/Ns. All fiducial parameters are listed in Table 1.

Setup
We injected an 8s BHNS merger signal, using the methods in Sect. 4.2, with our fiducial parameters into the LHV GW detector network starting at a frequency of 40 Hz and we used a cutoff frequency of 4096 Hz. We supposed the merger was sufficiently localised when the targeted radio follow-up was commenced. We then started observing with SKA1-Mid at = 1.43 GHz 11 days after the merger and took 20 observations geometrically spaced until 500 days post-merger. We set a constant SKA1-Mid measurement error at RMS = 2 Jy and assumed a data point was detected if it passed the 3 RMS threshold. Because we took multiple measurements of the same source instead of just a single observation, the detection threshold could be set lower for each measurement than in Sect. 3. All 20 observations at = 50 Mpc give a flux measurement above the threshold. At = 100 Mpc, the first and last two observations have flux measurements below this value and were not taken into account. At = 200 Mpc, only seven observations between approximately 55 and 183 days post-merger are above the detection threshold. We calculated the systematic errors as described in Sect. 4.5 for our fiducial parameter set. We show the resulting light curves in Fig. 5.
To perform the parameter estimation in ì , we chose the broad uninformative priors listed in Table 1. Specifically, we took uniform priors for M , , BH , NS , Λ NS , , , , and , while for 0 and we chose log-uniform priors. We assumed an isotropic prior for the orientation of the binary, . We initialised MultiNest with standard parameters through bilby and sampled using 1008 live points . For each , we did two runs: one incorporating both GW and EM data (log L joint ) and one using only GW data (log L GW ) for comparison.
For the two values of = 50 Mpc and = 100 Mpc, the GW signal network S/N is 117 and 59, respectively. These are well above the canonical threshold of S/N thresh ≥ 8. In those cases, we are focusing on the effects of an sGRB afterglow detection in the presence of an already loud GW signal. This is interesting for a few reasons. In the 3G detector network era, high GW S/N detections will become commonplace because of the increased sensitivity. There are other benefits to 3G detectors, such as their broader frequency range, that are not considered here. Still, a loud GW detection in our LHV network will serve as a reasonable testing ground for the extra information an sGRB afterglow detection with SKA1-Mid might give in this scenario. Furthermore, placing the merger relatively nearby allows all parts of the light curve to be well sampled at SKA1-Mid sensitivity. This enables a robust study of the radio afterglow in the limit of small relative measurement errors. Even so, as we describe in Sect. 3, most radio afterglows will be detected close to the detection limit. Here, the measurement errors become substantial relative to the observed flux. To investigate this scenario too, we also placed our fiducial binary at = 200 Mpc. At this distance, the light curve consists of fewer measurements, as shown in Fig. 5, and the relative error on those measurements is larger as well. The GW signal is also weaker but the network S/N of 29 is still above the detection threshold.
We show the results in the next section.

50 Mpc and 100 Mpc
In Fig. 6 we show the posteriors for ì GW for = 100 Mpc (left panel) and = 50 Mpc (right panel) when sampling from log L GW and when sampling from log L joint . The parameters ì EM are only inferred in the latter case. The posteriors of these parameters are shown in Fig. 7 for = 100 Mpc and for = 50 Mpc.
At 100 Mpc, including EM data of the afterglow provides a range of improvements in the estimates of the binary source  Fig. 7: Posterior distributions of the parameters ì EM . Purple contours give the inferred posterior distribution from the combination of the GW and radio data when the binary is placed at = 100 Mpc. Orange contours give the inferred posterior distribution from the combination of the GW and radio data when the binary is placed at = 50 Mpc. The fiducial injected values are indicated by the dashed black lines.
parameters over using the GW data alone. It is most notable in the marginalised posterior of , where the long tail towards more equal mass ratios can be ruled out using the SKA1-Mid observations. At high , corresponding to more equal BH and NS , the amount of disk mass likely becomes too large to still be consistent with the observed light curve. Both M and BH are well constrained from the GW data alone, and their estimates are slightly improved with the EM data added. Because the afterglow does not depend on either NS or in our model, the EM data provides no additional information on their already weak constraints by the GW data and we do not show their posteriors in Fig. 6.
Despite the high GW S/N at = 50 Mpc, the inclusion of EM data still gives notable improvements in the binary source parameter estimates at this distance too. Similar to the results at = 100 Mpc, the run with the EM data has less support for more equal mass ratios and gives some additional constraining power on M and BH compared to the run with only the GW data. At = 100 Mpc, the posterior range of Λ NS is somewhat reduced when including the EM data while at = 50 Mpc, the EM data gives no additional information on Λ NS . In both cases, the peak of the posteriors are shifted to higher values compared to the injected value. From our testing, we attribute this bias to the specific random noise realisation used for all inference runs (this is a known effect; see Nissanke et al. 2010).
An immediate advantage of knowing beforehand through a host galaxy association, is the resolved degeneracy between and in the GW inference ; The LIGO Scientific Collaboration and The Virgo Collaboration et al. 2017). The resulting GW data estimate on is essential to break various degeneracies in ì EM , specifically between and (Nakar & Piran 2021). Turning our attention to the posteriors of ì EM in Fig. 7, it is clear that the angles , and are quite well constrained at both distances. Still, at = 50 Mpc, the improved estimate of from the GW data helps improve the estimates on and further. The EM data gives almost no additional information on compared to only the GW data in Fig. 6. For the marginalised posterior of in Fig. 7, only the lower half of the prior range can  In the left panel, blue contours give the posterior distribution when only using the GW data, whereas the green contours give the posterior distribution when using both the GW and radio data. In the right panel, red contours give the inferred posterior distribution from the combination of the GW and radio data. The fiducial injected values are indicated by the dashed black lines. be confidently ruled out at either distance. We are able to infer 0 and to uncertainties of 0.5 − 2 dex at both distances.

200 Mpc
In Fig. 8 we show the posteriors at = 200 Mpc for ì GW (left panel) and ì EM (right panel). At this distance, it becomes considerably more difficult to fully characterise the binary. In regards to the binary source parameters, both M and BH are still inferred relatively well from the GW data, while the lower GW S/N impedes a precise estimate of , and Λ NS . Including the EM data gives some improvement in the inference of these parameters with the biggest improvement in akin to the results at = 50, 100 Mpc. Still, the large measurement and systematic errors on the radio observations allow the light curve to be only loosely constrained. This is clearly visible in the inference of the ì EM parameters. Most information on the geometry of the jet is now lost as the posterior of spans almost the entire prior range; for the estimate is similarly weak. The estimate on is weakened as well compared to the previous cases. From the posteriors of 0 and , we can only rule out the upper end of their respective prior ranges.

Discussion
In summary, we are able to recover most of the injected parameters of our fiducial binary with reasonable confidence at both 50 Mpc and 100 Mpc. The observed sGRB afterglow clearly provides additional constraints on the binary source parameters despite the high GW S/Ns. Conversely, crucial degeneracies are broken in the parameters of the afterglow light curve using the GW data. The GW observations thus play an essential role in fitting the corresponding EM data. At 200 Mpc, both the GW and EM data provide less information on the binary and the parameters of the sGRB afterglow in particular become fairly undetermined.
Inferring the parameters of both the GW data and the EM data simultaneously is not trivial. Various degeneracies come into play, arising from, for example, the conversion of four binary parameters to one jet energy parameter. These and other aspects warrant further discussion, below.
From the parameters in ì EM , and 0 are the most strongly correlated with the binary source parameters. We show the combined posterior of the binary source parameters together with 0 and in Appendix B for = 50, 100 Mpc. To improve our constraints on the binary parameters further, it would be helpful to thus constrain 0 and better too. This might be achieved, for in particular, by using additional EM data at different frequencies (e.g. Beniamini et al. 2015). If the observations at other frequencies are in the same cooling regime, however, such homothetic light curves give no additional information (e.g. Ryan et al. 2020). Better constraints on the binary parameters could in turn improve the estimates of 0 and .
For most of the data points on our light curves at = 50, 100 Mpc, the systematic errors described in Sect. 4.5 are substantially larger than the measurement errors. These systematics currently thus have a big effect on the inference power of well-observed EM data, where most observations have fluxes significantly above the sensitivity of a next-generation radio telescope such as SKA1-Mid. This is visible in Fig. 9, where we have computed the uncertainty in the fit of the light curve from the posterior samples of the combined GW and EM data inference at = 50 Mpc. The shaded area gives the 95% credible region of the resulting fit. The black error bars, corresponding to the SKA1-Mid measurement error on our 20 simulated observations, are much smaller than this uncertainty region for most data points. Thus, the systematic errors, and not the measurement errors, govern how well the light curve can be constrained at these distances. At 200 Mpc, the situation changes. While the systematic errors are not insignificant, the measurement errors are now the dominant source of uncertainty for all seven data points. Not only does this limit our ability to gain additional constraints on the binary parameters, it also hinders any attempt to resolve the structure of the jet. Besides, for example, a good estimate of the inclination angle, a well-measured afterglow peak flux and time are necessary to determine the jet opening angle (Nakar & Piran 2021). We conclude that this is not the case for our light curve at 200 Mpc and that this might not hold true for most afterglows that will be detected close the telescope sensitivity limit either.
For a low GW S/N, a well-measured afterglow light curve might provide better constraints on the binary source parameters than what we found with the poorly sampled light curve at 200 Mpc discussed in this work. However, the aforementioned systematics will still limit how much information can be gained. Even if those errors are brought down in the future, other systematics in the sGRB afterglow, such as interstellar scintillation at higher frequencies (Aksulu et al. 2020), can also have an effect on the EM data inference. A proper understanding of these and other systematic uncertainties and how they propagate through the Bayesian inference framework (see e.g. Carson et al. 2019 for a discussion on the systematics from different EOSs) will be essential to extract as much information as possible out of sGRB afterglow observations. Additional information on the inclination angle through, for example, Very Long Baseline Interferometry measurements (Mooley et al. 2018;Hotokezaka et al. 2019) will be required as well if this angle is not sufficiently inferred from the low S/N GW data. Furthermore, to fully characterise the sGRB afterglow, it might be necessary to infer , and too instead of setting them to often used canonical values (e.g. Aksulu et al. 2022). In this work, we have also ignored the uncertainties in the jet launching mechanism and focused on the systematics of Sect. 4.5. We leave a full discussion of the uncertainties in the launching mechanism for future work but give some brief remarks here. Salafia & Giacomazzo (2021, SG21) note that the fit of BZ in Eq. 11 pertains to simulations of BH-accretion disk systems hosted by active galactic nuclei. Such systems differ in terms of, for example, the accretion rate from BH merger remnants of compact object mergers. The estimated posterior of the accretionto-jet energy conversion efficiency of GW170817 in SG21 is still consistent with BZ . The spread in this posterior covers multiple decades in , however, and it is thus hard to make definitive conclusions. More EM detections accompanying BNS mergers in the future will bring down the uncertainty in for this class of compact object mergers. It remains to be seen whether this efficiency is similar for BHNS mergers given, for example, the differences in ejecta properties. Updated NR simulations will help in answering this question, even if an EM detection of a BHNS sGRB afterglow does not occur soon.
Given the large uncertainties in deriving the energy of the jet from the binary source parameters, it will also be necessary to better understand the dependence of jet features, such as the opening angle on the source properties. Lazzati et al. (2021), for instance, perform a hydrodynamic simulation of a relativistic jet propagating in BNS merger ejecta and characterise the jet properties. A (semi-)analytical description derived from the mapping of the output of many such simulations to a variety of associated input source parameters will be a valuable tool in the inference of a single observed source. It is worth bearing in mind that the influence of the central BH engine and surrounding ejecta on the jet propagation and opening angle is still under much debate, certainly in BHNS mergers (Kyutoku et al. 2021). Systematics will be a limiting factor when incorporating such descriptions too. We have only focused on sGRB afterglow observations in this work. Combining multiple sources of EM radiation (see e.g. Dietrich et al. 2020) will be key to exploit the full potential of multi-messenger astronomy. Incorporating kilonova physics and emission (Raaĳmakers et al. 2021a;Nicholl et al. 2021) could improve the estimates on the binary source parameters further, which we will explore in future work.
We end this discussion with some technical considerations. Perhaps the most straightforward and often utilised way of combining multi-messenger emission is taking the marginalised posteriors of the analysis of one emission type as prior input for the analysis of another type of emission. Nicholl et al. (2021), for example, took 1D GW parameter estimates of GW170817 as input for some of the priors of their kilonova model. In this way, a (re-)analysis of the GW data is avoided, making the inference clearly more efficient computationally. Furthermore, it is not always necessary to do a full re-analysis of the data if information on only one parameter (e.g. the inclination angle) is needed, which is readily available from the marginalised posterior. If there are strong correlations in the parameter space, however, such a sequential approach comes at the cost of losing that valuable information between sets of parameters. This information is preserved when the full likelihood of the GW data is used as in our work. To consistently combine not only GW and afterglow emission but also kilonova emission moving forwards, we argue that a joint likelihood Bayesian analysis is preferred (and perhaps necessary) over a sequential approach because of the many correlations between the parameters.
A method that sits somewhere in between a sequential and a joint approach is making a multi-dimensional density estimate of the sampled GW posterior to preserve the correlations. This estimate can then be used in place of the full likelihood. We have experimented with kernel density estimates, generated by software packages such as kalepy (Kelley 2021), but did not find these approximations to always converge to the same results as the full likelihood. We suspect this to be a result of insufficient accuracy in the estimates of the tails of the distributions. Other density estimates, such as methods based on Gaussian processes (D'Emilio et al. 2021), could fair better in this regard but exploring such alternatives was beyond the scope of this work.

Summary and conclusion
In this paper we have looked at the prospects for detecting and analysing BHNS mergers using both GW and radio emission. We modelled the radio emission from an accompanying sGRB jet afterglow launched from the ejecta disk surrounding the merger remnant. We used fits from the literature, connecting the binary source properties to the accreted disk mass and to the accretionto-jet-energy conversion efficiency of the jet.
We first performed a population synthesis study to derive the expected rates of GW BHNS observations with an sGRB afterglow detection at SKA1-Mid radio frequency and sensitivity. We explored the impact of the current generation of GW detectors (2G) versus the future generation (3G), the NS EOS, and the BH spin of the binary, BH . We find that rates of around one combined detection per year with a 2G detector network and SKA1-Mid are only possible for an unrealistically high BH of ∼ 0.8. This is a result of the strong dependence of the jet energy, and thus radio flux, on BH . The increased sensitivity of a 3G network will allow for similar rates near one combined detection per year for lower BH of ∼ 0.2. Such values of the BH spin are more in line with the expected values in BHNS mergers based on current GW detections. The probability of finding a combined GW and radio detection of a BHNS merger in the near future is thus low. This will increase substantially with forthcoming GW detectors that will be able to detect many more BHNS mergers, increasing the chances of an sGRB localisation and radio detection. Furthermore, transitioning from a soft to a hard NS EOS increases the expected rates two-to three-fold.
We then examined our ability to infer the properties of a fiducial BHNS merger using simulated aLIGO and Virgo GW observations in conjunction with SKA1-Mid radio observations of its sGRB afterglow. We employed a recent waveform model for the GW signal and modelled the full afterglow light curve using afterglowpy. We performed a joint Bayesian analysis, combining the likelihood of the GW data with the likelihood of the radio data. We also incorporated systematic errors from the conversion of the binary source parameters to the disk mass. We placed our fiducial binary at three distances to study the effect of different GW and sGRB afterglow detection S/Ns. We find that it is possible to simultaneously recover both the binary source parameters, such as the chirp mass, and parameters of the jet afterglow, such as the opening angle, with reasonable confidence for nearby binaries. When placing our fiducial binary at 100 Mpc, including the radio data provided various improvements in the inference of the binary source parameters compared to just using the GW data. This was most noticeable in the wellimproved estimate of the mass ratio. At 50 Mpc, despite the already high GW S/N, the estimates of the binary source parameters were also improved with the EM data used in the inference. At these distances, we find that the systematic errors were likely the limiting factor in how much information could be extracted from the radio light curve. On the other hand, the inclination angle estimate of the GW data broke various degeneracies in the afterglow jet parameters of the radio data. Specifically, parameters pertaining to the geometry of the jet were well determined from the combination of GW and radio data at both distances. Other parameters, such as the circumburst density, were recovered with larger relative uncertainties. Thus, in a close-by real life detection, we can be reasonably confident that we can determine the properties of the system. If the binary is farther away, it becomes much harder to get a complete picture of the system. For our fiducial binary at 200 Mpc, only the chirp mass and BH spin were easily inferred, while most of the information on the jet parameters was lost. If an afterglow is observed close to the telescope detection threshold, as for our binary at 200 Mpc, we conclude that the benefits of including radio data in the analysis of the binary are limited.
In the future, it will be necessary to get a more complete understanding of the connection between the binary source parameters and the resulting properties of the jet afterglow. Lowering the associated systematic uncertainties is important for really leveraging the depth of information encompassed in the afterglow signal. Numerical relativity simulations will continue to be instrumental in this regard, notably when combined with future GW detections of BHNS mergers to constrain the physical processes further.
To summarise, we conclude that the upcoming generation of GW detectors will provide a unique opportunity to study the population of BHNS mergers, even in the absence of an sGRB afterglow. For high S/N GW signals, which will be routine occurrences with future GW detectors, tight constraints on the inferred binary parameters are possible. If an associated afterglow is observed, which is, as we have argued, not unlikely in the 3G detector era, a combined analysis of the GW and radio data will be essential for characterising the jet afterglow properties. For a sufficiently bright afterglow, this joint analysis can also provide additional information on the binary parameters, but only if the modelling systematics are well understood.

Appendix A: Additional BHNS population figures
In Fig. A.1, we show the flux distribution of sGRB afterglows associated with BHNS mergers detected in GWs with a 3G detector network including ET. The different distributions correspond to different values for the average BH used in the simulations. Also displayed is the SKA1-Mid 5 RMS detection limit. Clearly, most of the sGRB afterglows are below the detection limit even with a radio telescope as sensitive as SKA1-Mid. As mentioned in Sect. 3, any change in the sensitivity of SKA1-Mid will have a big influence on the detection rates. This is especially evident for the distribution with BH = 0.8, which has a very steep slope below the 5 RMS detection limit. In Fig. A.2, we show the combined GW and radio detection rates as a function of 0 and . All rates are normalised to their respective total rate given in Fig. 3. We used the normalised rates of the ET and BH = {0.4, 0.8} combination and the aLIGO and BH = 0.8 combination. We assumed a hard NS EOS. Compared to the intrinsic distribution for 0 , the means of the detected distributions are a factor of 2.5 − 3.5 higher as larger circumburst densities lead to larger peak fluxes (Eq. 13). For , this selection effect is less prominent, with means ∼ 1.5 higher in the detected distributions. While the intrinsic distributions for 0 and are similar, it is truncated at 10 −2 for leading to lower means in the detected distributions.

Appendix B: Additional posterior of GW and EM data inference
In Fig. B.1, we show the combined posterior distribution of the binary source parameters together with 0 and at distances = 50, 100 Mpc. Both 0 and have varying degrees of correlation with the binary source parameters. We hypothesise that in the setup used in this work, improved estimates on 0 and , through different sources of EM emission for example, can help improve the estimates on the binary source parameters further. If we are instead able to improve our estimates on the binary source parameters, this will, in turn, also help with inferring 0 and .