Prospects for kilonova signals in the gravitational-wave era

The binary neutron star merger gravitational-wave signal GW170817 was followed by three electromagnetic counterparts, including a kilonova arising from the radioactivity of freshly synthesized $r$-process elements in ejecta from the merger. Finding kilonovae after gravitational-wave triggers is crucial for (i) the search for further counterparts, such as the afterglow, (ii) probing the diversity of kilonovae and their dependence on the system's inclination angle, and (iii) building a sample for multi-messenger cosmology. During the third observing run of the gravitational-wave interferometer network, no kilonova counterpart was found. We aim to predict the expected population of detectable kilonova signals for the upcoming O4 and O5 observing runs of the LIGO-Virgo-KAGRA instruments. Using a simplified criterion for gravitational-wave detection and a simple GW170817-calibrated model for the kilonova peak magnitude, we determine the rate of kilonovae in reach of follow-up campaigns and their distributions in magnitude for various bands. We briefly consider the case of GW190425, the only binary neutron star merger confirmed since GW170817, and obtain constraints on its inclination angle from the non-detection of its kilonova, assuming the source was below the follow-up thresholds. We also show that non-gravitational-wave-triggered kilonovae can be a numerous class of sources in future surveys and briefly discuss associations with short bright gamma-ray bursts. We finally discuss the detection of the jetted outflow afterglow in addition to the kilonova.


