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

# Where and When: Optimal Scheduling of the Electromagnetic Follow-up of Gravitational-wave Events Based on Counterpart Light-curve Models

Om Sharan Salafia¹,²,³✉, Monica Colpi¹,³, Marica Branchesi⁴,⁵, Eric Chassande-Mottin⁶, Giancarlo Ghirlanda¹,², Gabriele Ghisellini², and Susanna D. Vergani²,⁷

¹ Università degli Studi di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy; omsharan.salafia@brera.inaf.it, omsharan.salafia@gmail.com

² INAF—Osservatorio Astronomico di Brera Merate, via E. Bianchi 46, I-23807 Merate, Italy

³ INFN—Sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

⁴ Universit degli Studi di Urbino “Carlo Bo,” via A. Saffi 2, I-61029, Urbino, Italy

⁵ INFN, Sezione di Firenze, via G. Sansone 1, I-650049, Sesto Fiorentino, Italy

⁶ APC, Université Paris Diderot, CNRS/IN2P3, CEA/Irfu, Obs de Paris, Sorbonne Paris Cité, France

⁷ GEPI, Observatoire de Paris, PSL Research University, CNRS, Place Jules Janssen, F-92195 Meudon, France

Received 2017 April 19; revised 2017 August 1; accepted 2017 August 6; published 2017 August 31

## Abstract

The electromagnetic (EM) follow-up of a gravitational-wave (GW) event requires scanning a wide sky region, defined by the so-called “skymap,” to detect and identify a transient counterpart. We propose a novel method that exploits the information encoded in the GW signal to construct a “detectability map,” which represents the time-dependent (“when”) probability of detecting the transient at each position of the skymap (“where”). Focusing on the case of a neutron star binary inspiral, we model the associated short gamma-ray burst afterglow and macronova emission using the probability distributions of binary parameters (sky position, distance, orbit inclination, mass ratio) extracted from the GW signal as inputs. The resulting family of possible light curves is the basis for constructing the detectability map. As a practical example, we apply the method to a simulated GW signal produced by a neutron star merger at 75 Mpc whose localization uncertainty is very large (~1500 deg²). We construct observing strategies for optical, infrared, and radio facilities based on the detectability maps, taking VST, VISTA, and MeerKAT as prototypes. Assuming limiting fluxes of r ~ 24.5, J ~ 22.4 (AB magnitudes), and 500 μJy (1.4 GHz) for ~1000 s of exposure each, the afterglow and macronova emissions are successfully detected with a minimum observing time of 7, 15, and 5 hr respectively.

**Key words:** gamma-ray burst: general – gravitational waves – methods: statistical – stars: binaries – stars: neutron

## 1. Introduction

The first detection (The LIGO Scientific Collaboration & the Virgo Collaboration 2016c) of gravitational waves (GWs hereafter) from the inspiral and merger of a black hole binary, followed by a second (The LIGO Scientific Collaboration & the Virgo Collaboration 2016b) and a third one (The LIGO Scientific Collaboration et al. 2017), suddenly turned these fascinating, theoretical objects into real astronomical sources.

When such a compact binary coalescence is detected, the analysis of the GW signal and comparison with carefully constructed templates of the waveform (Abbott et al. 2016a; The LIGO Scientific Collaboration & the Virgo Collaboration 2016a) enable the extraction of precious information about the parameters of the binary and of the remnant. The identification of an electromagnetic (EM) counterpart would further increase the scientific outcome of the detection by, e.g., enabling the identification of the host galaxy, providing hints about the environment surrounding the merger, constraining the theoretical models of EM counterparts, and reducing degeneracies in the GW extrinsic parameter space (Pankow et al. 2016).

Several observatories, covering a large fraction of the EM spectrum, recently developed dedicated programs for the EM follow-up of GW events. The present main limitation for the detection of a possible EM counterpart is the large uncertainty on the sky localization of the GW source (see, e.g., Berry et al. 2014; Singer et al. 2014). The problem might be alleviated by targeting bright galaxies within the localization uncertainty region (Nuttall & Sutton 2010; Abadie et al. 2012), and the selection of target galaxies can also take into account the

sky-position-conditional posterior distribution of the source luminosity distance (Hanna et al. 2013; Nissanke et al. 2013; Gehrels et al. 2016; Singer et al. 2016a). The aim of this work is to propose an additional way to use the information encoded in the GW signal to optimize the follow-up strategy for each single event, namely, to combine posterior distributions of the compact binary parameters and available models of the EM emission to predict the best timing for the observation of different parts of the GW skymap. Such an approach can be applied in cases when a model of the expected EM counterpart is available, and it is especially useful when the light curve predicted by the model depends on (some of) the compact binary parameters.

### 1.1. The First Electromagnetic Follow-ups

The observation campaigns that followed up the first detections of GWs were very extensive. Hundreds of square degrees within the GW sky localization were covered by wide-field telescopes (Abbott et al. 2016b). Target areas were selected in order to maximize the contained GW source posterior sky-position probability, incorporating telescope visibility constraints (e.g., Kasliwal et al. 2016). In some cases, models of the expected EM counterpart emission were used to estimate the optimal search depth (e.g., Soares-Santos et al. 2016); other searches combined the posterior sky-position probability map with the areal density and luminosity of nearby galaxies to select the best target fields (e.g., Díaz et al. 2016; Evans et al. 2016). Observations were concentrated during the first days after the events and repeated weeks to months
later to search for both rapid and slowly evolving possible counterparts.

## 1.2. Candidate EM Counterparts

It is not clear whether an EM counterpart should be expected in the case of a binary of black holes (BH-BH), due to the unlikely presence of matter surrounding the binary (but see Loeb 2016; Perna et al. 2016; Yamazaki et al. 2016); on the other hand, if the merger involves a black hole and a neutron star (BH-NS) or two neutron stars (NS-NS), there are solid reasons to believe that EM emission should take place. The most popular mechanisms for such an emission in both BH-NS and NS-NS cases include prompt (gamma-ray) and afterglow (panchromatic) emission from a short gamma-ray burst (SGRB) jet, and “macronova” (optical/infrared) emission from ejecta launched during and after the merger, powered by the decay of unstable heavy nuclei resulting from r-process nucleosynthesis taking place within the neutron-rich ejecta during the early expansion phase.

Many other promising EM counterparts have been proposed, e.g., the long-lasting radio transient (Nakar & Piran 2011) arising from the deceleration of the dynamical ejecta due to interaction with the interstellar medium (ISM), the jet cocoon emission (Lazzati et al. 2017; Gottlieb et al. 2017), or the spin-down-powered emission described by Siegel & Ciolfi (2016) in the case when a (meta-)stable neutron star is left after the merger. To keep the discussion as simple as possible, in this paper we will only consider the (optical and radio) SGRB afterglow and the dynamical ejecta macronova as examples, leaving the possibility of applying the present approach to other EM counterparts to future works.

## 1.3. The SGRB Afterglow

The detectability of the SGRB prompt emission depends crucially on the jet viewing angle $\theta_v$, i.e., the angle between the jet axis and our line of sight. If the viewing angle is larger than the jet half-opening angle $\theta_{jet}$ (in other words, if the jet points away from the Earth), the prompt emission flux received by an observer on Earth is severely suppressed (e.g., Salafia et al. 2016) due to relativistic beaming (by the compactness argument, the bulk Lorentz factor in GRB jets must be comparable to or larger than one hundred—e.g., Lithwick & Sari 2001—and estimates based on observations are sometimes even larger than a thousand—as in the short burst GRB090510; see Ghirlanda et al. 2009; Ackermann et al. 2010). Since the typical half-opening angle $\theta_{jet}$ is somewhere between 5° and 15° (e.g., Berger 2014), the prompt emission goes undetected in the majority of cases (for an isotropic population, the probability that $\theta_v < 15°$ is less than 2%). Soon after producing the prompt emission, the jet starts interacting significantly with the ISM, and a shock develops (Meszaros & Rees 1996). Electrons in the shocked ISM produce synchrotron radiation, giving rise to a fading afterglow (observed for the first time by Beppo-SAX; Costa et al. 1997). Since the consequent deceleration of the jet reduces the relativistic beaming, an off-axis observer (who missed the prompt emission) could in principle detect the afterglow before it fades (Rhoads 1997): in this case, the afterglow is said to be an orphan. No convincing detection of such a transient has been claimed to date, consistent with predictions for current and past surveys (Ghirlanda et al. 2014,

2015), but future deep surveys (e.g., MeerKAT in the radio—Booth et al. 2009, LSST in the optical—Ivezic et al. 2008, eROSITA in the X-rays—Merloni et al. 2012) are anticipated to detect tens to thousands of such events per year.

Given the large uncertainty on the expected rate of NS-NS and BH-NS detections by the aLIGO and Advanced Virgo facilities in the near future (LIGO Scientific Collaboration et al. 2010; Dominik et al. 2015; Kim et al. 2015; de Mink & Mandel 2016; Abbott et al. 2016) and the rather low expected fraction of GW events with an associated SGRB jet pointing at the Earth (Metzger & Berger 2011; Wanderman & Piran 2014; Ghirlanda et al. 2016; Patricelli et al. 2016), the inclusion of orphan afterglows as potential counterparts is of primary importance to test the SGRB-compact binary coalescence connection.

## 1.4. The Dynamical Ejecta Macronova

Despite the idea dating back to almost 20 years ago (Li & Paczyński 1998), the understanding of a possible macronova emission following a compact binary merger has been expanded relatively recently, as a result of the combined effort of researchers with expertise in a wide range of areas. A non-exhaustive list of the main contributions should include

1. numerical simulations of the merger dynamics (relativistic simulations by many groups using different approaches—e.g., Rezzolla et al. 2010; Tanaka & Hotokezaka 2013; Kiuchi et al. 2014; Wanajo et al. 2014; Bauswein & Stergioulas 2015; East et al. 2015; Giacomazzo et al. 2015; Just et al. 2016; Radice et al. 2016; Ruiz et al. 2016; Sekiguchi et al. 2016, 2015; Ciolfi et al. 2017; Dietrich et al. 2017—and non-relativistic simulations, especially by Stephan Rosswog and collaborators—e.g., Rosswog et al. 2014);

2. studies to assess the efficiency of r-process nucleosynthesis and the consequent heating rate due to heavy-element decay in various ejecta (Freiburghaus et al. 1999; Rosswog et al. 2000; Korobkin et al. 2012; Wanajo et al. 2014; Hotokezaka et al. 2015; Lippuner & Roberts 2015; Eichler et al. 2016; Rosswog et al. 2016);

3. atomic structure modeling which revealed the role of lanthanides in the ejecta opacity evolution (Kasen et al. 2013);

