<figure><img src="image_1.png" /></figure>

HAL
open science

# Mergers of black hole-neutron star binaries and rates of associated electromagnetic counterparts

Mukul Bhattacharya, Pawan Kumar, George Smoot

► To cite this version:

Mukul Bhattacharya, Pawan Kumar, George Smoot. Mergers of black hole-neutron star binaries and rates of associated electromagnetic counterparts. Monthly Notices of the Royal Astronomical Society, 2019, 486 (4), pp.5289-5309. 10.1093/mnras/stz1147. hal-01880713

HAL Id: hal-01880713

https://hal.science/hal-01880713

Submitted on 29 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Mergers of black hole–neutron star binaries and rates of associated electromagnetic counterparts

Mukul Bhattacharya<sup>1,★</sup> Pawan Kumar<sup>2</sup> and George Smoot<sup>3,4,5</sup>

Accepted 2019 April 25. Received 2019 April 25; in original form 2018 August 31

## ABSTRACT

Black hole–neutron star (BHNS) binaries are amongst promising candidates for the joint detection of electromagnetic (EM) signals with gravitational waves (GWs) and are expected to be detected in the near future. Here we study the effect of the BHNS binary parameters on the merger ejecta properties and associated EM signals. We estimate the remnant disc and unbound ejecta masses for BH mass and spin distributions motivated from the observations of transient low-mass X-ray binaries and a specific NS equation of state (EoS). The amount of r-process elements synthesized in BHNS mergers is estimated to be a factor of ~$10^2–10^4$ smaller than BNS mergers, due to the smaller dynamical ejecta and merger rates for the former. We compute the EM luminosities and light curves for the early- and late-time emissions from the ultrarelativistic jet, sub-relativistic dynamical ejecta and wind, and the mildly relativistic cocoon for typical ejecta parameters. We then evaluate the low-latency EM follow-up rates of the GW triggers in terms of the GW detection rate $\dot{N}_{GW}$ for current telescope sensitivities and typical BHNS binary parameters to find that most of the EM counterparts are detectable for high BH spin, small BH mass, and a stiffer NS EoS when NS disruption is significant. Based on the relative detection rates for given binary parameters, we find the ease of EM follow-up to be: ejecta afterglow > cocoon afterglow ≳ jet prompt > ejecta macronova > cocoon prompt > jet afterglow >> wind macronova >> wind afterglow.

**Key words:** accretion, accretion discs – equation of state – gravitational waves – gamma-ray burst: general – stars: neutron.

## 1 INTRODUCTION

In addition to being one of the most likely candidates for the progenitor of short gamma-ray bursts (sGRBs; Eichler et al. 1989; Nakar 2007), the coalescence of BNS and black hole–neutron star (BHNS) compact binaries is also considered to be a promising source of gravitational waves (GWs) detectable by the current generation of ground-based interferometers such as Advanced LIGO (aLIGO; LIGO Scientific Collaboration 2010), Advanced Virgo (Acernese et al. 2015), and KAGRA (Somiya 2012). At present, six BBH and one BNS mergers have been detected and BHNS merger detections are anticipated in the upcoming aLIGO/Virgo runs at a rate of more than one per year. The observations from binaries with at least one neutron star will help us constrain the equation of state (EoS) of dense nuclear matter in these compact objects (Takami,

Rezzolla & Baiotti 2015; Clark et al. 2016). Furthermore, the detection of accompanying electromagnetic (EM) counterparts is a potent prospect (Metzger & Berger 2012), as it will improve source localization and redshift estimate, and provide more information about the merger process as well as the properties of the compact binary.

Many numerical relativity simulations of compact binary mergers have been carried out in the recent past to study the merger dynamics and properties of the ejected material as well as the post-merger remnant (Foucart 2012; Foucart et al. 2013, 2014; Kyutoku et al. 2015; Foucart, Hinderer & Nissanke 2018; Kawaguchi et al. 2015, 2016). The outcome of the BHNS merger critically depends on both the binary separation at tidal disruption of the NS ($R_{tid}$) and the location of the innermost stable circular orbit (ISCO) of the BH ($R_{ISCO}$). The most important parameters that determine the BHNS merger dynamics and ejecta properties are the mass ratio of the binary system (Foucart et al. 2012; Kyutoku et al.), the BH spin magnitude and orientation (Pannarale, Tonita & Rezzolla 2011;

<footnote n="1"><sup>1</sup>Department of Physics, University of Texas at Austin, Austin, TX 78712, USA</footnote>

<footnote n="2"><sup>2</sup>Department of Astronomy, University of Texas at Austin, Austin, TX 78712, USA</footnote>

<footnote n="3"><sup>3</sup>Institute for Advanced Study, Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, 999077 Hong Kong, China</footnote>

<footnote n="4"><sup>4</sup>PCCP; APC, Université Paris Diderot, Université Sorbonne Paris Citè, 75013 Paris, France</footnote>

<footnote n="5"><sup>5</sup>BCCP; LBNL & Physics Dept. University of California at Berkeley, Berkeley, CA 94720, USA</footnote>

<footnote>* E-mail: mukul.b@utexas.edu</footnote>
Foucart 2012), and the EoS of the NS (Duez et al. 2010; Kyutoku, Shibata & Taniguchi 2010).

If the NS is tidally disrupted, most of the infalling matter is accreted by the BH; however, some material remains outside for longer time-scales in the form of a hot accretion disc, eccentric bound orbits (tails), and/or unbound (dynamical) ejecta around the BH. A considerable fraction of the remnant accretion disc mass can also get unbound over long time-scales due to winds driven by magnetic fields (Kiuchi et al. 2015), viscous flows (Fernandez & Metzger 2013), and/or neutrinos (Dessart et al. 2009). The ejected material from a BHNS merger can be classified into four different components: a sGRB jet, dynamical ejected mass, magnetic/viscous/neutrino-driven winds, and a cocoon from jet interacting with other ejecta – see Fig. 1 for a schematic sketch and Table 1 for the properties of these ejecta components and their associated emission. The dynamical ejected mass is found to be much more anisotropic for BHNS mergers compared to BNS mergers with most of the matter being concentrated around the orbital plane (Kyutoku et al. 2015). Recent studies (Barbieri et al. 2019; Coughlin & Dietrich 2019) have shown that joint multimessenger analysis of the GW and EM signals, especially those arising from the sGRB jet and the dynamical ejected mass, can break degeneracies in the GW parameter space to provide better constraints on the BHNS binary parameters.

The BNS merger GW170817 is the only known astrophysical source at present with both GW and multiwavelength EM radiation observations (e.g. Abbott et al. 2017a,b; Mooley et al. 2018a). The steadily increasing afterglow luminosity observations (Ghirlanda et al. 2019; Margutti et al. 2018; Mooley et al. 2018a) are better explained with a jet–cocoon system whereby a highly relativistic narrow jet propagates through an environment of denser and mildly relativistic material (Ioka & Nakamura 2018; Lazzati et al. 2018), as opposed to a mildly relativistic isotropic fireball (Haggard et al. 2017) and/or top-hat jet seen off-axis (Alexander et al. 2017). The very long baseline interferometry (VLBI) observations of the rapid turnover around the radio light-curve peak and the very fast subsequent decline suggest that the early-time afterglow is most likely due to the emission from the wider-angle cocoon whereas the late-time emission is powered by a collimated relativistic jet with an opening angle <5° viewed from an angle of ~20° (Ghirlanda et al. 2019; Mooley et al. 2018b).

There can be several possible EM counterparts for BHNS mergers: optical–infrared (IR) transients known as macronova driven by radioactive decay of neutron-rich dynamical ejecta on time-scales of $t \sim$ few hours to days (Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010a), synchrotron radio afterglow from ejecta components decelerated by the surrounding interstellar medium (ISM) with $t \sim$ weeks to years (Nakar & Piran 2011; Metzger & Berger 2012; Piran, Nakar & Rosswog 2013), relativistic sGRB jet prompt emission on short time-scales $t \lesssim 1$ s (Eichler et al. 1989; Nakar 2007), and gamma/X-ray prompt emission from cocoon with $t \lesssim 1$ s (Lazzati et al. 2018). Even though all four ejecta components contribute to synchrotron afterglow emission directly, different ejecta component afterglows peak at different time-scales and with different intensities due to their separate masses and kinetic energies (KEs).

In this paper, we explore the dependence of BHNS merger ejecta properties and associated EM signals on the mass ratio of the binary, the magnitude of the BH spin, and the unknown NS EoS. For simplicity, we will only consider the cases where BH spin is completely aligned with the orbital angular momentum and the

binary orbit is of low eccentricity. We also evaluate the detection rates for low-latency EM follow-up of GW triggers restricted to the nearby Universe, determined by the aLIGO sensitivity at a given frequency. The structure of this paper is as follows. In Section 2, we evaluate the dynamical/unbound ejecta mass, the disc mass, and the dynamical ejecta velocity for a given BH mass, spin distribution, and NS EoS. In Section 3, we calculate the horizon distance from aLIGO sensitivity to find the GW event detection rate using BHNS merger rates from population synthesis models. In Section 4, we discuss all the possible EM counterparts for BHNS mergers and further estimate their rates for EM follow-up of GW triggers in Section 5. Finally, we summarize our results in Section 6.

## 2 BHNS MERGER EJECTA PROPERTIES FROM BINARY PARAMETERS

In this section, we estimate the mass of bound ejecta ($M_{ej,b}$) along with the dynamical ejecta mass ($M_{ej,dyn}$) and velocity ($v_{ej,dyn}$) for BHNS binaries with given BH mass ($M_{BH}$) and spin ($a_{BH}$) distributions, and an NS EoS. We consider different distributions for $M_{BH}$ and $a_{BH}$ in order to find the range of BHNS binary parameters for which the disruption of NS and subsequent mass ejection is more likely. The tidal disruption of NS by a BH is facilitated by a larger NS radius ($R_{NS}$), a smaller $M_{BH}$, and a larger $a_{BH}$. The BHNS merger ejecta calculations are done for a single NS mass $M_{NS}$ and therefore $R_{NS}$ is directly related to the NS EoS (as discussed in Section 2.2). For the rest of this paper, we use units in which $G = c = 1$ and $M_{BH}/M_{NS}$ is the Arnowitt–Deser–Misner (ADM) mass of the BH/NS at infinite separation.

### 2.1 Bound and dynamical ejecta mass

The outcome of any BHNS merger is two-fold: either the NS plunges directly into the BH before it can get tidally disrupted or the tidal forces on the NS are sufficiently strong to disrupt it before reaching the $R_{ISCO}$ of the BH. The binary parameters that essentially determine this outcome are the mass ratio of the compact objects $q = M_{BH}/M_{NS}$, the orientation and magnitude of $a_{BH}$, and $R_{NS}$ based on the EoS of the NS. The mass of the disrupted NS present outside the BH at late times ~10 ms is almost entirely determined by the relative positions of $R_{td}$ and $R_{ISCO}$. While $R_{td} \sim R_{NS}(M_{BH}/M_{NS})^{1/3}$ in Newtonian theory is obtained by equating the self-gravity of the NS with the tidal force from the BH, the ISCO radius is given by

$$ \frac{R_{ISCO}}{M_{BH}} = 3 + Z_2 - \operatorname{sign}(a_{BH}) \sqrt{(3 - Z_1)(3 + Z_1 + 2 Z_2)}, \quad (1) $$

where

$$ Z_1 = 1 + (1 - a_{BH}^2)^{1/3} [(1 + a_{BH})^{1/3} + (1 - a_{BH})^{1/3}] $$

$$ Z_2 = \sqrt{3a_{BH}^2 + Z_1^2}. $$

Here, $a_{BH} = a/M_{BH}$ is the dimensionless spin parameter of the BH.

Using the numerical simulation results obtained from three different evolutionary codes (Etienne et al. 2009; Kyutoku et al. 2011; Foucart et al. 2011, 2012, 2013; Lovelace et al. 2013; Foucart et al. 2014; Kyutoku et al. 2015; Brege et al. 2018), Foucart et al. (2018) proposed a model to determine the combined mass $M_{ej,tot}$ of the remnant disc, the tidal tails, and the unbound ejecta remaining
<figure><img src="image_4.png" /><figcaption>**Figure 1.** A schematic diagram for the various ejecta components from BHNS mergers: *ultrarelativistic* sGRB jet, *mildly relativistic* cocoon, and *sub-relativistic* dynamical ejecta and magnetic/viscous/neutrino-driven winds. The early-time emissions are (optical/IR) macronova, (gamma/X-ray) sGRB jet prompt, and (gamma/X-ray) cocoon prompt while the interaction of the ejecta components with the ambient ISM results in the late-time radio and optical afterglows from the dynamical ejecta, jet, wind, and cocoon.</figcaption></figure>

<table><caption>**Table 1.** Properties of different ejecta components and their associated EM signals.</caption><thead><tr><th>Component</th><th>Mass (M<sub>☉</sub>)</th><th>Average Γ</th><th>KE (erg)</th><th>EM signal</th><th>Band</th><th>L<sub>peak</sub> (erg s<sup>-1</sup>)</th><th>t<sub>peak</sub></th></tr></thead><tbody><tr><td rowspan="2">sGRB jet</td><td rowspan="2">10<sup>-4</sup>-10<sup>-3</sup></td><td rowspan="2">10-30</td><td rowspan="2">10<sup>51</sup>-10<sup>52</sup></td><td>prompt</td><td>gamma/X-ray</td><td>10<sup>48</sup>-10<sup>49</sup></td><td>&lt;second</td></tr><tr><td>afterglow</td><td>optical/radio</td><td>10<sup>36</sup>-10<sup>38</sup></td><td>~day-week</td></tr><tr><td rowspan="2">Cocoon</td><td rowspan="2">10<sup>-5</sup>-10<sup>-3</sup></td><td rowspan="2">3-10</td><td rowspan="2">10<sup>49</sup>-10<sup>51</sup></td><td>prompt</td><td>gamma/X-ray</td><td>10<sup>47</sup>-10<sup>48</sup></td><td>~second-minute</td></tr><tr><td>afterglow</td><td>optical/radio</td><td>10<sup>36</sup>-10<sup>39</sup></td><td>~week-month</td></tr><tr><td rowspan="2">Dynamical ejecta</td><td rowspan="2">10<sup>-4</sup>-10<sup>-1</sup></td><td rowspan="2">~1</td><td rowspan="2">10<sup>48</sup>-10<sup>51</sup></td><td>macronova</td><td>optical/IR</td><td>10<sup>40</sup>-10<sup>42</sup></td><td>~day-week</td></tr><tr><td>afterglow</td><td>optical/radio</td><td>10<sup>35</sup>-10<sup>38</sup></td><td>~year-decade</td></tr><tr><td rowspan="2">Wind</td><td rowspan="2">10<sup>-4</sup>-10<sup>-3</sup></td><td rowspan="2">~1</td><td rowspan="2">10<sup>50</sup>-10<sup>51</sup></td><td>macronova</td><td>optical/IR</td><td>10<sup>39</sup>-10<sup>40</sup></td><td>~hour-day</td></tr><tr><td>afterglow</td><td>optical/radio</td><td>10<sup>32</sup>-10<sup>34</sup></td><td>~year-decade</td></tr></tbody></table>
outside the BH $t \sim 10$ ms after the merger has occurred,

$$ \frac{M_{\mathrm{ej,tot}}}{M_{\mathrm{b,NS}}} = \mathrm{Max} \left[ \left( 0.406 q^{0.333}(1 - 2C_{\mathrm{NS}}) - 0.139 \frac{R_{\mathrm{ISCO}}}{R_{\mathrm{NS}}} + 0.255 \right)^{1.761}, 0 \right], \quad (2) $$

where $C_{NS} = M_{NS}/R_{NS}$ is the NS compactness parameter and $M_{b,NS}$ is the baryon mass of the NS. The error in $M_{ej,tot}$ is estimated to be $[(M_{ej,tot}/10)^2 + (1/100)^2]^{1/2}$. The NS mass is fixed to be $M_{NS} = 1.35 M_{\odot}$ whereas $M_{b,NS}$ depends on the EoS considered for the NS (see Section 2.2). Kawaguchi et al. (2016) used the numerical relativity simulation results of the Kyoto group (Kawaguchi et al. 2015; Kyutoku et al. 2015) to extend the formalism developed by Foucart (2012), and determine the mass and average velocity of the dynamical ejecta from BHNS mergers with the least-squares fit

$$ \frac{M_{\mathrm{ej,dyn}}}{M_{\mathrm{b,NS}}} = \mathrm{Max} \left[ 4.464 \times 10^{-2} q^{0.250} \left( \frac{1 - 2C_{\mathrm{NS}}}{C_{\mathrm{NS}}} \right) - 2.269 \times 10^{-3} q^{1.352} \frac{R_{\mathrm{ISCO}}}{M_{\mathrm{BH}}} + 2.431 \left( 1 - \frac{M_{\mathrm{NS}}}{M_{\mathrm{b,NS}}} \right) - 0.4159, 0 \right], \quad (3) $$

$$ v_{\mathrm{avg,dyn}} = (1.533 \times 10^{-2} q + 1.907 \times 10^{-1}) c. \quad (4) $$

