© O. M. Boersma and J. van Leeuwen 2022

1. 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 gravitational-wave (GW) observations (Abbott et al. 2021a). 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 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¹. 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 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. 2021; 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 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 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 multi-messenger 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 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.

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 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 ():

(1)

with best-fit values α = 0.406, β = 0.139, γ = 0.255, and δ = 1.761. Here, η = Q/(1 + Q)² and R̂_ISCO = R_ISCO/M_BH is the innermost stable circular orbit radius normalised by the BH mass:

(2)

(3)

(4)

To compute the compactness of the NS C_NS, we used the approximately universal C-Love relation for NSs of Eq. (78) in Yagi & Yunes (2017):

(5)

with best-fit coefficients a₀ = 0.360, a₁ = −0.0355 and a₂ = 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):

(6)

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). Best-fit 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 (Mrem) 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):

(7)

We again turned to fits of NR simulations (Krüger & Foucart 2020) for M_dyn:

(8)

with best-fit coefficients a₁ = 0.007116, a₂ = 0.001436, a₄ = −0.02762, n₁ = 0.8636, and n₂ = 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):

(9)

Here we assumed average values for the free parameters of ξ₁ = 0.18 and ξ₂ = 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 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 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 gamma-ray 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:

(10)

where the factor one-half accounts for an identical counter-jet. McKinney & Gammie (2004) describe general-relativistic magneto-hydrodynamical 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 accretion-to-jet energy conversion efficiency:

(11)

where is the spin of the final remnant BH and Ω_H is the dimensionless angular frequency at the BH horizon:

(12)

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énez-Forteza et al. 2017). Compared to the formulae stated in SG21, in the regime an additional factor of 10⁻² 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 Mdisk = 0 or that the minimal accreted disk mass for sGRB creation, Macc = 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 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 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 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; 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 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 gamma-ray-burst 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):

(13)

where θ_c is the opening angle of the jet core, n₀ 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 on-axis isotropic-equivalent jet energy E₀. We used a standard formula for top-hat jets to convert between the two: E₀ = 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 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 p = 2.2 and ϵ_e = 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 n₀ 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 log-normal distribution for both parameters with mean μ = 10⁻³ and standard deviation σ = 0.75. Furthermore, ϵ_B was constrained to the range [10⁻⁴, 10⁻²]. 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 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 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 (RNS = 11.5 km) or a hard NS EOS (RNS = 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 (RNS = 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 SKA1-Mid 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 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 (ι ≤ θ_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 SKA1-Mid 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 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). E₀ 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 log-linear 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 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.

4. 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.

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., Romero-Shaw 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:

(14)

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:

(15)

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ℒ_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 log-likelihood of this network takes the standard form (see e.g., Cutler & Flanagan 1994):

(16)

where the inner product for two GW strains a(t) and b(t) is defined as

(17)

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 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 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³. The chirp mass is defined as

(18)

(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 built-in 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 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, t_obs; (ii) the core opening angle of the jet, θ_c, which determines the effective width of the jet core; (iii) the on-axis isotropic-equivalent jet energy, E₀; (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₀. We treated n₀ 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:

(19)

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₀, 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₀ 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

(20)

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₀, 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.

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 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 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 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.

Fig. 5. sGRB afterglow radio light curves, computed with afterglowpy, of our fiducial binary at the nominal SKA1-Mid frequency, ν = 1.43 GHz. The curves are plotted on a log–log scale. Twenty simulated observations were taken between 11 and 500 days post-merger. The orange circles indicate the 20 observations above the 3σrms detection threshold when the binary is placed at dL = 50 Mpc. The purple circles indicate the 16 observations above the 3σrms detection threshold when the binary is placed at dL = 100 Mpc. The red circles indicate the seven observations above the 3σrms detection threshold when the binary is placed at dL = 200 Mpc. The solid black line shows the SKA1-Mid 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₀ and ϵ_B 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 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 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 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 dL = 100 Mpc (left panel) and dL = 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 dL = 100 Mpc. Orange contours give the inferred posterior distribution from the combination of the GW and radio data when the binary is placed at dL = 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 SKA1-Mid 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₀ 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₀ 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 dL = 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₀ are the most strongly correlated with the binary source parameters. We show the combined posterior of the binary source parameters together with n₀ 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₀ 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₀ 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 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 d_L = 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.

Fig. 9. Resulting fit to the light curve of our fiducial binary at dL = 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 SKA1-Mid measurement errors.

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 ϵ_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 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 accretion-to-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 (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 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.

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 accretion-to-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 well-improved 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.

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 NWA-ORC, financed by NWO; and from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 617199 (‘ALERT’).

References

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 n₀ 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₀, 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₀ and ϵ_B are similar, it is truncated at 10⁻² for ϵ_B leading to lower means in the detected distributions.

Fig. A.1. Flux distribution in log10 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 SKA1-Mid 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 (RNS = 13.0 km) is used in all cases. The rates are shown as a function of n0 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₀ and ϵ_B at distances d_L = 50,  100 Mpc. Both n₀ 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₀ 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₀ and ϵ_B.

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

All Tables

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). Best-fit 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 (Mrem) 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 Mdisk = 0 or that the minimal accreted disk mass for sGRB creation, Macc = 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 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 (RNS = 11.5 km) or a hard NS EOS (RNS = 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 (RNS = 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 SKA1-Mid frequency, ν = 1.43 GHz. The curves are plotted on a log–log scale. Twenty simulated observations were taken between 11 and 500 days post-merger. The orange circles indicate the 20 observations above the 3σrms detection threshold when the binary is placed at dL = 50 Mpc. The purple circles indicate the 16 observations above the 3σrms detection threshold when the binary is placed at dL = 100 Mpc. The red circles indicate the seven observations above the 3σrms detection threshold when the binary is placed at dL = 200 Mpc. The solid black line shows the SKA1-Mid 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 dL = 100 Mpc (left panel) and dL = 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 dL = 100 Mpc. Orange contours give the inferred posterior distribution from the combination of the GW and radio data when the binary is placed at dL = 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 dL = 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 dL = 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 SKA1-Mid measurement errors.
In the text

	Fig. A.1. Flux distribution in log10 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 SKA1-Mid 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 (RNS = 13.0 km) is used in all cases. The rates are shown as a function of n0 and ϵB of the detected sources.

In the text

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

In the text