4. simulations including neutrino physics to model the neutrino-driven wind and the associated macronova (Dessart et al. 2008; Martin et al. 2015; Perego et al. 2017).

Results (especially for the dynamical ejecta) from various research groups are beginning to converge, and the dependence of the emission features on the parameters of the binary is in the process of being understood. Both analytical and numerical models capable of predicting the light curve have been developed recently (Barnes & Kasen 2013; Grossman et al. 2013; Dietrich & Ujevic 2016; Kawaguchi et al. 2016; Barnes et al. 2016; Rosswog et al. 2017). The emission from the dynamical ejecta is generally thought to be isotropic, which is an advantage compared to the SGRB afterglow (which is instead beamed) from the point of view of EM follow-up. The energy reservoir is the ejected mass $M_{ej}$, which depends most prominently on the mass ratio $q = M_1/M_2$ of the binary and on the neutron star compactness, which in turn reflects the mass of
the neutron star and its equation of state (EoS). Exciting claims of the detection of possible macronova signatures in the afterglows of few short GRBs (Tanvir et al. 2013; Jin et al. 2015, 2016; Yang et al. 2015) are in the process of being tested by intensive observational campaigns. All of this makes the macronova emission an extremely interesting candidate EM counterpart.

## 1.5. *Outline of This Work*

In Section 2.1, we introduce the idea of a follow-up strategy as a collection of observations that partially fill a “search volume” (search sky area × typical transient duration), stressing that the GW “skymap” (the sky-position probability density) gives information about where to observe, but not about when. In Section 2.1.1, we show how a priori information about the EM counterpart can be used to quantitatively define how likely the detection of the EM emission is if the observation is performed at time t, thus providing some information about how to explore the temporal dimension of the “search volume.” In Section 2.2, we suggest that the same approach can be extended to use a posteriori information extracted from the GW signal, provided that we have a way to link the properties of the inspiral to those of the EM counterpart (as shown in Section 2.2.2). In Section 2.2.4, we go one step further by introducing the idea that the information on the inspiral parameters that we can extract from the GW signal is dependent on sky position, and thus the clues (that we obtain from the GW signal analysis) about when to observe can also depend on the sky position. In Section 3, we introduce a method to extract such information from the “posterior samples” obtained from the analysis of a GW signal, and in Section 4, we apply it to a synthetic example to show how it optimizes the follow-up strategy. Finally, we discuss the results in Section 5, and we draw our conclusions in Section 6.

## 2. Where and When to Look

### 2.1. *A Sketch of the Design of a Follow-up Observation Strategy*

A short time after the detection of a compact binary coalescence signal, the LIGO Scientific Collaboration and Virgo Collaboration share information about the event with a network of astronomical facilities interested in the EM follow-up. The most fundamental piece of information for the follow-up is the so-called “skymap,” i.e., the posterior sky-position probability density, which we denote as $P(\alpha|S)$. It represents the probability per unit solid angle that the source is at sky position $\alpha$, say, $\alpha = (\text{R.A., decl.})$, given the GW signal $S$ detected by the interferometers ($S$ here represents all information contained in the strain amplitudes measured by all interferometers in the network). In what follows, we will most often call this probability density the “skymap probability.” Imagine that the EM counterpart appears at the GW position right after the event and never turns off. Assume that it can be found by comparison with previously available images of the sky, and that it can be easily identified by its spectrum or by another method. An ultra-simplified sketch of the obvious follow-up strategy would then be the following:

1. find the smallest sky area $A_\omega$ containing a large fraction $\omega$ (say, $\omega = 90\%$) of the skymap probability $P(\alpha|S)$;

2. divide such area into patches of size $A_{\text{FoV}}$ corresponding to the field of view of the instrument;

3. observe the patches in decreasing order of skymap probability$^8$;

4. for each patch:

(a) identify the new sources by comparison with archive images;

(b) perform a set of operations, including e.g., cross-matching with catalogs and spectral characterization, to discard known variable sources and unrelated transients in order to identify the counterpart.

The expected EM counterparts are transients, and thus a first modification to the above sketch must take into account the time constraints coming from our a priori knowledge of the transient features. If we have a physical or phenomenological model of the transient and we have some hint about the distribution of the parameters of such a model, we can construct a prior probability $P(F(t) > F_{\text{lim}})$ that the transient flux (in a chosen band) $F$ is above some limiting flux $F_{\text{lim}}$ at a given time $t$ after the GW event. Hereafter, we will call such a quantity the “a priori detectability.” The probability of detecting the transient at time $t$ by observing a sky position $\alpha$ with an instrument with field of view $A_{\text{FoV}}$ and limiting flux $F_{\text{lim}}$ is then

$$P(\det|t, \alpha, \text{FoV}) \sim A_{\text{FoV}} P(\alpha|S) \times P(F(t) > F_{\text{lim}}), \quad (1)$$

where we are assuming a relatively small field of view in order to consider $P(\alpha|S)$ constant over its area. This is nothing more than saying that the best place to look for the transient is the point of maximum skymap probability, at the time of highest a priori detectability. The probability of detection decreases both moving away from the point of maximum skymap probability and observing at a time when $P(F(t) > F_{\text{lim}})$ is smaller.

Let us work in the simplifying assumption that all observations have the same exposure $T$ and the same limiting flux $F_{\text{lim}}$. Let us denote by $T_{\text{tr}}$ the most conservative (i.e., largest) estimate of the transient duration, and let us define the “search volume” $V_\omega = A_\omega \times T_{\text{tr}}$ (we refer to this set as a “volume” because it is three dimensional, even though the dimensions are solid angle × time—see Figure 1). The follow-up strategy can then be thought of as an optimization problem, where one wants to (partially) fill the search volume $V_\omega$ with $N$ observations $V_{\text{obs},i} = A_{\text{FoV}}(\alpha_i) \times (t_{\text{obs},i}, t_{\text{obs},i} + T)$ (where $\alpha_i$ and $t_{\text{obs},i}$ are, respectively, the sky coordinates of the center of the field of view and the starting time of the $i$th observation) in order to maximize the detection probability $P(\det|\text{strategy})$, which can be written as

$$P(\det|\text{strategy}) = \sum_{i=1}^{N} \int_{A_{\text{FoV}}(\alpha_i)} P(\alpha|S) d\alpha \times P(F(t_{\text{obs},i}) > F_{\text{lim}}), \quad (2)$$

with the constraint $NT \le T_{\text{tr}}$ (see Figure 1).

The above paragraphs are essentially a formal description of the most basic follow-up strategy one can think of, which can be reduced to the principle “try to arrange the observations in order to cover the largest possible fraction of the GW skymap around the time when the flux is expected to be high enough for

<footnote>$^8$To keep the discussion as simple as possible, we are neglecting the limitations due to observing conditions.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Graphical representation of the follow-up strategy as a “volume filling” problem. $A_\omega$ represents the sky region that contains a fraction $\omega$ of the sky-position probability. The “search volume” is defined as the set of points $V_\omega = \{(\alpha, t)|\alpha \in A_\omega \text{ and } t \in (0, T_{tr})\}$. Observations are sets that intersect the search volume, defined by a field of view $A_{FoV}(\alpha_i)$ centered about sky position $\alpha_i$, an exposure time $T$, and an observation time $t_{obs,i}$ such that $V_{obs,i} = \{(\alpha, t)|\alpha \in A_{FoV}(\alpha_i) \text{ and } t \in (t_{obs,i}, t_{obs,i} + T)\}$. Observations made by the same instrument cannot overlap on the time axis (unless the instrument can see more than one field at the same time). For the detection to be successful, the EM counterpart must be located within one of the $A_{FoV}(\alpha_i)$ and its light curve must be above the detection threshold during the corresponding exposure time.</figcaption></figure>

a detection.” In this approach, the proper construction of the a priori detectability $P(F(t) > F_{lim})$ is key: it defines the time span within which the observations are to be performed, while the posterior sky-position probability density defines the search area.

### 2.1.1. How to Construct the A Priori Detectability

In order to construct the a priori detectability $P(F(t) > F_{lim})$, one must assume some prior probability density of the model parameters. Let us consider a simple, illustrative example. First, we construct a synthetic population of NS-NS inspirals whose properties roughly reproduce those expected for the population detected by Advanced LIGO; we then associate to each of them a jet afterglow and a macronova, under some assumptions. The detectable fraction of light curves in a given band, at a given time, will then constitute our estimate of the a priori detectability for this particular case. For the jet afterglow, we assume<fnref n="9" /> that all SGRB jets have an isotropic kinetic energy $E_K = 10^{50}$ erg and a half-opening angle $\theta_{jet} = 0.2$ radians (11°:5), and that they are surrounded by a relatively tenuous ISM with constant number density $n_{ISM} = 0.01 \text{ cm}^{-3}$. We fix the microphysical parameters<fnref n="10" /> so that the only remaining parameters needed to predict the afterglow light curve of the SGRB are the distance $d_L$ and the viewing angle $\theta_v$. We will link the viewing angle to the binary orbit inclination, and the distance will be obviously set equal to that of the binary. Assuming two opposite jets launched

perpendicular to the binary orbital plane, we have

$$ \theta_v(\iota) = \begin{cases} \iota & 0 \le \iota < \pi/2 \\ \pi - \iota & \pi/2 \le \iota < \pi \end{cases} \qquad (3) $$

where $\iota$ is the angle between the normal to the orbital plane and the line of sight.

For the dynamical ejecta macronova, we evaluate the disk mass and the ejecta velocity using the fitting formulas of Dietrich & Ujevic (2017), and we use them as inputs to compute the light curve following Grossman et al. (2013; using a constant gray opacity $\kappa = 10 \text{ cm}^2 \text{ g}^{-1}$), assuming a black-body spectrum with effective temperature equal to that of the photosphere. The input compact binary parameters in this case are the masses $M_1$ and $M_2$. To determine the compactness and the baryon mass of the neutron stars, which are necessary to associate the dynamical ejecta mass $M_{ej}$ and velocity $v_{ej}$ to the merger through the fitting formulas of Dietrich & Ujevic (2017), we assume the H4 EoS (Lackey et al. 2006; Glendenning & Moszkowski 1991), which has a mid-range stiffness<fnref n="11" /> among those that are compatible with the observational constraints (Özel & Freire 2016).

First, we need to derive the proper distributions of the distance and orbital plane inclination of the inspiral population detected by our interferometer network. For simplicity, we neglect the dependence of the network sensitivity on sky position and on the binary polarization angle $\psi$, and we assume that the maximum luminosity distance $d_{L,max}$ out to which an NS-NS inspiral can be detected depends only on the binary plane inclination $\iota$ with respect to the line of sight, namely,

