Issue 
A&A
Volume 664, August 2022



Article Number  A160  
Number of page(s)  15  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202243267  
Published online  25 August 2022 
Investigating the detection rates and inference of gravitationalwave and radio emission from black hole neutron star mergers^{⋆}
^{1}
Anton Pannekoek Institute, University of Amsterdam, Postbus 94249, 1090 GE Amsterdam, The Netherlands
email: o.m.boersma@uva.nl
^{2}
ASTRON, The Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The Netherlands
Received:
4
February
2022
Accepted:
16
May
2022
Context. Black hole neutron star (BHNS) mergers have recently been detected through their gravitationalwave (GW) emission. While no electromagnetic emission has yet been confidently associated with these systems, observing any such emission could provide information on, for example, the neutron star equation of state. Black hole neutron star mergers could produce electromagnetic emission as a short gammaray burst (sGRB) and/or an sGRB afterglow upon interaction with the circummerger medium.
Aims. 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.
Methods. We simulated a population of BHNS mergers, making use of recent stellar population synthesis results, and estimated 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 spin, and the neutron star equation of state. We then performed a multimessenger Bayesian inference study on a fiducial BHNS merger. We simulated its sGRB afterglow and GW emission as input to this study, using recent models for both, and take systematic errors into account.
Results. The expected rates of a combined GW and radio detection with the currentgeneration GW detectors are likely low. Due to the much increased sensitivity of future GW detectors such as the Einstein Telescope, the chances of an sGRB localisation and radio detection increase substantially. The unknown distribution of the black hole spin has a big influence on the detection rates, however, and it is a large source of uncertainty. Furthermore, when placing our fiducial BHNS merger at 50 and 100 Mpc, we are able to infer both the binary source parameters and the parameters of the sGRB afterglow simultaneously if we combine the GW and radio data. The radio data provide useful extra information on the binary parameters, such as the mass ratio, but this is limited by the systematic errors involved. For our fiducial binary at 200 Mpc, it is considerably more difficult to adequately infer the parameters of the system.
Conclusions. The probability of finding an sGRB afterglow of a BHNS merger is low in the near future but will rise significantly when the nextgeneration GW detectors come online. Combining information from GW data with radio data is crucial for characterising the jet properties. 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.
Key words: gravitational waves / stars: neutron / stars: black holes / radio continuum: stars
Data used to plot the images have been uploaded at: http://doi.org/10.5281/zenodo.6573093
© O. M. Boersma and J. van Leeuwen 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1. Introduction
The detection of the first two black hole neutron star (BHNS) mergers, GW200105 and GW200115, by the advanced Laser Interferometer GravitationalWave Observatory (aLIGO) and Virgo detector network (Aasi et al. 2015; Acernese et al. 2014) completed the trifecta of compact binary coalescence gravitationalwave (GW) observations (Abbott et al. 2021a). The third GW catalogue (GWTC3) 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 2021b). The Kamioka Gravitational Wave Detector (KAGRA; Somiya 2012; The KAGRA Collaboration 2013) came online and was included in the network in early 2020. This rich source of data has brought forth numerous exciting results on its own in areas of research such as tests of general relativity (The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration 2021c) and the cosmic expansion history (The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration 2021a).
Of particular interest to astronomers are the GW observations of BNS and BHNS systems because of the possibility of a complementary electromagnetic (EM) signature^{1}. This possibility was first confirmed with the BNS source GW170817 (LIGO Scientific Collaboration and Virgo Collaboration 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 gammaray 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 followup campaigns being performed after their detections (e.g., Antier et al. 2020; Abbasi et al. 2021; Abe et al. 2021; Paterson et al. 2021; Ridnaia et al. 2020; Kasliwal et al. 2020; Ashkar et al. 2021; Anand et al. 2021)^{2}.
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 circummerger 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 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 multimessenger analyses have been performed extensively on GW170817 (e.g., Radice et al. 2018; Coughlin et al. 2018; Radice & Dai 2019; Coughlin et al. 2019; Raaijmakers et al. 2020; Capano et al. 2020; Dietrich et al. 2020; Breschi et al. 2021b), and various general Bayesian multimessenger frameworks have been developed with this goal in mind (e.g., Breschi et al. 2021a; Raaijmakers 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 multimessenger 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 multimessenger 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.
2. Ejecta outflows of BHNS mergers
2.1. 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 Q = M_{BH}/M_{NS}, with M_{BH} the mass of the BH and M_{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 nonprecessing 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 bestfit 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 ():
with bestfit values α = 0.406, β = 0.139, γ = 0.255, and δ = 1.761. Here, η = Q/(1 + Q)^{2} and R̂_{ISCO} = R_{ISCO}/M_{BH} is the innermost stable circular orbit radius normalised by the BH mass:
To compute the compactness of the NS C_{NS}, we used the approximately universal CLove relation for NSs of Eq. (78) in Yagi & Yunes (2017):
with bestfit coefficients a_{0} = 0.360, a_{1} = −0.0355 and a_{2} = 0.000705. While this relation is not a perfect substitute for integrating a NS EOS to obtain C_{NS}, it performs well across a wide range of EOSs and has much less computational overhead. From C_{NS} we could, again approximately, compute (Lattimer & Prakash 2001):
The dependence of M_{rem} on χ_{BH}, Q, and Λ_{NS} is visualised in Fig. 1.
Fig. 1. Visualisation of the fitting formula for the remnant mass in Foucart et al. (2018), showing the dependence on the BH spin (χ_{BH}), mass ratio (Q), and NS tidal deformability (Λ_{NS}). Bestfit values for the free parameters are assumed. In each of the three panels, χ_{BH}, Q, 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 (M_{rem}) in solar masses. 
Given M_{rem}, we can calculate the amount of disk mass (M_{disk}) if we have an estimate of the amount of unbound dynamical ejecta mass (M_{dyn}):
We again turned to fits of NR simulations (Krüger & Foucart 2020) for M_{dyn}:
with bestfit coefficients a_{1} = 0.007116, a_{2} = 0.001436, a_{4} = −0.02762, n_{1} = 0.8636, and n_{2} = 1.6840.
Parts of the disk itself can also become unbounded through disk wind outflows that are either thermally or magnetically driven. Raaijmakers et al. (2021a) derive a simple formula that broadly captures the dependence found in NR simulations of the disk wind ejecta mass M_{ej, disk} on Q (Fernández et al. 2020):
Here we assumed average values for the free parameters of ξ_{1} = 0.18 and ξ_{2} = 0.29.
2.2. 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 ultrarelativistic 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 ultrarelativistic outflow can overcome such ejecta with sufficient collimation. Just et al. (2016) still find that BHNS mergers are able to harbour ultrarelativistic 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 accretiontojet 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 accretiontojet energy conversion efficiency (see Barbieri et al. 2019 for a similar derivation). For a total accreted disk mass M_{acc} = M_{disk} − M_{ej, disk}, we took M_{acc} ≥ 0.03 M_{⊙} as a necessary condition to create an sGRB of ∼1 s duration (Stone et al. 2013; Pannarale & Ohme 2014).
Taking a fairly typical gammaray burst prompt emission efficiency f_{γ} 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 onehalf accounts for an identical counterjet. McKinney & Gammie (2004) describe generalrelativistic magnetohydrodynamical simulations of spinning Kerr BHs with thick accretion disks. SG21 fit the following functions, extending the fit of McKinney (2005), to those simulations for the BZ accretiontojet energy conversion efficiency:
where is the spin of the final remnant BH and Ω_{H} is the dimensionless angular frequency at the BH horizon:
To capture the dependence of 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énezForteza et al. 2017). Compared to the formulae stated in SG21, in the regime an additional factor of 10^{−2} is necessary (Salafia & Giacomazzo 2022).
Analogous to Fig. 1, the dependence of E_{k} on χ_{BH}, Q, 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 M_{rem}, comparable mass binaries with high Λ_{NS} produce the highest jet energies. There is a strong dependence of E_{k} on χ_{BH} as well. This is partly inherited from M_{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.
Fig. 2. Same as Fig. 1 except now the final kinetic energy of the jet in Eq. (10) is visualised. The regions of the parameter space with no colour indicate either that M_{disk} = 0 or that the minimal accreted disk mass for sGRB creation, M_{acc} = 0.03 M_{⊙}, is not reached. 
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.
3. 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.
3.1. Population synthesis of BHNS mergers
B21 provide both the final properties of the BHNS mergers after star formation and the merger rate per redshift (z). We followed their methods, which are derived from Neijssel et al. (2019), and integrated over 250 redshift shells from z = 0 to z = 0.5 to get the total amount of BHNS mergers per year in this volume of spacetime. 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 M_{BH} and M_{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 signaltonoise 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; Husa et al. 2016; 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 secondgeneration (2G) GW detectors at design sensitivity and a future network of thirdgeneration (3G) detectors including the Einstein Telescope (ET) with the ETD 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 followup, 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 gammarayburst jet E_{k} 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 R_{NS} = 11.5 km (consistent with GW observations; see e.g., The LIGO Scientific Collaboration and The Virgo Collaboration 2018) or R_{NS} = 13.0 km (consistent with NICER observations; see Miller et al. 2019; Raaijmakers et al. 2021b; Riley et al. 2021). As the distribution of BH spins in BHNS mergers is still highly uncertain, we used an equal (average) value for all the mergers in a single population and varied this value as χ_{BH} = {0, 0.2, 0.4, 0.6, 0.8} to study its influence. In the next section we use E_{k} to obtain the peak flux magnitude of the possible sGRB afterglow.
3.2. 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 θ_{c} is the opening angle of the jet core, n_{0} is the circumburst density, ϵ_{e}, ϵ_{B} and p are shock microphysics parameters, ν is the observing frequency, d_{L} 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 E_{k} but on the onaxis isotropicequivalent jet energy E_{0}. We used a standard formula for tophat jets to convert between the two: E_{0} = E_{k}/(1 − cos(θ_{c})), where the opening angle is fixed to a representative value of θ_{c} = 0.1 rad (≈5.7 deg) (Beniamini & van der Horst 2017). The literature on gammarayburst afterglow modelling often uses such a simplified tophat jet approximation instead of a structured jet model (see e.g., Aksulu et al. 2022 for a recent study). The tophat 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 offaxis. 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 tophat 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 p = 2.2 and ϵ_{e} = 0.1 to fiducial values. These values are consistent with the observed gammarayburst population (Beniamini & van der Horst 2017; Aksulu et al. 2022). Both n_{0} and ϵ_{B} have a much broader population distributions than ϵ_{e} or p so we did not fix these parameters. We took the same approach as Duque et al. (2019) and considered a lognormal distribution for both parameters with mean μ = 10^{−3} and standard deviation σ = 0.75. Furthermore, ϵ_{B} 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 SKA1Mid 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 F_{p, 1.43} passes the 5σ_{rms} threshold.
We have used 20 combinations of detector type, R_{NS} and χ_{BH} in total. We ran each combination for 1000 realisations of one year. We show the averaged detection rates in Fig. 3.
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 SKA1Mid 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 (R_{NS} = 11.5 km) or a hard NS EOS (R_{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. 
3.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 R_{NS} = 11.5 km and R_{NS} = 13.0 km, respectively. We provide some additional figures of the simulated populations in Appendix A. We focus on Figs. 3 and 4 below.
Fig. 4. Normalised GW and radio detection rates for the ET and χ_{BH} = {0.4, 0.8} combinations (left and middle panel) and for the aLIGO and χ_{BH} = 0.8 combination (right panel). The rates are normalised to the combined detection rates of Fig. 3, given by (from left to right): 3.3 yr^{−1}, 44 yr^{−1}, and 1.7 yr^{−1}. A hard NS EOS (R_{NS} = 13.0 km) is used in all cases. The rates are shown as a function of the distance and inclination angle of the detected sources. 
For a network of 2G detectors, we find that the SKA1Mid detection rates of sGRB afterglows associated with BHNS GW events are likely low. Only for a population of BHNSs with probably unrealistically high χ_{BH} ∼ 0.8 (Fragione 2021), the combined detection rate per year nears unity. Transitioning from a soft (R_{NS} = 11.5 km) to a hard (R_{NS} = 13.0 km) NS EOS does increase the rates two to threefold 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 nonzero 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 SKA1Mid 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 spacetime volume, which enables many more combined detections with SKA1Mid. 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 SKA1Mid 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 onaxis (ι ≤ θ_{c}).
A 3G network is expected to detect mergers at significantly larger redshifts than our maximum of z = 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 < z < 0.5. This implies SKA1Mid will be sensitive enough to detect afterglows with z > 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 SKA1Mid combination further but we note that a sufficient localisation for followup is not guaranteed at larger distances (Maggiore et al. 2020).
The combined detection rates are strongly dependent on the assumed χ_{BH} and an approximate loglinear relation is visible in Fig. 3. We attribute this partly to the formulation of the BZ accretiontojet energy conversion efficiency in Eq. (11). E_{0} is directly proportional to this efficiency and by extension F_{p, 1.43} as well. For our set of χ_{BH}, we are usually in the regime where (Zappa et al. 2019). As the final BH spin is approximately linear in χ_{BH} (Zappa et al. 2019), the loglinear dependence on χ_{BH} follows naturally from Eq. (11) in this regime. The dependence of M_{acc} on χ_{BH} is less intuitively understood, but we observe a roughly loglinear 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 SKA1Mid 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.
4. Multimessenger 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.
4.1. 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., RomeroShaw et al. 2020) and we proceeded similarly to Raaijmakers 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 d and a model of the signal h(θ), the posterior density function is proportional to:
where ℒ(dθ, h(θ)) is the likelihood function and p(θh(θ)) 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 N individual likelihoods assuming the noise streams in the various detectors are uncorrelated:
In practice, we sampled from the joint loglikelihood where the product over the individual likelihoods is replaced with a sum over their natural log counterpart. To sample from logℒ_{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.
4.2. 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 loglikelihood of this network takes the standard form (see e.g., Cutler & Flanagan 1994):
where the inner product for two GW strains a(t) and b(t) is defined as
Here the tilde superscript implies a Fourier transform, and S_{n}(f) is the power spectral density of the detector’s noise. In bilby, the likelihood of Eq. (16) is calculated using the GravitationalWaveTransient class (see Ashton et al. 2019 for more details).
For our GW waveform model h_{GW}(θ_{GW}), we chose a recent model, ‘SEOBNRv4_ROM_NRTidalv2_NSBH’ (Matas et al. 2020), which is specifically tailored to alignedspin BHNS mergers. It includes tidal effects and also accounts for NS tidal disruption in the late inspiral, merger and ringdown part of the waveform. The data d^{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, M_{BH}, and the mass of the NS, M_{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, ℳ_{c}, and the mass ratio, q^{3}. The chirp mass is defined as
(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, d_{L}. 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 t_{c} and the phase of coalescence ϕ_{c} 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 builtin functions in bilby.
4.3. EM likelihood
For our model h_{EM}(θ_{EM}) of the jet synchrotron afterglow at radio frequencies, we turned to the afterglowpy Python package (Ryan et al. 2020). It provides functionality to quickly calculate synthetic light curves at various frequencies for afterglows 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 tophat 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, t_{obs}; (ii) the core opening angle of the jet, θ_{c}, which determines the effective width of the jet core; (iii) the onaxis isotropicequivalent jet energy, E_{0}; (iv) the observing angle, θ_{obs}. As in Sect. 3.2, we assumed θ_{obs} = ι; (v) the truncation angle, θ_{w}, which determines how far the wings of the jet extend; (vi) the power law index, b, which determines how energetic the wings of the jet are; (vii) the circumburst density, n_{0}. We treated n_{0} as being uniform, though not fixed, in value; (viii) the index of the Lorentz factor distribution, p, of the accelerated electrons and the fraction, ϵ_{e}, of thermal energy in those electrons. As in Sect. 3.2, we set p = 2.2 and ϵ_{e} = 0.1; (ix) the fraction, ϵ_{B}, of magnetic energy relative to thermal energy; (x) the fraction of accelerated electrons, ξ_{N}. We fixed ξ_{N} = 1, in line with Ryan et al. (2020); and (xi) the luminosity distance, d_{L}, 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 k data points d^{EM} generated using afterglowpy. This leads to the following likelihood for our EM data:
4.4. Connecting the EM and GW likelihoods
Similar to Sect. 3, we expressed the kinetic jet energy, E_{k}, 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 E_{k}, or more precisely E_{0}, directly but used Eqs. (1) through (12) to calculate E_{k} from ℳ_{c}, q, χ_{BH} and Λ_{NS}. After converting the resulting E_{k} to E_{0} using the specific formula for structured jets in the afterglowpy code, the EM data points, d^{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ℒ_{joint} = logℒ_{GW} + logℒ_{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.
4.5. 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 M_{rem} and M_{dyn} each with means given by Eqs. (1) and (8), respectively. The standard deviation for each distribution is assumed to be a relative error of 15% in M_{rem} (Foucart et al. 2018) and an absolute error of σ_{dyn} = 0.0047 M_{⊙} in M_{dyn} (Krüger & Foucart 2020). For each generated combination of samples, we calculated the k data points of our light curve as detailed in Sect. 4.4. We set σ_{sys, i} equal to the standard deviation of the resulting flux distribution for each data point d_{i}^{EM}. We discuss more potential sources of systematic errors in Sect. 6.
Having specified our Bayesian framework, we apply it to a fiducial BHNS merger.
4.6. Fiducial BHNS merger
For our fiducial BHNS merger, we chose binary source parameters comparable to GW200105 (Abbott et al. 2021a) 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, θ_{c}, θ_{w}, b, n_{0}, and ϵ_{B}, 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 d_{L}, we varied the injected value to generate GW and sGRB afterglow detections with different S/Ns. All fiducial parameters are listed in Table 1.
Parameters of our fiducial BHNS merger and the chosen prior type and range for the parameter inference.
4.7. 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 followup was commenced. We then started observing with SKA1Mid at ν = 1.43 GHz 11 days after the merger and took 20 observations geometrically spaced until 500 days postmerger. We set a constant SKA1Mid 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 d_{L} = 50 Mpc give a flux measurement above the threshold. At d_{L} = 100 Mpc, the first and last two observations have flux measurements below this value and were not taken into account. At d_{L} = 200 Mpc, only seven observations between approximately 55 and 183 days postmerger 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.
Fig. 5. sGRB afterglow radio light curves, computed with afterglowpy, of our fiducial binary at the nominal SKA1Mid frequency, ν = 1.43 GHz. The curves are plotted on a log–log scale. Twenty simulated observations were taken between 11 and 500 days postmerger. The orange circles indicate the 20 observations above the 3σ_{rms} detection threshold when the binary is placed at d_{L} = 50 Mpc. The purple circles indicate the 16 observations above the 3σ_{rms} detection threshold when the binary is placed at d_{L} = 100 Mpc. The red circles indicate the seven observations above the 3σ_{rms} detection threshold when the binary is placed at d_{L} = 200 Mpc. The solid black line shows the SKA1Mid 3σ_{rms} detection threshold with σ_{rms} = 2 μJy. The error bars give the total error as a combination of the measurement error and the systematic error. 
To perform the parameter estimation in θ, we chose the broad uninformative priors listed in Table 1. Specifically, we took uniform priors for ℳ_{c}, q, χ_{BH}, χ_{NS}, Λ_{NS}, ψ, θ_{c}, θ_{w}, and b, while for n_{0} and ϵ_{B} we chose loguniform 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^{4}. For each d_{L}, we did two runs: one incorporating both GW and EM data (logℒ_{joint}) and one using only GW data (logℒ_{GW}) for comparison.
For the two values of d_{L} = 50 Mpc and d_{L} = 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 SKA1Mid might give in this scenario. Furthermore, placing the merger relatively nearby allows all parts of the light curve to be well sampled at SKA1Mid 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 d_{L} = 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.
5. Results
5.1. 50 Mpc and 100 Mpc
In Fig. 6 we show the posteriors for θ_{GW} for d_{L} = 100 Mpc (left panel) and d_{L} = 50 Mpc (right panel) when sampling from logℒ_{GW} and when sampling from logℒ_{joint}. The parameters θ_{EM} are only inferred in the latter case. The posteriors of these parameters are shown in Fig. 7 for d_{L} = 100 Mpc and for d_{L} = 50 Mpc.
Fig. 6. Posterior distributions of the parameters θ_{GW} as inferred for our fiducial binary at d_{L} = 100 Mpc (left panel) and d_{L} = 50 Mpc (right 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. The dark and light shaded regions give the 68% and 95% confidence levels, respectively, for all posteriors shown. The fiducial injected values are indicated by the dashed black lines. 
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 d_{L} = 100 Mpc. Orange contours give the inferred posterior distribution from the combination of the GW and radio data when the binary is placed at d_{L} = 50 Mpc. The fiducial injected values are indicated by the dashed black lines. 
At 100 Mpc, including EM data of the afterglow provides a range of improvements in the estimates of the binary source parameters over using the GW data alone. It is most notable in the marginalised posterior of q, where the long tail towards more equal mass ratios can be ruled out using the SKA1Mid observations. At high q, corresponding to more equal M_{BH} and M_{NS}, the amount of disk mass likely becomes too large to still be consistent with the observed light curve. Both ℳ_{c} 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 d_{L} = 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 d_{L} = 100 Mpc, the run with the EM data has less support for more equal mass ratios and gives some additional constraining power on ℳ_{c} and χ_{BH} compared to the run with only the GW data.
At d_{L} = 100 Mpc, the posterior range of Λ_{NS} is somewhat reduced when including the EM data while at d_{L} = 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 d_{L} beforehand through a host galaxy association, is the resolved degeneracy between d_{L} and ι in the GW inference (Abbott et al. 2017; The LIGO Scientific Collaboration and The Virgo Collaboration 2017). The resulting GW data estimate on ι is essential to break various degeneracies in θ_{EM}, specifically between ι and θ_{c} (Nakar & Piran 2021). Turning our attention to the posteriors of θ_{EM} in Fig. 7, it is clear that the angles ι, θ_{c} and θ_{w} are quite well constrained at both distances. Still, at d_{L} = 50 Mpc, the improved estimate of ι from the GW data helps improve the estimates on θ_{c} and θ_{w} further. The EM data gives almost no additional information on ι compared to only the GW data in Fig. 6. For the marginalised posterior of b in Fig. 7, only the lower half of the prior range can be confidently ruled out at either distance. We are able to infer n_{0} and ϵ_{B} to uncertainties of 0.5−2 dex at both distances.
5.2. 200 Mpc
In Fig. 8 we show the posteriors at d_{L} = 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 ℳ_{c} and χ_{BH} are still inferred relatively well from the GW data, while the lower GW S/N impedes a precise estimate of q, ι and Λ_{NS}. Including the EM data gives some improvement in the inference of these parameters with the biggest improvement in q akin to the results at d_{L} = 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 θ_{c} spans almost the entire prior range; for b the estimate is similarly weak. The estimate on θ_{w} is weakened as well compared to the previous cases. From the posteriors of n_{0} and ϵ_{B}, we can only rule out the upper end of their respective prior ranges.
Fig. 8. Posterior distributions of the parameters θ_{GW} (left panel) and θ_{EM} (right panel) as inferred for our fiducial binary at d_{L} = 200 Mpc. 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. 
6. 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}, ϵ_{B} and n_{0} are the most strongly correlated with the binary source parameters. We show the combined posterior of the binary source parameters together with n_{0} and ϵ_{B} in Appendix B for d_{L} = 50, 100 Mpc. To improve our constraints on the binary parameters further, it would be helpful to thus constrain n_{0} and ϵ_{B} better too. This might be achieved, for ϵ_{B} 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 n_{0} and ϵ_{B}.
For most of the data points on our light curves at d_{L} = 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 wellobserved EM data, where most observations have fluxes significantly above the sensitivity of a nextgeneration radio telescope such as SKA1Mid. 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 d_{L} = 50 Mpc. The shaded area gives the 95% credible region of the resulting fit. The black error bars, corresponding to the SKA1Mid 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 wellmeasured 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.
Fig. 9. Resulting fit to the light curve of our fiducial binary at d_{L} = 50 Mpc. The 95% credible region of the fit, indicated by the shaded region, is computed from the posterior samples of the combined GW and EM data inference in Sect. 5. The black circles give the 20 simulated observations, and the error bars correspond to the SKA1Mid measurement errors. 
For a low GW S/N, a wellmeasured 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 ϵ_{e}, p and ξ_{N} 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 BHaccretion 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 accretiontojet 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 multimessenger astronomy. Incorporating kilonova physics and emission (Raaijmakers 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 multimessenger 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 reanalysis 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 multidimensional 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.
7. 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 accretiontojetenergy 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 SKA1Mid 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 SKA1Mid 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 threefold.
We then examined our ability to infer the properties of a fiducial BHNS merger using simulated aLIGO and Virgo GW observations in conjunction with SKA1Mid 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 closeby 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.
It has been hypothesised that BBH mergers can also be a source of EM emission (see e.g., Liebling & Palenzuela 2016).
See also GCN Archive for S200105ae 2020, GCN Archive for S200115j 2020.
Acknowledgments
We thank Pikky Atri, Geert Raaijmakers and Samaya Nissanke for the helpful discussions and suggestions on this work. We thank the anonymous referee for their thoughtful comments which have improved this work. This research was supported by Vici research program ‘ARGO’ with project number 639.043.815, financed by the Dutch Research Council (NWO). J.V.L. further acknowledges funding through CORTEX (NWA.1160.18.316), under the research programme NWAORC, financed by NWO; and from the European Research Council under the European Union’s Seventh Framework Programme (FP/20072013)/ERC Grant Agreement No. 617199 (‘ALERT’).
References
 Aasi, A. J., Abbott, B. P., Abbott, R., et al. 2015, Class. Quant. Grav., 32, 074001 [Google Scholar]
 Abbasi, R., Ackermann, M., Adams, J., et al. 2021, ApJ, submitted [arXiv:2105.13160] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L13 [CrossRef] [Google Scholar]
 Abbott, R., Abbott, T. D., Abraham, S., et al. 2021a, ApJ, 915, L5 [NASA ADS] [CrossRef] [Google Scholar]
 Abbott, R., Abbott, T. D., Abraham, S., et al. 2021b, ApJ, 913, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Abe, S., Asami, S., Gando, A., et al. 2021, ApJ, 909, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Acernese, F., Agathos, M., Agatsuma, K., et al. 2014, Class. Quant. Grav., 32, 024001 [Google Scholar]
 Adhikari, R. X., Aguiar, O., Arai, K., et al. 2020, Class. Quant. Grav., 37, 165003 [NASA ADS] [CrossRef] [Google Scholar]
 Aksulu, M. D., Wijers, R. A. M. J., van Eerten, H. J., & van der Horst, A. J. 2020, MNRAS, 497, 4672 [NASA ADS] [CrossRef] [Google Scholar]
 Aksulu, M. D., Wijers, R. A. M. J., van Eerten, H. J., & van der Horst, A. J. 2022, MNRAS, 511, 2848 [NASA ADS] [CrossRef] [Google Scholar]
 Alexander, K. D., Berger, E., Fong, W., et al. 2017, ApJ, 848, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Anand, S., Coughlin, M. W., Kasliwal, M. M., et al. 2021, Nat. Astron., 5, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Antier, S., Agayeva, S., Almualla, M., et al. 2020, MNRAS, 497, 5518 [NASA ADS] [CrossRef] [Google Scholar]
 Ascenzi, S., Lillo, N. D., Haster, C.J., Ohme, F., & Pannarale, F. 2019, ApJ, 877, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Ashkar, H., Brun, F., Füßling, M., et al. 2021, J. Cosmol. Astropart. Phys., 2021, 045 [CrossRef] [Google Scholar]
 Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, ApJS, 241, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Babak, S., Taracchini, A., & Buonanno, A. 2017, Phys. Rev. D, 95, 024010 [NASA ADS] [CrossRef] [Google Scholar]
 Barbieri, C., Salafia, O. S., Perego, A., Colpi, M., & Ghirlanda, G. 2019, A&A, 625, A152 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barrett, J. W., Gaebel, S. M., Neijssel, C. J., et al. 2018, MNRAS, 477, 4685 [Google Scholar]
 Beniamini, P., & van der Horst, A. J. 2017, MNRAS, 472, 3161 [NASA ADS] [CrossRef] [Google Scholar]
 Beniamini, P., Nava, L., Duran, R. B., & Piran, T. 2015, MNRAS, 454, 1073 [NASA ADS] [CrossRef] [Google Scholar]
 Beniamini, P., Nava, L., & Piran, T. 2016, MNRAS, 461, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Boersma, O. M., Leeuwen, J. V., Adams, E. A. K., et al. 2021, A&A, 650, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Braun, R., Bonaldi, A., Bourke, T., Keane, E., & Wagg, J. 2019, ApJ, submitted [arXiv:1912.12699] [Google Scholar]
 Breschi, M., Gamba, R., & Bernuzzi, S. 2021a, Phys. Rev. D, 104, 042001 [NASA ADS] [CrossRef] [Google Scholar]
 Breschi, M., Perego, A., Bernuzzi, S., et al. 2021b, MNRAS, 505, 1661 [NASA ADS] [CrossRef] [Google Scholar]
 Broekgaarden, F. S., Berger, E., Neijssel, C. J., et al. 2021, MNRAS, 508, 5028 [NASA ADS] [CrossRef] [Google Scholar]
 Burgay, M., D’Amico, N., Possenti, A., et al. 2003, Nature, 426, 531 [NASA ADS] [CrossRef] [Google Scholar]
 Capano, C. D., Tews, I., Brown, S. M., et al. 2020, Nat. Astron., 4, 625 [NASA ADS] [CrossRef] [Google Scholar]
 Carson, Z., Steiner, A. W., & Yagi, K. 2019, Phys. Rev. D, 99, 043010 [NASA ADS] [CrossRef] [Google Scholar]
 Chornock, R., Berger, E., Kasen, D., et al. 2017, ApJ, 848, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Coughlin, M. W., Dietrich, T., Doctor, Z., et al. 2018, MNRAS, 480, 3871 [NASA ADS] [CrossRef] [Google Scholar]
 Coughlin, M. W., Dietrich, T., Margalit, B., & Metzger, B. D. 2019, MNRAS, 489, L91 [NASA ADS] [CrossRef] [Google Scholar]
 Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Science, 358, 1556 [NASA ADS] [CrossRef] [Google Scholar]
 Cutler, C., & Flanagan, E. E. 1994, Phys. Rev. D, 49, 2658 [Google Scholar]
 D’Emilio, V., Green, R., & Raymond, V. 2021, MNRAS, 508, 2090 [CrossRef] [Google Scholar]
 Dietrich, T., Coughlin, M. W., Pang, P. T. H., et al. 2020, Science, 370, 1450 [Google Scholar]
 Dobie, D., Murphy, T., Kaplan, D. L., et al. 2021, MNRAS, 505, 2647 [NASA ADS] [CrossRef] [Google Scholar]
 Duque, R., Daigne, F., & Mochkovitch, R. 2019, A&A, 631, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eerten, H. V., Horst, A. V. D., & MacFadyen, A. 2012, ApJ, 749, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Eichler, D., Livio, M., Piran, T., & Schramm, D. N. 1989, Nature, 340, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Fernández, R., Foucart, F., & Lippuner, J. 2020, MNRAS, 497, 3221 [CrossRef] [Google Scholar]
 Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
 Finn, L. S., & Chernoff, D. F. 1993, Phys. Rev. D, 47, 2198 [NASA ADS] [CrossRef] [Google Scholar]
 Foucart, F., Hinderer, T., & Nissanke, S. 2018, Phys. Rev. D, 98, 081501 [NASA ADS] [CrossRef] [Google Scholar]
 Fragione, G. 2021, ApJ, 923, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Fragione, G., & Loeb, A. 2021, MNRAS, 503, 2861 [NASA ADS] [CrossRef] [Google Scholar]
 GCN Archive for S200105ae 2020, https://gcn.gsfc.nasa.gov/other/S200105ae.gcn3 [Google Scholar]
 GCN Archive for S200115j 2020, https://gcn.gsfc.nasa.gov/other/S200115j.gcn3 [Google Scholar]
 Gompertz, B. P., Levan, A. J., & Tanvir, N. R. 2020, ApJ, 895, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Haggard, D., Nynka, M., Ruan, J. J., et al. 2017, ApJ, 848, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Hallinan, G., Corsi, A., Mooley, K. P., et al. 2017, Science, 358, 1579 [NASA ADS] [CrossRef] [Google Scholar]
 Hannam, M., Schmidt, P., Bohé, A., et al. 2014, Phys. Rev. Lett., 113, 151101 [Google Scholar]
 Hild, S., Abernathy, M., Acernese, F., et al. 2011, Class. Quant. Grav., 28, 094013 [Google Scholar]
 Hinderer, T., Nissanke, S., Foucart, F., et al. 2019, Phys. Rev. D, 100, 063021 [NASA ADS] [CrossRef] [Google Scholar]
 Hinshaw, G., Larson, D., Komatsu, E., et al. 2013, ApJS, 208, 19 [Google Scholar]
 Hotokezaka, K., Nissanke, S., Hallinan, G., et al. 2016, ApJ, 831, 190 [NASA ADS] [CrossRef] [Google Scholar]
 Hotokezaka, K., Nakar, E., Gottlieb, O., et al. 2019, Nat. Astron., 3, 940 [NASA ADS] [CrossRef] [Google Scholar]
 Husa, S., Khan, S., Hannam, M., et al. 2016, Phys. Rev. D, 93, 044006 [NASA ADS] [CrossRef] [Google Scholar]
 JiménezForteza, X., Keitel, D., Husa, S., et al. 2017, Phys. Rev. D, 95, 064024 [CrossRef] [Google Scholar]
 Just, O., Obergaulinger, M., Janka, H.T., Bauswein, A., & Schwarz, N. 2016, ApJ, 816, L30 [NASA ADS] [CrossRef] [Google Scholar]
 Kasliwal, M. M., Anand, S., Ahumada, T., et al. 2020, ApJ, 905, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Kelley, L. Z. 2021, J. Open Source Softw., 6, 2784 [NASA ADS] [CrossRef] [Google Scholar]
 Khan, S., Husa, S., Hannam, M., et al. 2016, Phys. Rev. D, 93 [Google Scholar]
 Krüger, C. J., & Foucart, F. 2020, Phys. Rev. D, 101, 103002 [CrossRef] [Google Scholar]
 Kyutoku, K., Ioka, K., & Shibata, M. 2013, Phys. Rev. D, 88, 041503 [NASA ADS] [CrossRef] [Google Scholar]
 Kyutoku, K., Shibata, M., & Taniguchi, K. 2021, Liv. Rev. Relat., 24, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Lattimer, J. M., & Prakash, M. 2001, ApJ, 550, 426 [NASA ADS] [CrossRef] [Google Scholar]
 Lazzati, D., Perna, R., Ciolfi, R., et al. 2021, ApJ, 918, L6 [NASA ADS] [CrossRef] [Google Scholar]
 Liebling, S. L., & Palenzuela, C. 2016, Phys. Rev. D, 94, 064046 [NASA ADS] [CrossRef] [Google Scholar]
 LIGO Scientific Collaboration and Virgo Collaboration (Abbott, B., et al.) 2017, Phys. Rev. Lett., 119, 161101 [CrossRef] [Google Scholar]
 Maggiore, M., Broeck, C. V. D., Bartolo, N., et al. 2020, J. Cosmol. Astropart. Phys., 2020, 050 [CrossRef] [Google Scholar]
 Matas, A., Dietrich, T., Buonanno, A., et al. 2020, Phys. Rev. D, 102, 043023 [NASA ADS] [CrossRef] [Google Scholar]
 Mészáros, P., & Rees, M. J. 1992, MNRAS, 257, 29P [CrossRef] [Google Scholar]
 McKinney, J. C. 2005, ApJ, 630, L5 [NASA ADS] [CrossRef] [Google Scholar]
 McKinney, J. C., & Gammie, C. F. 2004, ApJ, 611, 977 [Google Scholar]
 Metzger, B. D., & Berger, E. 2012, ApJ, 746, 48 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24 [Google Scholar]
 Mooley, K. P., Frail, D. A., Ofek, E. O., et al. 2013, ApJ, 768, 165 [NASA ADS] [CrossRef] [Google Scholar]
 Mooley, K. P., Deller, A. T., Gottlieb, O., et al. 2018, Nature, 561, 355 [Google Scholar]
 Nagakura, H., Hotokezaka, K., Sekiguchi, Y., Shibata, M., & Ioka, K. 2014, ApJ, 784, L28 [CrossRef] [Google Scholar]
 Nakar, E., & Piran, T. 2011, Nature, 478, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Nakar, E., & Piran, T. 2021, ApJ, 909, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Nakar, E., Piran, T., & Granot, J. 2002, ApJ, 579, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Neijssel, C. J., VignaGómez, A., Stevenson, S., et al. 2019, MNRAS, 490, 3740 [Google Scholar]
 Nicholl, M., Margalit, B., Schmidt, P., et al. 2021, MNRAS, 505, 3016 [NASA ADS] [CrossRef] [Google Scholar]
 Nissanke, S., Holz, D. E., Hughes, S. A., Dalal, N., & Sievers, J. L. 2010, ApJ, 725, 496 [Google Scholar]
 Pan, Y., Buonanno, A., Taracchini, A., et al. 2014, Phys. Rev. D, 89, 084006 [NASA ADS] [CrossRef] [Google Scholar]
 Pannarale, F., & Ohme, F. 2014, ApJ, 791, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Paschalidis, V., Ruiz, M., & Shapiro, S. L. 2015, ApJ, 806, L14 [NASA ADS] [CrossRef] [Google Scholar]
 Paterson, K., Lundquist, M. J., Rastinejad, J. C., et al. 2021, ApJ, 912, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Raaijmakers, G., Greif, S. K., Riley, T. E., et al. 2020, ApJ, 893, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Raaijmakers, G., Nissanke, S., Foucart, F., et al. 2021a, ApJ, 922, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021b, ApJ, 918, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Radice, D., & Dai, L. 2019, Eur. Phys. J. A, 55, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Radice, D., Perego, A., Zappa, F., & Bernuzzi, S. 2018, ApJ, 852, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Ridnaia, A., Svinkin, D., & Frederiks, D. 2020, J. Phys.: Conf. Ser., 1697, 012030 [NASA ADS] [CrossRef] [Google Scholar]
 Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, ApJ, 918, L27 [CrossRef] [Google Scholar]
 RomeroShaw, I. M., Talbot, C., Biscoveanu, S., et al. 2020, MNRAS, 499, 3295 [NASA ADS] [CrossRef] [Google Scholar]
 Ryan, G., van Eerten, H., Piro, L., & Troja, E. 2020, ApJ, 896, 166 [CrossRef] [Google Scholar]
 Salafia, O. S., & Giacomazzo, B. 2021, A&A, 645, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Salafia, O. S., & Giacomazzo, B. 2022, A&A, 660, C1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Somiya, K. 2012, Class. Quant. Grav., 29, 124007 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, N., Loeb, A., & Berger, E. 2013, Phys. Rev. D, 87, 084053 [NASA ADS] [CrossRef] [Google Scholar]
 The KAGRA Collaboration (Aso, Y., et al.) 2013, Phys. Rev. D, 88, 043007 [NASA ADS] [CrossRef] [Google Scholar]
 The LIGO Scientific Collaboration and The Virgo Collaboration (Abbott, B., et al.) 2018, Phys. Rev. Lett., 121, 161101 [NASA ADS] [CrossRef] [Google Scholar]
 The LIGO Scientific Collaboration and The Virgo Collaboration (Abbott, B. P., et al.) 2017, Nature, 551, 85 [NASA ADS] [CrossRef] [Google Scholar]
 The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration (Abbott, R., et al.) 2021a, ArXiv eprints [arXiv:2111.03606] [Google Scholar]
 The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration (Abbott, R., et al.) 2021b, ArXiv eprints [arXiv:2112.06861] [Google Scholar]
 The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration (Abbott, R., et al.) 2021c, ArXiv eprints [arXiv:2111.03604] [Google Scholar]
 Yagi, K., & Yunes, N. 2017, Phys. Rep., 681, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Zappa, F., Bernuzzi, S., Pannarale, F., Mapelli, M., & Giacobbo, N. 2019, Phys. Rev. Lett., 123, 041102 [NASA ADS] [CrossRef] [Google Scholar]
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 SKA1Mid 5σ_{rms} detection limit. Clearly, most of the sGRB afterglows are below the detection limit even with a radio telescope as sensitive as SKA1Mid. As mentioned in Sect. 3, any change in the sensitivity of SKA1Mid 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 n_{0} and ϵ_{B}. 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 n_{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 ϵ_{B}, this selection effect is less prominent, with means ∼1.5 higher in the detected distributions. While the intrinsic distributions for n_{0} and ϵ_{B} are similar, it is truncated at 10^{−2} for ϵ_{B} leading to lower means in the detected distributions.
Fig. A.1. Flux distribution in log_{10} of sGRB afterglows associated with BHNS mergers detected in GWs with a 3G detector network that includes ET. The distributions are shown for different average BH spins: χ_{BH} = 0 (purple), 0.2 (red), 0.4 (green), 0.6 (orange), and 0.8 (blue). The vertical dashed line shows the SKA1Mid 5σ_{rms} detection limit. 
Fig. A.2. Normalised GW and radio detection rates for the ET and χ_{BH} = {0.4, 0.8} combinations (left and middle panel) and for the aLIGO and χ_{BH} = 0.8 combination (right panel). The rates are normalised to the combined detection rates of Fig. 3, given by (from left to right): 3.3 yr^{−1}, 44 yr^{−1}, and 1.7 yr^{−1}. A hard NS EOS (R_{NS} = 13.0 km) is used in all cases. The rates are shown as a function of n_{0} and ϵ_{B} of the detected sources. 
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 n_{0} and ϵ_{B} at distances d_{L} = 50, 100 Mpc. Both n_{0} and ϵ_{B} have varying degrees of correlation with the binary source parameters. We hypothesise that in the setup used in this work, improved estimates on n_{0} and ϵ_{B}, 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 n_{0} and ϵ_{B}.
Fig. B.1. Posterior distributions of the binary source parameters together with n_{0} and ϵ_{B}. Purple contours give the inferred posterior distribution from the combination of the GW and radio data when the binary is placed at d_{L} = 100 Mpc. Orange contours give the inferred posterior distribution from the combination of the GW and radio data when the binary is placed at d_{L} = 50 Mpc. The fiducial injected values are indicated by the dashed black lines. 
All Tables
Parameters of our fiducial BHNS merger and the chosen prior type and range for the parameter inference.
All Figures
Fig. 1. Visualisation of the fitting formula for the remnant mass in Foucart et al. (2018), showing the dependence on the BH spin (χ_{BH}), mass ratio (Q), and NS tidal deformability (Λ_{NS}). Bestfit values for the free parameters are assumed. In each of the three panels, χ_{BH}, Q, 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 (M_{rem}) in solar masses. 

In the text 
Fig. 2. Same as Fig. 1 except now the final kinetic energy of the jet in Eq. (10) is visualised. The regions of the parameter space with no colour indicate either that M_{disk} = 0 or that the minimal accreted disk mass for sGRB creation, M_{acc} = 0.03 M_{⊙}, is not reached. 

In the text 
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 SKA1Mid 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 (R_{NS} = 11.5 km) or a hard NS EOS (R_{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. 

In the text 
Fig. 4. Normalised GW and radio detection rates for the ET and χ_{BH} = {0.4, 0.8} combinations (left and middle panel) and for the aLIGO and χ_{BH} = 0.8 combination (right panel). The rates are normalised to the combined detection rates of Fig. 3, given by (from left to right): 3.3 yr^{−1}, 44 yr^{−1}, and 1.7 yr^{−1}. A hard NS EOS (R_{NS} = 13.0 km) is used in all cases. The rates are shown as a function of the distance and inclination angle of the detected sources. 

In the text 
Fig. 5. sGRB afterglow radio light curves, computed with afterglowpy, of our fiducial binary at the nominal SKA1Mid frequency, ν = 1.43 GHz. The curves are plotted on a log–log scale. Twenty simulated observations were taken between 11 and 500 days postmerger. The orange circles indicate the 20 observations above the 3σ_{rms} detection threshold when the binary is placed at d_{L} = 50 Mpc. The purple circles indicate the 16 observations above the 3σ_{rms} detection threshold when the binary is placed at d_{L} = 100 Mpc. The red circles indicate the seven observations above the 3σ_{rms} detection threshold when the binary is placed at d_{L} = 200 Mpc. The solid black line shows the SKA1Mid 3σ_{rms} detection threshold with σ_{rms} = 2 μJy. The error bars give the total error as a combination of the measurement error and the systematic error. 

In the text 
Fig. 6. Posterior distributions of the parameters θ_{GW} as inferred for our fiducial binary at d_{L} = 100 Mpc (left panel) and d_{L} = 50 Mpc (right 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. The dark and light shaded regions give the 68% and 95% confidence levels, respectively, for all posteriors shown. The fiducial injected values are indicated by the dashed black lines. 

In the text 
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 d_{L} = 100 Mpc. Orange contours give the inferred posterior distribution from the combination of the GW and radio data when the binary is placed at d_{L} = 50 Mpc. The fiducial injected values are indicated by the dashed black lines. 

In the text 
Fig. 8. Posterior distributions of the parameters θ_{GW} (left panel) and θ_{EM} (right panel) as inferred for our fiducial binary at d_{L} = 200 Mpc. 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. 

In the text 
Fig. 9. Resulting fit to the light curve of our fiducial binary at d_{L} = 50 Mpc. The 95% credible region of the fit, indicated by the shaded region, is computed from the posterior samples of the combined GW and EM data inference in Sect. 5. The black circles give the 20 simulated observations, and the error bars correspond to the SKA1Mid measurement errors. 

In the text 
Fig. A.1. Flux distribution in log_{10} of sGRB afterglows associated with BHNS mergers detected in GWs with a 3G detector network that includes ET. The distributions are shown for different average BH spins: χ_{BH} = 0 (purple), 0.2 (red), 0.4 (green), 0.6 (orange), and 0.8 (blue). The vertical dashed line shows the SKA1Mid 5σ_{rms} detection limit. 

In the text 
Fig. A.2. Normalised GW and radio detection rates for the ET and χ_{BH} = {0.4, 0.8} combinations (left and middle panel) and for the aLIGO and χ_{BH} = 0.8 combination (right panel). The rates are normalised to the combined detection rates of Fig. 3, given by (from left to right): 3.3 yr^{−1}, 44 yr^{−1}, and 1.7 yr^{−1}. A hard NS EOS (R_{NS} = 13.0 km) is used in all cases. The rates are shown as a function of n_{0} and ϵ_{B} of the detected sources. 

In the text 
Fig. B.1. Posterior distributions of the binary source parameters together with n_{0} and ϵ_{B}. Purple contours give the inferred posterior distribution from the combination of the GW and radio data when the binary is placed at d_{L} = 100 Mpc. Orange contours give the inferred posterior distribution from the combination of the GW and radio data when the binary is placed at d_{L} = 50 Mpc. The fiducial injected values are indicated by the dashed black lines. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.