Introduction
The first detection of electromagnetic counterparts to a gravitational-wave (GW) event was a truly historic event (Abbott et al. 2017b). The coalescence of two neutron stars detected by the LIGO-Virgo instruments on August 17, 2017 (GW170817; Abbott et al. 2017d) was followed 1.7 s later by a weak short gamma-ray burst (GRB) observed by Fermi and Integral (Goldstein et al. 2017;Savchenko et al. 2017;Abbott et al. 2017c). A search of the error box with optical telescopes led, after 11 hours, to the discovery of a kilonova in the spheroidal galaxy NGC 4993 at ∼40 Mpc (Abbott et al. 2017b, and references therein). Then, after respectively 9 and 16 days the afterglow was detected in X-rays with Chandra (Troja et al. 2017) and in radio with the Very Large Array (VLA; Hallinan et al. 2017). The afterglow light curve was atypical, with a steady rise to a maximum at about 170 days post-merger D'Avanzo et al. 2018;Nynka et al. 2018;Resmi et al. 2018a;Mooley et al. 2018b). While such a behavior could result from either a radial or an angular structure of the ejecta (e.g., Gill & Granot 2018), very long baseline interferometry observations showing a displacement of the unresolved source by about 2.5 mas in 5 months (Mooley et al. 2018a;Ghirlanda et al. 2019) provided firm evidence for the latter. Joint fits to the afterglow photometry and The kilonova dubbed AT 2017gfo reached a peak magnitude of ∼ 17 in the g, r, i, and z bands and was followed for two to three weeks in the infrared (z to K bands), where the decline is shallower than in the visible. The data are well fit by the combination of a blue, mostly polar component, which declines faster, and a more isotropic red component. The red component results from the high opacity lanthanide-rich material tidally ejected during coalescence or blown off from the accretion disk (Cowperthwaite et al. 2017;Tanvir et al. 2017;Villar et al. 2017a;Tanaka et al. 2017). Different possible origins have been considered for the blue, lanthanide-poor component, which is expected to be present only if the central core does not directly collapse into a black hole (Metzger 2019) 2 . In most cases, the core is probably short lived and collapses into a black hole. In 1 Throughout this paper we systematically state median 90% credible intervals. When primary sources state results with other credible intervals, we scaled such intervals assuming a Gaussian distribution. 2 In AT 2017gfo, a third "purple component" with an intermediate lanthanide fraction can be added to improve the fit (Cowperthwaite et al. 2017;Villar et al. 2017a); however, a thorough preference for a three-component model has not been established so far.
Article number, page 1 of 10 arXiv:2103.00943v2 [astro-ph.HE] 19 May 2021 the case of GW170817, the eventual GW signal of the black hole ring-down was too weak to confirm if or when this collapse had occurred (Abbott et al. 2017d). This single event represented a breakthrough in the understanding of merger physics and the study of the origin of heavy elements in the Universe (for a review, see Ciolfi 2020). Also, it allowed the first standard siren measurement of the Hubble constant (Abbott et al. 2017a;Hotokezaka et al. 2019).
During the third LIGO-Virgo-KAGRA Collaboration (LVKC) observing run, O3, the only confirmed binary neutron star merger so far (i.e., GW190425) was located at 159 +69 −71 Mpc (The LIGO Scientific Collaboration & The Virgo Collaboration 2019; Abbott et al. 2020a). With only one of the LIGO interferometers taking data at the time of the event, the GW sky map was very large, nearly 7500 deg 2 . It could only be partially explored by the various follow-up efforts. The Pan-STARRS and Zwicky Transient Facility (ZTF) telescopes achieved the largest coverage of the event , and no kilonova was found. It remains unknown if the kilonova was weaker than the detection limits or simply outside the area covered by the searches.
Following the premature end of O3 at the end of March 2020, the GW detectors are expected to resume operations in mid-2022 with a binary neutron star merger range increased by about 50% (Abbott et al. 2020b). The GW discovery rate of binary neutron star mergers should then reach 10 +52 −10 per calendar year, a factor of ten larger than during O3. The participation of KAGRA in this run will reduce the average sky surface and volume of searches for kilonova counterparts. However, these kilonovae will be located at larger distances on average, possibly impeding their detection.
The goal of this paper is to obtain the expected distributions of various physical parameters for kilonovae that should be detected in association with GW signals of merging binary neutron stars during O4 and beyond. Using a simple parametrization of the kilonova magnitude as a function of viewing angle in different spectral bands, we obtain the distributions (i) in magnitude for visible and nearinfrared bands and (ii) in viewing angle for different limiting magnitudes of the kilonova follow-up search. We also estimate the corresponding discovery rate, again for different limiting magnitudes.
In the case of GW190425, we show how a constraint on the viewing angle can be obtained from the lack of counterparts in three bands. We derive this constraint by assuming that the source was indeed located in the areas searched during follow-up efforts, but below detection threshold.
When a kilonova is found, the sky location is known with an arcsecond accuracy, which allows searching for the afterglow in X-rays or radio. We calculate the fraction of sources that can be detected in radio in addition to the kilonova with the VLA as well as their distribution in viewing angle. In the visible, the afterglow is likely to be initially outshone by the kilonova. This is expected as long as the viewing angle is larger than the opening angle of the jet core, which is most probable when the alert is given by the GW detectors.
This paper is organized as follows: In Sect. 2 we first obtain the distributions in distance and viewing angle of the neutron star merging events detected in GWs. Our simplified parametrization of the kilonova magnitude as a function of viewing angle is presented in Sect. 3. The resulting distributions in magnitude and viewing angle of the detectable kilonovae are shown in Sect. 4 together with the constraints that can be obtained on the viewing angle for GW190425. The possibility to observe the afterglow when the kilonova has been found is discussed in Sect. 5. Finally, our results are discussed in Sect. 6, and Sect. 7 is the conclusion.

Distribution in distance
Since we are interested in kilonovae in association with a GW signal, their distance and viewing angle distributions simply follow those of GW-triggered neutron star merger events. For a given sky-position-averaged horizon D H , corresponding to a binary neutron star having its rotation axis pointing toward the observer, detection at viewing angle θ v is possible only at distances D such that (Schutz 2011) Then, for D ≤ D 0 = D H / √ 8, all sources are detected, while for D 0 < D ≤ D H , they are progressively lost until, for D H , only those pointing directly at the observer remain. The resulting distribution in distance is represented in Fig. 1a, which shows a maximum at D/D H = 0.63. Figure 1b gives the corresponding cumulative distribution.
In Table 1 we state our assumed sky-position-averaged horizons for past and upcoming GW observing runs. These were taken from Abbott et al. (2020b) Table 1: Horizon distances assumed for the various GW observing runs, as used in the detection criterion in Eq. 1, and parameters for the kilonova peak absolute AB magnitude dependence on the viewing angle, as given in Eq. 4.