Here, $M_{ej,dyn}$ is also estimated ~10 ms after the merger and does not include the possible ejecta component from the remnant BH accretion disc. Similar to $M_{ej,tot}$, the error in $M_{ej,dyn}$ is estimated as $[(M_{ej,dyn}/10)^2 + (M_{\odot}/200)^2]^{1/2}$. The mass of the bound ejecta component (in the remnant disc and tidal tails) is calculated as $M_{ej,b} = M_{ej,tot} - M_{ej,dyn}$. In general, the majority of the bound ejecta mass is contained in the accretion disc ($M_{disc} \approx M_{ej,b}$). It should also be noted that these estimates are limited to low-eccentricity BHNS binary orbits for BH spins aligned with the orbital angular momentum and in the range of parameters $M_{BH} = 1-7 M_{NS}$, $a_{BH} = -0.5-0.97$, and $C_{NS} = 0.13-0.18$ ($R_{NS} = 11-16$ km).

From equations 2 and 3, the ejecta masses $M_{ej,tot}$ and $M_{ej,dyn}$ are both larger for a smaller $M_{BH}$, a larger $a_{BH}$, and a stiffer NS EoS (smaller $C_{NS}$), which is expected from the physics of tidal disruption. While a smaller $M_{BH}$ means a larger tidal force from the BH on the infalling NS between $R_{td}$ and $R_{ISCO}$, resulting in a larger NS disruption, the increase in $a_{BH}$ reduces $R_{ISCO}$ – thereby causing a larger amount of NS matter to be stripped. For a stiffer NS EoS, the increase in NS tidal deformability is expected to give off larger ejecta material during its disruption.

The ratio of the r-process elements synthesized from BHNS mergers to BNS mergers is approximately the product of the relative ratios of the dynamical ejecta produced and the merger rate. While $M_{ej,dyn} \sim 10^{-4}-10^{-3} M_{\odot}$ for BHNS mergers (Kawaguchi et al. 2016) is about an order of magnitude smaller than that for typical BNS merger events such as GW170817 (Abbott et al. 2017c), the relative merger rate $R_{BHNS}/R_{BNS}$ is a fairly uncertain quantity and ranges from ~$10^{-3}$ to ~$10^{-1}$ with a most probable value of ~$10^{-2}$ (Abadie et al. 2010). Therefore, BHNS mergers are expected to synthesize only about 0.01–1 per cent of the r-process elements found in the Universe.

## 2.2 Distributions for BHNS binary parameters

Here, we consider different sets of binary parameters in order to estimate the ejecta properties for typical BHNS mergers. We perform Monte Carlo simulations to estimate the ejecta masses $M_{ej,tot}$ and $M_{ej,b} \approx M_{disc}$ – see Appendix A for the results. We consider BHNS mergers with a given NS EoS as well as specific $M_{BH}$ and $a_{BH}$ distributions as described below.

(i) **BH mass distribution:** We consider Gaussian, power-law, and exponential $M_{BH}$ distributions inferred from transient low-mass X-ray binary (LMXB) observations. Ozel et al. (2010) compiled the dynamical BH mass measurements in transient LMXBs (McClintock & Remillard 2006; Remillard & McClintock 2006) to find that most of the BH masses were clustered within ~6–10 $M_{\odot}$, whereas there were no observed BHs with $M_{BH} \sim 2–5 M_{\odot}$. The narrow mass distribution of the 16 observed BHs used in their sample was found to be consistent with a Gaussian at $(7.8 \pm 1.2) M_{\odot}$. Farr et al. (2011) then obtained a best-fitting power-law distribution, $P(M_{BH}) \propto M_{BH}^{-6.4}$, for $6.1 M_{\odot} < M_{BH} < 23 M_{\odot}$, with the same data as Ozel et al. (2010). The lower-end cut-off of this distribution accounts for the paucity of BHs within the mass range ~2–5 $M_{\odot}$. The BH masses in transient LMXBs have also been modelled using an exponentially decaying distribution with a sharp cut-off (Ozel et al. 2010; Steiner, McClintock & Reid 2012), $P(M_{BH}) = M_{\mathrm{scale}}^{-1} \exp[-(M_{BH} - M_c)/M_{\mathrm{scale}}]$ for $M_{BH} > M_c$, where $M_c = 6.30 M_{\odot}$ and $M_{\mathrm{scale}} = 1.57 M_{\odot}$. This parametric mass distribution is motivated by the mass distributions of pre-SNe stars as well as the energetics of SNe explosions. While we also consider $M_{BH}$ distributions with peak around 10–25 $M_{\odot}$ for low-metallicity simulated BHNS systems (Kruckow et al. 2018), such large BH masses do not produce almost any post-merger ejecta material (both bound and dynamical ejecta from equations 2 and 3) for the associated EM counterparts to be detected.

(ii) **BH spin distribution:** The magnitude of the BH spin in a BHNS binary can be modified by mass transfer via accretion from its NS companion either through an accretion disc or a common-envelope phase (O'Shaughnessy et al. 2008). A significant fraction of the intrinsic BH spin is expected to be attained during its formation via the collapse of the progenitor star. At later times, the BH can further spin up due to repeated accretion episodes. At present, we do not have any direct observations of the BH spin, which makes it very difficult to confidently constrain the BH spins in merging BHNS systems. Estimates of $a_{BH}$ from quasi-periodic oscillations and/or iron line profiles are largely dependent on how the disc accretion is modelled for XRBs (Reynolds et al. 1999; Miller et al. 2004; Narayan 2005; O'Shaughnessy et al. 2005; Torok et al. 2005) and can vary within a large range based on the assumptions used. Due to these uncertainties in the determination of $a_{BH}$, here we only consider two simplistic distributions for the BH spin: (a) $a_{BH} = \mathrm{Uniform}(0, 0.97)$, and (b) $a_{BH} = A$, where $0 \le A \le 0.97$ is a constant. Note that the upper limit of $a_{BH} = 0.97$ is assumed as the ejecta mass estimate fits are only valid up to that spin magnitude.

(iii) **NS equation of state:** Along with the mass and spin magnitude of the BH, the size of the NS specified by its EoS plays an important role in determining the amount of matter ejected in any BHNS merger event. The disruption of the NS by the BH results in larger ejected mass for a stiffer NS EoS, that is a larger NS radius. The NS parameters $M_{b,NS}$, $R_{NS}$, and $C_{NS} = M_{NS}/R_{NS}$ are closely associated with the NS EoS adopted. In this work, we consider two piecewise polytropic EoSs (Read et al. 2009) that are commonly used in the numerical relativity simulations:
(a) APR4 EoS model: soft NS EoS with $R_{NS} = 11.1$ km, $M_{b,NS} = 1.50 M_\odot$, and $C_{NS} = M_{NS}/R_{NS} = 0.180$ for $M_{NS} = 1.35 M_\odot$.

(b) H4 EoS model: stiff NS EoS with $R_{NS} = 13.6$ km, $M_{b,NS} = 1.47 M_\odot$, and $C_{NS} = M_{NS}/R_{NS} = 0.147$ for $M_{NS} = 1.35 M_\odot$.

### 3 GW EVENT DETECTION RATE FROM LIGO SENSITIVITY

In this section, we first calculate the detection horizon distance ($D_{hor}$) for GW observations based on the aLIGO strain sensitivity curve and signal-to-noise ratio (SNR). We further evaluate the event detection rate of BHNS mergers within the observable volume using the plausible estimates for the BHNS coalescence rates from population synthesis models. The design strain sensitivity of the initial LIGO detectors was improved by a factor of ~10 for aLIGO in the entire operating frequency range ~10 Hz to 10 kHz (LIGO Scientific Collaboration 2010), thereby increasing the observable volume of the Universe significantly by ~1000 times for the first observing run (O1) in 2015 September. The source orientation and sky location averaged range for observing a binary coalescence event with aLIGO design sensitivity and at an SNR of 8 will be ~300 Mpc for ~1.4 M_\odot BNS inspirals, ~650 Mpc for ~1.4–8 M_\odot BHNS inspirals, and ~1.6 Gpc for ~10 M_\odot BBH inspirals.

Here we consider a well-separated BHNS binary with point mass constituents moving in a circular orbit to estimate the GW strain amplitude at a distance $r$ from the source. Even though the orbit of the BHNS system is elliptical initially, it can be shown that it gets quickly circularized due to GW emission (Peters 1964). The lowest order quadrupolar contribution gives the GW strain amplitude (Postnov & Yungelson 2014)

$$ 
\begin{aligned}
h_{GW} &= \left(\frac{32}{5}\right)^{1/2} \frac{(GM_{ch})^{5/3}}{rc^4} (\pi f_{GW})^{2/3} \\
&= (7.56 \times 10^{-23}) M_{\text{tot,M}\odot}^{2/3} \mu_{\text{red,M}\odot} f_{\text{GW,Hz}}^{2/3} r_{\text{Mpc}}^{-1}
\end{aligned} 
\quad (5) $$

where $M_{ch} = \mu_{\text{red}}^{3/5} M_{\text{tot}}^{2/5}$ is the chirp mass, $\mu_{\text{red}} = \mu_{\text{red,M}\odot} \times M_\odot$ is the reduced mass, $M_{\text{tot}} = M_{\text{tot,M}\odot} \times M_\odot$ is the total mass of the binary, $f_{GW} = f_{GW,Hz} \times 1$ Hz = $(1/\pi)\sqrt{GM_{\text{tot}}/a^3}$ is the GW frequency that is twice the orbital frequency, $a$ is the orbital separation, and $r_{\text{Mpc}} = r/(1 \text{ Mpc})$. The GW frequency is maximum, $f_{GW,max} \sim f_{\text{ISCO}} = (1/\pi)\sqrt{GM_{\text{tot}}/R_{\text{ISCO}}^3}$, at $a=R_{\text{ISCO}}$. It should be noted that we have excluded geometrical factors of order unity in the expression of $h_{GW}$ and equation (5) corresponds to the case where the observer line of sight is along the orbital angular momentum axis of the binary.

We use the aLIGO design sensitivity curve (LIGO Scientific Collaboration 2010) to estimate the GW strain amplitude spectral density (ASD) as a function of the detection frequency. As aLIGO is expected to reach its design sensitivity (lowest achievable noise level with current detector capabilities) by 2019, we fitted the design strain ASD in the entire detector frequency range $f \sim 10$ Hz to 10 kHz to obtain

$$ \sqrt{S_n(f)} = 2.6 \times 10^{-23} f^{-0.44} + 1.3 \times 10^{-19} f^{-3.06} + 3.8 \times 10^{-27} f, \quad (6) $$

which is better by a factor of ~3–4 in the most sensitive frequency band ~100–300 Hz compared to the aLIGO sensitivity during the detection of GW170814. In equation (6), $S_n(f)$ is the strain power spectral density (PSD) of the interferometer. The noise level of

the interferometer is mainly determined by the thermal noise and quantum noise.

The SNR $\rho$ for a matched-filter search is determined by the strain PSD $S_n(f)$, Fourier transform of the GW strain amplitude $\tilde{h}(f)$, and $f_{\text{ISCO}}$,

$$ \rho = \sqrt{4Z(i) \int_{0}^{f_{\text{ISCO}}} \frac{|\tilde{h}(f)|^2}{S_n(f)} df}, \quad (7) $$

where $Z(i) = (1 + 6\cos^2 i + \cos^4 i)/8$ accounts for the variation of the radiation power emitted for a binary inclination angle $i$ (Sathyaprakash & Schutz 2009; Schutz 2011). The Fourier transform of the GW strain amplitude is (Abadie et al. 2010)

$$ |\tilde{h}(f)| = \frac{2c}{D} \left( \frac{5G\mu_{\text{red}}}{96c^3} \right)^{1/2} \left( \frac{GM_{\text{tot}}}{\pi^2 c^3} \right)^{1/3} f^{-7/6} \quad (8) $$

with $D$ being the luminosity distance to the source. The detection horizon distance is the maximum distance at which a GW source can be detected with SNR $\rho$,

$$ D_{\text{hor}} = \sqrt{\frac{5}{6} \frac{(GM_{\text{ch}})^{5/6}}{c^{3/2}\pi^{2/3}} \frac{1}{\rho}} \sqrt{Z(i) \int_{0}^{f_{\text{ISCO}}} \frac{1}{S_n(f)} f^{-7/3} df}. \quad (9) $$

Substituting $\rho = 8$ and $S_n(f)$ from equation (6) gives

$$ D_{\text{hor}} = \frac{1.214 \text{ Mpc}}{8} \left( \frac{M_{\text{ch}}}{M_{\odot}} \right)^{5/6} \sqrt{Z(i) \int_{0}^{f_{\text{ISCO}}} \frac{f^{-7/3}}{[S_n(f)]^2} df}, \quad (10) $$

where $\overline{S_n(f)} = 1.986 \times 10^{-4} f^{-0.436} + f^{-3.058} + 2.942 \times 10^{-8} f$. As most of the radiation power is concentrated around the binary rotation axis ($i=0$), the GW detectors preferentially select more face-on binaries with smaller $D_{\text{hor}}$ (Nissanke et al. 2010). The expression for $D_{\text{hor}}$ is derived assuming that the GW source is at a small redshift $z$.

The detection rate $\dot{N}_{GW}$ for a binary system is given by the product of the binary coalescence rate $R_V$ and the detector's horizon volume $V_{\text{hor}}$, $\dot{N}_{GW} = R_V \times V_{\text{hor}} = R_V \times (4/3)\pi D_{\text{hor}}^3$. While the upper limit for BNS coalescence rates is determined by the luminosity distribution of the observed binary pulsars in the Milky Way (Kalogera et al. 2004), it is much more difficult to estimate the BHNS/BBH coalescence rates due to the paucity of observations. The BHNS coalescence rates are based on relatively poorly constrained population synthesis models. The plausible BHNS coalescence rates range from $R_{V,low} \sim 5.80 \times 10^{-4} \text{ Mpc}^{-3} \text{ Myr}^{-1}$ to $R_{V,high} \sim 1.16 \text{ Mpc}^{-3} \text{ Myr}^{-1}$ with a most probable estimate of $R_{V,mp} \sim 3.48 \times 10^{-2} \text{ Mpc}^{-3} \text{ Myr}^{-1}$ (O'Shaughnessy et al. 2008; Abadie et al. 2010). For large distances $D_{\text{hor}} \gtrsim 30 \text{ Mpc}$, the GW detection rate from BHNS sources can then be written as (Abadie et al. 2010)

$$ 
\begin{aligned}
\dot{N}_{GW}(R_V) &= (0.0348 \text{ Myr}^{-1}) \frac{4}{3} \pi D_{\text{hor}}^3 R_{V,0} (2.26)^{-3} \\
&= (0.0126 \text{ year}^{-1}) R_{V,0} D_{\text{hor},2}^3,
\end{aligned} 
\quad (11) $$

where, $R_{V,0} = R_V/R_{V,mp}$, $D_{\text{hor},2} = D_{\text{hor}}/(100 \text{ Mpc})$, and the $(2.26)^{-3}$ factor is included to account for the average over all source orientations and sky locations (Finn & Chernoff 1993). For typical BHNS binaries containing a $1.4 M_\odot$ NS and a $7.8 M_\odot$ BH with source orientation averaged $D_{\text{hor}} \sim 650 \text{ Mpc}$, the most likely GW detection rate with aLIGO design sensitivity is $\dot{N}_{GW}(R_{V,mp}) \sim 3.46 \text{ yr}^{-1}$. Equation (11) is obtained for a Milky-Way-type galaxy and by assuming that $R_V$ is independent of the redshift for distances $D \lesssim D_{\text{hor}} \sim 650 \text{ Mpc}$. Depending on the value of $R_V$, $\dot{N}_{GW}$
can vary by almost three orders of magnitude from $\dot{N}_{\text{GW,low}} \sim 5.77 \times 10^{-2}$ yr$^{-1}$ to $\dot{N}_{\text{GW,high}} \sim 1.16 \times 10^2$ yr$^{-1}$.

## 4 EM SIGNALS ASSOCIATED WITH BHNS MERGERS

In Section 2, we estimated the bound disc mass $M_{\text{disc}}$ and the dynamical ejecta mass (velocity) $M_{\text{ej,dyn}}$ (v$_{\text{ej,dyn}}$) for BHNS coalescence events with a given set of binary parameters. While most of the disrupted NS mass remaining outside the BH at longer time-scales is contained in the remnant accretion disc, a significant fraction is present in the ejecta surrounding the remnant BH. The different ejecta components can be distinguished based on their masses and KEs: an ultrarelativistic sGRB jet, mildly relativistic cocoon, and sub-relativistic dynamical ejecta and magnetic/viscous/neutrino-driven winds. These ejecta components are also directly associated with the EM counterparts of the BHNS merger events. Fig. 1 shows the approximate geometrical distribution of the ejecta components and their associated early- and late-time EM emission in various observing bands. In the following, we briefly discuss the major possible EM counterparts accompanying GW signals from BHNS mergers.

### 4.1 Macronova