$$ d_{L,\max}(\iota) = d_{L,\max}(0) \sqrt{\frac{1}{8}(1 + 6 \cos^2 \iota + \cos^4 \iota)}, \quad (4) $$

where $d_{L,max}(0)$ is the maximum luminosity distance out to which our network can detect a face-on inspiral. This expression accounts for the fact that gravitational radiation from a compact binary inspiral is anisotropic (Schutz 2011). Assuming that NS-NS mergers are uniformly distributed in space and have isotropic orientations, their distance and inclination distributions are $P(d_L) \propto d_L^2$ and $P(\iota) \propto \sin^4 \iota$. By the above assumptions, the probability that a binary with luminosity distance $d_L$ and inclination $\iota$ is detected is

$$ P(\det|d_L, \iota) = \begin{cases} 1 & \text{if } d_L < d_{L,\max}(\iota) \\ 0 & \text{otherwise.} \end{cases} \qquad (5) $$

By Bayes' theorem, the probability distribution of the distance and inclination of a detected NS-NS inspiral is then $P(d_L, \iota|\det) \propto P(\det|d_L, \iota) \times P(d_L) \times P(\iota)$, which gives

$$ P(d_L, \iota|\det) \propto \begin{cases} d_L^2 \sin \iota & \text{if } d_L < d_{L,\max}(\iota) \\ 0 & \text{otherwise.} \end{cases} \qquad (6) $$

The corresponding probability distribution of the inclination for detected inspirals (which is obtained by marginalization of Equation (6) over $d_L$) is then the well-known

$$ P(\iota|\det) = 7.6 \times 10^{-2}(1 + 6 \cos^2 \iota + \cos^4 \iota)^{3/2} \sin \iota. \qquad (7) $$

<footnote n="9">⁹ These are typical reference values for short GRBs, though they suffer from the still limited number of reliable measurements available.</footnote>

<footnote n="10">¹⁰ We refer here to the standard synchrotron afterglow model, and we set the microphysical parameters $p = 2.5$, $\epsilon_e = 0.1$, and $\epsilon_B = 0.01$. Such values, typical of long GRBs (e.g., Panaitescu & Kumar 2002; Ghisellini et al. 2009; Ghirlanda et al. 2015), seem to be representative for SGRBs as well (Fong et al. 2015), despite the much smaller sample of broadband light curves available.</footnote>

<footnote n="11">¹¹ This translates into mid-range values of the corresponding $M_{ej}$ and $v_{ej}$.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** "A priori detectability," i.e., the a priori probability that the EM counterpart of a compact binary inspiral is detected if the observation is performed at a time *t* after the merger, for observations in radio at 1.4 GHz, in infrared (IR) in the *J* band, and in optical in the *r* band, with limiting fluxes of 0.5 mJy in radio and 22.4 AB magnitude in IR and optical. The radio and optical probabilities account only for the jet afterglow, while the IR probability accounts only for the dynamical ejecta macronova. Based on a series of simplifying assumptions—see Section 2.1.1.</figcaption></figure>

For the distribution of the masses $M_1$ and $M_2$, we simply assume a normal distribution with mean $1.35 M_\odot$ and sigma $0.1 M_\odot$ for both of them, which reproduces the mass distribution of known galactic NS-NS binaries (Özel & Freire 2016).

Now, to construct the a priori detection probability $P(F(t) > F_{\text{lim}})$, we adopt the following Monte Carlo approach:

1. we construct our synthetic population of $N$ inspiral sampling distances and inclinations from $P(d_L, t|\det)$ and the masses $M_1$ and $M_2$ from the assumed normal distribution;

2. we compute the flux $F_i = F_A(t, d_{L,i}, \theta_v(\nu_i))$ of the jet afterglow or $F_i = F_M(t, d_{L,i}, M_{1,i}, M_{2,i})$ of the macronova in the chosen band for each sample;

3. we estimate $P(F(t) > F_{\text{lim}})$ as the fraction of $F_i$'s that exceed $F_{\text{lim}}$.

Figure 2 shows the a priori detectability computed with the above method for $d_{L,max}(0) = 100$ Mpc (which corresponds roughly to the sky-position-averaged aLIGO range for an optimally oriented NS-NS inspiral withO1 sensitivity; see Abbott et al. 2016c) for radio, infrared, and optical observations (see the caption for details). The sensitivity of the optical observations was chosen to match the limiting magnitude of the VST follow-up of GW 150914; see Abbott et al. (2016b). The flux of the jet afterglow has been computed using BOXFIT v. 1.0 (van Eerten et al. 2011).

A more accurate a priori detectability would require us to use astrophysically motivated priors on the other model parameters, such as the kinetic energy $E_K$, the ISM number density $n_{\text{ISM}}$, etc. Moreover, the actual intrinsic mass distribution of neutron stars that merge within the frequency band of GW detectors might differ significantly from the assumed one. The curves shown in Figure 2, thus, must be taken as illustrative.

The definitely higher detectability of the macronova is due to the fact that its emission is assumed to be isotropic, while the jet is fainter for off-axis observers.

## 2.2. Two Steps Further: How to Improve the Strategy Using Posterior Information on the Other Parameters of the Binary

### 2.2.1. The Full Posterior Probability Density in Parameter Space

Parameter estimation techniques applied to a compact binary coalescence signal $S$ result in a posterior probability density $P(\xi|S)$, where $\xi \in \mathbb{R}^n$ is a point in the $n$-dimensional parameter space. The “skymap probability,” i.e., the posterior sky-position probability density $P(\alpha|S)$, is essentially the $P(\xi|S)$ marginalized over all parameters but the sky position. Much more information is contained in the full posterior probability density, though, and some of it can be used to improve the design of the EM follow-up strategy.

### 2.2.2. Relevant Parameters in Our Case

In Section 2.1.1, we already made use of two extrinsic parameters of the compact binary inspiral that are relevant for the SGRB afterglow and the macronova, namely, the luminosity distance $d_L$ and the binary inclination$^{12}$ $\ell$.

Recent works based on numerical simulations of NS-NS and BH-NS mergers (e.g., Foucart 2012; Giacomazzo et al. 2013; Hotokezaka & Piran 2015; Dietrich & Ujevic 2017; Kawaguchi et al. 2016) seem to indicate that the amount of matter in the remnant disk and in the dynamical ejecta, plus some other properties of the latter such as the velocity profile, depend in a quite simple way (once an EoS is assumed) on the parameters of the binary prior to the merger, especially the masses $M_1$ and $M_2$ and the effective spin $\chi_{\text{eff}}$ of the black hole in the BH-NS case. Such information can be used to predict the observed light curve of the associated macronovas (e.g., Dietrich & Ujevic 2017; Kawaguchi et al. 2016) and, with greater uncertainty, the energy in the GRB jet (as in Giacomazzo et al. 2013).

Summarizing, at least the following compact binary coalescence parameters are relevant in order to predict the light curve of the SGRB and/or of the dynamical ejecta macronova associated with the merger:

1. the luminosity distance $d_L$ and the associated redshift $z$;

2. the orbital plane inclination $\ell$ with respect to the line of sight;

3. the component masses $M_1$ and $M_2$;

4. the effective spin $\chi_{\text{eff}}$ of the black hole in the BH-NS case.

Figure 3 represents a sketch of how the above parameters influence the properties of the SGRB afterglow and the dynamical ejecta macronova associated to the merger. The same approach can be adopted to link the properties of other EM counterparts (such as the long-lasting radio transient described by Nakar & Piran 2011 or the X-ray spin-down-powered transient described by Siegel & Cioffi 2016) to those of the binary, whose distributions can be constrained by the GW signal.

In Section 4, we will use the fitting formulas provided in Dietrich & Ujevic (2017) to compute the posterior ejecta mass distribution associated with the example NS-NS inspiral treated in that section. We will refrain from deriving the SGRB jet energy from the disk mass as suggested in Figure 3,

<footnote n="12">12 The inclination of the orbital plane with respect to the line of sight can also be relevant for the neutrino-driven wind (Martin et al. 2015) and the disk wind (Kasen et al. 2014) macronovas, due to their axial geometry and to the possibility that the dynamical ejecta act as a “lanthanide curtain” obscuring their optical emission (Rosswog et al. 2016) if the binary is observed edge-on.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Sketch of the dependence of some parameters of the SGRB and dynamical ejecta macronova on the progenitor binary parameters. The masses and spin of the binary components, through the dynamics and assuming an equation of state (EoS), determine the masses $M_{ej}$ (Kawaguchi et al. 2016; Dietrich & Ujevic 2017) and $M_{disk}$ (Foucart 2012; Giacomazzo et al. 2013) of the dynamical ejecta and the accretion disk on the remnant compact object, respectively. The dynamics also determine the velocity profile $dm/dv_{ej}$ of the ejecta. Accretion on the remnant converts disk rest mass (with some unknown efficiency) into jet kinetic energy $E_K$ (e.g., Giacomazzo et al. 2013), which constitutes the energy reservoir of the SGRB afterglow. By making assumptions on the remaining parameters, the light curves of the macronova and SGRB afterglow (possibly) associated with the merger can be predicted, taking into account the luminosity distance $d_L$ and the inclination $i$ of the binary.</figcaption></figure>

however, because that would require a detailed discussion about the proper disk mass energy conversion efficiency to be used, which is outside the scope of this work. We leave such a discussion to a future work.

### 2.2.3. A Posteriori Detectability

In Section 2.1, we introduced the idea of the “a priori detectability,” $P(F(t) > F_{lim})$, which can be regarded as the basic tool to set the timing of observations for the EM follow-up if no specific information about the source is available. Once the GW signal $S$ is observed, the information it carries can be used to construct a posterior probability $P(F(t) > F_{lim}|S)$ to better plan such observations. We will call it the “a posteriori detectability.” If the a priori detectability $P(F(t) > F_{lim})$ is constructed using the prior distributions of the parameters of the EM transient model, then the a posteriori detectability $P(F(t) > F_{lim}|S)$ is obtained exactly the same way (as exemplified in Section 2.1.1), but using the posterior distributions of the relevant parameters.

### 2.2.4. Detectability Maps

Several parameters of a compact binary inspiral are degenerate to some degree, i.e., the same signal $S$ can be produced by different combinations of the parameter values. These combinations, though, are not just uniformly distributed in some subset of the parameter space, but rather, they follow fundamental relations that depend both on the nature of the source (the binary inspiral) and on the properties of the detector network (the locations and orientations of the interferometers, their antenna patterns, the noise power spectrum). In particular, the distance, inclination, polarization angle, chirp mass, and sky position of the binary share a certain degree of degeneracy: the same signal $S$ can be produced by different combinations of values of these parameters and different realizations of the detector noises, which is the obvious reason why the sky-position uncertainty is so large. For this reason, if we restrict the posterior probability density in parameter space to a certain point of the skymap, i.e.,

