© R. Duque et al. 2019

1. Introduction

The first detection of the gravitational waves (GW) from the inspiral phase of a binary neutron star merger (Abbott et al. 2017, 2019) was followed by all three electromagnetic counterparts expected after the coalescence: a short gamma-ray burst (GRB; Goldstein et al. 2017; Savchenko et al. 2017) and its multi-wavelength afterglow (AG) in the X-ray (Haggard et al. 2017; Margutti et al. 2017; Troja et al. 2017; D’Avanzo et al. 2018), radio (Hallinan et al. 2017) and optical bands (Lyman et al. 2018), and an optical-IR thermal transient source (Evans et al. 2017; Valenti et al. 2017; Coulter et al. 2018). The latter transient allowed for the localization of the event with sub-arcsecond precision in the S0-type galaxy NGC 4993 at a distance of 40.7 ± 3.3 Mpc (Palmese et al. 2017; Cantiello et al. 2018). This thermal transient showed evidence of heating as a result of r-process nucleosynthesis (Tanvir et al. 2017; Cowperthwaite et al. 2017; Pian et al. 2017; Smartt et al. 2017), classifying it as a “kilonova” (KN), the first to feature such detailed observations.

According to estimates on joint GW-GRB events by Beniamini et al. (2019), such a combined detection of GW-GRB-AG-KN should remain rare. Indeed, GRB170817A would not have been detected at a distance larger than 50 Mpc (Goldstein et al. 2017) or from a viewing angle larger than ∼25°.

In the context of the present O3 and future observing runs of the LIGO-Virgo Collaboration (LVC), we would like to know what to expect regarding the afterglows of future binary neutron star (BNS) mergers: the rates of GW and joint GW-AG events; the distributions of event observables, such as distance D, viewing angle θ_v, radio afterglow time of peak t_p and peak flux F_p, and remnant proper motion μ; and the sensitivity of these to the population’s characteristics and to the multi-messenger detector configuration.

In this paper, we approach these questions from the angle of a population study. More precisely, we start from the basis of prior knowledge on binary neutron star mergers from short GRB science and observations of GW170817 and its counterparts. This allows us to prescribe some intrinsic variability in mergers of the population, such as their kinetic energy, external medium density, shock microphysical conditions, etc. We then calculated the radio afterglows arising from this population of mergers in the hypothesis of a jet-dominated afterglow emission. Finally, we applied the thresholds of current and future GW and radio detectors to determine the events within the population which would be detectable in both the GW and radio domains. By describing and studying this jointly-detectable population, we make predictions for observing runs in the future.

Saleem et al. (2018a,b) consider such a population model for short GRB afterglows, with detailed numerical models for the multi-band afterglow and GW signals, and based on uninformative (e.g., uniform or log-uniform) event parameter distributions. Moreover, the discussion therein is centered around the parameter-space constraints that can be derived from detections or non-detections of the afterglow in various electromagnetic bands, and on the conditions necessary to observe afterglows from off-axis lines of sight. Here, we discuss, rather, the effect of the population parameters and detector configurations on the expected population of afterglows.

More recently, Gottlieb et al. (2019) also follows this approach. Once again, our study differs in that it accounts for some prior knowledge of short GRBs (such as their luminosity function) and observations from GW170817 to sharpen our predictions for future events. This allows for a detailed study of the impact of the population parameters on these predicted observations.

Our study relies on the following assumptions.

• (i)

  In most cases the merger produces a successful central jet. This is clearly a strong assumption that will have to be validated by future observations. Recent studies based on the population of short GRBs (Beniamini et al. 2019) and the dynamics of jets interacting with merger ejecta (Duffell et al. 2018) suggest that successful jets could be common for these phenomena.

• (ii)

  Around its peak, the afterglow is dominated by the contribution of the jet’s core. In this case, the afterglow peak flux depends on the jet’s core parameters (core isotropic kinetic energy E_iso, c and opening angle θ_j, shock microphysics parameters ϵ_e, ϵ_B, and p), on its environment (external medium density n), and on its viewing conditions (the distance and viewing angle). The population’s distributions of most of these quantities remain quite uncertain. Our knowledge of these relies on the afterglow fitting of a limited sample of short GRB afterglows with known distance (e.g., Berger 2014; Fong et al. 2015), on the expectation that short GRBs generally occur in a low-density environment, and on the results obtained from the (up to now) single multi-messenger event GRB170817A (e.g., van Eerten et al. 2019).

This publication is organized as follows. In Sect. 2 we describe our methods to calculate the afterglow emission and detail the input parameters of our population model. In Sect. 3 we provide the criteria for GW and radio afterglow detection that we applied to deduce the jointly-detectable population. In Sect. 4 we report general results on this population, describe its observable characteristics, and study their sensitivity to the detector configurations and to the population model parameters. In Sect. 5, we discuss our results in the context of future multi-messenger campaigns and how they may help to interpret observations therein.

2. A population model for binary neutron star merger afterglows

2.1. Radio afterglows from structured jets

The multi-wavelength afterglow of GW170817 is associated with the deceleration of a structured relativistic jet emitted by the central source formed in the merger, and more precisely, to the synchrotron emission of electrons accelerated by the forward shock propagating in the external medium. It was observed for more than 300 days. Its time evolution is similar at all wavelengths, indicating that the emission was produced in the same spectral regime of the slow-cooling synchrotron process, that is, ν_m <  ν <  ν_c, and that the non-thermal distribution of shock-accelerated electrons assumed a non-evolving slope p  ∼  2.2 (Nakar et al. 2018; Mooley et al. 2018a; van Eerten et al. 2019). The observed slow rise of the afterglow is clear evidence of the lateral structure of the outflow, which was observed off-axis with a viewing angle of θ_v  ∼  15−25°. Initially, the observed flux is dominated by the beamed emission from mildly-relativistic, mildly-energetic material coming towards the observer (Mooley et al. 2018b). As the jet decelerates, relativistic beaming becomes less efficient and regions closer to the jet axis start to contribute to the flux. The peak of the afterglow is reached when the jet’s core, which is open to a few degrees at the angle θ_j, is revealed. This occurs when the jet’s core’s Lorentz factor has decreased to Γ ∼ 1/θ_v.

An alternative possibility is a quasi-spherical, mildly relativistic ejecta with a steep radial structure (Kasliwal et al. 2017; Gottlieb et al. 2018; Hotokezaka et al. 2018; Mooley et al. 2018b; Gill & Granot 2018). This is now strongly disfavored by radio very long baseline interferometry (VLBI) observations, which demonstrate an apparent superluminal motion of the remnant with β_app ∼ 4 (Mooley et al. 2018c) and constrain the angular size of the radio source to below 2 mas (Ghirlanda et al. 2019), thus confirming the emergence of a relativistic core jet at the peak of the afterglow.

