On August 17th, 2017, during the second run of Advanced LIGO (LIGO Scientific Collaboration et al. 2015) and Advanced Virgo (Acernese et al. 2015) collaborations, the first gravitational-wave (GW) signal from the merger of two neutron stars (BNS) was detected, GW170817. Multiple independent detections associated with this event were made throughout the entire electromagnetic (EM) spectrum, including a short gamma-ray burst (Abbott et al. 2017b; Goldstein et al. 2017; Savchenko et al. 2017; Alexander et al. 2017; Haggard et al. 2017; Hallinan et al. 2017; Margutti et al. 2017; Coulter et al. 2017; Troja et al. 2017 ) and the kilonova (kN) AT2017gfo found at optical and IR wavelengths (Cowperthwaite et al. 2017; Kasliwal et al. 2017; Tanvir et al. 2017; Soares-Santos et al. 2017; Andreoni et al. 2017; Kasliwal et al. 2019) . A record number of studies were carried out on this single event, and many open questions remain with regards to the rates of neutron mergers, the nucleosynthesis of heavy elements, the mass range and equation of state of neutron stars, the diversity in the kN explosion population and their use as BNS mergers as distance yardsticks E-mail: ana.sagues-carracedo@fysik.su.se in cosmology. At only ∼ 40 Mpc distance, in the near-by galaxy NGC 4993, and close to face-on towards us (Hotokezaka et al. 2019) , GW170817/AT2017gfo was extremely favourable for detection.
On April 1st, 2019, LIGO/Virgo started the third observing run (O3) on the search of GWs. Since then, astronomers work closely with the LIGO and Virgo collaboration (LVC), who provide search maps with their public alerts used by the astronomy community to carry out target of opportunity (ToO) follow-up programs. One year later, by March 27 2020, the O3 run was prematurely suspended due to the COVID-19 pandemic. In spite of 14 merger alerts been issued by the LVC including at least one neutron star, i.e., either BNS or neutron star -black hole (NS-BH), no kN was securely identified. In this work we explore possible explanations for the lack of optical counterparts, in spite of rapid responses by the optical astronomy community to the LVC alerts. We quantify the kN detection probability for the typical specifications of the follow-up surveys adopted, and suggest potential improvements. Furthermore, we investigate the feasibility of serendipitous kN detection, independently of LVC triggers.
Significant volumes of space are explored with high cadence thanks to the advent of optical sky surveys such as PanSTARRS (Kaiser et al. 2010) , the AllSky Automated Survey for Supernovae (ASASSN, Shappee et al. 2014) , the Dark Energy Survey (DES, Dark Energy Survey Collaboration et al. 2016) , the Asteroid Terrestrial impact Last Alert System (ATLAS, Tonry et al. 2018) , the Zwicky Transient Facility (ZTF, Graham et al. 2019) , the Gravitational-Wave Optical Transient Observer (GOTO, Dyer et al. 2018) , the Mobile Astronomical System of TElescope Robots (MAS-TER, Lipunov et al. 2010) , the deca-degree optical transient imager (DDOTI, Watson et al. 2016) , MeerLICHT (Bloemen et al. 2016) , KMT-Net (Kim et al. 2016) , the Rapid Action Telescope for Transient Objects (TAROT, Boër 2001) or the Global Rapid Advanced Network Devoted to the Multi-messenger Addicts (GRANDMA, Antier et al. 2020a,b) . These surveys have also been carrying out regular targeted observations searching for potential EM counterparts of GW alerts in the regions posted by the LVC during O3, whenever the weather and the observable region of the sky made it possible.
In this work, we combine theoretical kN models with survey simulation tools to explore the detection probability. We focus on optical wavelengths, making use of the ZTF gri filter system, at effective wavelengths 4814, 6422 and 7883Å, as a representative set of photometric observations. Viewing-angle dependent models created with the radiative transfer code possis (Bulla 2019) were generated, and we choose a fiducial model that best fits the photometric data of AT2017gfo as our bench mark, also when comparing other models.
In Section 2, we describe the properties of synthetic kN lightcurves from possis. In Section 3, we present the details of the observing strategies explored, along with the distribution of transients generated and the detection criterion considered. Section 4 focuses on the detection probabilities for the LVC alerts in O3. In Section 5, we describe the population of synthetic kNe selected under various scenarios. Section 6 considers attainable constraints on kN rates from serendipitous searches, unrelated to LVC alerts. Finally in Section 7, we discuss and summarize the results.
In this work, we consider kN models synthesized using the Monte Carlo radiative transfer code possis (Bulla et al. 2015; Bulla 2019) . Assuming homologous expansion and analytic functions for the wavelength-and time-dependence of the opacity, possis performs time-dependent radiative transfers simulations for multi-dimensional models. Timedependent spectral energy distributions (SEDs) are computed for different viewing angles and used to construct multi-band lightcurves.
Compared to the version of possis in Bulla (2019), here we use an improved version, where the temperature is no longer parameterized and uniform throughout the ejecta, but estimated self-consistently from the mean intensity of the radiation field at each time and in each zone. Furthermore, we adopt thermalization efficiencies th from Barnes et al. (2016) rather than assuming th = 0.5 as in Bulla (2019) . For more information about the code, we refer the reader to Bulla (2019) .
A two-dimensional geometry tailored to BNS mergers is used, comprising two ejecta components as in Bulla (2019): a lanthanide-rich component distributed within an angle ±φ around the merger/equatorial plane, and a lanthanide-free component at higher latitudes (see fig. 1 in Bulla 2019). A density profile ∝ −3 , where is the velocity, is adopted for both components between min = 0.025c (instead of min = 0 as in Bulla 2019) and max = 0.3c (where c is the speed of light). The three free model parameters explored in this work are:
• the total ejecta mass, m ej • the half-opening angle, φ • the viewing angle, Θ
For m ej , we consider 0.01 to 0.10 M with a mass step of 0.01 M , and compute half-opening angles, φ, from 15 to 75 • with steps of 15 • . A total number of 50 simulations are generated, each comprising predictions for N obs = 11 viewing angles Θ from pole (cos Θ = 1) to equator (cos Θ = 0). SEDs used in this work are made available at https://github.com/mbulla/kilonova_models.
The ugrizyJH photometric data of AT2017gfo (Andreoni et al. 2017; Arcavi et al. 2017; Chornock et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Evans et al. 2017; Kasliwal et al. 2017; Pian et al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Troja et al. 2017; Utsumi et al. 2017; Valenti et al. 2017) can be reproduced by a fiducial model with m ej = 0.06 M and φ = 60 • for viewing angles close to polar, as shown in Fig. 1 for the three filters of interest in this work, gri. The figure also shows that kNe are fast-evolving transients with absolute peak magnitudes around -16 mag, i.e., relatively faint compared to many types of supernovae.
The absolute peak magnitudes as a function of ejecta mass ranges from -12.6 to -17.4 mag, with a decay rate for 3 days after peak (∆m 3 ) varying from 0.4 to 2.2 mag, depending on the inclination angle and the broadband filter used, as shown in Fig. 2 . For the fiducial model (m ej = 0.06 M and φ = 60 • ), the optical peak magnitudes vary from -13.8 to -16.1 mag, and ∆m 3 from 0.6 to 1.5 mag for different viewing angles.
Finally, Fig. 3 shows the observed r-band magnitude as a function of redshift, corresponding to the range of absolute magnitudes at peak, emphasizing the faintness of kNe, and thus the limited redshift range that can be probed with small telescopes.
Brighter lightcurves result from larger m ej , smaller φ, and Θ close to viewing angles towards the pole. The difference in brightness and shape of the synthetic kN lightcurves is less pronounced at redder wavelengths. ∆m 3 is higher for shorter wavelengths and smaller m ej , being the difference small when varying φ.
We perform survey simulations using simsurvey 1 (Feindt et al. 2019 ), a python based simulation tool for astronom- ical surveys. simsurvey uses Monte Carlo methods to generate synthetic lightcurves of transients, e.g., supernovae or kNe. It uses a timeseries of SEDs for the transient model, from which the photometry is calculated for a given filter, including cosmological redshift. The input for the simulation is a survey schedule containing the information about cadence, field of view, pointings, filter choice, sky brightness and zero points of each exposure. simsurvey is a tool developed within the ZTF collaboration, and it has been used successfully for forecasting rates of supernovae (Feindt et al. 2019) .
Besides including the range of kN models with random orientations described above, we also account for dimming by dust in the interstellar medium of the host galaxy, as well as Galactic extinction based on the Schlegel, Finkbeiner, Davis (SFD98) reddening maps. In both cases, we use the standard Milky-Way extinction law for the wavelength dependence (Cardelli et al. 1989) . For this work, we have picked a region of the sky with a mean Milky Way extinction of E(B − V) = 0.08 mag. The output of simsurvey is a python dictionary that contains all the parameters of lightcurves for both detected and undetected transients.
We simulate ToO triggers under optimal conditions, i.e., shortly after the merger of two neutron stars. The synthetic observations start between 2.4 to 72 hours after the merger, and continue for seven days, with observations twice per night. Although we use ZTF gri filters for the simulations, we expect little dependence on the specifics of these filters, i.e., our findings should be applicable to most optical followup campaigns during O3. We simulate lightcurves down to a hypothetical m lim = 25 mag in all filters. Magnitudes are in AB system throughout this work.
We generate a fixed number of synthetic kN lightcurves randomly distributed in the field of view of the instrument. Transients are placed at distances sampled from a distribution based on the volumetric rate as a function of redshift up to z = 0.15.
The kN viewing angle is randomized, i.e., has constant probability for cos Θ. Thus, because of the assumed axial symmetry in the models, viewing angles closer to the merger plane are more probable to happen than those closer to the pole.
Next, we describe the adopted detection criterion. In dedicated searches of kNe following a GW alert, a basic requirement is non-detection prior to the merger time, as well as non-detection after sufficient time has elapsed, beyond which the lightcurves are expected to have faded. To have a nondetection we require a prior visit with no signal passing the threshold. However, in this study we simplify the situation in that we do not simulate data before the merger, nor after the lightcurves fade. Furthermore, we separate the synthetic observations by at least one hour, a common strategy to discard potential false positives from asteroids. We require at least two photometric detections with a SNR ≥ 5.
LVC alerts involving neutron stars from LVC in O3, 6 BNS and 8 NS-BH mergers, were followed-up by ground-based optical/NIR telescopes around the globe, following the ability of each telescope to cover the merger location and weather conditions permitting. The searches were performed trying to optimize coverage of the candidate skymap provided by the LVC. Unfortunately, no identification of a kN was made, nor any other EM detection that could be securely associated to the GW signal. In this section, we explore the detection probability for kNe at the distances reported in the LVC alerts, with the purpose to explore possible explanations for the non-detections. Since the modelled grid are tailored to BNS mergers (Section 2), we will focus on the six BNS triggers listed in Table 1 . The efficiency for detection of NS-BH counterparts have to be treated separately, as their lightcurves are predicted to be brighter and longer lasting compared to those from BNS mergers (Rosswog et al. 2017; Barbieri et al. 2019; Kawaguchi et al. 2020) . S190425z was the first BNS candidate in O3. The LAL-Inference localization pipeline (Veitch et al. 2015) provided 7461 deg 2 in the 90% credible region and an estimated luminosity distance of 156 ± 41 Mpc. The day after, a new alert for a possible BNS was sent (S190426c), with a better localization (1131 deg 2 ) but more than twice further away (377 ± 100 Mpc). These first two candidates were followed-up within the GROWTH collaboration (Coughlin et al. 2019 ,Goldstein et al. 2019 , probing the fast response to the low latency skymap and coverage capacity. No alerts with BNS candidates were released in the following few months until the beginning of September with S190901ap and S190910h. The skymap for both of these events was poorly localized,14753 deg 2 and 24264 deg 2 , and the luminosity distances inferred to be 241 ± 79 Mpc and 233 ± 88 Mpc, respectively. These were the last events detected before the break in October and the end of O3a. During O3b, two non-retracted BNS candidates were sent, S191213g and S200213t, both at a mean distance of 201 Mpc (± 81 and 80 Mpc, respectively). The closest BNS candidate in O3, at 156 Mpc, is four times more distant than GW170817 (40 Mpc), and the rest of the events happened more than five times further away. Another remarkable factor is the poor localization comparing with GW170817, 28 deg 2 90% confident area (Abbott et al. 2017a) . We seek to understand which is the limiting magnitude required to detect a kN. We study the detectability at different luminosity distances. We focus on the results at 200 Mpc as a representative distance for the BNS candidates in O3.
The efficiency is calculated from the fraction of the synthetic lightcurve population that fulfill the adopted detection criterion. As expected, the kN detection probability improves for short distances and viewing angles closer to the pole, i.e., where the line of sight points to the lanthanidefree region. Fig. 4 shows the detection efficiency for limiting magnitudes 21, 22 and 23 mag, assuming the fiducial model matching AT2017gfo, where the only free parameters are distance and viewing angle. We also indicate the estimated mean distances from the GW signals of the six BNS trigger by LIGO/Virgo during O3. The shape of the 90% probability contours shows the dependency with the viewing angle for the fiducial model, the probability drops significantly towards equatorial viewing angles in which the photons are heavily absorbed by the lanthanides. As can be appreciated, unlike GW170817, none of the much more distant six LVC BNS events during O3 had a sizeable probability of detection at a limiting magnitude of 21 mag, which was very typical for the O3 follow-up campaigns. For 22 mag, only favourable viewing angles can be probed, whereas surveys reaching 23 mag would have had a good chance to detect a counterpart. We also find that redder wavelength coverage, like i-band reaching the same limiting magnitude, significantly increases the volume in which a kN is likely observable. Table 1 , shows the minimum limiting magnitude at which the integrated detection probability over all viewing angles is more than 90% at distance of the BNS merger candidates from LVC in O3 assuming the fiducial model.
The probability of kN detection depends on the assumed physical parameters of the model. Fig. 5 indicates the estimated probabilities at 200 Mpc for observations with gr and gri, all at an assumed limiting magnitude of 22 mag. Results for other limiting magnitudes can be found in the appendix A. The probabilities are calculated for two subsets of detections: events with viewing angles pointing to the lanthanide-free component Table 2 shows probabilities at 40 Mpc (GW170817), 100 Mpc and 200 Mpc, for the fiducial model and two extremes cases, i.e., models that provide the brightest and the faintest kN.
The detection probability strongly depends on the physical parameters of the kN and the orientation, i.e., large m ej , small φ, and Θ LF are associated with brighter kNe and thus provide higher probability values.
The observations are constraining on the parameter space of the model that provides high detection probability. We aim to be sensitive to a wide range of parameters with our observations. The inclination of the merger is unknown when the ToO is triggered. Viewing angles close to face-on provide the brighter lightcurves, but they are less likely to happen (see section 3.2). Therefore, it is preferable to plan for less desirable circumstances. 23 mag depth in gri is needed to constraint most of the parameter space at 200 Mpc.
Next, we explore the model properties for the population of kNe that pass the detection criterion for a given magnitude limit. We group the viewing angles of the detected sample in two sets: one with orientation such that the line of sight points to the lanthanide-free component (Θ LF ) and the other pointing to the lanthanide-rich component (Θ LR ), defined in Section 4.2. We also separate the detections for synthetic observations into those that include i-band and those that do not. To better assess the selected samples under different considerations, we scale the number of detections to a BNS rate of 1000 Gpc −1 yr −1 (Scolnic et al. 2018) for 20 000 sq. degrees. The summary of our findings is shown in Table 3 . Given the red nature of kNe, including i-band to the same depth as g and r in the survey leads to more than 50% percent improvement on the number of detections. kN lightcurves show a slower decay rate in redder bands (smaller ∆m 3 , see Fig. 2 ), they are less affected by the flux suppression of the lanthanide-rich component and are less affected by dust extinction.
The increase in detections with limiting magnitude is exponential, as more volume becomes available, V ∝ 10 0.6·m lim . In our simulations, we find that g-band observations, although useful for background rejection and characterization of the kN properties, do not add to the detection efficiency. All the kNe detected with SNR ≥ 5 do also pass the detection criterion in r and i-band.
The distribution of detections in viewing angle and distance is shown in Fig. 6 . The detected population is affected by Malmquist bias (Malmquist 1922) , i.e., given a magnitude limit, only the brightest objects in the population are detected. For the synthetic population, as volume increases with distance, most kNe are generated at farther distances. Furthermore, equatorial viewing angles are more common (see Section 3.2). Thus, the brightest events correspond to lower distances and polar viewing angles. The combination of the two effects gives a median viewing angle of about 40 • for the population of selected kNe, assuming the fiducial model.
When searching for transients in the survey alerts, astronomers look for at least two data points with a sufficient signal to noise and time separation to consider a reliable discovery. The chances to observe a lightcurve that can be identified as a kN increase when the transient is observable for a more extended period. The duration of the observable lightcurve depends on the properties of the kN and vary with distance. Fig. 7 shows the range from the minimum to the maximum length of the detected lightcurves per distance for the fiducial model, being length defined as the post-merger epoch of the last 5σ data point.
Next, we explore the feasibility for kN detections in timedomain surveys, irrespective of external GW triggers, for a volumetric BNS rate of 1000 Gpc −3 yr −1 , consistent with Scolnic et al. (2018) . We make several simplifying assumptions, including losses due to candidate selection inefficiencies, bad weather, or other observing conditions. Every kNe observable within the covered volume with lightcurves of at least two 5σ data points with a minimum separation of one hour is detected. Furthermore, a constant magnitude limit is assumed for the entire duration of an all-sky survey, i.e., 20 000 sq.degrees. We assume that the distribution of BNS is homogeneous on the sky and rescale the number of detections using efficiencies calculated in Section 4. With the fiducial model, we find that for m lim = 20.5 mag, 3 kNe are observable per year, around 20 for limiting magnitude of 22 mag and 71 for 25 mag, assuming constant depth with gri wavelength coverage (see Fig. 8 and Table 3) .
The detection of GW-EM counterparts has proven to be very challenging. During the LIGO/Virgo O3 run, six BNS triggers were issued, and in spite of large efforts from the entire time-domain astronomy community, no transient resembling AT2017gfo was found. Furthermore, fast-evolving and faint kN transients are too rare to facilitate serendipitous discoveries. The fast evolution in bolometric luminosity and color typically leads to lightcurves with few data points, hard to identify photometrically, and even more to get a good quality spectrum. In some cases, the lightcurve would fade below threshold within the first day from the merger. S190425z S191213g
S200213t S190910h S190901ap S190426c S190425z S191213g
S200213t S190910h S190901ap S190426c 0. 90 0 0.9 00 S190425z S191213g
S200213t S190910h S190901ap S190426c S190425z S191213g
S200213t S190910h S190901ap S190426c 0. 90 0 0 .9 0 0 S190425z S191213g
S200213t S190910h S190901ap S190426c S190425z S191213g
S200213t S190910h S190901ap S190426c 0 . 9 0 0 0 . 9 0 0 0 1 gri gr 0 1 Table 3 . Summary of simulation statistics. n KNe, gr and n KNe, gr i correspond to the number of detections scaled to a given BNS rate of 1000 Gpc −3 yr −1 (Scolnic et al. 2018) for 20 000 sq. degrees. Θ LF, gr and Θ LF, gr i correspond to the fraction of detections with viewing angles towards the lanthanide free component (Θ LF ) observing with gr and gri.Θ gr andΘ gr i are the mean viewing angle in degrees, andD L, gr andD L, gr i , the mean luminosity distance in Mpc of the detected kNe. lightcurve coverage demanding SNR ≥ 5, for limiting magnitudes ranging from 21 to 25 mag. For a given limiting magnitude, the efficiency is reduced when requiring longer lightcurves. E.g., for m lim = 23 mag, there is more than 90% probability for detecting a kilonova at 200 Mpc over 3 days since the merger. Requiring a duration of 5 days of the ob-servable lightcurve, the efficiency drops down to about 10%. For viewing angles towards the lanthanide-free component, reaching 22 mag is required in gri, and 23 mag for the viewing angles towards the lanthanide-rich component, to reach 90% probability with a 3 days lightcurve.
In addition to the intrinsic observational challenges posed by the rapid evolution of kN lightcurves, the distances to the BNS triggers during the third run of LIGO and Virgo have been around 200 Mpc, about five times larger than for GW170817. Some of the events were detected with a single interferometer, which translates into poor localization. Large search areas, in combination with long distances, poses great challenges for the detection of EM counterparts. Survey scheduling under such conditions is even harder, considering the lack of knowledge of the underlying kN population. During O3, the effort to cover these large search areas on the sky was remarkable (Gompertz et al. 2020; Antier et al. 2020b; Andreoni et al. 2020; Ackley et al. 2020; Vieira et al. 2020; Watson et al. 2020; Gompertz et al. 2020; Goldstein et al. 2019 ), but the adopted depth was generally insufficient to identify the counterpart, according to the models we have studied. Several collaborations such as DECam-GROWTH 2 (Anand et al. 2020) , ENGRAVE (Ackley et al. 2020) or GRANDMA (Antier et al. 2020a,b) had triggered ToO observations with a global multi-telescope network when it was feasible. It has been shown that current capabilities can afford to cover large areas on the sky, but the reached depth was insufficient to detect the EM.
Clearly, GW170817 appears to have been an exceptionally lucky occurrence, both with respect to distance, localisation and orientation, well outside the bulk of the BNS mergers in O3 (see Fig. 6 ). Due to their low rate, e.g., com-pared with other types of transients like supernovae, the probability of another nearby kN with viewing angle close to the pole, is very low. The asymmetry of the kN emission and the higher likelihood of viewing angles close to the merger plane means that most of the kNe are going to be produced under conditions much less favourable for detection than GW170817.
In agreement with previous studies of follow up campaigns during O3 (Coughlin et al. 2020 ,Hosseinzadeh et al. 2019 , we conclude that deeper observations are required to ensure high detection probability. At a distance of 200 Mpc, the required magnitude to observe a kN similar to AT2017gfo at lanthanide-rich viewing angles with three days lasting lightcurve is at least 23 mag. For an AT2017gfo-like event at that distance, pointing towards us with a viewing angle within the lanthanide-free component, m lim = 22 mag could be sufficient. This values change for different physical parameters of the kN model (see Section 2.1).
The optical lightcurves of kNe become very red in just a few days. The nucleosynthesis of heavy nuclei translates into longer-lived emission at longer wavelengths. The blue emission fades rapidly, and it suffers more from absorption at viewing angles close to the merger plane, where the dynamical ejecta leaves most of the lanthanides. Thus, observations in red filters tend to be much more efficient. In this work we show that adding i-band observations in an optical survey yields, on average, more than two times more detections than e.g., observing with only g-and r-band. ing luminosities, continued deep observations over multiple days is required, something that is hard to achieve with current instruments due to unpredictable weather patterns and moon phase among other factors. The Vera C. Rubin Observatory (VRO/LSST) with the Legacy Survey of Space and Time (LSST, LSST Science Collaboration et al. 2009) with its 8-meter-class telescope, planned to be online by 2023, will discover transients that are too far or too faint for the current available surveys, including serendipitous discoveries of kNe (Setzer et al. 2019) .
In this work, we have used models tailored to BNS. However, we expect that results for NS-BH lightcurves would change quantitatively but not qualitatively. In which case, the lack of EM counterpart detections for the distant NS-BH LVC candidates during O3 would also be justified.
In conclusion, detection and identification of GW electromagnetic counterparts from BNS mergers is very challenging. The survey magnitude depths during follow-up campaigns in O3 have been generally insufficient, considering the asymmetry of the kN emission, the sometimes poor localization and distances, typically 5 times further than GW170817. Despite the lack of kN detections during O3, the multimessenger astronomy community have demonstrated a remarkable ability to respond promptly to LVC alerts, even those with wide search areas. Thus, the future is promising, also as more detectors join the gravitational waves observatory network, leading to better localization. Thus, we expect that with improved search strategies described here, future campaigns will be successful. 