we take the sky-position-conditional posterior distribution of the physical parameters of the binary, in principle it will depend on the chosen sky position. Knowing the sky-position-conditional posterior probability distribution $P(d_L, i, M_1, M_2, ..., \alpha, S)$ of the relevant binary parameters at sky position $\alpha$, we can thus derive the corresponding distribution of the properties of the EM counterpart at that particular sky position, which means that we can construct a sky-position-conditional posterior detectability $P(F(t) > F_{lim}|\alpha, S)$ that can be used as the basis of the EM follow-up strategy. We call this quantity the “detectability map.”

## 2.3. Recap

It is useful to summarize here the steps of increasing complexity that led us to the definition of the detectability maps:

1. We started by assuming an unrealistic model of the EM counterpart: a source that turns on at the GW time and never turns off. In this case, no timing information is needed for the follow-up strategy, which simply consists of scanning the localization uncertainty area, starting from the most probable sky location, until the source is found;

2. If a model of the counterpart is available and prior distributions of the model parameters can be assumed (based on available astrophysical data or on an educated guess), the “a priori detectability” $P(F(t) > F_{lim})$ can be constructed, as shown in Section 2.1.1. This is the best follow-up timing information that can be constructed based on a priori knowledge only;

3. After a GW signal $S$ is detected and parameter estimation has been performed, prior distributions of the model parameters can be (partly) replaced with posterior distributions derived from the signal: the “a posteriori detectability” $P(F(t) > F_{lim}|S)$ can be constructed. This exploits information contained in the GW signal, but it is still independent of the sky position;
4. If the counterpart is assumed to be located at a certain sky position, the corresponding sky-position-conditional posterior distributions can be used in place of the full posterior distributions. Indeed, given a signal S, compact binary inspiral parameters compatible with S and a particular sky position are in general different from those compatible with S and another sky position. By varying the assumed sky position on a grid that covers the whole skymap, one can then construct the “detectability map” P(F(t) > F<sub>lim</sub>|α, S).

Let us now introduce a method to compute the detectability maps and apply it to a practical example.

### 3. How to Construct and How to Use the Detectability maps

#### 3.1. Extraction of the Sky-position-conditional Posterior Distributions Using a Simple Method Based on “Inverse-Distance Weighting”

The extraction of the sky-position-conditional posterior distribution requires some multidimensional kernel density estimation (KDE) technique to be applied to the posterior samples obtained from a parameter estimation pipeline run on the GW signals recorded by the detectors. Since the aim of this work is to propose a new approach in the design of the EM follow-up rather than to discuss the technical subtleties of such multidimensional KDE, we adopt the following simple and intuitive method, which can be replaced with a more accurate one in a possible application of our approach to a real case.

Our simplified method to extract the sky-position-conditional posterior distribution P(q|α, S) of a quantity q at sky position α is based on the concept of “inverse-distance weighting” (Shepard 1968): we assume that each posterior sample {αᵢ, qᵢ, d<sub>L,i</sub>, vᵢ, ...} contributes to the P(q|α, S) with a weight that is a decreasing function of the angular distance δ(α, αᵢ) between the posterior sample and the sky position α. In particular, we assign a Gaussian weight to each posterior sample:

$$w_i \propto \exp\left[-\frac{1}{2}\left(\frac{\delta(\alpha, \alpha_i)}{\sigma(\alpha)}\right)^2\right], \quad (8)$$

where the bandwidth σ(α) is taken as

$$\sigma(\alpha) = \sqrt{\frac{\sum_{i=1}^{N} [\delta(\alpha, \alpha_i) - \langle\delta\rangle]^2}{N-1}} \times N^{-1/5}, \quad (9)$$

where ⟨δ⟩ is the arithmetic mean of the δ(α, αᵢ). The normalization of the weights is given by $\sum_{i=1}^{N} w_i = 1$.

The ideas behind this method are simply that the closer the posterior sample is to the sky position α, the more it contributes to the conditional posterior distribution at that sky position, and that the influence of the posterior sample decreases as a Gaussian with increasing angular distance. The choice of the bandwidth (Equation (9)) is just “Silverman’s rule of thumb” (Silverman 1982) for a Gaussian KDE in a one-dimensional parameter space (namely, the angular distance space).

The mean of q at sky position α is thus computed as

$$\langle q \rangle_{\alpha} = \sum_{i=1}^{N} w_i q_i, \quad (10)$$

and similarly, the variance

$$\mathrm{Var}_{\alpha}(q) = \sum_{i=1}^{N} w_i q_i^2 - \left( \sum_{i=1}^{N} w_i q_i \right)^2. \quad (11)$$

More generally, the sky-position-conditional posterior distribution of q at sky position α is approximated as

$$P(q|\alpha, S) \sim \sum_{i=1}^{N} w_i K\left(\frac{q - q_i}{\sigma_q}\right), \quad (12)$$

where K(x) is some kernel function and σq is its bandwidth.

We performed tests to show that the above method yields consistent results (see the Appendix). As one might expect, the results are accurate in sky regions where the distribution of posterior samples is sufficiently dense.

#### 3.2. The Detectability Map

By the above method, we can thus define the sky-position-conditional posterior detectability estimate (i.e., the detectability map) as

$$P(F(t) > F_{\mathrm{lim}}|\alpha, S) = \sum_{i=1}^{N} w_i \int_{F_{\mathrm{lim}}}^{\infty} K\left(\frac{F - F_i(t)}{\sigma_F}\right) dF, \quad (13)$$

where Fᵢ(t) represents the flux (in the chosen band) at time t of the light curve computed using the ith posterior sample parameter values, Fᵢ(t) = F(t, d<sub>L,i</sub>, vᵢ, ...). If we approximate the kernel functions with delta functions K(x) ~ δ(x), the expression becomes

$$P(F(t) > F_{\mathrm{lim}}|\alpha, S) \sim \sum_{i=1}^{N} w_i H(F_i(t) - F_{\mathrm{lim}}), \quad (14)$$

where H(x) is the Heaviside function.

#### 3.3. The Earliest, Best, and Latest Detection Time Maps

By the above method, the information encoded in the GW signal is used to estimate the detectability of the EM counterpart at a given time t, if it is at a certain sky position α. By setting a minimum required detectability λ, one can define (for each sky position α) a time interval during which P(F(t) > F<sub>lim</sub>|α, S) ≥ λ. If the detectability map P(F(t) > F<sub>lim</sub>|α, S) never reaches λ for a certain sky position, that position is “hopeless,” i.e., there is too little chance of detecting the EM counterpart if it is located there. The earliest and latest detection time maps are then defined respectively as

$$\begin{cases} t_{E,\lambda}(\alpha) = \inf \{t | p(t, \alpha) \ge \lambda\} \\ t_{L,\lambda}(\alpha) = \sup \{t | p(t, \alpha) \ge \lambda\} \end{cases} \quad (15)$$

where we set p(t, α) = P(F(t) > F<sub>lim</sub>|α, S) for ease of reading. Irrespective of λ, the best detection time map can be defined as

$$t_B(\alpha) = \operatorname{argmax} \{p(t, \alpha)\}. \quad (16)$$

These maps are the simplest piece of information about “where and when to observe” that can be constructed using the detectability map. A follow-up strategy should then try to arrange observations so that a field centered at α is observed at a time as close as possible to t<sub>B</sub>(α), and in any case not earlier than t<sub>E,λ</sub>(α) or later than t<sub>L,λ</sub>(α). Let us now introduce a simple
algorithm to construct such a follow-up strategy, after which
we will be able to show a practical example (Section 4).

3.4. *A Follow-up Strategy Construction Algorithm*

In order to perform a first test of the approach outlined in the
preceding sections, we will apply it to a simulated event, and
we will construct a “simulated follow-up strategy” based on it.
To this end, we use an unambiguous algorithm to define the
strategy for a given event and a given observing facility. To
keep the discussion as simple as possible, we work in an
idealized setting where all points of the skymap are observable
by our facility during some pre-defined time windows. We
assume that each observation covers an area A<sub>FoV</sub> at observing
frequency ν<sub>obs</sub> and that the limiting flux F<sub>lim</sub> for detection is
independent of sky position and is always reached after an
integration time T<sub>int</sub>. The outline of the algorithm is the
following:

1. we divide the skymap into patches, each representing a potential field to be observed;

2. we define a list of available time slots (i.e., possible observing time windows) on the time axis, starting 24 hr after the event (posterior samples are typically obtained after several hours or even days, so an earlier start would be unrealistic. Note that this does not mean that we discourage an earlier follow-up—which is of great importance especially for the optical and X-ray afterglows—but only that this method might not be applicable with very low latency<fnref n="13" />);

3. starting from the patch with the highest sky-position probability, we check if the detectability map $P(F(t) > F_{\text{lim}}|\alpha, S)$ within that patch exceeds $\lambda$ at some time within the available time slots: if it does not, we discard the patch (we choose not to observe it); if it does, we schedule the observations that cover that patch at the time when the detectability is highest, and we mark the corresponding time slots as no longer available;

4. we proceed to the next patch in descending order of skymap probability until the available time is over, or until all patches have been processed.

To keep the implementation of the above steps as simple as
possible, we use a HEALPix tessellation of the sky (Gorski
et al. 2005) to define the observable fields: it is a way to divide
the sky (i.e., a sphere) into equal area patches, called “pixels.”
The “order” N<sub>side</sub> of the HEALPix tessellation defines the
number of pixels the sky is divided into, namely,
N<sub>pixels</sub> = 12N<sub>side</sub><sup>2</sup>. The pixel area is then A<sub>pixel</sub> = 3438N<sub>side</sub><sup>-2</sup> deg<sup>2</sup>. Thus, we replace the actual observations of duration T<sub>int</sub>
and field of view A<sub>FoV</sub> with “pseudo observations” of area
A<sub>pixel</sub> and effective duration T = T<sub>int</sub> × A<sub>pixel</sub>/A<sub>FoV</sub>, choosing
N<sub>side</sub> in order to minimize the difference between A<sub>pixel</sub> and
A<sub>FoV</sub>. The algorithm is thus implemented as follows:

1. We consider the posterior samples produced by a parameter estimation sampler (multiple sampling algorithms are implemented in LIGO's LALINFERENCE-parameter estimation tool<fnref n="14" />) applied to a simulated signal

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Schematic representation of the strategy construction algorithm. At a given step, some of the available time slots have been already assigned, while others are still available. Observations that cover the current pixel are assigned to the available time slot corresponding to the maximum detectability.</figcaption></figure>

$S$. These are points in the compact binary inspiral
parameter space distributed according to the posterior
probability density;

2. We find the 90% sky-position confidence region of the source based on the posterior samples;

3. We divide the sky into pixels according to a HEALPix tessellation of order Nside, and we consider only those pixels that fall inside the 90% sky-position confidence region of the simulated signal;

4. We associate to each pixel p the integral of the skymap probability density $P_p = \int_{A_{\text{pixel}}(p)} P(\alpha|S) d\Omega$ over the pixel area;

5. We associate to each pixel p the earliest and latest detection times $t_{E,\lambda}(p)$ and $t_{L,\lambda}(p)$ averaged over that pixel; pixels for which the detectability never reaches $\lambda$ are excluded from the list of possible observations;

6. We divide the time axis into contiguous intervals (slots) of duration *T* (the effective time needed to cover one pixel), and we mark some of the slots as “available” for the follow-up (based on the characteristics of the instrument);

7. We sort the pixels in order of decreasing $P_p$, and starting from the first ($p=1$) pixel, we do the following:

(a) we check that at least one available time slot is comprised between $t_{E,\lambda}(p)$ and $t_{L,\lambda}(p)$: if not, no observation of the pixel is scheduled; otherwise, the available time slot where the detectability is maximum is assigned to the observation of the pixel;

(b) we proceed to the next pixel, until all available time slots are assigned, or until all pixels have been processed.

Figure 4 shows a schematic representation of the algorithm
described above.

The output of the above algorithm is thus a list of
observation times $t_{\text{obs}}(p)$ that cover (part of) the skymap,
giving priority to pixels with high skymap probability and
trying to use the available time in a way that maximizes the
probability of detecting the transient. The inputs of the
algorithm are:

1. the posterior sample list based on $S$;

2. the detectability threshold $\lambda$;

3. the available time windows;

<footnote n="13">¹³ In the example of Section 4, the only parameters actually needed to predict the possible SGRB afterglow light curves are $d_L$ and $\nu$. Posterior distributions of these parameters (so-called “extrinsic”) can be obtained with a low latency analysis tool such as BAYESTAR (Singer & Price 2016).</footnote>

<footnote n="14">¹⁴ http://software.ligo.org/docs/lalsuite/lalinference/</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Light curves of the EM counterparts associated with injection 28840 of F2Y. Jet afterglow: radio (1.4 GHz—blue line) and optical (r filter—red line); macronova: infrared (J filter—yellow line). The same assumptions as in Section 2.1.1 were adopted. The limiting fluxes for detection adopted in the test example (Section 4) are shown with horizontal black dashed lines.</figcaption></figure>

4. the instrument observing frequency $\nu_{\text{obs}}$, the field-of-view area $A_{\text{FoV}}$, the limiting flux $F_{\text{lim}}$, and the corresponding integration time $T_{\text{int}}$ (which should also include the slew time);

5. The HEALPix tessellation order $N_{\text{side}}$, which should give a pixel area $A_{\text{pixel}}$ close to $A_{\text{FoV}}$.

The observations are given in order of decreasing “importance” (skymap probability). Once the list of observations is produced, one can decide to perform only the first *N* observations that fit into the available telescope time. If all observations suggested by the algorithm can be performed within the available time, the excess time can be used, e.g., to take comparison images of some fields at different times (for the identification of transients or uncatalogued variable sources, in the absence of previously available images) or to perform deeper observations for the characterization of the candidates.

We are now ready to construct an example of how the use of sky-position-conditional posterior probabilities can help in the definition of an EM follow-up strategy.

## 4. Test Example: Injection 28840, an NS-NS Merger with an Associated Orphan Afterglow

As a test example, we consider injection number 28840 from the “First Two Years of Electromagnetic Follow-up with Advanced LIGO and Virgo” study (Singer et al. 2014, F2Y hereafter). The injection simulates the inspiral of a neutron star binary with $M_1 = 1.59 M_\odot$ and $M_2 = 1.53 M_\odot$ at a luminosity distance $d_L = 75$ Mpc, with orbital plane inclination $i = 14^\circ$, at sky position (R.A., decl.) $= (23^\circ 27' 12''$, $-10^\circ 30' 0''$), detected by the two-detector Advanced LIGO network on MJD 55483.27839 (i.e., at 06:40:53 of 2010 October 14—this is just a simulated event) adopting an early sensitivity curve corresponding to a binary NS range of 55 Mpc (Barsotti & Fritschel 2012). Assuming that a relativistic jet with isotropic equivalent kinetic energy $E_K = 10^{50}$ erg and half-opening angle $\theta_{\text{jet}} = 0.2$ rad (which is less than the viewing angle, thus the afterglow is an orphan) is launched perpendicular to the orbital plane right after the merger, we computed its afterglow light curve assuming an ISM number density $n_{\text{ISM}} = 0.01 \text{ cm}^{-3}$, adopting the same microphysical parameters as in Section 2.1.1, using BOXFIT v. 1.0 (van Eerten

et al. 2011). The light curves at $\nu_{\text{obs}} = 1.4$ GHz and in the r filter are shown in Figure 5.

To produce our sky-position-conditional posterior distributions, we use 7962 posterior samples produced by one of the LALINFERENCE parameter estimation samplers. Since we fixed the jet isotropic equivalent kinetic energy $E_K = 10^{50}$ erg, in this example the only binary parameters that are relevant to the posterior detectability of the jet afterglow are the luminosity distance and the binary plane inclination. Figure 6 shows the best detection time map at $\nu_{\text{obs}} = 1.4$ GHz produced using these samples (see the caption for additional information). The 90% sky-position confidence area (represented by the larger red contours) covers approximately 1500 deg$^2$.

### 4.1. Optical Search

For our virtual optical EM follow-up, we adopt parameters inspired by the VST follow-up of GW 150914 (Abbott et al. 2016b). We consider observations in the r band ($\nu_{\text{obs}} = 4.8 \times 10^{14}$ Hz) with a detection limit $F_{\text{lim}} = 22.4$ AB magnitude, reached after $T_{\text{int}} = 100$ seconds (slew + integration). The field of view is $A_{\text{FoV}} = 1$ deg$^2$. With the adopted flux limit, the optical light curve (Figure 5) becomes too faint to be detected after a few hours. Consistently, the detectability $P(F(t) > F_{\text{lim}}|\alpha, S)$ is below 1% at all times $t > 1$ day over the whole skymap except for a subregion of the skymap with a total area of 185 deg$^2$. We run the follow-up strategy construction algorithm allowing 3 hr of available time per night from day 1 to day 15. We set $N_{\text{side}} = 64$, which gives $A_{\text{pixel}} = 0.84$ deg$^2$ and $T = 84$ s. Even adopting the very low detectability limit $\lambda = 0.01$ (i.e., allowing for observations with a detectability as low as 1%), only 24 “pseudo observations” (corresponding to 22 observations) are scheduled (all during the first available night), totalling 36 minutes of telescope time. The position of the EM counterpart is not contained in any of the observed fields.

The short integration time and the relatively shallow detection limit of the VST follow-up of GW 150914 are good in order to cover the largest possible area during the next few nights after the event; if we wish to have some chance of detecting a relatively dim optical afterglow like that in Figure 5, though, we need to go deeper, i.e., we need longer integration times. We thus repeat the optical search with the same parameters as above, but with $T_{\text{int}} = 1000$ s and $F_{\text{lim}} = 24.5$ m$_{AB}$. We also raise the minimum detectability to $\lambda = 0.05$, to avoid pointing fields with a very low detectability. The algorithm outputs 54 “pseudo observations,” corresponding to 46 observations, totalling 15 hr of telescope time. The $p = 25$ observation (i.e., the 25th observation in descending order of $P_p$), scheduled $2^d 2^h 34^m 45^s$ after the event, contains the EM counterpart. The flux at that time is $24.58$ m$_{AB}$, slightly dimmer than the required detection threshold, but detectable with a threshold S/N of 5 under optimal observing conditions.

As explained in Section 3.4, the follow-up observations suggested by the algorithm are given in descending order of sky-position probability, thus astronomers can choose to perform only the first *N* observations if not enough telescope time is available to complete them all. The field containing the counterpart is at $p = 25$ in this case. To complete the first 25 observations, 7 hr of telescope time are needed: this represents the minimum amount of telescope time for the EM counterpart to be detected by this facility in this case.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Mercator projection of the best detection time map $t_B(\alpha)$ at $\nu_{\text{obs}} = 1.4$ GHz for injection 28840 of F2Y. Red lines represent the contours of the 50% and 90% sky-position confidence areas; the black solid line is the locus of sky positions that yield a GW signal arrival time delay of $\Delta t = 0.359$ ms (the “true” arrival time difference associated with the injection position) between the two LIGO detectors, and the black dashed lines correspond to $\Delta t \pm 1$ ms. The white contours are the isolines of the average antenna pattern of the aLIGO network (the four smaller closed circular contours enclose the local and absolute minima, while the two large closed contours near the lower-left and the upper-right corners enclose the absolute maxima). The injection true position is marked with a star symbol. Regions marked by blue contours are those where the posterior detectability of the EM counterpart of our test example reaches the required 5% threshold.</figcaption></figure>

Figure 7 shows the positions and times of the observations scheduled by the algorithm in this case. During the first night, essentially all points of the detectability map are above the limit $\lambda$, thus observations are concentrated around the centers of the two large uncertainty regions (see Figure 6 for an all-sky view), which are the points of largest skymap probability. During the second night, the detectability has fallen below the limit in most of the central parts of the two uncertainty regions, thus the algorithm moves toward points of lower skymap probability, but higher detectability. The evolution of the detectability map proceeds in a similar fashion until the fifth night, when the detectability at all points of the skymap eventually falls below $\lambda$.

## 4.2. Radio Search

The parameters of our virtual radio follow-up are inspired by MeerKAT, the South African SKA precursor. Sixteen (of the eventual 64) 13.5 m dishes have already been integrated into a working radio telescope and produced their “first light” image<fnref n="15" /> in 2016 July. We assume a field of view $A_{\text{Fov}} = 1.7$ deg² at $\nu_{\text{obs}} = 1.4$ GHz. We conservatively estimate a 50 µJy rms noise for a $T_{\text{int}} = 1000$ s observation (slew + integration), assuming 16 working dishes. In a large-area survey, a $10\sigma$ detection is usually required to avoid a large number of false

alarms, thus we set $F_{\text{lim}} = 0.5$ mJy, i.e., we require the flux to be 10 times the rms noise for the detection to be considered confident. We allow a maximum of 20% of the available time from day 1 to day 100 to be dedicated to the follow-up (in practice, we allow an available time window of 4.8 hr each day), and we adopt a $\lambda = 0.05$ detectability limit. The best detection time map with the chosen parameters is shown in Figure 6, where the blue contours represent regions where the detectability reaches the required detectability limit at some time $t > 1$ day (382 deg² in total). Setting $N_{\text{side}} = 64$, the follow-up construction algorithm outputs 337 “pseudo observations” (corresponding to 169 pointings, totalling 47 hr of telescope time), which are represented in Figure 8. The $p = 35$ observation contains the counterpart. It is scheduled $5^\text{d}0^\text{hr}9^\text{m}27^\text{s}$ after the event, when the flux of the EM counterpart is 0.52 mJy (see Figure 5), which means that the afterglow is detected at better than $10\sigma$. The first 35 (pseudo) observations make up 4.7 hr of telescope time: only this amount needs to be actually allocated for the follow-up to successfully detect the radio counterpart.

## 4.3. Infrared Search of the Associated Macronova

For the same event, we computed the infrared (J band) light curve (see Figure 5) of the associated macronova with the same assumptions as in Section 2.1.1. Due to the rather large masses of the binary components and to their similar mass, the

<footnote n="15">¹⁵ Media release at http://www.ska.ac.za/media-releases/.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Time and position of the optical follow-up observations of our test example (Section 4). The star marks the injection position.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Time and position of the radio follow-up observations of our test example (Section 4). The star marks the injection position.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Time and position of the infrared follow-up observations of our test example (Section 4). The star marks the injection position.</figcaption></figure>

dynamically ejected mass is small ($M_{ej} \approx 5.6 \times 10^{-3} M_{\odot}$ according to the Dietrich & Ujevic 2017 fitting formula assuming the H4 EoS). The light curve is thus quite dim. It peaks between the second and the third day, slightly brighter than 22.4 $m_{AB}$ in the J band. The effective temperature of the photosphere at peak is $T_{peak} \sim 2900$ K. To detect such a transient with a telescope like VISTA, an integration time of the order of 1000 s is needed. We thus perform our virtual follow-up strategy with the following parameters inspired by VISTA: we choose a limiting flux $F_{lim} = 22.4 m_{AB}$ in the J band with $T_{int} = 1200$ s, and we set $A_{Fov} = 1.5 \text{ deg}^2$. Again, we assume that 3 hr per night are dedicated to the follow-up. The a priori detectability (see Figure 2) of the macronova for this limiting flux is high, meaning that most of the possible light curves exceed 22.4 $m_{AB}$, thus we set $\lambda = 0.5$ to limit the search to points of the skymap that reach a detectability at least as good as the a priori one.

Adopting the above parameters, the algorithm outputs 123 “pseudo observations” (corresponding to 69 pointings, totalling 23 hr of telescope time), which are represented in Figure 9. The 79th observation in descending order of $P_p$ contains the counterpart. It is scheduled $6\text{d}0\text{h}18\text{m}27\text{s}$ after the event, when the flux of the macronova in the J band is just below the 22.9 $m_{AB}$ and the effective temperature of the spectrum is around 2600 K. The flux is lower than $F_{lim}$ (which was chosen to represent an indicative limit for achieving a S/N ~ 10), but according to ESO’s exposure time calculator<fnref n="16" />, such emission would be detected through the VISTA telescope at ESO in the J band with an S/N of 6, assuming a 1000 s integration, in

<footnote n="16">¹⁶ We queried the VIRCAM ETC at http://www.eso.org assuming a 1.2 airmass and a seeing of 0.8 arcsec.</footnote>
<table><caption>**Table 1**

Algorithm Input Parameters Used to Construct the Example Follow-up Strategies</caption><thead><tr><th>Instrument</th><th>ν<sub>obs</sub> (Hz)</th><th>F<sub>lim</sub> (μJy)</th><th>A<sub>Fov</sub> (deg<sup>2</sup>)</th><th>T<sub>int</sub> (s)</th><th>λ</th><th>N<sub>side</sub></th><th>Available Time<sup>a</sup></th><th>Det. Time<sup>b</sup> (hr)</th></tr></thead><tbody><tr><td>VST-like (shallow)</td><td>4.8 × 10<sup>14</sup> (r filter)</td><td>4 (22.4 m<sub>AB</sub>)</td><td>1</td><td>100</td><td>0.01</td><td>64</td><td>3 hr/night</td><td>...</td></tr><tr><td>VST-like (deep)</td><td>4.8 × 10<sup>14</sup> (r filter)</td><td>0.58 (24.5 m<sub>AB</sub>)</td><td>1</td><td>1000</td><td>0.05</td><td>64</td><td>3 hr/night</td><td>7</td></tr><tr><td>MeerKAT-like</td><td>1.4 × 10<sup>9</sup></td><td>500</td><td>1.7</td><td>1000</td><td>0.05</td><td>64</td><td>20%</td><td>4.7</td></tr><tr><td>VISTA-like</td><td>2.4 × 10<sup>14</sup> (J filter)</td><td>4 (22.4 m<sub>AB</sub>)</td><td>1.5</td><td>1200</td><td>0.5</td><td>64</td><td>3 hr/night</td><td>14.5</td></tr></tbody></table>

**Notes.**

a Available observing time starts 24 hr after the event.

b Minimum observation time needed for detection.

optimal observing conditions. Our virtual infrared follow-up would thus again result in a detection in a search with a threshold at S/N ~ 5. To achieve it, 14.5 hr of telescope time are needed (the amount of time for the first 79 most important observations to be performed).

### 4.4. Comparison with Follow-up Strategies Based on the A Priori Detectability Only

For comparison, we performed additional optical, radio, and infrared searches using the same parameters as before (listed in Table 1), but replacing the a posteriori detectability $P(F(t) > F_{lim}|\alpha, S)$ with the a priori detectability $P(F(t) > F_{lim})$ computed in Section 2.1.1. This should simulate a search based on a priori information only. The results are the following:

1. Optical search: the counterpart is in the field of view $5^{\circ}22^{\circ}53^{\circ}32^{\circ}$ after the event, when the flux of the jet afterglow in the r band is as low as 26 $m_{AB}$, which is definitely too faint for a detection;

2. Radio search: the counterpart is in the field of view $20^{\circ}22^{\circ}18^{\circ}0^{\circ}$ after the event, when the flux of the jet afterglow at 1.4 GHz is 120 μJy, which is significantly below our required limiting flux. Even assuming that the sensitivity was good enough for a detection, the facility should have allocated at least 77.6 hr of telescope time for this single follow-up in order to include the observation that contains the counterpart;

3. Infrared search: none of the 122 fields of view whose observation is scheduled by the algorithm contains the counterpart.

We conclude that the use of posterior information from the GW signal has a decisive impact on the EM follow-up in our example test case.

## 5. Discussion

With this work, we proposed the new idea that information on the compact binary inspiral parameters extracted from a GW signal can be used to predict (to some extent) the best timing for observation of the possible EM counterpart. In practice, the probability distributions of the binary parameters inferred from the GW signal are fed to a model of the candidate EM counterpart in order to define a family of possible light curves. The possible light curves are then used to construct the “detectability maps,” which represent an estimate of how likely is the detection of the EM counterpart with a given instrument, if the observation is performed at time $t$ looking at sky position $\alpha$ = (R.A., decl.). In order to apply the idea to a practical example, we introduced an explicit method to construct the

detectability maps (Section 3.1) and an algorithm that uses these maps to define an EM follow-up observing strategy (Section 3.4). We then applied the method to a synthetic example, showing that it significantly improves the effectiveness of the EM follow-up (Section 4).

In order to keep the treatment as simple as possible, we adopted many simplifications at various stages of the discussion, and we intentionally avoided mentioning some secondary details or including too much complexity. Let us briefly address some of the points that were not discussed in the preceding sections.

### 5.1. Model Dependence and Inclusion of Priors on Unknown Parameters

The approach clearly relies on the availability of models of the EM counterparts, and on our confidence in the predictions of these models. On the other hand, since the light curves are treated in a statistical sense, the models only need to represent correctly the peak flux and light-curve general evolution. Fine details are lost in the processing of the light curves, and are thus unnecessary. Moreover, it is very straightforward to include our uncertainty on the, model parameters unrelated to the GW signal. In all examples discussed in this work we fixed the values of such parameters (e.g., the kinetic energy $E_K$ and ISM density $n_{ISM}$ of the SGRB afterglow, or the NS EoS). A better approach (at the cost of a higher computational cost) would be to assume priors for these parameters, i.e., to assign a probability distribution to the values of these parameters based on some prior information (e.g., available astrophysical data, if any) or on theoretical arguments. In this case, multiple light curves of the counterpart must be computed for each posterior sample, using different values of the unknown parameters sampled from the assumed priors. This should be the most effective way of incorporating the uncertainty on these parameters in the computation of the detectability maps (it applies as well to the a priori detectability and to the a posteriori detectability). We will explore the effect of the inclusion of such priors in the construction of a priori detectabilities, a posteriori detectabilities, and detectability maps for a range of potential EM counterpart models in a future work.

### 5.2. Sky-position Dependence of the Parameters

In Section 2.2.4, we stated that, in general, the sky-position-conditional posterior distribution of the inspiral parameters depends on the assumed sky position. The main driver of this dependence is the sky projection of the antenna patterns (i.e., the sensitivity to different polarizations) of the interferometers of the network: the distribution of distances, inclinations, and
mass ratios compatible with a given signal, assuming that the source is at a particular sky position, is especially constrained by what the interferometers can or cannot detect if the source is at that sky position. As an intuitive example, say that a particular sky position corresponds to the maximum sensitivity of one of the detectors with respect to a particular polarization, and say that the incident GW that yields the signal picked up by that detector contains no such polarization, then only combinations of parameters for which the corresponding component of the strain is smaller than the limit set by the sensitivity are admissible if that sky location is assumed. If another sky location is assumed, the constraints change accordingly. As a general trend, points of the sky where the network has a higher sensitivity will correspond to a larger average distance of the source, thus implying a lower detectability for both the macronova and the SGRB afterglow. If the sky-position uncertainty region is smaller than the typical angular scale over which the antenna pattern varies, the dependence of the posterior distributions of the parameters on sky position becomes less important: the a posteriori detectability contains most of the relevant information in that case. With Advanced Virgo joining the network, in many cases the sky-position uncertainty region will still extend over a few hundreds of square degrees (Abbott et al. 2016c); on the other hand, better information on the two polarization states of the GW signal will be available (the two interferometers of the aLIGO network are almost anti-aligned, and thus the ability of the network to distinguish between the two polarization states is rather poor as of now). In the next decades, third-generation interferometers will again face the same issues with sky localization. The use of detectability maps instead of the a posteriori detectability alone is thus likely to remain useful with more advanced networks as well.

### 5.3. The Choice of Injection 28840

The injection event used to construct the example presented in the last section was selected among those in the F2Y study. We considered a two-detector case (LIGO only), as it leads in general to larger localization uncertainties. We looked through the list for an event that was quite distant and whose orbit inclination was sufficiently inclined for the jet to be off-axis. The 28840 injection event luminosity distance is indeed rather large ($d_L = 75$ Mpc), the jet is slightly off-axis ($\iota = 14^\circ$), the sky-position uncertainty is large (more than $1500 \text{ deg}^2$), and the injection position is rather far away from the maximum of the skymap probability. The latter condition makes a search based only on the a priori information particularly ineffective, because the exploration of the skymap proceeds slowly (using relatively small field instruments) from the center of the skymap (where the skymap probability is high) to the periphery (where the source is actually located). In cases when the sky location is better reconstructed (i.e., the source is closer to the point of maximum skymap probability), the improvement in the EM follow-up effectiveness thanks to the detectability maps (with respect to a strategy based only on the skymap probability and on a priori information on the EM counterpart characteristics) could be less striking. We plan to study systematically the relative improvement in a future work.

### 5.4. A Better Strategy Construction Algorithm

The algorithm (Section 3.4) used to construct the virtual EM follow-up strategy of our example is admittedly oversimplified. An algorithm suited for real application should be able to:

1. take into account the actual observability constraints on all points of the skymap at a given time, e.g., the setting of tiles below the horizon;

2. use a better (instrument specific) way of dividing the sky into potentially observable fields;

3. consider the impact of airmass, expected seeing, dust extinction, stellar density in the field, and other variables on the detectability;

4. potentially use a different integration time for each tile, in order to maximize the detection probability;

5. avoid sequences of widely separated pointings, which would result in a waste of time in slewing.

Some recent works have already addressed, at least in part, some of the above points. Ghosh et al. (2015) discussed an algorithm for the optimization of the tiling, which is totally compatible with our approach, since the detectability maps do not set a preferred tiling. Rana et al. (2017) developed and compared some ingenious algorithms that aim to maximize the sky-position probability in the search, taking into account per-tile setting and rising times. Their approach does not account for the time evolution of the EM counterpart luminosity, though, and it can result in the paradoxical situation in which the highest probability tile is observed *before* the time at which the flux is high enough for a detection. Incorporating the information from the detectability maps in their method could be the starting point for a realistic automated strategy construction algorithm based on the ideas presented in this work. Concerning point (4) above, both Coughlin & Stubbs (2016) and Chan et al. (2017) found that an equal integration time in all observations could be sub-optimal with respect to the EM counterpart detection probability. Their assumptions, though, are significantly different from ours: both assume a constant luminosity of the EM counterpart (i.e., they ignore the time variation of the flux), and they do not use astrophysically motivated priors on luminosity (the former use a flat prior, while the latter use Jeffrey's prior $L^{-1/2}$ limited to the range spanned by the peak luminosities of the macronova model in Barnes & Kasen 2013). It is not straightforward to figure out if their results are applicable to our more general case as well.

As an additional caveat, we only considered the case where one epoch of observation per field is enough to identify transient sources. Realistically, this is only feasible when previous images can be used as reference, or where source catalogs are complete up to the survey limit. The identification of interesting transients and the removal of those unrelated to the GW source require at least two epochs of observations, which are not taken into account in our example algorithm.

### 5.5. Use in Conjunction with the “Galaxy Targeting” Approach

The use of detectability maps is entirely compatible with a search based on targeting candidate host galaxies. The observation of each target galaxy would simply need to be performed as close to the corresponding best detection time (as defined in Section 3.3) as possible. Since choosing a target galaxy corresponds to assuming a known distance to the
source, an intriguing further refinement of the present method could be to consider the posterior distribution of binary parameters conditioned on both sky position and distance. This would have a great impact especially on the binary orbit inclination. In other words, it would be possible to associate a fairly well-defined binary orbit inclination to each galaxy. The consequence would be that some galaxies (typically the most distant ones) would be better candidate hosts for an SGRB afterglow with respect to others, depending on the associated binary orbit inclination.

## 5.6. Computational Feasibility of the Approach

The Monte Carlo approach adopted in this work, in which “all possible” light curves of an event must be computed (at least one per posterior sample, which means around $10^4$ light curves per observing band—which must be increased by one or two orders of magnitude if priors on unknown parameters are included) requires a computationally effective way to produce the light curves. Since the aim of the approach is to assess the detectability, rather than to fit the model to observational data, simple analytic models that capture the main features of the expected light curves are better than complex numerical models suited for parameter estimation. The macronova model by Grossman et al. (2013) is a good example: a few minutes are sufficient to compute $10^4$ light curves on a laptop computer using this model. Parallelization is instead unavoidable<fnref n="17" /> when such a large number of off-axis SGRB afterglow light curves are to be produced using BOXFIT (van Eerten et al. 2011).

## 5.7. Reverse Engineering: Tuning the Models Using Information from the EM Counterparts

A fascinating possible future application could be to use the detectability maps to test the underlying EM counterpart models and the assumptions on the priors: when a relatively large number of inspirals involving at least one neutron star has been detected, it will be possible to use the detectability maps to estimate how likely the detection (or non-detection) of the corresponding EM counterparts would have been with a particular choice of priors and adopting a particular model. This could help tune the models in order to better predict the detectability of subsequent events and may be an alternative way to gain insights into the population properties of the EM counterparts.

# 6. Conclusions

The electromagnetic follow-up of a GW event is one of the major challenges that transient astronomy will face in subsequent years. The large localization uncertainty regions and the relatively low expected luminosity of the candidate counterparts call for highly optimized observation strategies. The results of this work show that information from the GW signal can be used to make event-specific adjustments to the

EM follow-up strategy, and that such adjustments can significantly improve the effectiveness of the search at least in some cases, as shown in the example in Section 4. Advances in the theoretical modeling of the EM counterparts (e.g., in our ability to predict the amount of mass ejected during the merger) will increase the effectiveness of this approach, and are thus of great importance for the astronomy community as well.

We wish to thank Leo Singer and Christopher Berry for making the posterior samples of the F2Y study available to us. We thank Man Leong Chan, Christopher Messenger, and Zhiping Jin for pointing out some relevant references that were missing in the first draft. M.B. acknowledges financial support from the Italian Ministry of Education, University and Research (MIUR) through grant FIRB 2012 RBFR12PM1F. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 653477. We also acknowledge financial support from the UnivEarthS Labex programme at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02).

# Appendix
## Testing our Inverse-distance-weighting-based Method of Sky-position-conditional Density Estimation

### A.1. Test 1: Reconstruction of the Position-conditional Mean and Standard Deviation of a Known Distribution

As a first test, we constructed a set of 10,000 mock posterior samples whose underlying probability distribution in parameter space has a known analytic form. The parameter space is four-dimensional, the parameters being right ascension, R.A.; declination, decl.; luminosity distance, $d_L$; and a fourth quantity, $q$, with no physical meaning. The R.A., decl., and $d_L$ parameters are independent. The R.A. is normally distributed with mean 12 hr and sigma 30 minute, the decl. is uniformly distributed between -7°30' and +7°30', the luminosity distance is normally distributed with mean 100 Mpc and sigma 20 Mpc, and the distribution of the quantity $q$ depends on R.A.: its position-conditional distribution is normal, with mean and sigma given by

$$ 
\begin{aligned}
\mu_q(\text{R.A.}) &= \sin^2 \left[ 8\pi \left( \frac{\text{R.A.}}{12 \text{ hr}} - 1 \right) \right] \\
\sigma_q(\text{R.A.}) &= \frac{1}{5} + \frac{2}{3} \left( \frac{\text{R.A.}}{12 \text{ hr}} - 1 \right).
\end{aligned} 
\quad (17) $$

The 90% position probability area of the posterior samples is about 435 deg² wide and its shape is approximately a rectangle extending in R.A. from 11 to 13 hr and in decl. from -7°5' to +7°5'. Figure 10 shows the reconstructed position-conditional means and standard deviations computed in a set of points along the decl. = 0° axis. Both moments are reconstructed with an acceptable accuracy within the 90% sky-position probability area.

### A.2. Test 2: Comparison with the “Going the Distance” Study

The “Going the Distance” study (Singer et al. 2016a, GTD hereafter) and especially the related supplement (Singer et al. 2016b) represent an important practical step toward the use of posterior distributions of parameters other than the sky position to inform and improve the electromagnetic follow-up. In their

<footnote n="17">¹⁷ In the specific case adopted in this paper, actually, the only two free parameters were the luminosity distance $d_L$ and the viewing angle $\theta_v$. Since the involved redshifts are very low, it was sufficient to compute a small number of base light curves for each observing frequency, at a fixed distance, with varying viewing angles. The actual light curves were then computed by interpolation of the base light curves over the viewing angle, after rescaling to the correct distance, ignoring the negligible change in the rest-frame frequency corresponding to the observer frame frequency considered. Indeed, interpolation over a table of precomputed light curves can be a general way of reducing the computational cost of the evaluation of the EM counterpart models.</footnote>
<figure><img src="image_11.png" /></figure>

**Figure 10.** Test to assess the capability of our method to recover the mean and the standard deviation of the sky-position-conditional posterior distribution of a quantity. Left panel: the blue line represents the position-conditional mean of the true underlying distribution of quantity *q*, while the red dots are the position-conditional means derived with our method. The black dotted line shows the density of samples around each point normalized to the maximum density, while the black dashed vertical lines show the approximate R.A. limits of the 90% position probability area. The red crosses in the lower panel are the residuals of the computed means with respect to the true values. The accuracy in the reconstruction of the means clearly depends on the density of samples in the surrounding area. Right panel: same as the left panel, but for the position-conditional standard deviation. The reconstruction accuracy is clearly lower than in the case of the mean, but it remains acceptable in the region of high sample density.

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

**Figure 11.** Comparison between the mean of the sky-position-conditional posterior distribution of luminosity distance of injection 18951 of Singer et al. (2014) as computed with our method and that given in Singer et al. (2016b). Left panel: the color coding shows the fractional deviation of the mean luminosity distance of our method compared to that of Singer et al. (2016b; “GTD” stands for “Going the Distance,” i.e., the title of Singer and collaborators’ work). The outer (inner) red boundary represents the contour of the sky area containing 90% (50%) of the posterior sky-position probability. The star marks the actual position of the injection. Right panel: same as the left panel, but the comparison is on the standard deviation.

approach, distance information encoded in the signal is used in conjunction with galaxy catalogs as the basis for a follow-up strategy based on pointing candidate host galaxies to maximize the counterpart detection probability. In the supplement, Singer and collaborators show a step-by-step procedure to download and visualize the sky-position-conditional posterior distribution of the luminosity distance of injection 18951 from the F2Y study. We took that procedure as a starting point and used it to compare the sky-position-conditional mean and standard deviation of luminosity distance derived with our method to those of the GTD study. Figure 11 shows the relative difference between the quantities computed with the two methods. Again, the difference is very small except for regions where the density of the posterior samples is small, i.e., at the borders and outside the 90% sky-position confidence region.

## ORCID iDs

Om Sharan Salafia https://orcid.org/0000-0003-4924-7322

## References

* Abadie, J., Abbott, B. P., Abbott, R., et al. 2012, A&A, 539, A124

* Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, ApJ, 832, L21

* Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, PhRvL, 116 arXiv:1602.03840

* Abbott, B. P., Abbott, R., Abbott, T. D. M. C., et al. 2016b, ApJL, 826 arXiv:1602.08492

* Abbott, B. P., Abbott, R., Abernathy, M. R., et al. 2016c, LRR, 19, 1

* Ackermann, M., Atwood, W. B., Axelsson, M., et al. 2010, ApJ, 716, 1178

* Barnes, J., & Kasen, D. 2013, ApJ, 775, 18

* Barnes, J., Kasen, D., Wu, M.-R., & Martinez-Pinedo, G. 2016, ApJ, 829, 110

* Barsotti, L., & Fritschel, P. 2012, Early aLIGO Configurations: Example Scenarios Toward Design Sensitivity, Tech. Rep LIGO-T1200307-v4
Bauswein, A., & Stergioulas, N. 2015, PhysRevD, 91, 124056

Berger, E. 2014, ARA&A, 52, 43

Berry, C. P. L., Mandel, I., Middleton, H., et al. 2014, ApJ, 804, 114

Booth, R. S., de Blok, W. J. G., Jonas, J. L., & Fanaroff, B. 2009, arXiv:0910.2935

Chan, M. L., Hu, Y.-M., Messenger, C., Hendry, M., & Heng, I. S. 2017, ApJ, 834, 84

Ciolfi, R., Kastaun, W., Giacomazzo, B., et al. 2017, Phys Rev D, 95, 063016

Costa, E., Frontera, F., Heise, J., et al. 1997, Natur, 387, 783

Coughlin, M., & Stubbs, C. 2016, ExA, 42, 165

de Mink, S. E., & Mandel, I. 2016, MNRAS, 460, 3545

Dessart, L., Ott, C., Burrows, A., Rosswog, S., & Livne, E. 2008, ApJ, 690, 1681

Díaz, M. C., Beroiz, M., Peñuela, T., et al. 2016, ApJL, 828, L16

Dietrich, T., & Ujevic, M. 2017, CQG, 34, 105014

Dietrich, T., Ujevic, M., Tichy, W., Bernuzzi, S., & Brügmann, B. 2017, PhRvD, 95, 024029

Dominik, M., Berti, E., O'Shaughnessy, R., et al. 2015, ApJ, 806, 263

East, W. E., Paschalidis, V., Pretorius, F., & Shapiro, S. L. 2015, PhRvD, 93, 024011

Eichler, M., Arcones, A., Käppeli, R., et al. 2016, Journal of Physics: Conference Series, 665, 012054

Evans, P. A., Kennea, J. A., Palmer, D. M., et al. 2016, MNRAS, 462, 1591

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

Foucart, F. 2012, PhRvD, 86, 124007

Freiburghaus, C., Rosswog, S., & Thielemann, F.-K. 1999, ApJL, 525, L121

Gehrels, N., Cannizzo, J. K., Kanner, J., et al. 2016, ApJ, 820, 136

Ghirlanda, G., Burlon, D., Ghisellini, G., et al. 2014, PASA, 31, e022

Ghirlanda, G., Ghisellini, G., & Nava, L. 2009, A&A, 510, L7

Ghirlanda, G., Salafia, O. S., Pescalli, A., et al. 2016, A&A, 594, A84

Ghirlanda, G., Salvaterra, R., Campana, S., et al. 2015, A&A, 578, A71

Ghisellini, G., Nardini, M., Ghirlanda, G., & Celotti, A. 2009, MNRAS, 393, 253

Ghosh, S., Bloemen, S., Nelemans, G., Groot, P. J., & Price, L. R. 2015, A&A, 592, A82

Giacomazzo, B., Perna, R., Rezzolla, L., et al. 2013, ApJL, 762, L18

Giacomazzo, B., Zrake, J., Duffell, P. C., MacFadyen, A. I., & Perna, R. 2015, ApJ, 809, 39

Glendenning, N., & Moszkowski, S. 1991, PhRvL, 67, 2414

Gorski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Gottlieb, O., Nakar, E., & Piran, T. 2017, arXiv:1705.10797

Grossman, D., Korobkin, O., Rosswog, S., & Piran, T. 2013, MNRAS, 439, 757

Hanna, C., Mandel, I., & Vousden, W. 2013, ApJ, 784, 8

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

Hotokezaka, K., Wanajo, S., Tanaka, M., et al. 2015, MNRAS, 459, 35

Ivezic, Z., Tyson, J. A., Abel, B., et al. 2008, arXiv:0805.2366

Jin, Z.-P., Hotokezaka, K., Li, X., et al. 2016, NatCo, 7, 12898

Jin, Z.-P., Li, X., Cano, Z., et al. 2015, arXiv:1507.07206

Just, O., Obergaulinger, M., Janka, H.-T., Bauswein, A., & Schwarz, N. 2016, ApJL, 816, L30

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

Kasen, D., Fernandez, R., & Metzger, B. 2014, MNRAS, 450, 1777

Kasliwal, M. M., Cenko, S. B., Singer, L. P., et al. 2016, ApJL, 824, L24

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

Kim, H.-J., Yoon, S.-C., & Koo, B.-C. 2015, ApJ, 809, 131

Kiuchi, K., Kyutoku, K., Sekiguchi, Y., Shibata, M., & Wada, T. 2014, PhRvD, 90, 041502

Korobkin, O., Rosswog, S., Arcones, A., & Winteler, C. 2012, MNRAS, 426, 1940

Lackey, B. D., Nayyar, M., & Owen, B. J. 2006, PhRvD, 73, 024021

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

LIGO Scientific Collaboration, Virgo Collaboration, Abadie, J., et al. 2010, CQGra, 27, 173001

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

Lippuner, J., & Roberts, L. F. 2015, ApJ, 815, 82

Lithwick, Y., & Sari, R. 2001, ApJ, 555, 540

Loeb, A. 2016, ApJL, 819, L21

Martin, D., Perego, A., Arcones, A., et al. 2015, ApJ, 813, 2

Merloni, A., Predehl, P., Becker, W., et al. 2012, arXiv:1209.3114

Meszaros, P., & Rees, M. J. 1996, ApJ, 476, 232

Metzger, B. D., & Berger, E. 2011, ApJ, 746, 16

Nakar, E., & Piran, T. 2011, Natur, 478, 82

Nissanke, S., Kasliwal, M., & Georgieva, A. 2013, ApJ, 767, 124

Nuttall, L. K., & Sutton, P. J. 2010, PhRvD, 82, 102002

Özel, F., & Freire, P. 2016, ARA&A, 54, 401

Panaitescu, A., & Kumar, P. 2002, ApJ, 571, 779

Pankow, C., Sampson, L., Perri, L., et al. 2016, ApJ, 834, 154

Patricelli, B., Razzano, M., Cella, G., et al. 2016, JCAP, 2016, 056

Perego, A., Yasin, H., & Arcones, A. 2017, JPhG, 44, 84007

Perna, R., Lazzati, D., & Giacomazzo, B. 2016, ApJ, 821, L18

Radice, D., Galeazzi, F., Lippuner, J., et al. 2016, MNRAS, 460, 3255

Rana, J., Singhal, A., Gadre, B., Bhalerao, V., & Bose, S. 2017, ApJ, 838, 108

Rezzolla, L., Baiotti, L., Giacomazzo, B., Link, D., & Font, J. A. 2010, CQG, 27, 114105

Rhoads, J. 1997, ApJL, 487, L1

Rosswog, S., Feindt, U., Korobkin, O., et al. 2016, CQGra, 34, 104001

Rosswog, S., Freiburghaus, C., & Thielemann, F. K. 2000, NuPhA, 688, 344

Rosswog, S., Korobkin, O., Arcones, A., Thielemann, F. K., & Piran, T. 2014, MNRAS, 439, 744

Ruiz, M., Lang, R. N., Paschalidis, V., & Shapiro, S. L. 2016, ApJL, 824, L6

Salafia, O. S., Ghisellini, G., Pescalli, A., Ghirlanda, G., & Nappo, F. 2016, MNRAS, 461, 3607

Schutz, B. F. 2011, CQGra, 28, 125023

Sekiguchi, Y., Kiuchi, K., Kyutoku, K., & Shibata, M. 2015, PhRvD, 91, 064059

Sekiguchi, Y., Kiuchi, K., Kyutoku, K., Shibata, M., & Taniguchi, K. 2016, PhRvD, 93, 124046

Shepard, D. 1968, Proc. 1968 ACM National Conf. (New York: ACM Press), 517

Siegel, D., & Ciolfi, R. c. 2016, AJ, 819, 14

Silverman, B. W. 1982, AnSta, 10, 795

Singer, L. P., Chen, H.-Y., Holz, D. E., et al. 2016a, ApJL, 829, L15

Singer, L. P., Chen, H.-Y., Holz, D. E., et al. 2016b, ApJS, 226, 10

Singer, L. P., & Price, L. R. 2016, PhRvD, 93, 024013

Singer, L. P., Price, L. R., Farr, B., et al. 2014, ApJ, 795, 105

Soares-Santos, M., Kessler, R., Berger, E., et al. 2016, ApJL, 823, L33

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

Tanvir, N. R., Levan, A. J., Fruchter, A. S., et al. 2013, Natur, 500, 547

The LIGO Scientific Collaboration & the Virgo Collaboration 2016a, PhRvX, 6, 041015

The LIGO Scientific Collaboration & the Virgo Collaboration 2016b, PhRvL, 116, 241103

The LIGO Scientific Collaboration & the Virgo Collaboration 2016c, PhRvL, 116, 061102

The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, B. P., et al. 2017, PhRvL, 118, 221101

van Eerten, H. J., van der Horst, A. J., & MacFadyen, A. I. 2011, ApJ, 749, 44

Wanajo, S., Sekiguchi, Y., Nishimura, N., et al. 2014, ApJL, 789, L39

Wanderman, D., & Piran, T. 2010, MNRAS, 406, 1944

Yamazaki, R., Asano, K., & Ohira, Y. 2016, PTEP, 2016, 051E01

Yang, B., Jin, Z.-P., Li, X., et al. 2015, NatCo, 6, 7323