Macronovae are thermal transients powered by the radioactive decay of unstable r-process elements via beta-decay and fission (Kasen & Ramirez-Ruiz 2010; Metzger et al. 2010a). This radiation is emitted by the adiabatically expanding ejecta on time-scales of ~a few days to months and is typically peaked in the optical/IR bands (Metzger et al. 2010b; Roberts et al. 2011). The photons escape the bulk and contribute significantly to the EM luminosity only when the density decreases sufficiently such that the photon diffusion time-scale matches the expansion time-scale. Both the dynamical ejecta associated with the tidal tails and the magnetic/viscous/neutrino-driven winds from the remnant accretion disc contribute to the macronova light curve (Fernandez & Metzger 2013; Kasen, Fernandez & Metzger 2015; Fernandez et al. 2017).

The dynamical ejecta from BHNS mergers is expected to be considerably anisotropic from numerical relativity simulations and is characterized by the opening angle in the equatorial (meridional) plane $\phi_{\text{ej,dyn}}(\theta_{\text{ej,dyn}})$. The dynamical ejecta sweeps out about half of the equatorial plane with $\phi_{\text{ej,dyn}} \approx \pi$ and is mainly concentrated around the equatorial plane with $\theta_{\text{ej,dyn}} \approx 10^{\circ}-20^{\circ}$ (Kyutoku et al. 2015). Kyutoku, Kunihito & Shibata (2013) approximated the anisotropic ejecta geometry using an axisymmetric cylinder with radial/z-direction velocity $v_{\parallel}/v_{\perp}$. The peak time for the emission is determined by photons escaping the ejecta along the z-direction due to shorter distance

$$t_{\text{peak}} = (4 \text{ d}) \kappa_1^{1/2} M_{\text{dyn},0.03}^{1/2} v_{\text{dyn},0.3}^{-1/2} \theta_{\text{dyn},0}^{1/2} \phi_{\text{dyn},0}^{-1/2} \quad (12)$$

where $\kappa = \kappa_1 \times 10 \text{ cm}^2 \text{ g}^{-1}$ is the opacity with $M_{\text{dyn},0.03} = M_{\text{ej,dyn}}/0.03 M_{\odot}$, $v_{\text{dyn},0.3} = v_{\text{ej,dyn}}/0.3c$, $\theta_{\text{dyn},0} = \theta_{\text{ej,dyn}}/0.2$, and $\phi_{\text{dyn},0} = \phi_{\text{ej,dyn}}/\pi$ being short-hand notations for ejecta parameters. The peak luminosity is estimated using the radioactive heating rate,

$$L_{\text{peak}} \approx \frac{f M_{\text{ej,dyn}} c^2}{t_{\text{peak}}} = (1.6 \times 10^{41} \text{ erg s}^{-1}) f_{-6} \kappa_1^{-1/2} M_{\text{dyn},0.03}^{1/2} v_{\text{dyn},0.3}^{1/2} \theta_{\text{dyn},0}^{-1/2} \phi_{\text{dyn},0}^{1/2} \quad (13)$$

where $f = f_{-6} \times 10^{-6}$ is the fractional ejecta energy radiated in time $t_{\text{peak}}$.

As opposed to the dynamical ejecta, the wind from the accretion disc is spherically symmetric with a mass $M_{\text{wind}} = \eta_{\text{wind}} M_{\text{disc}}$, where $\eta_{\text{wind}}$ is the fraction of the disc mass that gets unbound and contributes to the macronova emission. As the wind is sub-relativistic with velocity $v_{\text{wind}} \approx 0.07c$ (Hotokezaka & Piran 2015), the peak time and luminosity for the emission are

$$t_{\text{peak}} = (0.3 \text{ d}) \kappa_1^{1/2} \eta_{\text{wind},-3}^{1/2} M_{\text{disc},-2}^{1/2} v_{\text{wind},0.07}^{-1/2} \quad (14)$$

$$L_{\text{peak}} = (7.1 \times 10^{38} \text{ erg s}^{-1}) f_{-6} \kappa_1^{-1/2} \eta_{\text{wind},-3}^{1/2} M_{\text{disc},-2}^{1/2} v_{\text{wind},0.07}^{1/2} \quad (15)$$

where $\eta_{\text{wind},-3} = \eta_{\text{wind}}/10^{-3}$, $M_{\text{disc},-2} = M_{\text{disc}}/0.01 M_{\odot}$, and $v_{\text{wind},0.07} = v_{\text{wind}}/0.07c$. The temporal evolution of the macronova optical/IR luminosity is expected to be $\propto (t/t_{\text{peak}})^{-1.3}$, that is the same as the energy injection rate from the radioactive decay of the unstable r-process elements (Metzger et al. 2010a; Tanaka & Hotokezaka 2013). While the low-opacity blue macronova light curve evolving rapidly over the time-scale of ~2d for GW170817 is expected to result from the isotropic wind component, the high-opacity red macronova light curve with longer time-scale ~10d most likely arises due to the anisotropic dynamical ejecta component (Abbott et al. 2017b; Hinderer et al. 2018).

Due to the anisotropic distribution of dynamical ejecta for BHNS mergers, the macronova emission peaks with larger luminosity $L_{\text{peak}}$ and at smaller $t_{\text{peak}}$ compared to the isotropic dynamical ejecta ($\theta_{\text{ej,dyn}} \approx \pi/2$ and $\phi_{\text{ej,dyn}} \approx 2\pi$) from BNS mergers for a similar $M_{\text{ej,dyn}}$. However, as $M_{\text{ej,dyn}}$ for BNS mergers is expected to be larger by a factor of ~10, the macronova emission $L_{\text{peak}}$ from BNS mergers will be about ~1.5 times larger even after including the geometrical factors. The heating efficiency $f$ depends on the electron fraction in the ejecta that is highly uncertain (Wanajo et al. 2014). In this study, we consider an effective value of $f \approx 3 \times 10^{-6}$ (Metzger et al. 2010a). Similarly, the opacity $\kappa$ is also very uncertain due to our limited knowledge of the bound-bound transition line features of r-process elements (Kasen, Badnell & Barnes 2013) and we will hereby assume a fiducial value of $\kappa = 10 \text{ cm}^2 \text{ g}^{-1}$.

### 4.2 sGRB prompt

One of the most widely studied EM counterparts of BHNS mergers is short-duration (<2 s), highly relativistic collimated outflows, otherwise known as sGRBs. These energetic jets are primarily powered by the gravitational and/or rotational energy of the central compact remnant via accretion. The burst duration for emission in gamma/X-ray bands from sGRB jets is determined by the characteristic viscous time-scale for accretion of the remnant disc/torus on to the BH and is given by (Fernandez et al. 2015)

$$t_{\text{peak}} \approx t_{\text{visc}} = \frac{R^2}{v_a} \approx (0.19 \text{ s}) \alpha_{0.03}^{-1} R_2^{3/2} M_{\text{BH,0}}^{-1/2} \left(\frac{H}{R}\right)^2 \quad (16)$$

where $R = R_2 \times 100$ km is the radial extent of the torus in which most of its mass and angular momentum is concentrated, $v_a = \alpha c_s H$ is the disc viscosity due to turbulence, $\alpha = \alpha_{0.03} \times 0.03$ is the disc viscosity parameter, $c_s$ is the sound speed in the medium, $M_{\text{BH}} = M_{\text{BH,0}} \times 8 M_{\odot}$ and $H \approx c_s(GM_{\text{BH}}/R^3)^{1/2} \sim R$ is the vertical scale height of the disc. The isotropic energy output of sGRB jets is

$$E_{\text{jet,iso}} = \epsilon_{\text{jet}} M_{\text{disc}} c^2 = (10^{51} \text{ erg}) \epsilon_{\text{jet},-2} M_{\text{disc},-1}$$
where $\epsilon_{jet} = \epsilon_{jet,-2} \times 10^{-2}$ is the jet energy conversion efficiency factor and $M_{disc,-1} = M_{disc}/0.1 M_{\odot}$. The jet energy is closely associated with the remnant torus mass $M_{disc}$, which itself is determined by the BHNS binary properties such as $M_{BH}$, $a_{BH}$, and NS EoS (see Section 2). As the jet is highly relativistic, the observed sGRB energy is obtained by including the relativistic beaming factor: $E_{jet} = E_{jet,iso} \times (1/2)(1 - \cos\theta_{jet}) = (1/4)\theta_{jet}^2 E_{jet,iso}$, where $\theta_{jet}$ is the jet-opening angle. The peak luminosity for sGRB prompt emission is then given by

$$L_{\text{peak}} \approx \frac{E_{\text{jet}}}{t_{\text{visc}}} = (1.32 \times 10^{51} \text{ erg s}^{-1}) \\ \times \epsilon_{\text{jet},-2} M_{\text{disc},-1} \alpha_{0.03} R_2^{-3/2} M_{\text{BH},0}^{1/2} \left(\frac{H}{R}\right)^{-2} \theta_{\text{jet}}^2, \quad (17)$$

where all the jet and disc parameters except $\theta_{jet}$ are normalized to their typical values. While the sGRB prompt gamma/X-ray luminosity is roughly constant for $t \lesssim \text{min}$, it is expected to reduce considerably as $\propto (t/\text{min})^{-3}$ at later times and until $t \sim 10^2 \text{ s}$ (Kumar & Zhang 2015).

## 4.3 Cocoon prompt

The observed increase in the afterglow luminosity $\propto t^{0.8}$ for GW170817 up to $t \sim 150 \text{ d}$ after the merger and the subsequent rapid decline can be robustly explained for a compact source with a structured jet including a cocoon (Lazzati et al. 2018; Ghirlanda et al. 2019). The VLBI observations of GW170817 also suggest that the early-time radio emission over the first few months was dominated by the mildly relativistic cocoon and the late-time emission around $t \sim 150 \text{ d}$ was primarily due to the narrowly collimated relativistic jet with an opening angle $\theta_{jet} \lesssim 5^{\circ}$ seen from a viewing angle $\theta_v \sim 14^{\circ}-28^{\circ}$ (Abbott et al. 2017b; Mooley et al. 2018b). The interaction of the lighter outgoing ultrarelativistic sGRB jet with the denser material ejected previously (magnetic/viscous/neutrino-driven wind from the disc and/or dynamical ejecta) along the rotation axis of the remnant BH generates a hot mildly relativistic cocoon surrounding the jet. The cocoon expands approximately spherically after breaking out from the ejecta surface and produces prompt emission in gamma/X-ray bands with energy $E_c \sim E_{jet,iso}$.

The cocoon energy is written in terms of the jet luminosity $L_{jet}$ and the time $t_p$ that the jet needs to propagate across the denser ejecta as $E_c = L_{jet}t_p$. As the jet head traverses a distance $R_a$ with velocity $\beta_{jh}$ in time $t_p$, $E_c = L_{jet}R_a/(c\beta_{jh})$. The jet head velocity $\beta_{jh}$ is estimated by assuming that the lateral width of the cocoon can be determined by balancing the pressure of the cocoon at equilibrium with the ram pressure of the jet (Matzner 2003), $\beta_{jh} \approx (L_{jet}\pi^2/\rho_a t_p^2 \Omega_{jet}^2 c^5)^{1/5}$. Here, $\rho_a$ is the density of the ambient ejecta and $\Omega_{jet}$ is the solid angle occupied by the jet. Substituting $\beta_{jh}$ in $E_c = L_{jet}R_a/(c\beta_{jh})$ further gives $E_c = (L_{jet}R_a^5 \rho \Omega_{jet}^2 \pi^{-2})^{1/3} \sim 4 \times 10^{48} \text{ erg for typical parameters: } L_{jet} \sim 2 \times 10^{50} \text{ erg s}^{-1}$, $R_a \sim 10^8 \text{ cm}$, $\rho_a \sim 10^7 \text{ g cm}^{-3}$, and $\theta_{jet} \sim 16^{\circ}$ (Lazzati et al. 2017).

After breaking out from the surrounding ejecta, the cocoon fireball accelerates outwards adiabatically until it attains a saturation Lorentz factor $\Gamma_c$ at radial distance $R_{sat}$ and finally releases the advected radiation at the photosphere with radius (Meszaros & Rees 2000)

$$R_c = \left( \frac{E_c \sigma_T}{8 \pi m_p \Gamma_c^3 c^2} \right)^{1/2} = (4.2 \times 10^{11} \text{ cm}) E_{c,49}^{1/2} \Gamma_{c,1}^{-3/2}, \quad (18)$$

where $E_{c,49} = E_c/10^{49}$ erg and $\Gamma_{c,1} = \Gamma_c/10$. As the cocoon fireball cools adiabatically, its temperature drops as $T_c = T_{\text{in,c}}(R_c/R_{\text{sat}})^{-2/3} \approx 10 \text{ keV for typical parameters (Lazzati et al. 2017), where } T_{\text{in,c}} \text{ is the initial cocoon temperature. The peak cocoon prompt luminosity can then be written assuming a blackbody spectrum and including the cocoon relativistic beaming factor}$

$$L_{\text{peak}} = (1.0 \times 10^{49} \text{ erg s}^{-1}) E_{c,49}^{2/3} \Gamma_{c,1}^{5/3} R_{a,8}^{-1/3} \theta_c^2, \quad (19)$$

where $R_{a,8} = R_a/10^8 \text{ cm and } \theta_c$ is the cocoon opening angle. For the typical jet and cocoon parameters, the prompt emission lasts for the angular time-scale $t_{\text{ang}}$ given by

$$t_{\text{peak}} \sim t_{\text{ang}} \approx (0.14 \text{ s}) R_{c,0} \Gamma_{c,1}^{-2} = (0.14 \text{ s}) E_{c,49}^{1/2} \Gamma_{c,1}^{-7/2}, \quad (20)$$

where $R_{c,0} = R_c/4.2 \times 10^{11} \text{ cm}$. The gamma/X-ray luminosity for the prompt emission from the mildly relativistic cocoon is expected to be approximately constant until $t \sim \text{a few minutes and then fall off at later times as } \propto (t/\text{min})^{-1}$ (Lazzati et al. 2017).

As the relativistic jet propagates through the ejecta, the surrounding material forms a cocoon that accelerates outwards adiabatically and with mildly relativistic velocities. For a jet with a given opening angle, the condition for a cocoon to develop but not overcome the jet is (Matzner 2003)

$$\left(\frac{\theta_{\text{jet}}}{90^{\circ}}\right)^4 < \tilde{L} < \theta_{\text{jet}}^{-4}, \quad (21)$$

where $\tilde{L} = \beta_{jh}^2 \approx 5.31 \times 10^{-3} \theta_{jet}^{-8/5}$, for $t_p \sim R_a/c \approx 0.01 \text{ s}$ and $\Omega_{jet} = 4\pi(1 - \cos\theta_{jet}) = 2\pi\theta_{jet}^2$. Substituting $\tilde{L}$ in equation (21) further gives $\theta_{jet} < 31.05^{\circ}$ and thereby two possible outcomes: (1) a jet with a wide opening angle $\theta_{jet} \gtrsim 30^{\circ}$ that cannot overcome the cocoon fireball and gets choked, and (2) a narrow opening angle jet with $\theta_{jet} \lesssim 30^{\circ}$ that can penetrate the surrounding ejecta to look like an on-axis sGRB. While the cocoon geometry is still somewhat uncertain, the cocoon opening angle is expected to be directly correlated with the jet opening angle. We discuss the effect of jet and cocoon opening angles on their corresponding EM counterparts $L_{peak}$ and $t_{peak}$ along with the detection rate in the next section.

Although BHNS mergers do not produce polar-shocked material unlike BNS mergers and the dynamical ejecta is mostly concentrated around the equatorial plane, the presence of strong magnetic field/viscous flow/neutrino-driven disc winds and sufficient jet-launching time delay $t_{jet,d}$ can result in the accumulation of the required ejecta mass around the polar regions. The cocoon mass can be written as $M_{\text{cocoon}} = \eta_{\text{dyn}}M_{\text{dyn}} + \eta_{\text{disc}}M_{\text{disc}}$, where $\eta_{\text{dyn/disc}}$ is the contribution to the cocoon mass from the dynamical ejecta/disc wind component. The relative contributions from the dynamical ejecta and the wind components depend on $t_{jet,d}$, with the ratio $\eta_{\text{disc}}/\eta_{\text{dyn}} \propto t_{jet,d}$ as the sub-relativistic winds progressively drive more unbound disc mass towards the existing ejecta material along the BH rotation axis.

## 4.4 Synchrotron afterglow

As the energetic ejecta launched from BHNS mergers expands outwards, it interacts with the surrounding ISM and its KE gets converted to the internal energy of a propagating shock. A significant portion of this energy is transferred to the electrons that are collected from the ambient ISM. Synchrotron afterglow radiation is then emitted by these accelerated non-thermal electrons in the presence of B on time-scales of $\sim \text{a few weeks to a year, and is detectable in radio and (possibly) optical bands. While all ejecta components (sGRB}$
jet, dynamical ejecta, magnetic/viscous/neutrino-driven winds and cocoon) contribute to this emission, different components have different $L_{\text{peak}}$ and $t_{\text{peak}}$ due to their distinct masses, KEs, and geometrical distributions (Hotokezaka & Piran 2015). It should be noted that the *ultrarelativistic jet*, mildly relativistic cocoon, and *sub-relativistic* dynamical ejecta are spherically asymmetric while the *sub-relativistic* wind ejecta component is spherically symmetric. Here we assume that the dynamical ejecta/jet velocity $v_{\text{ej,dyn}}/v_{\text{jet}}$ is practically unchanged until $t \sim t_{\text{peak}}$, resulting in fixed ejecta/jet geometry.