The lateral structure of the jet is constrained by observations: the best fits are obtained for a sharp decay of the kinetic energy and Lorentz factor, with the angle θ from the jet axis (Lamb & Kobayashi 2017; Lazzati et al. 2017; Troja et al. 2018; Resmi et al. 2018; Margutti et al. 2018). Typically, for a top-hat core jet with an opening angle θ_j, an isotropic equivalent kinetic energy E_iso, c = 4πϵ_c and Lorentz factor Γ_c, surrounded by a power-law-structured sheath with the following kinetic energy per solid angle ϵ(θ) and Lorentz factor Γ₀(θ),

$$\begin{array}{c}\hfill ϵ(\theta )={ϵ}_{\mathrm{c}}\{\begin{array}{cc}1\hfill & \mathrm{if}\theta \le {\theta }_{\mathrm{j}}\\ {(\theta /{\theta }_{\mathrm{j}})}^{-a}\hfill & \mathrm{if}\theta \ge {\theta }_{\mathrm{j}}\end{array}\end{array}$$(1)

and

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_0(\theta )=1+({\mathrm{\Gamma }}_{\mathrm{c}}-1)\{\begin{array}{cc}1\hfill & \mathrm{if}\theta \le {\theta }_{\mathrm{j}}\\ {(\theta /{\theta }_{\mathrm{j}})}^{-b}\hfill & \mathrm{if}\theta \ge {\theta }_{\mathrm{j}}\end{array},\end{array}$$(2)

observations require steep slopes of a, b ≳ 2 (Gill & Granot 2018; Ghirlanda et al. 2019).

An example of the radio light curves obtained for such a structured jet is plotted in Fig. 1 (left) using parameters typical of the various fits to the data that have been published (see e.g., Gill & Granot 2018; van Eerten et al. 2019; Ghirlanda et al. 2019). To obtain this light curve, the dynamics of the material at different latitudes are computed independently. At a latitude of θ, the deceleration radius is $R_{\mathrm{dec}}(\theta )={\left(\frac{3ϵ(\theta )}{{\mathrm{\Gamma }}_0^2(\theta )n{m}_{\mathrm{p}}{c}^2}\right)}^{1/3}$, which is constant at the core and slowly increases with θ outside the core (R_dec(θ)∝θ^{(2b − a)/3}). For the values of a and b chosen in Fig. 1, $\frac{2b-a}{3}=0.17$.

Subsequently, the synchrotron emission of shock-accelerated electrons in the shock comoving frame is assumed to follow the standard synchrotron slow-cooling spectrum (Sari et al. 1998) including self-absorption. Finally, the observed light curve is computed by summing the contributions of all latitudes on equal-arrival time surfaces, taking relativistic beaming and Doppler boosting into account.

The separate contributions of the core jet and the sheath are plotted in Fig. 1 (left) and the emergence of the core at the peak is clearly visible. The evolution of the peak flux of the same structured jet is plotted in Fig. 1 (right) as a function of the viewing angle, as well as the ratio of the peak flux to the peak flux from the core jet only. Interestingly, this ratio is almost constant (∼1.5), except for the on-axis cases (θ_v ≤ θ_j) when the core is more dominant.

Fig. 1. Left: 3 GHz light curve of afterglow from structured jet at distance D = 42 Mpc, with sharp power-law structure (a = 4.5, b = 2.5; see Eqs. (1)–(2)), viewing angle θv = 22°, external density n = 3 × 10−3 cm−3, core jet with θj = 4°, Γc = 100, Eiso, c = 2 × 1052 erg, and with microphysics parameters p = 2.2, ϵe = 0.1 and ϵB = 10−4. The radio observations of GW170817 are also plotted and are compiled from Hallinan et al. (2017), Alexander et al. (2018), Margutti et al. (2018), Mooley et al. (2018a,b), Dobie et al. (2018). The respective contributions of the core jet (θ ≤ θj), the sheath and the total are plotted in dashed, dotted and solid red lines. The flux of the core jet computed with the simplified treatment described by Eq. (4) is also plotted in dashed blue line for comparison. Bottom right: peak flux of same structured jet as function of viewing angle θv (solid red line), peak flux of light curve from core jet only (dashed red line) or from sheath only (dotted red line), and peak flux of core jet as computed using scaling law in Eq. (6) (thin black line). Upper right: ratio of peak flux from whole outflow (core jet+sheath) to that of core jet only.

As long as the kinetic energy and the Lorentz factor decay steeply with θ, the core jet is expected to dominate the flux at the peak whatever the viewing angle, even if the precise value of the total-to-core ratio at the peak may slightly vary depending on the details of the assumed lateral structure. Therefore, as our population model takes only the properties of the afterglow at the peak into account, that is, when it can be detected more easily, it is enough in the following to compute only the contribution from the core jet, while keeping in mind that it may slightly underestimate the total peak flux by a factor of¹ ∼1.5.

The contribution from a top-hat core jet can efficiently be computed using the approximation suggested by Granot et al. (2002) to avoid the full integration over equal-arrival time surfaces. More precisely, we computed the flux $F_{\nu }^{\text{iso}}(t_{\text{obs}})$ of a spherical ejecta with initial Lorentz factor Γ_c and kinetic energy E_iso, c as in Panaitescu & Kumar (2000), and we corrected it by defining the on-axis afterglow flux including the jet break:

$$\begin{array}{c}\hfill F_{\nu }^{\mathrm{on}}(t_{\mathrm{obs}})=F_{\nu }^{\mathrm{iso}}(t_{\mathrm{obs}})\times \{\begin{array}{cc}1\hfill & \mathrm{if}\mathrm{\Gamma }{\theta }_{\mathrm{j}}\le 1\hfill \\ {\left(\mathrm{\Gamma }{\theta }_{\mathrm{j}}\right)}^{-2}\hfill & \mathrm{if}\mathrm{\Gamma }{\theta }_{\mathrm{j}}\ge 1\hfill \end{array}\end{array}$$(3)

and by calculating the afterglow from any viewing angle as:

$$\begin{array}{c}\hfill F_{\nu }(t_{\mathrm{obs}})=\{\begin{array}{cc}{F}_{\nu }^{\mathrm{on}}(t_{\mathrm{obs}})\hfill & \mathrm{if}{\theta }_v<{\theta }_j(\mathrm{on}\mathrm{axis})\hfill \\ a^3{F}_{\nu /a}^{\mathrm{on}}(b{t}_{\mathrm{obs}})\hfill & \mathrm{if}{\theta }_v>{\theta }_j(\mathrm{off}\mathrm{axis})\hfill \end{array}\end{array}$$(4)

with $a=\frac{1-\beta }{1-\beta cos({\theta }_{\mathrm{v}}-{\theta }_{\mathrm{j}})}$ and $b=\frac{1-\frac{R}{\mathit{ct}}}{1-\frac{R}{\mathit{ct}}cos({\theta }_{\mathrm{v}}-{\theta }_{\mathrm{j}})}$ where t is the source frame time. The latter correction, corresponding to the ratio of on-axis and off-axis arrival times is not included in Granot et al. (2002) and is found to improve the quality of the approximation. The corresponding light curve is plotted in Fig. 1 (left, blue dashed line) and corresponds well with the exact calculation.

2.2. An analytical model for afterglow peak properties

This afterglow calculation is standard for top-hat jets. Therefore analytical expressions for the afterglow properties at the peak are available. They depend slightly on the assumptions for a possible late jet lateral expansion. If lateral expansion is taken into account like in standard GRB afterglow theory, the peak flux of the radio light curve scales as (Nakar et al. 2002):

$$\begin{array}{c}\hfill F_{\mathrm{p},\nu }\propto E_{\mathrm{iso},\mathrm{c}}{\theta }_{\mathrm{j}}^2{n}^{\frac{p+1}{4}}{ϵ}_{\mathrm{e}}^{p-1}{ϵ}_{\mathrm{B}}^{\frac{p+1}{4}}{\nu }^{\frac{1-p}{2}}{D}^{-2}max{({\theta }_{\mathrm{j}},{\theta }_{\mathrm{v}})}^{-2p},\end{array}$$(5)

as long as the spectral regime remains ν_m <  ν <  ν_c. This leads to the following expression² of the peak flux at 3 GHz:

$$\begin{array}{cc}\hfill F_{\mathrm{p},3\mathrm{GHz}}=& 8.6E_{52}{\theta }_{\mathrm{j},-1}^2{n}_{-3}^{4/5}{ϵ}_{\mathrm{e},-1}^{6/5}{ϵ}_{\mathrm{B},-3}^{4/5}\hfill \\ \hfill & \times D_{100}^{-2}max{({\theta }_{\mathrm{j},-1},{\theta }_{\mathrm{v},-1})}^{-4.4}\mathrm{mJy},\hfill \end{array}$$(6)

where E₅₂ = E_iso, c/10⁵² erg, θ_j, −1 = θ_j/0.1 rad, n₋₃ = n/(10⁻³ cm⁻³), ϵ_e, −1 = ϵ_e/10⁻¹, ϵ_B, −3 = ϵ_B/10⁻³, θ_v, −1 = θ_v/0.1 rad and D₁₀₀ = D/(100 Mpc), and where we assume p = 2.2. This scaling law is plotted in Fig. 1 (right, solid black line) and corresponds well with the detailed calculation. The expression of the corresponding peak time is also given by Nakar et al. (2002), assuming lateral expansion of the jet. However, even if late observations of GRB170817A may give some evidence for the lateral expansion of the core jet (for example the temporal decay index being ∼ − p, as found by Mooley et al. 2018a), it is not clear if this expansion should be as strong for a core jet embedded in a sheath as for a “naked” top-hat jet. Therefore, we also consider a limit case without lateral expansion. This only affects the peak flux slightly and Eq. (6) above remains a good approximation in both cases. On the other hand, it modifies the peak time, leading to:

$$\begin{array}{cc}\hfill t_{\mathrm{p},\mathrm{no}\mathrm{ex}}& \sim t_{\mathrm{b}}\times {\left(\frac{{\theta }_{\mathrm{v}}}{{\theta }_{\mathrm{j}}}\right)}^{8/3}\hfill \\ \hfill & =137{\left(\frac{E_{52}}{n_{-3}}\right)}^{1/3}{\theta }_{v,.35}^{8/3}\mathrm{days}\hfill \end{array}$$(7)

and

$$\begin{array}{cc}\hfill t_{\mathrm{p},\mathrm{ex}}& \sim t_{\mathrm{b}}\times {\left(\frac{{\theta }_{\mathrm{v}}}{{\theta }_{\mathrm{j}}}\right)}^2\hfill \\ \hfill & =60{\left(\frac{E_{52}}{n_{-3}}\right)}^{1/3}{\theta }_{j,-1}^{2/3}{\theta }_{v,.35}^2\mathrm{days}\hfill \end{array}$$(8)

where t_b is the standard jet break time,

$$\begin{array}{c}\hfill t_{\mathrm{b}}=4.9{(E_{52}/n_{-3})}^{1/3}{\theta }_{j,-1}^{8/3}\mathrm{days}.\end{array}$$(9)

The ratio of the peak times without and with lateral expansion is simply

$$\begin{array}{c}\hfill \frac{t_{\mathrm{p},\mathrm{without}}}{t_{\mathrm{p},\mathrm{with}}}\sim {\left(\frac{{\theta }_{\mathrm{v}}}{{\theta }_{\mathrm{j}}}\right)}^{2/3}·\end{array}$$(10)

In the following, we use, unless otherwise noted, Eqs. (6)–(8) to estimate the peak flux and peak time of the radio afterglow from a binary neutron star merger.

2.3. Physical parameters of the jet and their intrinsic distribution

To generate a population of jet afterglows we had to fix the various parameters appearing in Eqs. (6)–(8). We first defined a set of “fiducial distributions” for these parameters. These are summarized in Table 1 and have been chosen as follows.

Core jet isotropic equivalent kinetic energy E_iso, c: we adopted a broken power-law distribution, which we deduced directly from the gamma-ray luminosity function of cosmological short GRBs, assuming a standard rest frame duration ⟨τ⟩ = 0.2 s and a standard efficiency f_γ = 0.2 (Beniamini et al. 2016) so that E_iso, c ∼ L_γ ⟨τ⟩/f_γ, leading to a density of probability

$$\begin{array}{cc}\hfill \varphi \left(E_{\mathrm{iso},\mathrm{c}}\right)& =\frac{1}{N}\frac{\mathrm{d}N}{\mathrm{d}{E}_{\mathrm{iso},\mathrm{c}}}\hfill \\ \hfill & \propto \{\begin{array}{cc}{E}_{\mathrm{iso},\mathrm{c}}^{-{\alpha }_1}\hfill & \mathrm{for}{E}_{\min }\le E_{\mathrm{iso},\mathrm{c}}\le E_{\mathrm{b}}\hfill \\ E_{\mathrm{iso},\mathrm{c}}^{-{\alpha }_2}\hfill & \mathrm{for}{E}_{\mathrm{b}}\le E_{\mathrm{iso},\mathrm{c}}\le E_{\max }\hfill \end{array}\hfill \end{array}$$(11)

where ϕ is normalized to unity. We fixed E_min = 10⁵⁰ erg and E_max = 10⁵³ erg. Luminosity functions for short GRBs were deduced based on observations from several teams (e.g., D’Avanzo et al. 2014; Wanderman & Piran 2015; Ghirlanda et al. 2016). For our fiducial model, we adopted that of Ghirlanda et al. (2016) with α₁ = 0.53, α₂ = 3.4 and a break luminosity of L_b = 2 × 10⁵² erg s⁻¹ leading to E_b = 2 × 10⁵² erg.

Jet opening angle θ_j: for simplicity we assumed a single value θ_j = 0.1 rad, which appears to be representative of the typical value found by fitting the afterglow light curve of GW170817 (e.g., Gill & Granot 2018; van Eerten et al. 2019), by the VLBI observations of the remnant (Mooley et al. 2018c) and also consistent with the results of prior short GRB afterglow fitting (Fong et al. 2015) and short GRB-BNS merger rate comparisons (Beniamini et al. 2019).

External density n: most BNS mergers are expected to occur in low density environments due to a long merger time, as assessed by the significant offsets of short GRB sources from their host galaxies (e.g., Nakar 2007; Troja et al. 2008; Fong et al. 2010; Church et al. 2011; Berger 2014). In the case of GW170817, afterglow fitting typically leads to n ∼ 10⁻³ cm⁻³, which is in agreement with the estimation of the HI content of the host galaxy NGC 4996 (Hallinan et al. 2017). The external densities observed in short GRBs cover the interval 10⁻³ − 1 cm⁻³ (Berger 2014), but this is probably biased towards high densities, which favor afterglow detection. In this paper we consider n ∼ 10⁻³ cm⁻³ as typical and set, therefore, a log-normal of mean 10⁻³ cm⁻³ with a standard deviation of 0.75 as the fiducial density distribution.

Microphysics parameters: we set the value of ϵ_e at 0.1. This is generally done in GRB afterglow modeling (Wijers & Galama 1999; Panaitescu & Kumar 2001; Santana et al. 2014; Beniamini & van der Horst 2017; D’Avanzo et al. 2018). Also, this value has been adopted or found in the fitting of the afterglow of GW170817 (e.g., van Eerten et al. 2019; Margutti et al. 2018; Nakar et al. 2018; Xie et al. 2018), and it is consistent with particle shock-acceleration simulations (Sironi & Spitkovsky 2011). Similarly, we set p at 2.2, which is also a common value used in GRB afterglow models and which is in agreement with shock-acceleration theory in the relativistic regime (e.g., Sironi et al. 2015). In the case of GW170817 this value can be directly measured from multi-wavelength observations (Mooley et al. 2018a,b; van Eerten et al. 2019). Finally, more diversity is generally expected for ϵ_B. We assumed a log-normal distribution similar to the one adopted for n with the additional constraint that ϵ_B is restricted to the interval [10⁻⁴, 10⁻²].

In Sect. 4.2 we consider possible alternatives to these fiducial distributions. In particular, we discuss the impact of substantial uncertainties in the distribution of kinetic energy, ϕ(E_iso, c), either due to the difficulty in measuring the luminosity function of short GRBs, or in simplifying assumptions regarding the duration and the efficiency of the prompt GRB emission. We also explore different values of the jet opening angle θ_j. As for the external density, we examine the effect of a change in the mean circum-merger density of the population of mergers.

3. From the intrinsic to the observed population: GW and radio joint detection

3.1. GW detection criterion

The signal-to-noise ratio (S/N, denoted by ρ) of a GW inspiral signal in a single LIGO-Virgo-type interferometer can be written as ${\rho }^2=\frac{{\mathrm{\Theta }}^2}{D^2}{\mathcal{M}}^{5/3}{S}_{\mathrm{I}}$ (Finn & Chernoff 1993), where D is the luminosity distance to the binary, ℳ is the chirp mass of the binary, S_I is a quantity depending only on the sensitivity profile of the interferometer, and Θ² represents the dependence of the signal to the binary sky position (θ, ϕ) and orientation (θ_v, ξ) with respect to the plane of the interferometer. Again, θ_v is the angle of the line of sight to the binary polar direction, which we also assumed is the jet-axis direction. The function Θ² admits a global maximum of ${\Theta }_M^2=16$ corresponding to an optimally positioned and oriented source (binary at the zenith and polar axis orthogonal to the instrument’s arms). This optimal binary can be detected out to a distance known as the horizon H of the instrument, which depends on the S/N threshold for detection, taken to be ρ₀ = 8 for the LVC network. We may thus rewrite the S/N in the following manner:

$$\begin{array}{c}\hfill {\rho }^2={\rho }_0^2\frac{{\mathrm{\Theta }}^2}{{\mathrm{\Theta }}_M^2}\frac{H^2}{D^2}·\end{array}$$(12)

In a full population study, one would have to draw all four angles and evaluate this criterion in every case. We chose to reduce the number of parameters to two: D and θ_v, as these are the ones relevant to the afterglow. Thus, we were compelled to review this criterion by averaging Θ² on the sky-position (θ, ϕ) and polarization angle ξ. This is done readily from the expression given in Finn & Chernoff (1993; Eq. (3.31)) and it is found analytically that:

$$\begin{array}{c}\hfill {⟨{\mathrm{\Theta }}^2⟩}_{\xi ,\theta ,\varphi }=\frac{1}{4\pi 2\pi }{\displaystyle \int }\mathrm{d}\xi \mathrm{d}\mathrm{\Omega }{\mathrm{\Theta }}^2=\frac{4}{5}(1+6{cos}^2{\theta }_{\mathrm{v}}+{cos}^4{\theta }_{\mathrm{v}})·\end{array}$$(13)

Hence, we chose to use the following criterion to determine those binaries of inclination θ_v and distance D which are detectable in GW on average in the sky:

$$\begin{array}{c}\hfill {⟨{\rho }^2⟩}_{\xi ,\theta ,\varphi }>{\rho }_0^2,\end{array}$$(14)

which is:

$$\begin{array}{c}\hfill \sqrt{\frac{1+6{cos}^2{\theta }_{\mathrm{v}}+{cos}^4{\theta }_{\mathrm{v}}}{8}}>\frac{D}{\overline{H}},\end{array}$$(15)

where we have denoted $\overline{H}=\sqrt{\frac{2}{5}}H\sim H/1.58$ for the sky-position-averaged horizon.

Averaging Θ² on inclination angle θ_v and imposing an S/N threshold as in Eq. (14) would lead to the simple detection criterion of D <  R, where R is the range of the instrument, that is, the maximum distance at which a binary can be detected on average in sky-position and in orientation. The range is linked to the horizon by H = 2.26R.

The criterion of Eq. (15) is valid for a detection using a single instrument. Instead, GW detection by the interferometer network is based on multi-instrument analysis, and is, thus, more complicated than the one described by our criterion. Furthermore, as was illustrated in the case of GW170817, true joint GW-AG detection requires pinpointing the source through the detection of the kilonova counterpart³, which in turn involves a small enough localization map from the GW data, and our criterion should incorporate this. We chose a simple localization criterion by supposing that a source is localized if it is detectable (once again on average) by the two most sensitive interferometers of the network. Thus our detection-localization criterion is that of Eq. (15), with the horizon of the second most sensitive instrument of the network, that is, the LIGO-Hanford instrument. Table 2 reports the corresponding values taken for our study for the LVC O2, O3 runs, and for the design interferometers.

In this framework, it can be shown that assuming homogeneity of the sources within the horizon and isotropy of the binary polar direction, the mean viewing angle of GW-detected events is ∼38° and the fraction of GW-detected events among all mergers is ∼29%, both regardless of the horizon value. This latter fraction is, thus, an absolute maximum for the fraction of mergers to be jointly detected in both radio and GW channels.

3.2. Radio detection criterion

The detection criterion in the radio band is simply for the 3 GHz peak flux, as determined by our models, to be larger than the sensitivity of the radio array available for follow-up at the time of consideration, which we denote by s. We note that this is a criterion aimed rather at detectability than detection. True detection only occurs if observations are pursued for long enough after the GW merger signal, as the flux rises to its maximum. We also note that events which are marginally detectable would only provide poor astrophysical output because the inference of jet parameters requires the fitting of an extended portion of the radio afterglow light curve.

We took three typical radio sensitivities reflecting the present and future capabilities of radio arrays: (i) the Karl Jansky Very Large Array (VLA) at s = 15 μJy, (ii) the phase I of the Square Kilometer Array (SKA1, Braun 2014), at s = 3 μJy, and (iii) phase II of the SKA (SKA2) and the Next Generation VLA (ngVLA, Selina et al. 2018), both at s = 0.3 μJy. The combination of this radio detection criterion with that of the GW domain (Eq. (15)) constituted our criterion for the joint detection of mergers.

In the sequel, we referred to a jointly GW and radio detectable population by using the names of the corresponding LVC run and limiting radio facility, as in “the joint population of the O3 and VLA configuration”. The corresponding GW interferometer horizons and radio sensitivities can be found respectively in Table 2 and in the paragraph above.

4. Results: detection rates and properties of the detected population

4.1. Fiducial model

In this section, we describe the population of events detected jointly in the GW and radio domains. For this purpose, we used a Monte Carlo approach in which we simulated N >  10⁶ binary systems within the sky-averaged horizon $\overline{H}$ using the parameter distributions of the fiducial model given in Sect. 2.3. By applying the detection criteria described in Sect. 3, we obtained N_GW systems detected by the GW interferometer network and N_joint systems which are also detectable by radio telescopes by their afterglows. Then we studied the properties of the sample of joint GW and radio detectable events.

4.1.1. Rate of joint GW and radio detectable events

Starting from the fraction f_GW = N_GW/N ∼ 29% of binary systems detected by the LVC (Sect. 3.1), we then estimated which fraction N_joint/N_GW of those events would also produce a detectable radio afterglow.

This fraction is shown in Fig. 2 (left) as a function of the GW horizon and for the VLA, SKA1 and SKA2/ngVLA limiting sensitivities. In the design configuration of the LVC network, 75% (resp. 2 times) more events are detectable by SKA1 (resp. SKA2/ngVLA) than by the VLA.

Fig. 2. Left: detectable fraction of radio afterglows among gravitational wave events as function of the horizon distance $H=1.58\overline{H}$. Right: expected number of joint detections normalized to case of O3-VLA configuration. In both panels the full (resp. dashed, resp. dotted) lines correspond to the VLA (resp. SKA1, resp. SKA2/ngVLA) being the limiting radio facility.

The number of joint event detections normalized to the O3-VLA configuration is represented in Fig. 2 (right). The number of joint detections behaves approximately as H^α, with α <  3 because of the reduction in radio detection efficiency when the distance is increased (Fig. 2, left). We find α ≃ 2.4 (resp. 2.6, resp. 2.8) for the VLA (resp. SKA1, resp. SKA2/ngVLA).

As illustrated in Fig. 2 (left), we find that the fraction of detected events decreases from 43% (O2) to 31% (O3) with the VLA limiting sensitivity. However the absolute number of detections increases (Fig. 2, right). With the design-SKA1 configuration, a fraction of 35% of GW events can lead to a radio detection. With regard to the much higher sensitivities that may be reached after 2030, more than a half of GW events may become detectable in the radio band. Indeed we found a percentage of 63% of detectable joint events in the design-SKA2/ngVLA configuration.

For the purposes of a comprehensive survey, the number of GW and jointly detected events per continuous year of GW network operation for future detector configurations can be found in Table 3, where we have used the horizons of Table 2 and our fiducial model.

4.1.2. Distance and viewing angle

The distributions of distances are shown in Fig. 3 for the O2-VLA and O3-VLA configurations, as well as for the entire GW population. In this figure, the distributions are not normalized to unity so as to show the decrease in the fraction of joint event as the horizon increases. It can be observed that as a result of the GW and radio detection thresholds, sources are progressively lost when the distance increases so that this distribution exhibits a near linear increase (as opposed to the dN/dD ∝ D² of the intrinsic homogeneous population). Up to $D=\overline{H}/\sqrt{8}\simeq 0.35\overline{H}$, all events are detected in GW, regardless of their orientation⁴ and the only selection is due to the radio sensitivity. Near the horizon, the maximum viewing angle allowing GW-detection is ${\theta }_{\max }\propto \sqrt{1-D/\overline{H}}$, which produces a strict decrease of event density.

Figure 4 shows the distribution of the viewing angles. For the joint events, the peak of the distribution takes place at small viewing angles θ_v ∼ 15−20° as a result of the rapid decline of the peak flux with angle ($F_{\mathrm{p}}\propto {\theta }_{\mathrm{v}}^{-2p}$), and about 39% (resp. 47%) of the events would be seen with a viewing angle θ_v <  20° for O2-VLA (resp. O3-VLA), the angle at which the GW170817 event was likely to have been seen. As shown in Fig. 5 (left), the mean viewing angle of the jointly observed events strongly depends on the GW horizon and the radio sensitivity. It appears that even if the radio sensitivity is not substantially improved while the GW interferometers reach their design horizons, a significant number of events would still be observable off-axis. The increasing and decreasing phases of the afterglows of these events would allow for an improved study of the jet structure, which is better revealed by off-axis events. Figure 5 (right) shows the fraction of on-axis events (θ_v ≤ θ_j) within the joint detections. It increases from 4.0% (O2-VLA) to 5.3% (O3-VLA).

Fig. 3. Left: differential distribution of distances (normalized to the sky-position-averaged horizon) of events in total GW population (full line), and detectable jointly in O2-VLA (dotted line) and O3-VLA (dashed line) configurations. Right: cumulative distribution of same events. Distributions have been normalized to the fraction of jointly detected events among all GW-detected events, see text for details.

4.1.3. Peak times, peak fluxes and synchrotron spectral regime

The distributions of peak times and peak fluxes in the O3-VLA configuration are represented in Fig. 6. The fraction of events with a peak time of fewer than 150 days (as observed in GRB 170817A) is 55% without jet expansion and 81% with lateral expansion. The distribution in peak flux is shown for all sources which are detected in gravitational waves within the sky-position-averaged horizon of O2, O3 and design instruments. It appears clear that for the present VLA sensitivity, most radio afterglows cannot be detected. This explains why improving the sensitivity of the future ngVLA would have such an impact on the joint detection rates, as shown in Fig. 2.

Fig. 4. Same as Fig. 3, for viewing angle. Distributions are normalized.

Fig. 5. Left: mean viewing angle of jointly detected events, as function of horizon and radio sensitivity s. Values of H for O2, O3 and design instruments, as well as s for VLA, SKA1 and SKA2/ngVLA, are indicated as dashed lines. Right: fraction of on-axis events (defined by θv ≤ θj) among jointly detected events in same H–s diagram.

Fig. 6. Left: cumulative distributions of peak times of radio afterglows for jointly detected events in O3-VLA configuration, assuming jet with (solid line) or without (dashed line) lateral expansion. Middle: differential distribution of peak fluxes of GW-detected events assuming horizon for O2 (dotted line), O3 (dashed line) and design (full line) configurations of GW interferometers. The radio sensitivities for VLA, SKA1 and SKA2/ngVLA are indicated for comparison. Right: differential distribution of maximum proper motion for jointly detected events assuming O2-VLA (dotted line), O3-VLA (dashed line), design-SKA1 (full line) configurations. The dotted vertical line shows the lower limit for detection of the proper motion.

The scaling law used to compute the radio emission at the peak assumes that the observing frequency remains in the same spectral regime of the slow-cooling synchrotron spectrum, that is, ν_m <  ν <  ν_c, and above the absorption frequency. Using our more detailed calculation of the radio afterglow from the core jet (Eq. (4) and accompanying text), we can ascertain if this condition is fulfilled for the bulk of the population. However, for mergers in high density environments (larger than 10 cm⁻³), this condition is no longer met as the absorption frequency ν_a and the injection frequency ν_i may be larger than the radio frequency at early times. For an event with n = 10 cm⁻³ and θ_v ∼ 40° (the mean viewing angle of the GW-detected population), the radio frequency typically meets the injection break at around 30 days (resp. before 10 days). In this case, chromatic light curves are expected.

4.1.4. Proper motion

Another independent observable is the proper motion of the remnant, which can be readily accessible to VLBI measurements as illustrated in the case of GRB170817A (Mooley et al. 2018c; Ghirlanda et al. 2019). However, a VLBI network can only measure the proper motion of a source if the total displacement of the source over the time it is detectable by the array is larger than the angular resolution of the network. For our purposes, we chose to simplify this criterion and adopt a dual condition on the flux (F_p > s_VLBI), as well as on the instantaneous proper motion μ_max (requiring that μ_max > μ_lim) of the remnant at the peak of the afterglow. Here, we adopted a limiting proper motion of μ_lim = 3 μas/day, in coherence with the uncertainties quoted in Mooley et al. (2018c), Ghirlanda et al. (2019).

An upper limit of the proper motion was obtained based on the assumption that the core jet contribution dominates at the peak of the radio light curve, and that the shocked region at the front of this core jet constitutes the visible remnant of the merger. In first approximation, afterglows from jets peak when the observer enters the focalization cone of the radiation, that is, Γ ∼ 1/θ_v. For an off-axis observer, we estimated, therefore, the maximum proper motion to be:

$$\begin{array}{c}\hfill {\mu }_{\max }=\frac{c{\beta }_{\mathrm{app}}}{D}=\frac{c}{D}f({\theta }_{\mathrm{v}})\end{array}$$(16)

with

$$\begin{array}{c}\hfill f({\theta }_{\mathrm{v}})=\frac{\sqrt{1-{\theta }_{\mathrm{v}}^2}sin{\theta }_{\mathrm{v}}}{1-\sqrt{1-{\theta }_{\mathrm{v}}^2}cos{\theta }_{\mathrm{v}}}·\end{array}$$(17)

Thus, for a given value of the viewing angle, there is a maximum distance beyond which the measurement of the proper motion is not possible. It is given by

$$\begin{array}{c}\hfill D_{\max }({\theta }_{\mathrm{v}})=\frac{c}{{\mu }_{\lim }}f({\theta }_{\mathrm{v}})\sim 50f({\theta }_{\mathrm{v}})\mathrm{Mpc}.\end{array}$$(18)

The distribution of μ_max is shown in Fig. 6 (right) for O2, O3 and design-level detectable GW events respecting the VLBI flux criterion. Here, we have take s_VLBI = 15 μJy, which represents a limiting S/N of ∼3 with the level of noise reported in Mooley et al. (2018c), Ghirlanda et al. (2019) for the e-MERLIN and High Sensitivity Array networks. Within the O2-VLA (resp. O3-VLA, resp. design-VLA) populations, the measurement of the proper motion is possible in only 79% (resp. 64%, resp. 42%) of the jointly-detectable events.

Finally, we represent in Fig. 7 various plots connecting two observable quantities: (θ_v, D), (θ_v, F_p, 3 GHz), (θ_v, t_p) and (D, μ_max). The grey dots correspond to events detectable in GW and the red crosses correspond to those detectable in both GW and radio. In the (θ_v, D) diagram, the limiting distance for the measurement of the proper motion (Eq. (18)) is shown. This figure illustrates the various observational biases discussed in this section, and especially the bias towards small viewing angles, mainly due to the limiting sensitivity of radio telescopes.

Fig. 7. Predicted samples of GW-detected (dots) and jointly detected (crosses) mergers in planes of various pairs of observables, for O3-VLA configuration. Upper left: distance and viewing angle, and maximum distance at which the proper motion can be measured (Dmax(θv), see Eq. (18), green dotted line). Upper right: peak flux and viewing angle. Lower left: time of radio afterglow peak (assuming lateral expansion of jet) and viewing angle. Lower right: maximum remnant proper motion and distance, and limiting proper motion μlim (green dotted line).

4.1.5. Kinetic energy, external density and microphysics parameters

The distributions of kinetic energy E_iso, c and external density n, which are not directly observable but may result from fits of the afterglow, are shown in Fig. 8. In each case, the distributions of jointly detected events have been compared to total GW detections. As the GW selection is independent of the kinetic energy and density, these distributions are the intrinsic distributions of our population model (Table 1). As could be expected, the detection of mergers leading to a relativistic core jet with a large E_iso, c is favored. In jointly detectable mergers, the energy distribution approximately remains a broken power-law, but the respective indices below and above the break decrease from 0.53 to 0.2 and 3.4 to 3.2, respectively. After applying the joint detectability criterion, the external density distribution is simply shifted to higher densities, as the mean value is increased by a factor of ∼2.2. The distribution of the microphysics parameter ϵ_B is very similar to that of the density, because both parameters appear with the same power in the peak flux (Eq. (6)).

Fig. 8. Differential distribution of core jet kinetic energy (left) and external medium density (right) of entire GW population (full) and jointly detected events (dashed), in O3-VLA configuration. Distributions are normalized to fraction of jointly detected events among all GW-detected events, see text for details.

We find that the Lorentz factor of the jet at the time of peak flux has a median value of ∼4. This is consistent with the approximate jet core de-beaming relation Γ ∼ (θ_v−θ_j)⁻¹ for a jet with θ_j = 0.1 rad seen with the mean viewing angle θ_v ∼ 24° for the O3-VLA configuration. This justifies keeping microphysics parameters representative of shock-acceleration in the ultra-relativistic regime for the whole population.

4.2. Impact of the population model uncertainties on the predicted population

4.2.1. Jet energy distribution and opening angle

Here we consider alternatives to our fiducial population model. With a distribution of jet energy that increases below the break (i.e., taking α₁ <  0 in Eq. (11)) as would be obtained from the short GRB luminosity function of Wanderman & Piran (2015) (α₁ = −0.9) and still assuming the same standard values for the duration and efficiency, the fraction of events detectable with the VLA would typically be three times smaller: 15% for O2 and 10% for O3 compared to 43% and 31% with α₁ = 0.5. The viewing angle and peak times would also be strongly affected, showing mean θ_v’s for O2 (resp. O3) at only 17° (resp. 14°), 98% (resp. 99%) of afterglows peaking before 150 days in the jet expansion hypothesis, and 90% (resp. 93%) without lateral expansion.

These results were obtained assuming a typical value f_γ = 20% for gamma-ray efficiency but one cannot exclude that f_γ may vary from from less than 1% to more than 20%. Estimates of f_γ from afterglow fitting have indeed been found to cover a large interval (e.g., Lloyd-Ronning & Zhang 2004; Zhang et al. 2007). Thus, the three orders of magnitude in isotropic gamma-ray luminosity could partially result from differences in efficiency, the jet isotropic kinetic energy being restricted to a smaller interval. Figure 9 shows the effect of increasing the minimum kinetic energy E_min from 10⁵⁰ to 10⁵² erg in the energy function deduced from the short GRB luminosity function of Wanderman & Piran (2015) on the detected fraction, the mean viewing angle and the fraction of afterglows to peak before 150 days. As could be expected, when the minimum kinetic energy was increased, the results came closer to those of the fiducial model in which the energy distribution function decreases below the peak.

Fig. 9. Left: fraction of joint events among GW events assuming short GRB luminosity function from Wanderman & Piran (2015) for different values of minimum kinetic energy Em (full lines), and for our fiducial model (dashed lines), in O3-VLA configuration. Results have been plotted for different values of the jet opening angle θj = 0.05, 0.1 and 0.15 rad (red, black and blue lines respectively). Middle: mean viewing angle (same color coding). Right: fraction of joint events with peaks earlier than 150 days post-merger, with hypothesis of jet lateral expansion, same color coding.

We also show in Fig. 9 the effect of changing the typical value of the jet opening angle from 0.1 to 0.05 or 0.15 rad. When θ_j = 0.05 (resp. 0.15), the detected fraction decreases (resp. increases) to 19% (resp. 40%), the mean viewing angle decreases (resp. increases) by 4°, and the peak of the light curve generally occurs later (resp. earlier).

4.2.2. External density distribution and microphysics

The afterglow peak flux depends greatly on the density of the circum-merger medium, as can be seen from Eq. (6). Therefore, the ability to detect the afterglows of a population of mergers depends strongly on the density of the media hosting the binaries upon merger. This is illustrated in Table 4, where we give this fraction for different external density distribution central values. These figures were calculated using the semi-analytical model (Eq. (4)), taking into account the full synchrotron spectrum including self-absorption.

If there is a significant population of fast-merging neutron star binaries, as proposed by Matteucci et al. (2014), Hotokezaka et al. (2015), Vangioni et al. (2016), Beniamini & Piran (2016), for example, these should merge in higher density environments, producing brighter afterglows and offering the possibility to notably contribute to the observed population even if they are intrinsically fewer in number. With a density distribution centered at n = 10 cm⁻³, and the other parameters taking their fiducial values, 99% of GW events produce detectable radio afterglows in the O3-VLA configuration.

The microphysics parameter ϵ_B enters in Eq. (6) with the same power as the external density. Moreover, our fiducial choices for the distribution of these two quantities are nearly the same (with simply a limit in range for ϵ_B) so that changing its central value affects the detected fraction in the same way.

5. Discussion and conclusion

We studied the population of binary neutron star mergers to be observed jointly through GW and their radio afterglows in future multi-messenger observing campaigns. To do so, we have assumed a likely population of mergers inspired by prior short GRB science, and simulated the GW and jet-dominated radio afterglow emissions from these. We have made predictions on the rates of such future events, and on how the observables from these events are expected to be distributed among themselves.

In the case of the ongoing O3 run of the LVC, and assuming our fiducial set of parameters (see Table 1), we predict that ∼30% of GW events should have a radio afterglow detectable by the VLA. These joint events should have mean viewing angles of ∼24°, and should peak earlier than 150 days post-merger in ∼80% of cases, assuming lateral expansion of the jet.

As illustrated in the 170817 event, the lateral material may also contribute to the peak flux of the radio afterglow. Figure 1 (right) shows that this may increase the peak flux by a factor of up to 1.5. In this case, our O3 predictions are revised to a fraction of 36% of events with detectable afterglow and a similar mean viewing angle of 25°.

5.1. Detecting the detectable events

Throughout this work, we have applied a threshold on the GW S/N and radio afterglow peak flux to label an event as jointly detected. Thus our study concerns the class of detectable events. The fraction of these detectable events that can actually be detected heavily depends on the accuracy of the localization, which is related to the size of the localization skymap provided by the GW network, the possible detection of an associated GRB, and the efficiency of the search for a kilonova in terms of covered sky area and limiting magnitude.

Assuming that all neutron star mergers produce a kilonova, and even taking into account a probable anisotropy of the emission, it is expected that 100% of kilonovae are detectable during O3 for a limiting magnitude of 21 (Duque et al. 2019; Mochkovitch et al. in prep.). However, the search for the kilonova becomes complicated due to the large localization skymaps of GW events. These error boxes can be improved in the case when an associated GRB is being detected, which should, however, tend to remain rare (see Sect. 4.1.2 and Beniamini et al. 2019). During O3, existing facilities with large fields of view in the optical band are experiencing difficulty in covering 100% of ≳1000 deg² skymaps with a deep limiting magnitude, especially in the southern sky. Instruments such as the Zwicky Transient Facility (Bellm et al. 2019) cover the northern sky down to magnitude r ∼ 20.5.

On the other hand, when the LIGO/Virgo instruments will have achieved their design sensitivities, it is expected that one or two new interferometers will join the GW network: KAGRA (Kagra Collaboration 2019) and Ligo-India (Iyer 2015). This should lead to a significant improvement in localization. In addition, the deep coverage of even large error boxes will be made possible by new facilities such as the Large Synoptic Survey Telescope (Ivezic et al. 2008). Nevertheless, there are other difficulties related to the search for the kilonova: kilonova/host galaxy contrast at large distances, possible large offsets, availability of a photometric and spectroscopic follow-up of the candidates, etc. In this context, it is clear that the efficiency of the search for kilonovae will always remain below 100%, but should increase with time and may potentially be high in the design configuration of the GW network.

Even assuming the localization has been acquired thanks to the detection of a kilonova, a continuous monitoring of the remnant up to ∼150 days may be necessary to detect the peak of the radio afterglow, in the case of marginally detectable events. Even then, only events with peaks somewhat larger than the radio threshold should yield observations of astrophysical interest because extended observations of the rise, peak and decay of the afterglow are necessary to resolve the structure and dynamics (e.g., expansion) of the ultra-relativistic jet, as illustrated by the case of GRB170817A.

As also illustrated in the case of GRB170817A, VLBI measurements are instrumental in assessing the presence and resolution of the jet structure. Therefore, an even more restrictive criterion for full event characterization might be the detectability of the remnant proper motion. As we demonstrated in Sect. 4.1.4, this would primarily decrease the fraction of events, especially in the context of design GW detectors.

5.2. Constraining the BNS population and short GRBs with radio afterglows

We have studied the sensitivity of the expected distributions of afterglow observables to the population parameters, as well as the radio and GW detector configurations. In particular, we have shown the uncertainties on the rates and typical viewing angles of these events which have been inherited from the uncertainties on the luminosity function of short GRBs and of the density of the media hosting the mergers.

In Sects. 4.1.1 and 4.1.2 we discussed the combined influence of increasing the radio and GW detector sensitivity on the rate and viewing angle of the events. At this stage, both improvements would be beneficial as the predicted rate in the current O3-VLA configuration is well below 50% of joint events among GW events.

As the GW horizon increases, more events will be seen on-axis. In the case of an opening angle of 0.1 rad and in the likely short-term evolution of the detector configuration from O3-VLA (current) to design-VLA (early 2020s), we predict that 5% to 10% of events should be seen on-axis, increasing the probability of a GRB counterpart in addition to the GW and afterglow. This figure is fully consistent with that announced in Beniamini et al. (2019) for events with a GW+GRB association. In fact, as a GW+classical GRB association is equivalent to an on-axis event, one may measure the opening angle of typical GRB jets by considering the ratio of GW events with afterglow counterpart only to events with both afterglow and GRB, thus exploiting the new possibility of observing afterglows without detecting the GRB.

More generally, a population of observed events corresponds to an intrinsic population of mergers. Thus, a comparison of a statistical number of joint event observations to our predicted distributions is a means of measuring fundamental parameters of the population of mergers, and of constraining GRB quantities, such as their energy function.

Finally, future facilities such as the SKA may bring detections of orphan radio afterglows through deep radio surveys. These would not be subject to the GW criterion and may probe another sub-population of mergers, bringing yet another class of constraints upon these phenomena.

In conclusion, regardless of the evolution of GW and radio detectors, multi-wavelength afterglows will remain instrumental in the study of binary neutron star mergers as windows on both their environment and their role as progenitors of short GRBs. They will bring invaluable insight at the level of the population of mergers and on an event-to-event basis, provided the sources may be accurately localized in the sky.

Acknowledgments

The authors acknowledge financial support from the Centre National d’Études Spatiales (CNES), and thank Paz Beniamini for useful discussions.

References

All Tables

All Figures

Fig. 1. Left: 3 GHz light curve of afterglow from structured jet at distance D = 42 Mpc, with sharp power-law structure (a = 4.5, b = 2.5; see Eqs. (1)–(2)), viewing angle θv = 22°, external density n = 3 × 10−3 cm−3, core jet with θj = 4°, Γc = 100, Eiso, c = 2 × 1052 erg, and with microphysics parameters p = 2.2, ϵe = 0.1 and ϵB = 10−4. The radio observations of GW170817 are also plotted and are compiled from Hallinan et al. (2017), Alexander et al. (2018), Margutti et al. (2018), Mooley et al. (2018a,b), Dobie et al. (2018). The respective contributions of the core jet (θ ≤ θj), the sheath and the total are plotted in dashed, dotted and solid red lines. The flux of the core jet computed with the simplified treatment described by Eq. (4) is also plotted in dashed blue line for comparison. Bottom right: peak flux of same structured jet as function of viewing angle θv (solid red line), peak flux of light curve from core jet only (dashed red line) or from sheath only (dotted red line), and peak flux of core jet as computed using scaling law in Eq. (6) (thin black line). Upper right: ratio of peak flux from whole outflow (core jet+sheath) to that of core jet only.

In the text

	Fig. 2. Left: detectable fraction of radio afterglows among gravitational wave events as function of the horizon distance $H=1.58\overline{H}$. Right: expected number of joint detections normalized to case of O3-VLA configuration. In both panels the full (resp. dashed, resp. dotted) lines correspond to the VLA (resp. SKA1, resp. SKA2/ngVLA) being the limiting radio facility.
In the text

Fig. 3. Left: differential distribution of distances (normalized to the sky-position-averaged horizon) of events in total GW population (full line), and detectable jointly in O2-VLA (dotted line) and O3-VLA (dashed line) configurations. Right: cumulative distribution of same events. Distributions have been normalized to the fraction of jointly detected events among all GW-detected events, see text for details.

In the text

	Fig. 4. Same as Fig. 3, for viewing angle. Distributions are normalized.
In the text

	Fig. 5. Left: mean viewing angle of jointly detected events, as function of horizon and radio sensitivity s. Values of H for O2, O3 and design instruments, as well as s for VLA, SKA1 and SKA2/ngVLA, are indicated as dashed lines. Right: fraction of on-axis events (defined by θv ≤ θj) among jointly detected events in same H–s diagram.
In the text

Fig. 6. Left: cumulative distributions of peak times of radio afterglows for jointly detected events in O3-VLA configuration, assuming jet with (solid line) or without (dashed line) lateral expansion. Middle: differential distribution of peak fluxes of GW-detected events assuming horizon for O2 (dotted line), O3 (dashed line) and design (full line) configurations of GW interferometers. The radio sensitivities for VLA, SKA1 and SKA2/ngVLA are indicated for comparison. Right: differential distribution of maximum proper motion for jointly detected events assuming O2-VLA (dotted line), O3-VLA (dashed line), design-SKA1 (full line) configurations. The dotted vertical line shows the lower limit for detection of the proper motion.

In the text

Fig. 7. Predicted samples of GW-detected (dots) and jointly detected (crosses) mergers in planes of various pairs of observables, for O3-VLA configuration. Upper left: distance and viewing angle, and maximum distance at which the proper motion can be measured (Dmax(θv), see Eq. (18), green dotted line). Upper right: peak flux and viewing angle. Lower left: time of radio afterglow peak (assuming lateral expansion of jet) and viewing angle. Lower right: maximum remnant proper motion and distance, and limiting proper motion μlim (green dotted line).

In the text

	Fig. 8. Differential distribution of core jet kinetic energy (left) and external medium density (right) of entire GW population (full) and jointly detected events (dashed), in O3-VLA configuration. Distributions are normalized to fraction of jointly detected events among all GW-detected events, see text for details.
In the text

Fig. 9. Left: fraction of joint events among GW events assuming short GRB luminosity function from Wanderman & Piran (2015) for different values of minimum kinetic energy Em (full lines), and for our fiducial model (dashed lines), in O3-VLA configuration. Results have been plotted for different values of the jet opening angle θj = 0.05, 0.1 and 0.15 rad (red, black and blue lines respectively). Middle: mean viewing angle (same color coding). Right: fraction of joint events with peaks earlier than 150 days post-merger, with hypothesis of jet lateral expansion, same color coding.

In the text