Distribution in viewing angle
The distribution in viewing angle of the gravitationally detected sources is represented in Figs. 1c and 1d. It peaks at θ v ∼ 30 deg, with an average value θ v = 38 deg. The fraction of sources in the conservative interval 10 deg < θ v < 20 deg corresponding to GW170817 is 14% (increased to 27% for 5 deg < θ v < 25 deg). Had GW170817 occurred during the O3 run, with D H = 157 Mpc, the distance to the source would have verified that D 170817 ∼ 40 Mpc < D 0 . Therefore, any merger at this distance would have been detectable, regardless of the inclination angle. In this case, the expected rate of binary We also indicate the median values and D 0 , the distance under which no selection in viewing angle occurs. The cusp at D 0 in the differential distribution in distance is nonphysical and a consequence of the simplified nature of the adopted GW detection criterion.
neutron star mergers up to the distance of GW170817 is simply given by where τ BNS = 320 +490 −240 Gpc −3 yr −1 is the local binary neutron star merger rate (Abbott et al. 2020d). This leads to an average rate of one event every R −1 = 12 +36 −7 yr. GW170817 not only was a nearby event but had a low inclination angle, θ 170817 v < θ 170817,max v ∼ 18 deg, according to the very long baseline interferometry observations (Mooley et al. 2018a;Ghirlanda et al. 2019). The detection of the radio afterglow and source proper motion was possible only up to a viewing angle of θ AG,max v ∼ 40 deg (e.g., Duque et al. 2019). Requiring θ v ≤ θ AG,max v -to get a rich multi-messenger data set with an inspiral signal as well as kilonova and afterglow photometry and imagery data -therefore leads to a rate of approximately that is, an average rate of one event every R −1 = 50 +149 −31 yr. The detection of the short GRB may require a even smaller viewing angle, , as GRB170817A was detected only at the ∼ 5σ level by the Gamma-Ray Burst Monitor aboard Fermi (Goldstein et al. 2017). Requiring θ v ≤ θ GRB,max v to get a full GW170817-like multi-messenger data set, including the short GRB, leads to an even lower rate, Abbott et al. 2020b), when D 170817 was not smaller than D 0 ∼ 30 Mpc, these rates were even lower. These numbers illustrate how lucky we were to detect GW170817 so early and how long we may have to wait to observe another equivalent event.

Kilonova magnitude dependence to viewing angle
The kilonova magnitude at the peak depends on the distributions of the mass, velocity, and composition of the ejected material as well as on the viewing conditions: distance and viewing angle. The ejection is anisotropic with neutron-rich, dynamical ejecta in the equatorial plane, where the formation of lanthanides leads to a large opacity while a relatively neutron-poor wind of lower opacity is blown in the polar direction (Fernández & Metzger 2016;Metzger 2019;Barnes 2020). This wind is expected to be present when a short-lived massive neutron star is formed before collapsing to a black hole, but probably not in the case of a direct collapse. The lanthanide-rich ejecta produces the "red kilonova," which peaks in the near-infrared, while the neutronpoor wind is responsible for the "blue kilonova" at optical wavelengths. The blue kilonova declines on a timescale of one day, whereas the timescale is one week for the red component. For our population model, our default scenario assumes that all kilonovae have a quasi-isotropic red component and a polar blue component. We obtain the peak absolute AB magnitude at a given wavelength and viewing angle from the following simple parametrization: where M λ,0 is the peak absolute magnitude for a polar viewer and ∆M λ is the amplitude of the polar effect.
The δM λ represents the intrinsic (i.e., non-viewing-anglerelated) variability in kilonova magnitudes linked to the abovementioned ejecta properties and, in turn, to the progenitor component masses and spins. For θ 0 = 60 deg, we find that the linear-in-cos θ v form of Eq. 4 reproduces the trends of sophisticated kilonova modeling work (e.g., Wollaeger et al. 2018, Kawaguchi et al. 2020, and the "asymmetric model" of Villar et al. 2017b). We chose δM λ to be uniformly distributed in [−1, 1], reproducing the expected variability in kilonova magnitude stemming from variability in ejecta mass, velocity, and opacity (Wollaeger et al. 2018, Eq. 33). The difference in magnitude between equatorial and polar views is moderate in the infrared but increases rapidly in the visible, already reaching about 4 mag in the r band. This is mainly due to the stronger anisotropy of the blue component. To calibrate this expression, we used AT 2017gfo, assuming θ v = 15 deg, as mentioned in the introduction. Corresponding values can be found in Table 1.
Calibrating Eq. 4 with AT 2017gfo supposes that this transient was representative of the kilonova population. This is the minimal hypothesis one can make while waiting for the number of kilonovae with robust angle measurements to increase in the future. We note that AT 2017gfo could have been brighter or dimmer than the average of the population our model seeks to encapsulate. We briefly indicate below how our results might change if this is indeed the case.
The polar ejecta may not be produced in all mergers, depending for instance on the post-merger formation of a massive neutron star before its collapse into a black hole (see, e.g., Metzger 2019). As such, we also consider below the possibility that a fraction of the kilonova population lacks the blue component, which would affect kilonova brightnesses in the bluer bands (see a preliminary luminosity function in Ascenzi et al. 2019 and related discussions in Gompertz et al. 2018 andKasliwal et al. 2020).