Most of the synchrotron emission is radiated once the ejecta decelerates considerably by accumulating a mass comparable to its own at a radius $R_{\text{dec,s}} = (3M_{\text{ej}}/4\pi m_{\text{p}}n)^{1/3} = (3.2 \text{ pc})M_{\text{ej,0.03}}^{1/3}n^{-1/3}$. Here $M_{\text{ej}} = M_{\text{ej,0.03}} \times 0.03 M_{\odot}$ is the spherically symmetric equivalent of ejecta mass, $m_{\text{p}}$ is the proton mass, and $n = n_{-2} \times 0.01 \text{ cm}^{-3}$ is the constant density of the ambient ISM. The value of $n$ can vary by orders of magnitude based on the location of the BHNS coalescence in the host galaxy and here we assume an average density $n \approx 0.01 \text{ cm}^{-3}$ that is similar to the typical sGRB median densities and slightly larger than the GW170817 expected circumburst density $\sim 10^{-4}-5 \times 10^{-3} \text{ cm}^{-3}$ (Fong et al. 2015; Hotokezaka et al. 2016; Mooley et al. 2018b). While the ejecta expands at a constant velocity $v_{\text{ej}}$ upto $R = R_{\text{dec,s}}$, it decelerates at a larger radius according to the Sedov–Taylor solution, $v_{\text{ej}}(R/R_{\text{dec,s}})^{-3/2}$. The deceleration time for spherical ejecta is then given by

$$t_{\text{dec,s}} = \frac{R_{\text{dec,s}}}{\Gamma_{\text{ej}}^2 v_{\text{ej}}} = (32.5 \text{ yr}) M_{\text{ej,0.03}}^{1/3} n_{-2}^{-1/3} v_{\text{ej,0.3}}^{-1} \Gamma_{\text{ej,0}}^{-2}, \quad (22)$$

where $\Gamma_{\text{ej}} = \Gamma_{\text{ej,0}} \times 1.0$ is the Lorentz factor of the ejecta component. For the anisotropic ejecta components (jet, cocoon, and dynamical ejecta), the deceleration radius includes angular factors, $R_{\text{dec}} = R_{\text{dec,s}}(\pi/2\theta_{\text{ej}})^{1/3}(2\pi/\phi_{\text{ej}})^{1/3} = (7.9 \text{ pc})M_{\text{ej,0.03}}^{1/3}n_{-2}^{-1/3}\theta_{\text{ej,0}}^{-1/3}\phi_{\text{ej,0}}^{-1/3}$ and the deceleration time is (Kyutoku et al. 2013)

$$t_{\text{dec}} = \frac{R_{\text{dec}}}{\Gamma_{\text{ej}}^2 v_{\parallel}} = (83.6 \text{ yr}) M_{\text{ej,0.03}}^{1/3} n_{-2}^{-1/3} v_{\text{ej,0.3}}^{-1} \theta_{\text{ej,0}}^{-1/3} \phi_{\text{ej,0}}^{-1/3} \Gamma_{\text{ej,0}}^{-2}, \quad (23)$$

where we have used the relation $v_{\parallel} \approx v_{\text{ej}}$. The peak emission time $t_{\text{peak}} \sim t_{\text{dec}}$ for synchrotron afterglow varies within a significant range based on the properties of the ejecta components: ~a few days to a week for the *ultrarelativistic* jet, ~a few weeks to a month for the *mildly relativistic* cocoon, ~a few years to a decade for the *sub-relativistic* dynamical ejecta and winds.

A considerable fraction $\epsilon_c/\epsilon_B$ of the post-shock internal energy is transferred to the non-thermal electrons/magnetic fields. For typical blast waves with $v_{\min} < v_a < v_c$ (see Appendix B for definitions), the synchrotron spectrum at a given time is

$$F_{\nu} = F_{\nu, \min} \begin{cases} \left( \frac{\nu_a}{\nu_{\min}} \right)^{(p+4)/2} \left( \frac{\nu}{\nu_{\min}} \right)^2, & \nu < \nu_{\min} \\ \left( \frac{\nu_a}{\nu_{\min}} \right)^{(p-1)/2} \left( \frac{\nu}{\nu_a} \right)^{5/2}, & \nu_{\min} \le \nu < \nu_a \\ \left( \frac{\nu}{\nu_{\min}} \right)^{(p-1)/2}, & \nu_a \le \nu < \nu_c \\ \left( \frac{\nu_c}{\nu_{\min}} \right)^{(p-1)/2} \left( \frac{\nu}{\nu_c} \right)^{-p/2}, & \nu \ge \nu_c \end{cases} \quad (24)$$

where $p \approx 2.5$ is the power-law index for the energy distribution of the accelerated electrons.

For most sensitive radio (optical) observations, $v_a \le \nu < v_c$ ($\nu \ge \nu_c$), and the peak specific flux/luminosity $F_{\nu,\text{peak}}/L_{\text{peak}}$ are computed using the respective $F_{\nu}$ expressions from equation (24). The flux $F_{\nu}$ across the whole spectrum increases for $t < t_{\text{dec}}$ and peaks at $t = t_{\text{dec}}$ for typical radio/optical observing frequencies

$v_{\text{obsr/o}} > v_a(t_{\text{dec}})$, $v_{\min}(t_{\text{dec}})$. The peak of the observed radio/optical specific flux $F_{\text{obsr/o}}$ at frequency $v_{\text{obsr/o}}$

$$F_{\text{vobsr}} = (1.1 \text{ mJy}) \epsilon_{e, -1}^{p-1} \epsilon_{B, -1}^{(p+1)/4} n_{-2}^{(p+1)/4} M_{\text{ej,0.03}} \times v_{\text{ej,0.3}}^{(5p-3)/2} D_2^{-2} \left( \frac{v_{\text{obsr}}}{1 \text{ GHz}} \right)^{-(p-1)/2}, \quad (25)$$

$$F_{\text{vobs}} = (0.2 \text{ Jy}) \epsilon_{e, -1}^{p-1} \epsilon_{B, -1}^{(p-2)/4} n_{-2}^{(3p-2)/12} M_{\text{ej,0.03}}^{2/3} \times v_{\text{ej,0.3}}^{(5p-4)/2} \theta_{\text{ej,0}}^{1/3} \phi_{\text{ej,0}}^{1/3} D_2^{-2} \left( \frac{v_{\text{obsr}}}{1 \text{ GHz}} \right)^{-p/2}, \quad (26)$$

where $D = D_2 \times 100 \text{ Mpc}$ is the distance to the binary source ignoring any cosmological effects. We have normalized $\epsilon_e \approx \epsilon_B \sim 0.1$ to their typical values with $\epsilon_{e, -1} = \epsilon_e/0.1$ and $\epsilon_{B, -1} = \epsilon_B/0.1$. The peak radio/optical luminosity $L_{\text{peak/o}}$ is then estimated from the peak specific flux $F_{\text{vobsr/o}}$ as $L_{\text{peak/o}} = 4\pi D^2(v_{\text{obsr/o}}/F_{\text{vobsr/o}})$.

As $v_a(t_{\text{dec}}), v_{\min}(t_{\text{dec}}) < v_{\text{obsr/o}}$, the synchrotron radio/optical afterglow luminosity for the *sub-relativistic* dynamical ejecta and wind increases as $\propto (t/t_{\text{peak}})^3$ until $t \approx t_{\text{peak}}$ and subsequently falls off as $\propto (t/t_{\text{peak}})^{-(15p-21)/10} \propto t^{-1.65}$ at later times (Piran et al. 2013). For the *ultrarelativistic* jet, the afterglow luminosity increases as $\propto (t/t_{\text{peak}})^3$ before $t = t_{\text{peak}}$ and then drops as $\propto (t/t_{\text{peak}})^{-1}$ for $t \gtrsim t_{\text{peak}}$ (Kumar & Zhang 2015), whereas for the *mildly relativistic* cocoon, the afterglow luminosity increases as $\propto (t/t_{\text{peak}})^3$ for $t \lesssim t_{\text{peak}}$ and drops beyond $t \sim t_{\text{peak}}$ as $\propto (t/t_{\text{peak}})^{-1.5}$ (Rezzolla & Kumar 2015; Lyman et al. 2018). The synchrotron radio afterglow peak luminosity from the dynamical ejecta for BHNS mergers will be about an order of magnitude smaller compared to BNS mergers due to the smaller dynamical ejecta mass for the former (see equation 25).

# 5 EM FOLLOW-UP OF GW TRIGGERS

In the previous section, we discussed the EM signals associated with GW emission from BHNS mergers. Here we estimate the detection rates corresponding to each of those EM counterparts for low-latency follow-up of GW triggers. In addition to enhancing the confidence on the astrophysical origin of the GW signal, the detection of EM counterparts aids precise source localization, redshift estimation and further provides information about the jet and the ejecta physics. The EM localization for a GW event via follow-up also makes it considerably easier to search for the transient in a relatively smaller sky area as compared to blind surveys across large sky areas.

## 5.1 Event rates for EM counterparts

For a given EM counterpart with peak luminosity $L_{\text{peak}}$, the maximum distance at which it can be detected using a telescope/survey with bolometric flux detection threshold $S_{\text{th}}$ in a given observing band is $D_{\text{max}} = \sqrt{L_{\text{peak}}/4\pi S_{\text{th}}}$. The $S_{\text{th}}$ values for some of the current gamma/X-ray, optical, and radio surveys that can be used to observe the EM signals accompanying BHNS GW triggers are listed in Table 2. The detection rate for the EM follow-up of BHNS GW events can then be written similarly to equation (11) as

$$ \begin{aligned} \dot{N}_{\text{EM}} &= (0.0348 \text{ yr}^{-1}) \times (2.26)^{-3} \int_{0}^{D_{\text{EM}}} R_{\nu,0} \cdot 4\pi D_{2}^{2} dD_{2} \\ &= (0.0126 \text{ yr}^{-1}) R_{\nu,0} D_{\text{EM},2}^{3} \end{aligned} \quad (27) $$
<table><caption>**Table 2.** EM telescope sensitivities in different observing bands.</caption><thead><tr><th>Observing band</th><th>Survey</th><th>Frequency range (Hz)</th><th>Nominal v<sub>obs</sub> (Hz)</th><th>Sensitivity (Jy-Hz)</th><th>Cadence (days)</th><th>Reference</th></tr></thead><tbody><tr><td>Gamma-ray</td><td>Fermi GBM</td><td>(1.93–9.67) × 10<sup>21</sup></td><td>4.84 × 10<sup>21</sup></td><td>10<sup>-6</sup></td><td>1</td><td>Meegan et al. (2009)</td></tr><tr><td>X-ray</td><td>Swift BAT</td><td>(0.36–3.63) × 10<sup>19</sup></td><td>1.21 × 10<sup>19</sup></td><td>10<sup>-8</sup></td><td>2</td><td>Cusumano et al. (2010)</td></tr><tr><td rowspan="3">Optical</td><td>LSST</td><td>(4.28–5.44) × 10<sup>14</sup></td><td>4.84 × 10<sup>14</sup></td><td>2.81 × 10<sup>-15</sup></td><td>3</td><td>Ivezic et al. (2019)</td></tr><tr><td>PTF</td><td>(4.28–5.44) × 10<sup>14</sup></td><td>4.84 × 10<sup>14</sup></td><td>6.78 × 10<sup>-14</sup></td><td>5</td><td>Law et al. (2009)</td></tr><tr><td>ZTF</td><td>(4.28–5.44) × 10<sup>14</sup></td><td>4.84 × 10<sup>14</sup></td><td>2.76 × 10<sup>-14</sup></td><td>1</td><td>Cao, Nugent & Kasliwal (2016)</td></tr><tr><td rowspan="3">Radio</td><td>Apertif</td><td>(1.00–1.75) × 10<sup>9</sup></td><td>1.32 × 10<sup>9</sup></td><td>1.32 × 10<sup>-21</sup></td><td>1</td><td>Oosterloo et al. (2010)</td></tr><tr><td>ASKAP</td><td>(0.70–1.80) × 10<sup>9</sup></td><td>1.12 × 10<sup>9</sup></td><td>1.12 × 10<sup>-18</sup></td><td>1</td><td>Johnston et al. (2009)</td></tr><tr><td>LOFAR</td><td>(0.10–2.00) × 10<sup>8</sup></td><td>4.50 × 10<sup>7</sup></td><td>4.50 × 10<sup>-19</sup></td><td>1</td><td>Fender et al. (2006)</td></tr></tbody></table>

where $D_{EM} = \min(D_{hor}, D_{max})$ and $D_2 = D/(100 \text{ Mpc})$. For obtaining equation (27), we have assumed that the BHNS coalescence rate $R_V$ is broadly independent of the distance, which is expected for relatively small distances, $D_{EM} \lesssim D_{hor} \sim 650 \text{ Mpc}$.

We use the following subscripts to denote the EM counterparts: $i = A1/2$ for a macronova from dynamical ejecta/wind, $B/C$ for sGRB/cocoon prompt emission, and $D1/2/3/4$ for dynamical ejecta/wind/sGRB jet/cocoon afterglow. We denote the radio/optical afterglow emission with the subscript r/o, and estimate the masses of all the ejecta components in terms of $M_{disc}$ and $M_{ej,dyn}$ (see Section 2). In Table 3, we list $t_{peak}(M_{ej}, v_{ej})$ and $L_{peak}(M_{ej}, v_{ej})$ used to evaluate $D_{max}(M_{ej}, v_{ej}, S_{th})$ and $\dot{N}_{EM}$ (from equation 27) for each EM counterpart.

(i) *Early-time emission*: The early-time EM emission consists of a macronova from the dynamical ejecta/wind, sGRB jet prompt, and cocoon prompt.