Apparent magnitude
From the known distance and viewing angle distributions and our adopted parametrization for the magnitude (Eq. 4), we can readily obtain the distribution of apparent AB magnitudes for kilonovae associated with GW detections. It is shown in Fig. 2a for the g, r, i, and z bands for the O4 observing run. If AT 2017gfo was in fact brighter than the average population, all the curves will have to be shifted by the corresponding difference, δmag = m − m 170817 . Changing the GW horizon implies an interplay between the maximum detection distance for the kilonova and the GW and thus does not result in a simple shifting of the magnitude distribution. However, we found that, to a good approximation, changing from O4 to O5, the magnitude distribution is shifted by about 5 Log(D H,O5 /D H,O4 ) = 1.6 mag.
The distribution of kilonovae in different magnitude ranges is summarized in Table 2 for three GW sensitivity hypotheses: O3, O4, and O5. It can be seen that there are very few kilonovae with m < 18 in all cases beyond O3. We note that recalibrating Eq. 4 assuming that AT 2017gfo was one magnitude brighter than average leads to dividing the expected fractions in the < 20 magnitude ranges for O4 by approximately three.
The kilonova magnitudes depend on the different merger ejecta and their physical conditions. The blue kilonova component is likely linked to neutrino-or magnetohydrodynamic-viscosity-driven winds from the transient remnant product and the accretion disk around such a product (Gill et al. 2019, and references therein).
It is possible that in some systems blue-enhancing ejection episodes are less effective (e.g., because of a short-lived merger remnant), leading to a lack of a blue kilonova component. We briefly consider this possibility, without seeking to know the fraction of these cases in the population.
If a fraction f red of the kilonovae lack the blue component, a simple approximation consists in stating that such kilonovae will be dimmer and thus transferred from the two brightest magnitude groups to the m > 20 group. This leads to the following for all bands: where the f 0 fractions are those listed in Table 2. We tested these approximations with our kilonova population model by emulating the absence of the blue component in a fraction f red of the synthetic kilonovae. We did this by adopting the θ v > θ 0 case of Eq. 4 for all viewing angles, as if only the red component were present. As the blue component affects the g, r, and i bands more, these expressions represent reasonable approximations of the exact results, while in the z band they somewhat overestimate the number of events that change from the < 20 to the > 20 magnitude groups.
The expected rates of kilonovae brighter than a given limiting r magnitude are shown in Fig. 2b for O3, O4, and O5, normalized to a GW neutron star coalescence detection rate of τ BNS,GW = 10 yr −1 for O4 (Abbott et al. 2020b). At the bright end of the distribution (i.e., r < 19), a fit to Fig. 2b shows that the rate approximately follows: where r lim is the limiting magnitude in the r band.
As an illustration, with f red = 0 and f red = 0.2 we expect one kilonova brighter than r = 19 every 1.6 years and every 2.0 years, respectively, and one brighter than r = 18 every 6.3 years and 7.9 years, independently of any future improvement in the sensitivity of GW detectors. On the other hand, the rate of kilonovae detectable by a followup with a limit magnitude r lim = 21 is increased by a factor of ∼ 3 between O3 and O4.
In Fig. 2b we also show the maximum r-band magnitude of any kilonova associated with a GW trigger for O3, O4, and O5, denoted by r max . These magnitudes are the search depths required to recover 100% of the kilonovae. Because our peak magnitude dependence with viewing angle saturates at θ 0 = 60 deg, these maximum-magnitude events have θ v = θ 0 and are placed at the largest distance to which the GW signal can be detected at this angle (i.e., at D/D H = (1 + 6 cos 2 θ 0 + cos 4 θ 0 )/8 ∼ 0.55).

Distribution in viewing angle for different limiting magnitudes
The distribution in viewing angle of the kilonovae associated with binary neutron star merger triggers that are brighter than a given limiting magnitude is shown in Fig. 3 for O4. As the limiting magnitude decreases, the median kilonova viewing angle -close to 36 deg in the entire population of GW triggers -significantly decreases: 26 deg for an r-band limiting magnitude of 21 and 21 deg for all r lim smaller than 20.
In Fig. 4 we study the distributions in distance and viewing angles for events detected in the gravitational or optical domains. For this figure only, we remove the intrinsic variability of kilonovae introduced in Eq. 4 (i.e., we set δM λ to 0) to clearly separate the different observing scenarios in the distance-viewing angle plane.
For limiting r-band magnitudes equal to or smaller than 20, practically all kilonovae that can be detected will follow a GW event if the interferometers are taking data at the corresponding time. Conversely, for deeper searches reaching r lim = 21 or 22, the fraction of "orphan kilonovae" without a GW alert increases and becomes dominant.
Recently, an archival study searching for kilonovae in 23 months of ZTF data was carried out. Down to a limiting magnitude of r lim ∼ 20.5 for source detection, no transient consistent with being a kilonova was identified (Andreoni et al. 2020). Considering that a kilonova can be safely detected and characterized only if its peak magnitude is at least one magnitude brighter than the limit of the survey, Fig. 4 allows us to estimate the number of expected kilonova detections over the 23 month period. Assuming perfect identification and a sky coverage of ∼ 50%, as appropriate for ZTF, we find between 0.4 and 2.6 detections, taking into account the uncertainty on the binary neutron star merger rate but not that linked to the kilonova model. Beyond the kilonova model uncertainties, an overestimated rate of mergers or the limitations of the kilonova identification algorithm as discussed in Andreoni et al. 2020 could also contribute to the non-detection. Future surveys and archival studies by other optical facilities (Almualla et al. 2020;Setzer et al. 2019) should clarify which of these options is the most likely.
GRB170817A, associated with GW170817, was very weak considering the distance (Goldstein et al. 2017) and cannot be considered as the off-axis view of a bright short GRB (Matsumoto et al. 2019). It was not produced by the core ultra-relativistic jet revealed by very long baseline interferometry observations (Mooley et al. 2018a;Ghirlanda et al. 2019); it was instead emitted from the material located at a higher latitude that is propagating toward us (Matsumoto et al. 2019), though not necessarily via the same mechanism as that of other short GRB that are observed on-axis at large distances. Therefore, the question of the short GRB-merger connection remains open, even if the evidence for the production of an ultra-relativistic jet in GRB170817 is clearly in line with this hypothesis.
The detection of a bright short GRB seen on-axis following GWs from a binary neutron star merger would represent direct evidence for this connection. Figure 4 allows us to discuss the probability of such events in the future. With a minimum peak luminosity of ∼ 10 50 erg/s (Wanderman & Piran 2015;Ghirlanda et al. 2016) and a peak energy on the order of 1 MeV (Nakar 2007), short GRBs seen on-axis are bright at any distance below 600 Mpc (peak flux on the order of 1 ph/cm 2 /s), and the main limitation for detection by gamma-ray satellites is their sky coverage. Assuming a typical jet opening angle θ j = 0.1 rad (see e.g., Fong et al. 2015;Beniamini et al. 2019), Fig. 4 clearly indicates that triple associations of GWs with kilonovae and bright short GRBs seen on-axis should remain especially rare: one event every 5-20 years in the whole sky, according to our calculations.  Fig. 2: Distribution of peak kilonova magnitudes and kilonova detection rate as a function of threshold. Left: Distribution of the peak AB magnitude in the g, r, i, and z bands predicted for kilonovae associated with GW triggers during O4. Right: Rate of kilonovae brighter than a given r magnitude associated with GW detections during O3, O4, and O5, assuming a GW neutron star coalescence detection rate of τ BNS,GW = 10 yr −1 for O4 (Abbott et al. 2020b The association of a bright short GRB with a kilonova even without a GW detection is also a solid argument in favor of the merger connection. Figure 4 shows that the rate of such a double association is more optimistic if the limiting magnitude in the r band is at least 21, with ∼ 2 such events per year according to our calculations. However, for such bright GRB associations, it has been noted that the optical kilonova signal should only outshine the afterglow flux in dense circum-merger media and with less energetic jets, allowing for an early-breaking or dimmer afterglow (Guessoum et al. 2018).
GRB130603B and GRB050709 were well within the parameter region allowing for the kilonova to appear (Fong et al. 2015), and yet the associated claimed kilonova components (Tanvir et al. 2013;Jin et al. 2016) were only marginally brighter than the afterglow and required a follow-up duration of longer than a week to be detected. Still, the potential of such sources for studying binary neutron star merger physics and the larger distances to which these can be detected encourages deep photometric followups of short GRBs.

Viewing angle-magnitude diagram for GW190425
GW190425, the only confirmed binary neutron star merger during LVKC observing run O3, was located at 159 +69 −71 Mpc (Abbott et al. 2020c). No kilonova was found during the follow-up that was conducted by several facilities. The deepest searches were led by ZTF to mag. 21 in the g and r bands, covering 21% of the probability enclosed in the very large final GW sky map of nearly 7500 deg 2 (Kasliwal et al. 2019;Coughlin et al. 2019), and by Pan-STARRS to mag. 21.5 in the i band, covering 28% of the initial GW sky map (Smith et al. 2019;Coughlin et al. 2020).
This non-detection can have different origins, the most obvious being simply that the kilonova was not located in the search area. But it is also interesting to explore the possibility that the kilonova was there but below the detection limit. This possibility is supported by recent modeling predicting that the kilonova of GW190425 was dimmer and faster decaying than AT 2017gfo because of its larger mass (Nicholl et al. 2021). Figure 5 illustrates the resulting constraints in viewing angle-magnitude diagrams, adopting Eq. 4 for the kilonova magnitude. For this particularly massive event, we increased the horizon accordingly with the larger chirp mass of GW190425 with respect to a 1.4+1.4 M system, following the relation D H ∝ M 5/6 (Schutz 2011; see the value in our Table 1).
The viewing angle to GW190425 is bounded below by the non-detection of the kilonova and bounded above by the detection of the GW signal. This results in the individual constraints from the g, r, and i bands that can be read off Relative occurrence rates of signals in the distance-viewing angle plane predicted for the O4 run. Colors indicate different detection scenarios: events detectable (i) both as GW triggers and kilonovae (green), (ii) as GW signals alone (blue), and (iii) only as kilonovae ("orphan kilonovae," red). For the kilonova detection, the four diagrams correspond to the limiting r-band magnitudes from 19 to 22. We also indicate the total occurrence rates in each detection scenario, assuming a GW detection rate of 10 yr −1 for O4.
thresholds for the largest viewing angle range, according to our model. The combined three-band constraint is θ 190425 v = 53.3 +9.8 −12.4 deg, to which a systematic uncertainty due to the kilonova model should be added (see Sect. 6). Finally, in the case where no blue kilonova was produced in that event -possibly because the central core of the merged object directly collapsed into a black hole -no useful constraint can likely be obtained. This last possibility is indeed worth considering because of the large masses of the two neutron stars inferred for GW190425.

Detecting the radio afterglow
When the kilonova is found, the location of the source is known with an arcsecond accuracy, which allows an efficient search for the afterglow. Assuming an index p = 2.2 for the power-law distribution of the shock-accelerated electrons, the peak flux of the radio afterglow at 3 GHz is given by Nakar et al. (2002) and Gottlieb et al. (2019): where ϕ = E 52 θ 2 j,−1 n 0.8 −3 1.2 e,−1 0.8 B,−3 collects the flux dependences on the parameters not related to the observing conditions. Here E 52 and θ j,−1 are the isotropic energy and opening angle of the jet core (in units of 10 52 erg and 0.1 rad), n −3 is the density of the external medium (in units of 10 −3 cm −3 ), e,−1 and B,−3 are the usual microphysics parameters (in units of 0.1 and 10 −3 ), and D 100 is the distance of the source (in units of 100 Mpc). The normalization of ϕ was chosen such that ϕ = 1 for GW170817-like afterglows (e.g., Lamb & Kobayashi 2017;Resmi et al. 2018b;Lazzati et al. 2018;Troja et al. 2019).
We find that for O4 and with ϕ = 1, the afterglow can be detected at the VLA 3 GHz sensitivity of 15 µJy for 37%, 56%, and 76% of kilonovae, assuming r-band limiting magnitudes of 21, 20, and 19, respectively, for the kilonova search. To construct an exploitable light curve, we can impose a radio flux threshold of three times the detection limit. In this case, these fractions become 23%, 36%, and 53%. In terms of absolute numbers, they respectively correspond to 1.1, 0.7, and 0.3 joint GW-kilonova-afterglow detections per year. For particularly energetic jets or dense circum-merger environments (i.e., with ϕ = 10), the fractions of kilonovae with radio afterglows at three times the VLA limit rise to 59%, 81%, and 97%. This corresponds to 2.9, 1.6, and 0.5 three-signal detections per year.
We represent in Fig. 6 the distribution in viewing angle of the afterglows that can be detected with the VLA at three times the detection limit for ϕ = 1 and 10 and different kilonova search limiting magnitudes. Due to the very strong dependence of the afterglow flux on the viewing angle (F ∝ θ −4.4 v for p = 2.2; see Eq. 7), the detection is only possible at small viewing angles, < 15 − 20 deg.

Discussion
We have presented a population model for kilonova counterparts to GW binary neutron star merger signals. In Sect. 2 we obtained the distributions in distance and viewing angle of the GW-trigger events (see Fig. 1). Then, using a simple parametrization of the kilonova peak AB magnitude in various spectral bands as a function of viewing angle, we computed the distributions of sources in three magnitude intervals (<18, 18-20, and >20) for O3 and the future O4 and O5 observing runs. We also considered the possibility that some kilonovae lack the blue component due to the lack of some mass ejection episodes during the merger.
The rate of kilonovae brighter than a given limiting magnitude was obtained (see Eq. 6 and Fig. 2b). This confirms the extraordinary chance of having observed the August 17 event so early. In Fig. 3 we studied the viewing angles of detectable kilonovae for different limiting magnitudes. The median of this distribution, about 36 deg for the GW trigger population, decreases as the search becomes shallower, reaching 26 deg and then 21 deg for r-band limiting magnitudes of 21 and lower than 20.
We then studied the regions in distance-viewing angle space where the gravitational signal or the kilonova can be detected (see Fig. 4). In each zone, we estimated the event rate, normalized to a GW trigger rate of 10 yr −1 , as expected for O4. For deep surveys reaching mag. 21-22, the rate of orphan kilonovae without a detectable GW signal becomes dominant, opening the way to detecting kilonova counterparts to short GRBs, which should become much more common. The progress in understanding the merger phenomenon by leveraging such signals motivates the effort to carry out these surveys and the optical follow-up of short GRBs.
As the sensitivity of GW detectors increases, more distant events will be found, but the rate of bright kilonovae (r < 19) will not change. Conversely, follow-ups reaching magnitudes 20-21 have the potential to find five to ten times more kilonovae beyond magnitude 19, leading to a potential discovery rate of ten events per year during O4 to more than twenty for O5 (see Fig. 2b). Obviously, going from potential to effective discoveries will require a deep inspection of the GW sky map by target-of-opportunity endeavors or the analysis of untargeted searches by high cadence facilities, such as the Vera C. Rubin facility (Margutti et al. 2018;Cowperthwaite et al. 2019 on-axis GW events r < 21 r < 20 r < 19 r < 18 Fig. 6: Distributions in viewing angle of the afterglows detectable in the radio band at three times the VLA threshold (45 µJy) following a GW-triggered kilonova detection with a limiting r-band magnitude of 18 to 21. The dotted lines represent the corresponding distributions for the kilonova sources, as in Fig. 3. In the left (resp. right) panel, a value ϕ = 1 (resp. ϕ = 10, corresponding to particularly energetic jets or dense circum-merger media) was adopted in Eq. 7.
We studied the case of radio afterglow counterparts to kilonova signals. We note that our results in this respect are consistent with our prior study of afterglow counterparts (Duque et al. 2019), during which the kilonova was not considered. That study and the present paper find consistent results in the limit of very deep kilonova search limiting magnitudes. However, the present paper shows that, as of O4 and with current follow-up capacities, the kilonova detection will become a limiting factor in the search for afterglow counterparts.
We also used our kilonova model to constrain the viewing angle of GW190425, the only confirmed binary neutron star merger event since GW170817. No associated kilonova was detected; this could simply be a consequence of the poor localization, which limited the search to less than 30% of the GW sky map. If, however, a kilonova was in the search area but was too weak to be detected, a constraint on the viewing angle can be obtained. We find that the viewing angle must have been 53 ± 10 deg, assuming there was an AT 2017gfo-like blue component. In our kilonova model there are two sources of uncertainty: one linked to the polar-to-equatorial view contrast (∆M λ in Eq. 4) and one linked to the calibration of the polar magnitude (M λ ). The former was fit to theoretical expectations from sophisticated modeling and the latter from calibration on GW170817. All of our results are sensitive to these procedures, though we underline that our general conclusions (see Sect. 7) should remain valid and our constraints on GW190425 should be seen as proof of concept. This method will be most useful in the case of genuine nondetections when the future, smaller GW sky maps will effectively be fully covered by follow-up. For future GW triggers, one could use the GW-measured progenitor properties (such as component masses and tidal deformability) to tailor the kilonova modeling to the specific events (Nicholl et al. 2021). In case of a non-detection, the viewing angle constraints thus obtained would be more robust because they are informed with the complete multi-messenger data set.
Both aspects of this uncertainty should improve in the coming years with the detection and observation of even a limited sample of kilonovae following GW signals, allowing us to explore both their intrinsic diversity and their properties under different viewing angles . When the burst afterglow is also detected, information on the external density and a better estimate of the viewing angle can be obtained, which might be completed, in the long term, with the potential observation and leveraging of the kilonova afterglow (Hotokezaka et al. 2018;Nakar et al. 2018;Kathirgamaraju et al. 2019).
Finally, as we have noted, kilonovae and even mild associated viewing angle measurements seem to be the only means for electromagnetic modeling to contribute to multimessenger cosmology and the resolution of the Hubble tension (Mastrogiovanni et al. 2020). The other counterparts are ruled out for their rareness. The effort to collect a kilonova sample and study source variability and viewing angle properties thus appears even more desirable in this regard.

Conclusion
We have presented a population study of kilonova counterparts to GW binary neutron star merger signals based on a simple viewing-angle-dependent model deduced from state-of-the-art modeling and calibrated on AT 2017gfo.
For shallow searches, the rate and viewing angle properties of the detected kilonova population are independent of the GW sensitivity. However, deep searches by target-ofopportunity endeavors and high-cadence surveys can probe a high-inclination population, detecting tens of events per year with design-type GW observing runs. Deep surveys will, however, be dominated by non-GW-triggered (orphan) kilonovae with possible short GRB associations.
We have proven the concept of constraining the inclination angle of systems in the case of a non-detection of the kilonova counterpart. Our method will become more effective in the case of a genuine non-detection when future, smaller GW sky maps are fully covered by follow-up.
Our results would be refined with a better understanding of kilonova emissions, and dedicated efforts to collect a sample are warranted. Such efforts are further motivated given the potential role of kilonovae in precision multimessenger cosmology.