(a) *Macronova optical/IR follow-up*: The quasi-thermal emission from the decay of unstable r-process elements in the wind/dynamical ejecta peaks in the optical/IR bands. We compute $\dot{N}_{EM,A1/2}$ for typical optical/IR sensitivities $S_{th,OIR} \sim 10^{-15}–10^{-13}$ Jy Hz (see Table 2) and ejecta parameters: $\kappa \sim 10 \text{ cm}^2 \text{ g}^{-1}$, $\theta_{ej,dyn} \sim 0.1745 (10')$, $\phi_{ej,dyn} \sim \pi$, and $f \sim 3 \times 10^{-6}$ (see Table 3).

(b) *sGRB and cocoon prompt gamma/X-ray follow-up*: The short-duration sGRB prompt emission and cocoon prompt emission peak in the gamma/X-ray bands. We estimate the sGRB/cocoon prompt follow-up rate $\dot{N}_{EM,B/C}$ for typical gamma/X-ray sensitivities $S_{th,GX} \sim 10^{-8}–10^{-6}$ Jy Hz (see Table 2). We use characteristic accretion disc parameters: $\alpha \sim 0.03$, $R \sim 100 \text{ km}$, $M_{BH} \sim 8 \text{ M}_\odot$, and $H \sim R \approx 100 \text{ km}$ with a jet-opening angle $\theta_{jet} = 1/\Gamma_{jet} \approx 0.1$. The hot mildly relativistic cocoon surrounding the jet has energy $E_c \sim \eta_c(\Gamma_{jet}/\Gamma_{coc})(\eta_{jet}M_{disc}c^2) \approx 3\eta_{jet}M_{disc}c^2 = (5.4 \times 10^{48} \text{ erg})\eta_{jet}, -4M_{disc},-2$ and opening angle $\theta_{coc} = 1/\Gamma_{coc} = 3\theta_{jet}$, where $\eta_c \sim 1$ is the fraction of the jet energy that is transferred to the cocoon. We consider $\Gamma_{coc} = (1/3)\Gamma_{jet}$ and $R_a = 10^8 \text{ cm}$ as the typical cocoon parameters.

(ii) *Late-time emission*: The late-time emission consists of synchrotron radio and optical afterglows from the interaction of the dynamical ejecta, wind, jet, and cocoon components with the ambient ISM. The $t_{peak}$ and $L_{peak}$ values for the afterglow emission from a particular ejecta component critically depend on its mass and velocity along with the geometry (see Table 1). Here we assume $n = 0.01 \text{ cm}^{-3}$, $\epsilon_e = \epsilon_B = 0.1$, $p = 2.5$, and $v_{obsr} = 1 \text{ GHz}$ to estimate the radio emission detection rates for each ejecta component (see Appendix C for optical emission detection rates). The typical radio band sensitivity varies within $S_{th,R} \sim 10^{-21}–10^{-18}$ Jy Hz (see Table 2).

We consider $\theta_{ej,dyn} \sim 0.1745$ and $\phi_{ej,dyn} \sim \pi$ for the anisotropic sub-relativistic dynamical ejecta. A significant fraction $\eta_{wind}$ of the disc mass can also become unbound due to the spherically symmetric sub-relativistic magnetic field/viscous flow/neutrino-driven winds and contribute to the afterglow emission. The wind afterglow follow-up detection rate in radio band is significantly smaller compared to the other ejecta components for typical ejecta parameters and telescope sensitivities, thereby making it very challenging to detect.

The mass of the ultrarelativistic jet can be written as $M_{jet} = \eta_{jet}M_{disc}$, where $\eta_{jet}$ is the fraction of the remnant accretion disc mass contributing to the jet kinetic energy. As the jet propagates outwards with highly relativistic velocities $v_{jet} \approx c$, the emission is significantly anisotropic with $\theta_{jet} = 1/\Gamma_{jet}$ and $\phi_{jet} = 2\pi$. We assume $\Gamma_{jet} \sim 10$ to be the typical jet Lorentz factor to compute the EM follow-up rate. We consider $\eta_{disc}/\eta_{dyn} \sim 0.1$ for the mildly relativistic anisotropic cocoon generated by the interaction of the outgoing ultrarelativistic jet with the previously ejected material (see Section 4.3). The cocoon Lorentz factor $\Gamma_{coc}$ is expected to be directly correlated with $\Gamma_{jet}$ and here we assume a typical value of $\Gamma_{coc} \sim (1/3)\Gamma_{jet}$. For $\Gamma_{jet} \sim 10–30$ (see Table 1), $\Gamma_{jet} \sim 3–10$ with a cocoon velocity $v_{coc} \approx c$. The emission from the relativistic cocoon is anisotropic with $\theta_{coc} = 1/\Gamma_{coc} = 3\theta_{jet}$ and $\phi_{coc} = 2\pi$.

## 5.2 Effect of binary parameters on the EM follow-up rates

Here we study the effect of the BHNS binary parameters ($M_{BH}$, $a_{BH}$, and NS EoS) on the $L_{peak}$ and $t_{peak}$ values of the EM signals accompanying the GW trigger. We determine the range of binary parameters for which NS tidal disruption is favoured and the accompanying EM signals are detectable with the current telescope sensitivities. Fig. 2 shows the effect of binary parameters on the $L_{peak}$ and $t_{peak}$ values of the EM counterparts for two different cases: (a) less NS disruption for APR4 $R_{NS} = 11.1$ km, $M_{BH} = 8.2 \text{ M}_\odot$, and $a_{BH} = 0.8$ and (b) more NS disruption for H4 $R_{NS} = 13.6$ km, $M_{BH} = 7.4 \text{ M}_\odot$, and $a_{BH} = 0.97$. For less (more) NS disruption, $M_{ej,tot}/M_{tot} = 1.736 \text{ per cent}$ (6.155 per cent) with $M_{ej,dyn}/M_{tot} = 0.003 \text{ per cent}$ (1.190 per cent), $M_{disc}/M_{tot} = 1.733 \text{ per cent}$ (4.965 per cent) and $v_{ej,dyn} = 0.284 \text{c} (0.275c)$. We normalize the values of the ejecta parameters and the observing band sensitivities as in Section 5.1 for early-time emission and radio afterglow components, and Appendix C for optical afterglow components.

The early-time sGRB prompt emission lasts for $t_{peak,B} = (1.88–1.98) \times 10^{-1} \text{ s} \lesssim \text{ second}$ with a peak luminosity $L_{peak,B} = (2.21–5.52) \times 10^{48} \text{ erg s}^{-1}$ in gamma/X-ray. The cocoon prompt emission peaks around $t_{peak,C} = (1.96–3.17) \times 10^1 \text{ s} \sim \text{second–minute}$ with a peak luminosity $L_{peak,C} = (0.62–$
<table><caption>**Table 3.** The EM follow-up detection parameters are listed for each component. We use the definitions: $v_{\text{obstr,GHz}} = v_{\text{obstr}}/(1 \text{ GHz})$, $M_{\text{dyn},-3} = M_{\text{ej,dyn}}/10^{-3} \text{ M}_{\odot}$, $v_{\text{dyn},0.2} = v_{\text{ej,dyn}}/0.2c$, $\epsilon_{\text{jet},-3} = \epsilon_{\text{jet}}/10^{-3}$, $\eta_{\text{jet},-4} = \eta_{\text{jet}}/10^{-4}$, $\Gamma_{\text{jet},1} = \Gamma_{\text{jet}}/10$, $\eta_{\text{dyn},-2} = \eta_{\text{dyn}}/10^{-2}$, and $\eta_{\text{disc},-3} = \eta_{\text{disc}}/10^{-3}$. The optical/IR, gamma/X-ray, and radio band sensitivities used to compute EM detection rates are $S_{\text{OIR}} = 2 \times 10^{-14} \text{ Jy Hz}$, $S_{\text{GX}} = 10^{-7} \text{ Jy Hz}$ and $S_{\text{R}} = 3 \times 10^{-17} \text{ Jy Hz}$, respectively. The EM follow-up rate for a given component is $N_{\text{EM}}/N_{\text{GW}} = \min(1, D_{\text{max}}^3/D_{\text{hor}}^3)$ for telescope sensitivity $S_{\text{th}}$.</caption><thead><tr><th>EM component</th><th>t<sub>peak</sub></th><th>L<sub>peak</sub></th><th>D<sub>max</sub></th><th>Relevant equations</th></tr></thead><tbody><tr><td>[A1]: Ejecta<br>macronova</td><td>(0.84 d) ×<br>M<sup>1/2</sup><sub>dyn,−3</sub>v<sup>−1/2</sup><sub>dyn,0.2</sub></td><td>(7.66 × 10<sup>40</sup> erg s<sup>−1</sup>)<br>×M<sup>1/2</sup><sub>dyn,−3</sub>v<sup>1/2</sup><sub>dyn,0.2</sub></td><td>(178.94 Mpc)M<sup>1/4</sup><sub>dyn,−3</sub> ×<br>v<sup>1/4</sup><sub>dyn,0.2</sub>(S<sub>th</sub>/S<sub>OIR</sub>)<sup>−1/2</sup></td><td>12, 13</td></tr><tr><td>[A2]: Wind<br>macronova</td><td>(0.3 d)η<sup>1/2</sup><sub>wind,−3</sub> ×<br>M<sup>1/2</sup><sub>disc,−2</sub>v<sup>−1/2</sup><sub>wind,0.07</sub></td><td>(2.13 × 10<sup>39</sup> erg s<sup>−1</sup>)η<sup>1/2</sup><sub>wind,−3</sub><br>×M<sup>1/2</sup><sub>disc,−2</sub>v<sup>1/2</sup><sub>wind,0.07</sub></td><td>(29.84 Mpc)η<sup>1/4</sup><sub>wind,−3</sub>M<sup>1/4</sup><sub>disc,−2</sub><br>×v<sup>1/4</sup><sub>wind,0.07</sub>(S<sub>th</sub>/S<sub>OIR</sub>)<sup>−1/2</sup></td><td>14, 15</td></tr><tr><td>[B]: sGRB<br>prompt</td><td>~0.19 s</td><td>(1.32 × 10<sup>47</sup> erg s<sup>−1</sup>) ×<br>ε<sub>jet,−3</sub>M<sup>1/2</sup><sub>disc,−2</sub>Γ<sup>−2</sup><sub>jet,1</sub></td><td>(105.05 Mpc)ε<sup>1/2</sup><sub>jet,−3</sub> ×<br>M<sup>1/2</sup><sub>disc,−2</sub>Γ<sup>−1</sup><sub>jet,1</sub>(S<sub>th</sub>/S<sub>GX</sub>)<sup>−1/2</sup></td><td>16, 17</td></tr><tr><td>[C]: Cocoon<br>prompt</td><td>(4.81 s)η<sup>1/2</sup><sub>jet,−4</sub> ×<br>M<sup>1/2</sup><sub>disc,−2</sub>Γ<sup>−7/2</sup><sub>jet,1</sub></td><td>(9.56 × 10<sup>46</sup> erg s<sup>−1</sup>) ×<br>η<sup>2/3</sup><sub>jet,−4</sub>M<sup>2/3</sup><sub>disc,−2</sub>Γ<sup>−1/3</sup><sub>jet,1</sub></td><td>(89.40 Mpc)η<sup>1/3</sup><sub>jet,−4</sub> ×<br>M<sup>1/3</sup><sub>disc,−2</sub>Γ<sup>−1/6</sup><sub>jet,1</sub>(S<sub>th</sub>/S<sub>GX</sub>)<sup>−1/2</sup></td><td>19, 20</td></tr><tr><td>[D1]: Ejecta<br>afterglow</td><td>(38.99 yr) ×<br>M<sup>1/3</sup><sub>dyn,−3</sub>v<sup>−1</sup><sub>dyn,0.2</sub></td><td>(6.31 × 10<sup>34</sup> erg s<sup>−1</sup>)M<sub>dyn,−3</sub><br>×v<sup>(5p−3)/2</sup><sub>dyn,0.2</sub>v<sup>(−p−3)/2</sup><sub>obsr,GHz</sub></td><td>(4.19 Mpc)M<sup>1/2</sup><sub>dyn,−3</sub> ×<br>v<sup>(5p−3)/4</sup><sub>dyn,0.2</sub>(S<sub>th</sub>/S<sub>R</sub>)<sup>−1/2</sup></td><td>23, 25</td></tr><tr><td>[D2]: Wind<br>afterglow</td><td>(9.61 yr)η<sup>1/3</sup><sub>wind,−3</sub><br>×M<sup>1/3</sup><sub>disc,−2</sub>v<sup>−1</sup><sub>wind,0.07</sub></td><td>(4.30 × 10<sup>30</sup> erg s<sup>−1</sup>)η<sub>wind,−3</sub><br>×M<sub>disc,−2</sub>v<sup>(5p−3)/2</sup><sub>wind,0.07</sub>v<sup>(−p−3)/2</sup><sub>obsr,GHz</sub></td><td>(0.035 Mpc)η<sup>1/2</sup><sub>wind,−3</sub>M<sup>1/2</sup><sub>disc,−2</sub><br>×v<sup>(5p−3)/4</sup><sub>wind,0.07</sub>(S<sub>th</sub>/S<sub>R</sub>)<sup>−1/2</sup></td><td>22, 25</td></tr><tr><td>[D3]: Jet<br>afterglow</td><td>(70.69 h)η<sup>1/3</sup><sub>jet,−4</sub><br>×M<sup>1/3</sup><sub>disc,−2</sub>Γ<sup>−5/3</sup><sub>jet,1</sub></td><td>(1.32 × 10<sup>35</sup> erg s<sup>−1</sup>)η<sub>jet,−4</sub><br>×M<sub>disc,−2</sub>v<sup>(−p−3)/2</sup><sub>obsr,GHz</sub></td><td>(6.06 Mpc)η<sup>1/2</sup><sub>jet,−4</sub> ×<br>M<sup>1/2</sup><sub>disc,−2</sub>(S<sub>th</sub>/S<sub>R</sub>)<sup>−1/2</sup></td><td>23, 25</td></tr><tr><td>[D4]: Cocoon<br>afterglow</td><td>(39.50 d)Γ<sup>−5/3</sup><sub>jet,1</sub><br>[η<sub>dyn,−2</sub>M<sub>dyn,−3</sub> +<br>η<sub>disc,−3</sub>M<sub>disc,−2</sub>]<sup>1/3</sup></td><td>(1.32 × 10<sup>36</sup> erg s<sup>−1</sup>) ×<br>[η<sub>dyn,−2</sub>M<sub>dyn,−3</sub> +<br>η<sub>disc,−3</sub>M<sub>disc,−2</sub>]v<sup>(−p−3)/2</sup><sub>obsr,GHz</sub></td><td>(19.18 Mpc)(S<sub>th</sub>/S<sub>R</sub>)<sup>−1/2</sup><br>[η<sub>dyn,−2</sub>M<sub>dyn,−3</sub> +<br>η<sub>disc,−3</sub>M<sub>disc,−2</sub>]<sup>1/2</sup></td><td>23, 25</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Effect of BHNS binary parameters on the EM counterpart $L_{\text{peak}}$ and $t_{\text{peak}}$ values: simulation results for two BHNS binaries with fixed $M_{\text{BH}}$, $a_{\text{BH}}$, and NS EoS distributions. Left-hand panel: less NS disruption for $R_{\text{NS}} = 11.1$ km (APR4 NS EoS), $M_{\text{BH}} = 8.2 \text{ M}_{\odot}$, and $a_{\text{BH}} = 0.8$. Right-hand panel: more NS disruption for $R_{\text{NS}} = 13.6$ km (H4 NS EoS), $M_{\text{BH}} = 7.4 \text{ M}_{\odot}$, and $a_{\text{BH}} = 0.97$. Blue and green symbols in both panels denote radio and optical afterglow components, respectively. The black dashed lines show the power-law fits for the canonical temporal evolution of the luminosity for each counterpart and the symbols represent their peak quantities. The representative bolometric band sensitivities are as listed in Table 2: gamma/X-ray $v_{\text{obs}} = 10^{20}$ Hz, optical/IR $v_{\text{obs}} = 5 \times 10^{14}$ Hz, and radio $v_{\text{obs}} = 10^{9}$ Hz.</figcaption></figure>

1.18) × 10<sup>48</sup> erg s<sup>−1</sup> in gamma/X-ray. The macronova emission from dynamical ejecta/wind peaks at a later time with $t_{\text{peak,A1/A2}} = (0.38 – 7.27)/(1.22 – 1.97)$ d ~ hours to a week and a peak optical/IR luminosity $L_{\text{peak,A1/A2}} = (0.49 – 9.16) \times 10^{41}/(0.87 – 1.41) \times 10^{40}$ erg s<sup>−1</sup>. Although the macronova emission peaks with significantly smaller luminosities compared to both sGRB

and cocoon prompt, it has comparable detection rates in the optical/IR bands and is a relatively easy follow-up target due to its peak times of ~ hours to a week. The late-time afterglow emission from the ultrarelativistic jet peaks at $t_{\text{peak,D3}} = (0.75–1.04) \times 10^1$ d with radio/optical luminosity $L_{\text{peak,D3r/o}} = (2.18 – 5.73) \times 10^{36}/(0.98 – 1.86) \times 10^{38}$ erg s<sup>−1</sup>. The mildly rel-
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Effect of binary parameters on the $\dot{N}_{EM}$ values for early-time macronova and late-time radio afterglow emission from anisotropic sub-relativistic dynamical ejecta: contour plots using simulation results for 300 BHNS binaries with $M_{BH} = (7.8 \pm 1.2) M_{\odot}$, $a_{BH} = \text{Uniform}(0, 0.97)$, and APR4/H4 NS EoS. For the ith EM counterpart, $\dot{N}_{EM,i}$ denotes $\dot{N}_{EM,i}/\text{max}(\dot{N}_{GW})$. Top left-hand panel: $\dot{N}_{A1}$ contour plot for APR4 NS EoS. Top right-hand panel: $\dot{N}_{A1}$ contour plot for H4 NS EoS. Bottom left-hand panel: $\dot{N}_{D1r}$ contour plot for APR4 NS EoS. Bottom right-hand panel: $\dot{N}_{D1r}$ contour plot for H4 NS EoS.</figcaption></figure>

ativistic cocoon afterglow peaks at longer time-scales with $t_{peak,D4} = (4.71-9.71) \times 10^4$ d ~ a month with a radio/optical luminosity $L_{peak,D4r/o} = (0.22-1.95) \times 10^{37}/(1.43-6.07) \times 10^{38}$ erg s$^{-1}$. The afterglow from sub-relativistic dynamical ejecta peaks around $t_{peak,D1} = (0.18-1.34) \times 10^2$ yr ~ a decade with luminosity in radio/optical band $L_{peak,D1r/o} = (0.01-2.97) \times 10^{37}/(0.03-1.31) \times 10^{38}$ erg s$^{-1}$. The late-time afterglow from the sub-relativistic wind component peaks with similar time-scales of $t_{peak,D2} = (2.45-3.38) \times 10^1$ yr ~ a decade and radio/optical luminosity $L_{peak,D2r/o} = (0.71-1.87) \times 10^{32}/(1.40-2.67) \times 10^{34}$ erg s$^{-1}$.

We propagate the error in $M_{ej,tot}$ and $M_{ej,dyn}$ obtained for both less and more NS disruption cases to estimate the error in $L_{peak}$ and $t_{peak}$ for each EM counterpart. We find that the gamma/X-ray $\Delta L_{peak,B} = (1.90-5.52) \times 10^{47}$ erg s$^{-1}$ for sGRB prompt emission whereas the error in peak time and gamma/X-ray luminosity for cocoon prompt emission is $\Delta t_{peak,C} = (0.84-1.58)$ s and $\Delta L_{peak,C} = (3.54-7.91) \times 10^{46}$ erg s$^{-1}$, respectively. While the error in the peak time for macronova emission from dynamical ejecta/wind is $\Delta t_{peak,A1/A2} = (0.40-3.24)/(0.05-0.10)$ d, the corresponding error in the optical/IR peak luminosity is $\Delta L_{peak,A1/A2} = (0.51-4.22) \times 10^{41}/(3.75-7.05) \times 10^{38}$ erg s$^{-1}$. The jet afterglow peak time and radio/optical peak luminosity error are relatively small with $\Delta t_{peak,D3} = (0.22-0.34)$ d and $\Delta L_{peak,D3r/o} = (1.88-5.73) \times 10^{35}/(0.84-1.86) \times 10^{37}$ erg s$^{-1}$, while the corresponding quantities for the late-time cocoon emission are marginally higher with $\Delta t_{peak,D4} = (1.84-3.30)$ d and $\Delta L_{peak,D4r/o} = (0.26-1.99) \times 10^{36}/(1.67-6.19) \times 10^{37}$ erg s$^{-1}$. The magnitude of the error in the ejecta masses implies that the error in dynamical ejecta afterglow peak time and radio/optical luminosity are significantly large with $\Delta t_{peak,D1} = (0.05-1.04) \times 10^2$ yr and $\Delta L_{peak,D1r/o} = (1.66-3.27) \times 10^{36}/(1.44-5.15) \times 10^{37}$ erg s$^{-1}$. Similarly, the error in the peak quantities for the wind afterglow emission are $\Delta t_{peak,D2} = (0.71-1.12)$ yr and $\Delta L_{peak,D2r/o} = (0.61-1.87) \times 10^{31}/(1.20-2.67) \times 10^{33}$ erg s$^{-1}$.

We further evaluate the EM follow-up rates $\dot{N}_{EM}$ from the $L_{peak}$ values of different counterparts and for given bolometric sensitivities $S_{th}$. Figs 3–6 show the effect of BHNS binary parameters $M_{BH}$, $a_{BH}$, and NS EoS on the EM follow-up rate for the early-time and late-time EM emissions from the anisotropic sub-relativistic dynamical ejecta, the isotropic sub-relativistic winds, the anisotropic ultrarelativistic sGRB jet, and the anisotropic mildly relativistic cocoon, respectively. The Monte Carlo simulation results are shown for 300 BHNS binaries with Gaussian $M_{BH} = (7.8 \pm 1.2) M_{\odot}$, $a_{BH} = \text{Uniform}(0, 0.97)$, and APR4/H4 NS EoS.

Fig. 3 shows the contour plots for the EM follow-up rates of early-time macronova and late-time radio afterglow emission from the dynamical ejecta for APR4 and H4 NS EoSs. For the soft APR4
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Effect of binary parameters on the $\dot{N}_{EM}$ values for early-time macronova and late-time radio afterglow emission from isotropic sub-relativistic wind: contour plots using simulation results for 300 BHNS binaries with $M_{BH} = (7.8 \pm 1.2) M_{\odot}$, $a_{BH} = \text{Uniform}(0, 0.97)$, and APR4/H4 NS EoS. Top left-hand panel: $\dot{N}_{A2}$ contour plot for APR4 NS EoS. Top right-hand panel: $\dot{N}_{A2}$ contour plot for H4 NS EoS. Bottom left-hand panel: $\dot{N}_{D2r}$ contour plot for APR4 NS EoS. Bottom right-hand panel: $\dot{N}_{D2r}$ contour plot for H4 NS EoS.</figcaption></figure>

NS EoS, the ejecta macronova emission has an optical/IR follow-up rate of $0.10 \le \dot{N}_{A1}/\dot{N}_{GW} \le 0.30$ for $0.75 \le a_{BH} \le 1.0$ and $6.0 \le M_{BH}/M_{\odot} \le 10.0$. However, the optical/IR follow-up rate is larger by a factor of $\sim 2$ for the stiffer H4 EoS, with $0.22 \le \dot{N}_{A1}/\dot{N}_{GW} \le 0.60$ for $0.7 \le a_{BH} \le 1.0$ and $5.0 \le M_{BH}/M_{\odot} \le 11.0$. The radio afterglow emission from the dynamical ejecta is considerable for $0.75 \le a_{BH} \le 1.0$ ($0.6 \le a_{BH} \le 1.0$) and $6.5 \le M_{BH}/M_{\odot} \le 10.0$ ($6.0 \le M_{BH}/M_{\odot} \le 11.0$) with $0.18 \le \dot{N}_{D1r}/\dot{N}_{GW} \le 0.50$ ($0.3 \le \dot{N}_{D1r}/\dot{N}_{GW} \le 0.9$) in the case of the APR4 (H4) EoS. The radio afterglow follow-up rates are larger by a factor of $\sim 1.5-2$ for the stiffer H4 EoS compared to the softer APR4 NS EoS due to more ejecta from NS disruption. Fig. 4 shows the contour plots for the follow-up rates of early-time wind macronova and late-time wind radio afterglow for the APR4 and H4 NS EoSs. The wind macronova has a significantly smaller optical/IR follow-up rate compared to the ejecta macronova with $4.0 \times 10^{-4} \le \dot{N}_{A2}/\dot{N}_{GW} \le 10^{-3}$ for $0.5 \le a_{BH} \le 0.8$ and $5.0 \le M_{BH}/M_{\odot} \le 10.0$ in the region where $15.0a_{BH} - M_{BH} \ge 2.5$ and for the APR4 NS EoS. The optical/IR follow-up rate is similar for the H4 NS EoS with $4.0 \times 10^{-4} \le \dot{N}_{A2}/\dot{N}_{GW} \le 1.1 \times 10^{-3}$ for $0.2 \le a_{BH} \le 0.7$ and $5.0 \le M_{BH}/M_{\odot} \le 10.5$ in the region where $M_{BH} - 12.0a_{BH} \le 2.0$. The radio afterglow emission from the wind ejecta component has a negligibly small follow-up detection rate:

for the APR4 NS EoS, $0.85 \times 10^{-8} \le \dot{N}_{D2r}/\dot{N}_{GW} \le 2.35 \times 10^{-8}$ when $0.7 \le a_{BH} \le 0.9$ and $5.0 \le M_{BH}/M_{\odot} \le 10.0$ whereas for the H4 NS EoS, $1.0 \times 10^{-8} \le \dot{N}_{D2r}/\dot{N}_{GW} \le 3.0 \times 10^{-8}$ when $0.55 \le a_{BH} \le 0.85$ and $4.5 \le M_{BH}/M_{\odot} \le 10.5$ in the region where $22.0a_{BH} - M_{BH} \ge 8.7$. Due to the significantly small $L_{\text{peak,D2r}}$ and consequently $\dot{N}_{D2r}$, any subsequent radio follow-up of GW triggers using wind afterglow emission is extremely difficult with the current telescope sensitivities (see Table 2).

The contour plots for the EM follow-up rates of early-time prompt and late-time radio afterglow emission from the sGRB jet for the APR4 and H4 EoS are shown in Fig. 5. The jet prompt emission has a follow-up rate of $0.18 \le \dot{N}_B/\dot{N}_{GW} \le 0.50$ for $0.8 \le a_{BH} \le 1.0$ and $6.0 \le M_{BH}/M_{\odot} \le 10.0$ in case of the APR4 EoS. The corresponding rate is larger by a factor of $\sim 1.5$ for the H4 NS EoS with $0.28 \le \dot{N}_B/\dot{N}_{GW} \le 0.78$ for $0.8 \le a_{BH} \le 1.0$ and $6.0 \le M_{BH}/M_{\odot} \le 11.0$. For the APR4 EoS, the jet radio afterglow has the follow-up rate $0.045 \le \dot{N}_{D3r}/\dot{N}_{GW} \le 0.125$ for $0.7 \le a_{BH} \le 1.0$ and $5.0 \le M_{BH}/M_{\odot} \le 10.0$. However, the radio afterglow follow-up rate for the H4 EoS is slightly larger with $0.06 \le \dot{N}_{D3r}/\dot{N}_{GW} \le 0.16$ for $0.7 \le a_{BH} \le 1.0$ and $4.5 \le M_{BH}/M_{\odot} \le 10.5$. Fig. 6 shows the EM follow-up rate contour plots for the early-time cocoon prompt and the late-time cocoon radio afterglow for the APR4 and H4 NS EoSs. For the softer APR4
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Effect of binary parameters on the $\dot{N}_{EM}$ values for early-time prompt and late-time radio afterglow emission from anisotropic ultrarelativistic sGRB jet: contour plots using simulation results for 300 BHNS binaries with $M_{BH} = (7.8 \pm 1.2) M_{\odot}$, $a_{BH} = \text{Uniform}(0, 0.97)$, and APR4/H4 NS EoS. Top left-hand panel: $\dot{N}_B$ contour plot for APR4 NS EoS. Top right-hand panel: $\dot{N}_B$ contour plot for H4 NS EoS. Bottom left-hand panel: $\dot{N}_{D3r}$ contour plot for APR4 NS EoS. Bottom right-hand panel: $\dot{N}_{D3r}$ contour plot for H4 NS EoS.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

EoS, the cocoon prompt emission has rate $0.14 \le \dot{N}_C/\dot{N}_{GW} \le 0.36$ for $0.7 \le a_{BH} \le 1.0$ and $6.0 \le M_{BH}/M_{\odot} \le 10.0$. The gamma/X-ray follow-up rate for the H4 EoS is similar with $0.15 \le \dot{N}_C/\dot{N}_{GW} \le 0.42$ for $0.6 \le a_{BH} \le 1.0$ and $4.5 \le M_{BH}/M_{\odot} \le 11.0$. The cocoon afterglow has a rate of $0.1 \le \dot{N}_{D4r}/\dot{N}_{GW} \le 0.3$ for $0.75 \le a_{BH} \le 1.0$, $6.0 \le M_{BH}/M_{\odot} \le 10.0$, and APR4 EoSs. The follow-up rate for the H4 EoS is larger by a factor of $\sim 3$ with $0.3 \le \dot{N}_{D4r}/\dot{N}_{GW} \le 0.82$ for $0.7 \le a_{BH} \le 1.0$ and $6.0 \le M_{BH}/M_{\odot} \le 10.5$.

## 6 DISCUSSIONS AND CONCLUSIONS

In this paper, we studied the effect of BHNS binary parameters on the properties of the ejected material as well as their associated EM emission. The tidal disruption of the NS and the outcome of the subsequent BHNS merger are primarily determined by $M_{BH}$, $a_{BH}$, and the EoS of the infalling NS. While the majority of the tidally disrupted NS mass gets directly accreted on to the BH on short viscous time-scales, a considerable fraction of the ejecta with mass $M_{\text{disc}}$ forms a remnant disc around the BH as the remaining ejecta material with mass $M_{\text{ej,dyn}}$ becomes gravitationally unbound from the compact remnant. We consider $M_{BH}$ and $a_{BH}$ distributions motivated from transient LMXB observations in order to estimate

the ejecta masses $M_{\text{disc}}$ and $M_{\text{ej,dyn}}$. We find that $M_{\text{disc}}$ and $M_{\text{ej,dyn}}$ are practically unaffected by the specific $M_{BH}$ distribution considered (see Appendix A), provided that the mean BH masses are similar, whereas both the ejecta masses increase with the increase in BH spin magnitude. The ejecta masses are expectedly larger for the stiffer H4 NS EoS due to larger disruption.

We obtain the aLIGO design strain sensitivity as a function of detection frequency to estimate the detection horizon volume and thereby the event detection rates for a BHNS binary with a given orientation using the coalescence rates from population synthesis models. Although the GW event detection rates can vary by about three orders of magnitude around the mean value depending on the poorly constrained BHNS coalescence rate, the mean rate for a fixed binary configuration is practically independent of the choice of $M_{BH}$ and $a_{BH}$ distributions and is unaffected by the NS EoS. In addition to being promising sources of GWs detectable by aLIGO, BHNS binaries are expected to have significant EM emission providing further information about the properties of the ejecta material as well as the post-merger remnant.

We discuss the possible EM counterparts for these BHNS mergers and estimate their peak emission parameters. The EM emission is primarily due to four ejecta components: the ultrarelativistic sGRB jet, the mildly relativistic cocoon, the sub-relativistic dynamical
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Effect of binary parameters on the $\dot{N}_{EM}$ values for early-time prompt and late-time radio afterglow emission from anisotropic mildly relativistic cocoon: contour plots using simulation results for 300 BHNS binaries with $M_{BH} = (7.8 \pm 1.2) M_{\odot}$, $a_{BH} = \text{Uniform}(0, 0.97)$, and APR4/H4 NS EoS. Top left-hand panel: $\dot{N}_C$ contour plot for APR4 NS EoS. Top right-hand panel: $\dot{N}_C$ contour plot for H4 NS EoS. Bottom left-hand panel: $\dot{N}_{D4r}$ contour plot for APR4 NS EoS. Bottom right-hand panel: $\dot{N}_{D4r}$ contour plot for H4 NS EoS.</figcaption></figure>

<figure><img src="image_13.png" /></figure>

ejecta, and the *sub-relativistic* magnetic/viscous/neutrino-driven winds. The early-time emission consists of the jet prompt, cocoon prompt, and macronova whereas the late-time emission comprises synchrotron afterglow from the jet, cocoon, dynamical ejecta, and wind. In contrast to the isotropic emission from the wind, the emissions from the jet, cocoon, and dynamical ejecta components are expected to be fairly anisotropic.

The low-latency EM follow-up observations of GW triggers helps us in constraining the EoS of neutron degenerate matter and jet + ejecta properties post-merger in addition to improving the source localization. The EM follow-up volume in all observing bands is restricted to the GW detection horizon volume dictated by the aLIGO sensitivity and makes it easier for the transient search in a smaller sky area compared to the blind surveys. We perform Monte Carlo simulations for 300 BHNS binaries with Gaussian $M_{BH} = (7.8 \pm 1.2) M_{\odot}$ and $a_{BH} = \text{Uniform}(0, 0.97)$ distributions, and APR4/H4 NS EoSs in order to estimate the follow-up rates of the EM counterparts in terms of the binary parameters. The EM follow-up rates for the stiffer H4 NS EoS are larger compared to the corresponding rates for the APR EoS by a factor of ~1.5–2 and ~1.5–3 for early-time and late-time emissions, respectively. Almost all EM counterparts are detectable for typically small $M_{BH}$ and large $a_{BH}$, leading to significant NS disruption and subsequent stronger EM emission from the ejecta material. The wind macronova and afterglow detectability are especially challenging with the current

radio telescope sensitivities due to the small isotropic luminosity of this sub-relativistic component.

Below we summarize the main results of this work.

(i) The mass ejected in a BHNS merger event is primarily determined by the binary parameters: mass ratio $q = M_{BH}/M_{NS}$, magnitude of $a_{BH}$, and $R_{NS}$ based on the NS EoS. We find that both $M_{ej,tot}$ and $M_{ej,dyn}$ increase as NS disruption increases with a decrease in $M_{BH}$ (regardless of the choice of the distribution, provided the mean mass is similar), increase in $a_{BH}$, and stiffer NS EoS (see Figs A1 and A2). Due to the smaller dynamical ejecta masses and merger rates for typical BHNS mergers compared to BNS mergers, the BHNS mergers are estimated to synthesize only about 0.1 per cent of the total r-process elements found in the Universe. The *ultrarelativistic* sGRB jet, the *mildly relativistic* cocoon, the *sub-relativistic* dynamical ejecta, and the *sub-relativistic* magnetic/viscous/neutrino-driven winds are the different BHNS merger ejecta components, each associated with some early- and/or late-time EM emission. Fig. 1 shows a schematic sketch of all the ejecta components along with their characteristic emission and Table 1 lists their properties.

(ii) For the typical ejecta parameters and range of BHNS binary parameters allowing NS disruption (see Fig. 2), we estimate the EM peak luminosities and emission time-scales for both early- and late-time emissions.
The early-time emissions are in chronological order:

(1) sGRB jet prompt emission in gamma/X-ray ($\nu_{\text{obs}} \sim 10^{20}$ Hz) with $t_{\text{peak,B}} \lesssim$ second and $L_{\text{peak,B}} \sim 10^{48-49}$ erg s$^{-1}$ until $t \lesssim$ min and $\propto (t/\text{min})^{-3}$ afterwards,

(2) cocoon prompt emission in gamma/X-ray with $t_{\text{peak,C}} \sim$ second–minute and $L_{\text{peak,C}} \sim 10^{47-48}$ erg s$^{-1}$ until $t \approx t_{\text{peak,C}}$ and $\propto (t/t_{\text{peak,C}})^{-1}$ at later times,

(3) ejecta macronova emission in optical/IR ($\nu_{\text{obs}} \sim 5 \times 10^5$ GHz) with $t_{\text{peak,A1}} \sim$ day–week and $L_{\text{peak,A1}} \sim 10^{40-42}$ erg s$^{-1}$,

(4) wind macronova emission in optical/IR with $t_{\text{peak,A2}} \sim$ hour–day and $L_{\text{peak,A2}} \sim 10^{39-40}$ erg s$^{-1}$.

The macronova luminosity for both dynamical ejecta and magnetic/viscous/neutrino-driven winds are roughly constant until $t \approx t_{\text{peak}}$ and fall off as $\propto (t/t_{\text{peak}})^{-1.3}$ at later times.

The corresponding order for the late-time radio synchrotron afterglow emissions with $\nu_{\text{obs}} \sim 1$ GHz is

(1) jet afterglow with $t_{\text{peak,D3}} \sim$ day–week and $L_{\text{peak,D3}} \sim 10^{36-38}$ erg s$^{-1}$. $L(t) \propto (t/t_{\text{peak,D3}})^3$ for $t \lesssim t_{\text{peak,D3}}$ whereas $L(t) \propto (t/t_{\text{peak,D3}})^{-1}$ for $t > t_{\text{peak,D3}}$,

(2) cocoon afterglow with $t_{\text{peak,D4}} \sim$ week–month and $L_{\text{peak,D4}} \sim 10^{36-39}$ erg s$^{-1}$. The cocoon luminosity steadily increases as $\propto (t/t_{\text{peak,D4}})^3$ for $t \lesssim t_{\text{peak,D4}}$ and drops as $\propto (t/t_{\text{peak,D4}})^{-1.5}$ at later times,

(3) dynamical ejecta afterglow with $t_{\text{peak,D1}} \sim$ year–decade and $L_{\text{peak,D1}} \sim 10^{35-38}$ erg s$^{-1}$,

(4) wind afterglow with $t_{\text{peak,D2}} \sim$ year–decade and $L_{\text{peak,D2}} \sim 10^{32-34}$ erg s$^{-1}$.

The afterglow luminosity for the sub-relativistic dynamical ejecta and wind components is expected to increase as $t^3$ upto $t \approx t_{\text{peak}}$ and then reduce as $t^{-1.65}$ at later times.

(iii) We estimate the follow-up detection rates of all early- and late-time EM counterparts for representative bolometric sensitivities: optical/IR $S_{\text{th}} \sim 2 \times 10^{-14}$ Jy Hz, gamma/X-ray $S_{\text{th}} \sim 10^{-7}$ Jy Hz, and radio $S_{\text{th}} \sim 3 \times 10^{-19}$ Jy Hz (see Figs 3–6, B1, and B2). We find that the early-time sGRB prompt gamma/X-ray, cocoon prompt gamma/X-ray, ejecta macronova optical/IR, and wind macronova optical/IR follow-up rates are $\dot{N}_B \sim (0.2 - 0.8) \dot{N}_{\text{GW}}$, $\dot{N}_C \sim (0.1 - 0.4) \dot{N}_{\text{GW}}$, $\dot{N}_{A1} \sim (0.1 - 0.6) \dot{N}_{\text{GW}}$, and $\dot{N}_{A2} \sim (0.4 - 1.1) \times 10^{-3} \dot{N}_{\text{GW}}$, respectively, for the BH mass range $6.0 M_{\odot} \lesssim M_{\text{BH}} \lesssim 10.0 M_{\odot}$ and spin range $0.7 \lesssim a_{\text{BH}} \lesssim 1.0$. The late-time synchrotron afterglow radio follow-up rates for the jet, cocoon, dynamical ejecta, and wind ejecta components are $\dot{N}_{D3r} \sim (0.04 - 0.16) \dot{N}_{\text{GW}}$, $\dot{N}_{D4r} \sim (0.10 - 0.82) \dot{N}_{\text{GW}}$, $\dot{N}_{D1r} \sim (0.18 - 0.90) \dot{N}_{\text{GW}}$, and $\dot{N}_{D2r} \sim (0.85 - 3.00) \times 10^{-8} \dot{N}_{\text{GW}}$, respectively, for $6.0 M_{\odot} \lesssim M_{\text{BH}} \lesssim 10.0 M_{\odot}$ and $0.7 \lesssim a_{\text{BH}} \lesssim 1.0$.

(iv) The low-latency EM follow-up rate $\dot{N}_{\text{EM}} \propto D_{\text{EM}}^{-1}$ (from equation 27) for each counterpart. Based on the follow-up rates for $6 M_{\odot} \lesssim M_{\text{BH}} \lesssim 10 M_{\odot}$ and $0.7 \lesssim a_{\text{BH}} \lesssim 1.0$, the relative ease of follow-up of BHNS EM counterparts for typical ejecta parameters and telescope band sensitivities is in the order: ejecta afterglow (radio) > cocoon afterglow (radio) $\gtrsim$ jet prompt (gamma/X-ray) > ejecta macronova (optical/IR) > cocoon prompt (gamma/X-ray) > jet afterglow (radio) $\gg$ wind macronova (optical/IR) $\gg$ wind afterglow (radio).

## REFERENCES

Abadie J. et al., 2010, Class. Quantum Gravity, 27, 173001

Abbott B. P. et al., 2017a, Phys. Rev. Lett., 119, 161101

Abbott B. P. et al., 2017b, ApJ, 848, L12

Abbott B. P. et al., 2017c, ApJ, 850, L39

Acernese F. et al., 2015, Class. Quantum Gravity, 32, 024001

Alexander K. D. et al., 2017, ApJ, 848, L21

Barbieri C., Salafia O. S., Perego A., Colpi M., Ghirlanda G., 2019, preprint (arXiv:1903.04543)

Brege W. et al., 2018, Phys. Rev. D, 98, 063009

Cao Y., Nugent P. E., Kasliwal M. M., 2016, PASP, 128, 114502

Clark J. A., Bauswein A., Stergioulas N., Shoemaker D., 2016, Class. Quantum Gravity, 33, 085003

Coughlin M. W., Dietrich T., 2019, preprint (arXiv:1901.06052)

Cusumano G. et al., 2010, A&A, 524, A64

Dessart L., Ott C. D., Burrows A., Rosswog S., Livne E., 2009, ApJ, 690, 1681

Duez M. D., Foucart F., Kidder L. E., Ott C. D., Teukolsky S. A., 2010, Class. Quantum Gravity, 27, 114106

Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126

Etienne Z. B., Liu Y. T., Shapiro S. L., Baumgarte T. W., 2009, Phys. Rev. D, 79, 044024

Farr W. M., Sravan N., Cantrell A., Kreidberg L., Bailyn C. D., Mandel I., Kalogera V., 2011, ApJ, 741, 103

Fender R. P. et al., 2006, in VI Microquasar Workshop: Microquasars, Beyond. SISSA, Trieste, Italy

Fernández R., Metzger B. D., 2013, MNRAS, 435, 502

Fernández R., Quataert E., Schwab J., Kasen D., Rosswog S., 2015, MNRAS, 449, 390

Fernández R., Foucart F., Kasen D., Lippuner J., Desai D., Roberts L. F., 2017, Class. Quantum Gravity, 34, 154001

Finn L. S., Chernoff D. F., 1993, Phys. Rev. D, 47, 2198

Fong W., Berger E., Margutti R., Zauderer B. A., 2015, ApJ, 815, 102

Foucart F., 2012, Phys. Rev. D, 86, 124007

Foucart F., Duez M. D., Kidder L. E., Teukolsky S. A., 2011, Phys. Rev. D, 83, 024005

Foucart F., Duez M. D., Kidder L. E., Scheel M. A., Szilagyi B., Teukolsky S. A., 2012, Phys. Rev. D, 85, 044015

Foucart F. et al., 2013, Phys. Rev. D, 87, 084006

Foucart F. et al., 2014, Phys. Rev. D, 90, 024026

Foucart F., Hinderer T., Nissanke S., 2018, Phys. Rev. D, 98, 081501

Ghirlanda G. et al., 2019, Science, 363, 968

Haggard D., Nynka M., Ruan J. J., Kalogera V., Cenko S. B., Evans P., Kennea J. A., 2017, ApJ, 848, L25

LIGO Scientific Collaboration, 2010, Class. Quantum Gravity, 27, 084006

Hinderer T. et al., 2018, preprint (arXiv:1808.03836)

Hotokezaka K., Piran T., 2015, MNRAS, 450, 1430

Hotokezaka K., Nissanke S., Hallinan G., Lazio T. J. W., Nakar E., Piran T., 2016, ApJ, 831, 190

Ioka K., Nakamura T., 2018, Prog. Theor. Exp. Phys., 2018, 043E02

Ivezic Z. et al., 2019, ApJ, 873, 111

Johnston S., Feain I. J., Gupta N., 2009, in Saikia D. J., Green D. A., Gupta Y., Venturi T., eds, ASP Conf. Ser. Vol. 407, The Low Frequency Radio Universe. Astron. Soc. Pac., San Francisco, p. 446.

Kalogera V. et al., 2004, ApJ, 601, L179

Kasen D., Ramirez-Ruiz E., 2010, ApJ, 714, 155

Kasen D., Badnell N. R., Barnes J., 2013, ApJ, 774, 25

Kasen D., Fernández R., Metzger B. D., 2015, MNRAS, 450, 1777

Kawaguchi K., Kyutoku K., Nakano H., Okawa H., Shibata M., Taniguchi K., 2015, Phys. Rev. D, 92, 024014

Kawaguchi K., Kyutoku K., Shibata M., Tanaka M., 2016, ApJ, 825, 52

Kiuchi K., Sekiguchi Y., Kyutoku K., Shibata M., Taniguchi K., Wada T., 2015, Phys. Rev. D, 92, 064034

Kruckow M. U., Tauris T. M., Langer N., Kramer M., Izzard R. G., 2018, MNRAS, 481, 1908

Kulkarni S. R., 2005, preprint (arXiv:astro-ph/0510256)

Kumar P., Zhang B., 2015, Phys. Rep., 561, 1

Kyutoku K., Shibata M., Taniguchi K., 2010, Phys. Rev. D, 82, 044049

Kyutoku K., Okawa H., Shibata M., Taniguchi K., 2011, Phys. Rev. D, 84, 064018

Kyutoku K., Kunihito I., Shibata M., 2013, Phys. Rev. D, 88, 041503(R)
Kyutoku K., Ioka K., Okawa H., Shibata M., Taniguchi K., 2015, Phys. Rev. D, 92, 044028

Law N. M. et al., 2009, PASP, 121, 1395

Lazzati D., Deich A., Morsony B. J., Workman J. C., 2017, MNRAS, 471, 1652

Lazzati D., Perna R., Morsony B. J., Lopez-Cámara D., Cantiello M., Ciolfi R., Giacomazzo B., Workman J. C., 2018, Phys. Rev. Lett., 120, 241103

Li L.-X., Paczyński B., 1998, ApJ, 507, L59

Lovelace G., Duez M. D., Foucart F., Kidder L. E., Pfeiffer H. P., Scheel M. A., Szilágyi B., 2013, Class. Quantum Gravity, 30, 135004

Lyman J. D. et al., 2018, preprint (arXiv:1801.02669)

Margutti R. et al., 2018, ApJL, 856, L18

Matzner C. D., 2003, MNRAS, 345, 575

McClintock J. E., Remillard R. A., 2006, in Lewin W., van der Klis M., eds, Compact Stellar X-ray Sources. Cambridge Univ. Press, Cambridge

Meegan C. et al., 2009, ApJ, 702, 791

Mészáros P., Rees M. J., 2000, ApJ, 530, 292

Metzger B. D. et al., 2010a, MNRAS, 406, 2650

Metzger B. D., Berger E., 2012, ApJ, 746, 48

Metzger B. D., Arcones A., Quataert E., Martinez-Pinedo G., 2010b, MNRAS, 402, 2771

Miller J. M. et al., 2004, ApJ, 601, 45

Mooley K. P. et al., 2018a, Nature, 554, 207

Mooley K. P. et al., 2018b, Nature, 561, 355

Nakar E., 2007, Phys. Rep., 442, 166

Nakar E., Piran T., 2011, Nature, 478, 82

Narayan R., 2005, Ap&SS, 300, 177

Nissanke S., Holz D. E., Hughes S. A., Dalal N., Sievers J. L., 2010, ApJ, 725, 496

O'Shaughnessy R., Kaplan J., Kalogera V., Belczynski K., 2005, ApJ, 632, 1035

O'Shaughnessy R., Kim C., Kalogera V., Belczynski K., 2008, ApJ, 672, 479

Oosterloo T., Verheijen M., van Cappellen W., 2010, in ISKAF, 2010 Science Meeting. SISSA, Trieste, Italy

Özel F., Psaltis D., Narayan R., McClintock J. E., 2010, ApJ, 725, 1918

Pannarale F., Tonita A., Rezzolla L., 2011, ApJ, 727, 95

Peters P. C., 1964, Phys. Rev. B, 136, 1224

Piran T., Nakar E., Rosswog S., 2013, MNRAS, 430, 2121

Postnov K. A., Yungelson L. R., 2014, Living Rev. Relativ., 17, 3

Read J. S., Lackey B. D., Owen B. J., Friedman J. L., 2009, Phys. Rev. D, 79, 124032

Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49

Reynolds C. S., Young A. J., Begelman M. C., Fabian A. C., 1999, ApJ, 514, 164

Rezzolla L., Kumar P., 2015, ApJ, 802, 95

Roberts L. F., Kasen D., Lee W. H., Ramirez-Ruiz E., 2011, ApJ, 736, L21

Sathyaprakash B. S., Schutz B. F., 2009, Living Rev. Relativ., 12, 2

Schutz B. F., 2011, Class. Quantum Gravity, 28, 125023

Somiya K., 2012, Class. Quantum Gravity, 29, 124007

Steiner J. F., McClintock J. E., Reid M. J., 2012, ApJ, 745, L7

Takami K., Rezzolla L., Baiotti L., 2015, Phys. Rev. D, 91, 064001

Tanaka M., Hotokezaka K., 2013, ApJ, 775, 113

Török G., Abramowicz M. A., Kluźniak W., Suchlík Z., 2005, A&A, 436, 1

Wanajo S., Sekiguchi Y., Nishimura N., Kiuchi K., Kyutoku K., Shibata M., 2014, ApJ, 789, L39

## APPENDIX A: EFFECT OF BINARY PARAMETERS ON EJECTA MASSES

Here we discuss the effect of the BHNS binary parameters ($M_{BH}$, $a_{BH}$, and NS EoS) on the remnant disc and the dynamical ejecta

masses. Fig. A1 shows the effect of the $M_{BH}$ distribution on $M_{disc}$ and $M_{ej,dyn}$ for $a_{BH} = 0.97$ and APR4/H4 NS EoSs. For the APR4 EoS, we find that $M_{disc} \sim 2.26 \times 10^{-1} - 4.15 \times 10^{-1} M_{\odot}$ with similar mean mass $\sim 3.5 \times 10^{-1} M_{\odot}$ for all $M_{BH}$ distributions. The dispersion in $M_{disc}$ is larger in case of a power-law distribution ($3.2 \times 10^{-2} M_{\odot}$) as compared to exponential ($2.8 \times 10^{-2} M_{\odot}$) and Gaussian ($2.6 \times 10^{-2} M_{\odot}$) distributions. In case of the H4 EoS, $M_{disc} \sim 2.5 \times 10^{-11} - 4.6 \times 10^{-1} M_{\odot}$ with mean mass $\sim 4.3 \times 10^{-1} M_{\odot}$ for all $M_{BH}$ distributions.

The dispersion in $M_{disc}$ for a power-law $M_{BH}$ distribution ($3.1 \times 10^{-2} M_{\odot}$) is significantly larger compared to exponential ($1.9 \times 10^{-2} M_{\odot}$) and Gaussian ($1.4 \times 10^{-2} M_{\odot}$) distributions. For the APR4 EoS, the dynamical ejecta mass $M_{ej,dyn}$ lies in the range $\sim 2.5 \times 10^{-2} - 4.6 \times 10^{-2} M_{\odot}$ with mean mass $\sim 4.4 \times 10^{-2} M_{\odot}$ for all $M_{BH}$ distributions. The dispersion in $M_{ej,dyn}$ is in the order: Gaussian ($3.4 \times 10^{-3} M_{\odot}$) > exponential ($2.9 \times 10^{-3} M_{\odot}$) > power-law ($2.5 \times 10^{-3} M_{\odot}$). For the stiffer H4 NS EoS, $M_{ej,dyn} \sim 7.0 \times 10^{-2} - 1.1 \times 10^{-1} M_{\odot}$ with an equal mean mass $\sim 1.0 \times 10^{-1} M_{\odot}$ for all three distributions. The $M_{ej,dyn}$ dispersion for Gaussian ($6.4 \times 10^{-3} M_{\odot}$) is larger compared to power-law ($6.1 \times 10^{-3} M_{\odot}$) and exponential ($5.0 \times 10^{-3} M_{\odot}$) distributions. As expected, the mean masses for both $M_{disc}$ and $M_{ej,dyn}$ are smaller for the softer APR4 NS EoS compared to the stiffer H4 NS EoS and are unaffected by the choice of the $M_{BH}$ distribution. The mass dispersion for $M_{disc}$ and $M_{ej,dyn}$ are relatively unaffected by the NS EoS.

Fig. A2 shows the effect of the $a_{BH}$ distribution on $M_{disc}$ and $M_{ej,dyn}$ for Gaussian $M_{BH}$ and the APR4/H4 NS EoS. For the APR4 NS EoS, the bound disc mass $M_{disc}$ lies within the range $\sim 0-4.0 \times 10^{-1} M_{\odot}$ while the mean mass for $a_{BH} = \text{Unif}(0, 0.97)$ ($7.8 \times 10^{-2} M_{\odot}$) is smaller compared to that for $a_{BH} = 0.75$ ($1.5 \times 10^{-1} M_{\odot}$) and $a_{BH} = 0.97$ ($3.5 \times 10^{-1} M_{\odot}$). The $M_{disc}$ dispersion is in the order: $a_{BH} = \text{Unif}(0, 0.97)$ ($9.5 \times 10^{-2} M_{\odot}$) > $a_{BH} = 0.75$ ($4.5 \times 10^{-2} M_{\odot}$) > $a_{BH} = 0.97$ ($2.4 \times 10^{-2} M_{\odot}$). For the stiffer H4 NS EoS, $M_{disc} \sim 0-4.6 \times 10^{-1} M_{\odot}$ with mean mass in the order: $a_{BH} = 0.97$ ($4.3 \times 10^{-1} M_{\odot}$) > $a_{BH} = 0.75$ ($2.4 \times 10^{-1} M_{\odot}$) > $a_{BH} = \text{Unif}(0, 0.97)$ ($1.5 \times 10^{-1} M_{\odot}$). The dispersion in $M_{disc}$ for $a_{BH} = \text{Unif}(0, 0.97)$ ($1.1 \times 10^{-1} M_{\odot}$) is larger than that of $a_{BH} = 0.75$ ($3.4 \times 10^{-2} M_{\odot}$) and $a_{BH} = 0.97$ ($1.5 \times 10^{-2} M_{\odot}$). The $M_{ej,dyn}$ for the APR4 EoS lies in the range $\sim 0-4.6 \times 10^{-2} M_{\odot}$ with mean $M_{ej,dyn}$ in the order: $a_{BH} = 0.97$ ($4.4 \times 10^{-2} M_{\odot}$) > $a_{BH} = \text{Unif}(0, 0.97)$ ($3.6 \times 10^{-3} M_{\odot}$) > $a_{BH} = 0.75$ ($1.1 \times 10^{-3} M_{\odot}$). The mass dispersion also depends on the $a_{BH}$ distribution with: $a_{BH} = \text{Unif}(0, 0.97)$ ($9.4 \times 10^{-3} M_{\odot}$) > $a_{BH} = 0.75$ ($2.2 \times 10^{-3} M_{\odot}$) > $a_{BH} = 0.97$ ($1.9 \times 10^{-3} M_{\odot}$). For the H4 NS EoS, $M_{ej,dyn} \sim 0-1.1 \times 10^{-1} M_{\odot}$ and the mean $M_{ej,dyn}$ for $a_{BH} = 0.97$ ($1.0 \times 10^{-1} M_{\odot}$) is larger than that for $a_{BH} = 0.75$ ($5.5 \times 10^{-2} M_{\odot}$) and $a_{BH} = \text{Unif}(0, 0.97)$ ($2.6 \times 10^{-2} M_{\odot}$). The $M_{ej,dyn}$ dispersion is in the order: $a_{BH} = \text{Unif}(0, 0.97)$ ($3.0 \times 10^{-2} M_{\odot}$) > $a_{BH} = 0.97$ ($6.2 \times 10^{-3} M_{\odot}$) > $a_{BH} = 0.75$ ($4.5 \times 10^{-3} M_{\odot}$). The mean masses for both $M_{disc}$ and $M_{ej,dyn}$ increase with the BH spin whereas the corresponding mass dispersions decrease with $a_{BH}$ magnitude. This is expected as a larger $a_{BH}$ results in more NS disruption and thereby larger $M_{disc}$ and $M_{ej,dyn}$ values. Similar to the varying $M_{BH}$ distribution case, the mean masses for $M_{disc}$ and $M_{ej,dyn}$ are smaller for the APR4 EoS compared to the H4 EoS while the mass dispersions are unaffected by the choice of NS EoS in general.
<figure><img src="image_14.png" /></figure>

Effect of different $M_{BH}$ distributions on $M_{\text{dyn}}$ for $a=0.97$ and APR4 NS EoS

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

Effect of different $M_{BH}$ distributions on $M_{\text{dyn}}$ for $a=0.97$ and H4 NS EoS

<figure><img src="image_17.png" /></figure>

**Figure A1.** Effect of $M_{BH}$ distribution on the bound disc mass $M_{\text{disc}}$ and the unbound dynamical ejecta mass $M_{\text{ej,dyn}}$: simulation results for 100 BHNS binaries with $a_{BH} = 0.97$, two different NS EoSs (APR4 and H4) and three different $M_{BH}$ distributions (Gaussian, power-law, and exponential). The vertical axes show the number of data points in each mass bin. Top left-hand panel: $M_{\text{disc}}$ distribution for APR4 EoS. Top right-hand panel: $M_{\text{disc}}$ distribution for H4 EoS. Bottom left-hand panel: $M_{\text{ej,dyn}}$ distribution for APR4 EoS. Bottom right-hand panel: $M_{\text{ej,dyn}}$ distribution for H4 EoS.
<figure><img src="image_18.png" /><figcaption>**Figure A2.** Effect of $a_{BH}$ distribution on the bound disc mass $M_{\text{disc}}$ and the unbound dynamical ejecta mass $M_{\text{ej,dyn}}$: simulation results for 100 BHNS binaries with Gaussian $M_{\text{BH}}$ distribution, two different NS EoSs (APR4 and H4), and three different $a_{BH} = 0.75, 0.97$, and Uniform(0, 0.97). Top left-hand panel: $M_{\text{disc}}$ distribution for APR4 EoS. Top right-hand panel: $M_{\text{disc}}$ distribution for H4 EoS. Bottom left-hand panel: $M_{\text{ej,dyn}}$ distribution for APR4 EoS. Bottom right-hand panel: $M_{\text{ej,dyn}}$ distribution for H4 EoS.</figcaption></figure>

## APPENDIX B: SYNCHROTRON EMISSION CHARACTERISTIC FREQUENCIES

The ejecta geometry does not influence either the magnitude of B or the energy distribution of the non-thermal electrons. The distribution of the Lorentz factor $\gamma_e$ of the accelerated electrons is expected to be a power law, $dN_e/d\gamma_e \propto \gamma_e^{-p}$ for $\gamma_e > \gamma_{e,\text{min}} = (4.5 \times 10^{-3})\epsilon_{e,-1}[(p-2)/(p-1)](m_p/m_e)v_{ej,0.3}^2$, where $m_e$ is the electron mass and $\gamma_{e,\text{min}}$ is the minimum electron Lorentz factor. The field strength B is written using the Rankine–Hugoniot relation, $B = \sqrt{9\pi\epsilon_B n m_p} v_{ej} = (2.0 \times 10^{-3} \text{ G})\epsilon_{B,-1}^{1/2} n_{-2}^{1/2} v_{ej,0.3}$.

The radio spectrum for synchrotron emission is characterized by the frequency of electrons $\nu_{\text{min}}$ with Lorentz factor $\gamma_{e,\text{min}}$ and the self-absorption frequency $\nu_a$. The frequency and specific flux in the absence of self-absorption corresponding to $\gamma_{e,\text{min}}$ are (Kyutoku

et al. 2015)

$$ \nu_{\text{min}} = (3.6 \times 10^5 \text{ Hz}) \left(\frac{p-2}{p-1}\right)^2 \epsilon_{e,-1}^2 \epsilon_{B,-1}^{1/2} n_{-2}^{1/2} v_{ej,0.3}^5 \quad (\text{B1}) $$

$$ F_{\nu, \text{min}} = (2.2 \text{ Jy}) \epsilon_{B,-1}^{1/2} n_{-2}^{1/2} M_{ej,0.03} v_{ej,0.3} D_2^{-2}. \quad (\text{B2}) $$

The self-absorption frequency $\nu_a$ is determined by equating the Rayleigh blackbody flux with the unabsorbed flux $F_{\nu,\text{min}}$. The cooling frequency

$$ \nu_c = (3.2 \times 10^{13} \text{ Hz}) \epsilon_{B,-1}^{-3/2} n_{-2}^{-5/6} M_{ej,0.03}^{-2/3} v_{ej,0.3}^{-1} D_2^{-2} \theta_{ej,0}^{2/3} \phi_{ej,0}^{2/3} \quad (\text{B3}) $$

corresponding to the electron Lorentz factor $\gamma_{e,c}$ above which radiative energy losses are significant is unimportant in the radio band and only affects the optical spectrum.
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure B1.** Effect of binary parameters on the $\dot{N}_{EM}$ values for late-time optical afterglow emission from anisotropic dynamical ejecta and isotropic wind: contour plots using simulation results for 300 BHNS binaries with $M_{BH} = (7.8 \pm 1.2) M_{\odot}$, $a_{BH} = \text{Uniform}(0, 0.97)$, and APR4/H4 NS EoS. Top left-hand panel: $\dot{N}_{D10}$ contour plot for APR4 NS EoS. Top right-hand panel: $\dot{N}_{D10}$ contour plot for H4 NS EoS. Bottom left-hand panel: $\dot{N}_{D20}$ contour plot for APR4 NS EoS. Bottom right-hand panel: $\dot{N}_{D20}$ contour plot for H4 NS EoS.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure B2.** Effect of binary parameters on the $\dot{N}_{EM}$ values for late-time optical afterglow emission from anisotropic sGRB jet and anisotropic cocoon: contour plots using simulation results for 300 BHNS binaries with $M_{BH} = (7.8 \pm 1.2) M_{\odot}$, $a_{BH} = \text{Uniform}(0, 0.97)$, and APR4/H4 NS EoS. Top left-hand panel: $N_{D3o}$ contour plot for APR4 NS EoS. Top right-hand panel: $\dot{N}_{D3o}$ contour plot for H4 NS EoS. Bottom left-hand panel: $\dot{N}_{D4o}$ contour plot for APR4 NS EoS. Bottom right-hand panel: $\dot{N}_{D4o}$ contour plot for H4 NS EoS.</figcaption></figure>

<figure><img src="image_22.png" /></figure>

## APPENDIX C: OPTICAL AFTERGLOW EMISSION DETECTION RATES

Here we estimate the optical afterglow emission follow-up rates for the ultrarelativistic jet, the mildly relativistic cocoon, the sub-relativistic dynamical ejecta, and the wind. As in Section 5.1, we assume $n = 0.01 \text{ cm}^{-3}$, $\epsilon_e = \epsilon_B = 0.1$, and $p = 2.5$ for all ejecta components. We use $v_{\text{obso}} = v_{\text{obso,GHz}} \times 1 \text{ GHz} = 5 \times 10^5 \text{ GHz}$ as the observing frequency in the optical band. The optical afterglow peak luminosity for the anisotropic sub-relativistic dynamical ejecta, isotropic sub-relativistic wind, anisotropic ultrarelativistic jet, and anisotropic mildly relativistic cocoon are

$$L_{\text{peak, D1o}} = (4.09 \times 10^{37} \text{ erg s}^{-1}) M_{\text{dyn, -3}}^{2/3} \times v_{\text{dyn, 0.2}}^{(5p-4)/2} v_{\text{obso, GHz}}^{-(p-2)/2} \quad \text{(C1)}$$

$$L_{\text{peak, D2o}} = (5.75 \times 10^{34} \text{ erg s}^{-1}) \eta_{\text{wind, -3}}^{2/3} M_{\text{disc, -2}}^{2/3} \times v_{\text{wind, 0.07}}^{(5p-4)/2} v_{\text{obso, GHz}}^{-(p-2)/2} \quad \text{(C2)}$$

$$L_{\text{peak, D3o}} = (4.01 \times 10^{38} \text{ erg s}^{-1}) \eta_{\text{jet, -4}}^{2/3} M_{\text{disc, -2}}^{2/3} \Gamma_{\text{jet, 1}}^{-1/3} \times v_{\text{obso, GHz}}^{-(p-2)/2} \quad \text{(C3)}$$

$$L_{\text{peak, D4o}} = (2.68 \times 10^{39} \text{ erg s}^{-1}) [\eta_{\text{dyn, -2}} M_{\text{dyn, -3}} + \eta_{\text{disc, -3}} M_{\text{disc, -2}}]^{2/3} \Gamma_{\text{jet, 1}}^{-1/3} \times v_{\text{obso, GHz}}^{-(p-2)/2} \quad \text{(C4)}$$

with the corresponding follow-up rates for $v_{\text{obso}} = 5 \times 10^5 \text{ GHz}$

$$\dot{N}_{\text{D1o}} = (1.51 \times 10^3 \text{ yr}^{-1}) M_{\text{dyn, -3}} v_{\text{dyn, 0.2}}^{3(5p-4)/4} \times (S_{\text{th, D1o}}/S_{\text{th, O}})^{-3/2} \quad \text{(C5)}$$

$$\dot{N}_{\text{D2o}} = (1.86 \times 10^{-4} \text{ yr}^{-1}) \eta_{\text{wind, -3}} M_{\text{disc, -2}} \times v_{\text{wind, 0.07}}^{3(5p-4)/4} (S_{\text{th, D2o}}/S_{\text{th, O}})^{-3/2} \quad \text{(C6)}$$

$$\dot{N}_{\text{D3o}} = (1.08 \times 10^2 \text{ yr}^{-1}) \eta_{\text{jet, -4}} M_{\text{disc, -2}} \Gamma_{\text{jet, 1}}^{-1/2} \times (S_{\text{th, D3o}}/S_{\text{th, O}})^{-3/2} \quad \text{(C7)}$$

$$\dot{N}_{\text{D4o}} = (1.88 \times 10^3 \text{ yr}^{-1}) [\eta_{\text{dyn, -2}} M_{\text{dyn, -3}} + \eta_{\text{disc, -3}} M_{\text{disc, -2}}] \Gamma_{\text{jet, 1}}^{-1/2} \times (S_{\text{th, D4o}}/S_{\text{th, O}})^{-3/2} \quad \text{(C8)}$$

where $S_{\text{th,O}} = 3 \times 10^{-20} \text{ Jy Hz}$.

As in Section 5.2, we study the effect of binary parameters on the optical follow-up rates of the afterglow emission from each
ejecta component using Monte Carlo simulations for 300 BHNS binaries with Gaussian $M_{BH} = (7.8 \pm 1.2) M_{\odot}$, $a_{BH} = \text{Uniform}(0, 0.97)$, and the APR4/H4 NS EoS. Fig. B1 shows the contour plots for the optical follow-up rates of the afterglow emission from the dynamical ejecta and wind for the APR4/H4 NS EoS. For the APR4 EoS, the dynamical ejecta optical afterglow has a follow-up rate $6.0 \times 10^{-4} \le \dot{N}_{D10}/\dot{N}_{GW} \le 1.8 \times 10^{-3}$ for $0.75 \le a_{BH} \le 1.0$ and $6.5 \le M_{BH}/M_{\odot} \le 10.0$. The optical follow-up rate for the H4 EoS is larger by a factor of $\sim 4$ compared to that for the APR4 EoS due to more ejecta material, with $2.2 \times 10^{-3} \le \dot{N}_{D10}/\dot{N}_{GW} \le 6.4 \times 10^{-3}$ for $0.65 \le a_{BH} \le 1.0$ and $7.0 \le M_{BH}/M_{\odot} \le 11.0$ in the region $6.7a_{BH} + M_{BH} \ge 15.0$. The afterglow emission from the wind has a very small follow-up rate $2.6 \times 10^{-9} \le \dot{N}_{D20}/\dot{N}_{GW} \le 7.2 \times 10^{-9}$ for $0.55 \le a_{BH} \le 0.85$, $5.0 \le M_{BH}/M_{\odot} \le 9.5$, and the APR4 NS EoS in the region $22.5a_{BH} - M_{BH} \ge 8.5$. For the stiffer H4 EoS, the corresponding rate is also considerably small with $2.8 \times 10^{-9} \le \dot{N}_{D20}/\dot{N}_{GW} \le 8.2 \times 10^{-9}$ for $0.5 \le a_{BH} \le 0.8$ and $4.5 \le M_{BH}/M_{\odot} \le 10.0$ in the region where $15.0a_{BH} - M_{BH} \ge 1.5$. As in the case of radio follow-up of wind afterglows, the significantly

small optical follow-up rates make the detection of wind afterglows especially challenging.

The contour plots for the follow-up rates of the optical afterglow from the sGRB jet and cocoon for the APR4/H4 NS EoSs are shown in Fig. B2. For the APR4 EoS, the relativistic jet afterglow has an optical follow-up rate $1.5 \times 10^{-3} \le \dot{N}_{D30}/\dot{N}_{GW} \le 4.2 \times 10^{-3}$ for $0.6 \le a_{BH} \le 0.9$ and $5.0 \le M_{BH}/M_{\odot} \le 10.0$ in the region $22.5a_{BH} - M_{BH} \ge 8.5$. The corresponding rate is similar for the H4 EoS with $1.8 \times 10^{-3} \le \dot{N}_{D30}/\dot{N}_{GW} \le 4.7 \times 10^{-3}$ for $0.5 \le a_{BH} \le 0.8$ and $4.5 \le M_{BH}/M_{\odot} \le 10.0$ in the region where $15.0a_{BH} - M_{BH} \ge 1.5$. The cocoon afterglow has an optical follow-up rate $6.0 \times 10^{-3} \le \dot{N}_{D40}/\dot{N}_{GW} \le 1.6 \times 10^{-2}$ with $0.7 \le a_{BH} \le 1.0$ and $5.5 \le M_{BH}/M_{\odot} \le 10.0$ for the APR4 EoS. For the H4 EoS, the optical follow-up rate is larger by a factor of $\sim 2$ with $1.0 \times 10^{-2} \le \dot{N}_{D40}/\dot{N}_{GW} \le 2.8 \times 10^{-2}$ for $0.6 \le a_{BH} \le 0.8$ and $4.5 \le M_{BH}/M_{\odot} \le 10.0$ in the region where $27.5a_{BH} - M_{BH} \ge 11.5$.

This paper has been typeset from a TeX/LA TeX file prepared by the author.