Where and when: optimal scheduling of the electromagnetic follow-up of gravitational-wave events based on counterpart lightcurve models

TL;DR

This paper introduces a method to optimize electromagnetic follow-up observations of gravitational-wave events by creating detectability maps based on lightcurve models and GW data, improving the chances of transient detection.

Contribution

It presents a novel approach that uses GW signal information to generate time-dependent detectability maps for efficient follow-up scheduling.

Findings

Successful detection of simulated afterglow and macronova emissions within hours.

Optimized observing strategies improve detection probability for large sky localizations.

Method applicable to optical, infrared, and radio follow-up facilities.

Abstract

The electromagnetic (EM) follow-up of a gravitational wave (GW) event requires to scan a wide sky region, defined by the so called "skymap", for the detection and identification of a transient counterpart. We propose a novel method that exploits information encoded in the GW signal to construct a "detectability map", which represents the time-dependent ("when") probability to detect 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 lightcurves is the basis to construct 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 (about 1500 square degrees). We construct observing strategies based on the detectability maps for optical, infrared and radio facilities, taking VST, VISTA and MeerKAT as prototypes. Assuming limiting fluxes of r ~ 24.5, J ~ 22.4 (AB magnitudes) and 500 uJy @ 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 hours respectively.

Tables

Table 1. Table 1: Algorithm input parameters used to construct the example follow-up strategies.

Instrument	${\nu }_{\mathrm{obs}}$ [Hz]	$F_{\lim }$ [$\mu$Jy]	$A_{\mathrm{FoV}}$ [deg2]	$T_{\mathrm{int}}$ [s]	$\lambda$	$N_{\mathrm{side}}$	Available timea	Det. timeb [h]
VST-like (shallow)	4.8$\times {10}^{14}$ (r filter)	4 (22.4 mAB)	1	100	0.01	64	3 h/night	-
VST-like (deep)	4.8$\times {10}^{14}$ (r filter)	0.58 (24.5 mAB)	1	1000	0.05	64	3 h/night	7
MeerKAT-like	1.4$\times {10}^9$	500	1.7	1000	0.05	64	20%	4.7
VISTA-like	2.4$\times {10}^{14}$ (J filter)	4 (22.4 mAB)	1.5	1200	0.5	64	3 h/night	14.5

Equations

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

$$P(\det |\mathrm{strategy})=\sum _{\mathrm{i}=1}^N{\int }_{A_{\mathrm{FoV}}({\alpha }_{\mathrm{i}})}P(\alpha |\mathcal{S})\mathrm{d}\alpha \times P(F(t_{\mathrm{obs},\mathrm{i}})>F_{\lim })$$

$$\theta_{\mathrm{v}}(\iota)=\left\{\begin{array}[]{lcr}\iota&&0\leq\iota<\pi/2\\ \pi-\iota&&\pi/2\leq\iota<\pi\\ \end{array}\right.$$

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

$$P(\mathrm{det}\,|\,d_{\mathrm{L}},\iota)=\left\{\begin{array}[]{lr}1&\mathrm{if}\;d_{\mathrm{L}}<d_{\mathrm{L,max}}(\iota)\\ 0&\mathrm{otherwise}\\ \end{array}\right.$$

$$P(d_{\mathrm{L}},\iota\,|\,\mathrm{det})\propto\left\{\begin{array}[]{lr}d_{\mathrm{L}}^{2}\sin\iota&\mathrm{if}\;d_{\mathrm{L}}<d_{\mathrm{L,max}}(\iota)\\ 0&\mathrm{otherwise}\\ \end{array}\right.$$

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

$$w_{\mathrm{i}}\propto \exp \left[-\frac{1}{2}{\left(\frac{\delta (\alpha ,{\alpha }_{\mathrm{i}})}{\sigma (\alpha )}\right)}^2\right]$$

$$\sigma (\alpha )=\sqrt{\frac{{\sum }_{\mathrm{i}=1}^N{\left[\delta (\alpha ,{\alpha }_{\mathrm{i}})-⟨\delta ⟩\right]}^2}{N-1}}\times N^{-1/5}$$

$${⟨q⟩}_{\alpha }=\sum _{\mathrm{i}=1}^N{w}_{\mathrm{i}}{q}_{\mathrm{i}}$$

$${\mathrm{Var}}_{\alpha }(q)={\left(\sum _{\mathrm{i}=1}^N{w}_{\mathrm{i}}{q}_{\mathrm{i}}\right)}^2-\sum _{\mathrm{i}=1}^N{w}_{\mathrm{i}}{q}_{\mathrm{i}}^2$$

$$P(q|\alpha ,\mathcal{S})\sim \sum _{\mathrm{i}=1}^N{w}_{\mathrm{i}}\mathrm{K}\left(\frac{q-q_{\mathrm{i}}}{{\sigma }_q}\right)$$

$$P(F(t)>F_{\lim }|\alpha ,\mathcal{S})=\sum _{\mathrm{i}=1}^N{w}_{\mathrm{i}}{\int }_{F_{\lim }}^{\infty }\mathrm{K}\left(\frac{F-F_{\mathrm{i}}(t)}{{\sigma }_F}\right)\mathrm{d}F$$

$$P(F(t)>F_{\lim }|\alpha ,\mathcal{S})\sim \sum _{\mathrm{i}=1}^N{w}_i\mathrm{H}(F_{\mathrm{i}}(t)-F_{\lim })$$

$$\left\{\begin{array}[]{l}\displaystyle t_{\mathrm{E,\lambda}}(\alpha)=\inf\left\{t\,|\,p(t,\alpha)\geq\lambda\right\}\\ \rule{0.0pt}{16.0pt}\displaystyle t_{\mathrm{L,\lambda}}(\alpha)=\sup\left\{t\,|\,p(t,\alpha)\geq\lambda\right\}\\ \end{array}\right.$$

$$t_{\mathrm{B}}(\alpha )=\mathrm{argmax}\left\{p(t,\alpha )\right\}$$

$$\begin{array}[]{l}\displaystyle\mu_{q}(\mathrm{RA})=\sin^{2}\left[8\pi\left(\frac{\mathrm{RA}}{12\mathrm{h}}-1\right)\right]\\ \rule{0.0pt}{20.0pt}\displaystyle\sigma_{q}(\mathrm{RA})=\frac{1}{5}+\frac{2}{3}\left(\frac{\mathrm{RA}}{12\mathrm{h}}-1\right)\\ \end{array}$$

Full text

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

††journal: ApJ

1 Introduction

The first detection (The LIGO Scientific Collaboration & the Virgo Collaboration, 2016c) of gravitational waves (GW 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, analysis of the GW signal and comparison with carefully constructed templates of the waveform (The LIGO Scientific Collaboration & the Virgo Collaboration, 2016a; Abbott et al., 2016b) enables the extraction of precious information about the parameters of the binary and of the remnant. The identification of an electromagnetic (EM) counterpart would increase further the scientific outcome of the detection, e.g. by enabling the identification of the host galaxy, by providing hints about the environment surrounding the merger, by constraining theoretical models of EM counterparts and by 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. Singer et al. 2014; Berry 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 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 lightcurve 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., 2016a). 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.  Evans et al., 2016; Díaz 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 Yamazaki et al. 2016; Perna et al. 2016; Loeb 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., 2016; Gottlieb et al., 2017) or the spindown-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 to apply 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_{\mathrm{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_{\mathrm{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_{\mathrm{jet}}$ is somewhere between $5^{\circ}$ and $15^{\circ}$ (e. g.  Berger, 2014), the prompt emission goes undetected in the majority of cases (for an isotropic population, the probability that $\theta_{\mathrm{v}}<15^{\circ}$ is less than 2 percent). 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 fo 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 orphan. No convincing detection of such a transient has been claimed to date, consistently with predictions for current and past surveys (Ghirlanda et al., 2015, 2014), 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; Kim et al., 2015; Dominik et al., 2015; de Mink & Mandel, 2016; The LIGO Scientific Collaboration et al., 2016) and the rather low expected fraction of GW events with an associated SGRB jet pointing at the Earth (Ghirlanda et al., 2016; Patricelli et al., 2016; Wanderman & Piran, 2014; Metzger & Berger, 2011), 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 dates back to almost twenty years ago (Li & Paczyński, 1998), the understanding of the 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:

• •

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

• •

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

• •

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

• •

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 begin 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 to predict the lightcurve have been developed recently (Barnes & Kasen, 2013; Grossman et al., 2013; Kawaguchi et al., 2016; Dietrich & Ujevic, 2016; Barnes et al., 2016; Rosswog et al., 2016). The emission from the dynamical ejecta is generally thought to be isotropic, which is an advantage with respect to the SGRB afterglow (which is instead beamed) from the point of view of the EM follow-up. The energy reservoir is the ejected mass $M_{\mathrm{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; Yang et al., 2015; Jin et al., 2015, 2016) are in the process of being tested by intensive observational campaigns. All this makes the macronova emission an extremely interesting candidate EM counterpart.

1.5 Outline of this work

In §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 $\times$ typical transient duration), stressing that the GW “skymap” (the sky position probability density) gives information about where to observe, but not about when. In §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 §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 §2.2.2). In §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 has a dependence 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 §3 we introduce a method to extract such information from the “posterior samples” obtained from the analysis of a GW signal, and in §4 we apply it to a synthetic example to show how it optimizes the follow-up strategy. Finally, we discuss the results in §5 and we draw our conclusions in §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\,|\,\mathcal{S})$. It represents the probability per unit solid angle that the source is at sky position $\alpha$, say $\alpha=(\rm RA,\,Dec)$, given the GW signal $\mathcal{S}$ detected by the interferometers ($\mathcal{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 “skymap probability”. Imagine 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:

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

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

observe the patches in decreasing order of skymap probability111to keep the discussion as simple as possible, we are neglecting the limitations due to observing conditions.; 4. 4.

for each patch:

1. (a)

identify the new sources by comparison with archive images; 2. (b)

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

The expected EM counterparts are transients, 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 model, we can construct a prior probability $P(F(t)>F_{\mathrm{lim}})$ that the transient flux (in a chosen band) $F$ is above some limiting flux $F_{\mathrm{lim}}$ at a given time $t$ after the GW event. Hereafter, we will call such quantity “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_{\mathrm{FoV}}$ and limiting flux $F_{\mathrm{lim}}$ is then

[TABLE]

where we are assuming a relatively small field of view in order to consider $P(\alpha\,|\,\mathcal{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_{\mathrm{lim}})$ is smaller.

Let us work in the simplifying assumption that all observations have the same exposure $T$ and the same limiting flux $F_{\mathrm{lim}}$. Let us denote by $T_{\mathrm{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_{\mathrm{tr}}$ (we refer to this set as a “volume” because it is 3-dimensional, even though the dimensions are $\mathrm{solid}\,\mathrm{angle}\times\mathrm{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_{\mathrm{obs,i}}=A_{\mathrm{FoV}}(\alpha_{\mathrm{i}})\times(t_{\mathrm{obs,i}},\,t_{\mathrm{obs,i}}+T)$ (where $\alpha_{\mathrm{i}}$ and $t_{\mathrm{obs,i}}$ are respectively the sky coordinates of the center of the field of view and the starting time of the $\mathrm{i}$-th observation) in order to maximise the detection probability $P(\mathrm{det}\,|\,\mathrm{strategy})$, which can be written as

[TABLE]

with the constraint $NT\leq T_{\mathrm{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 a detection”. In this approach, the proper construction of the a priori detectability $P(F(t)>F_{\mathrm{lim}})$ is the 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_{\mathrm{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; then we associate to each of them a jet afterglow and a macronova, under some assumptions. The detectable fraction of lightcurves 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 assume222these are typical reference values for short GRBs, though they suffer from the still limited number of reliable measurements available. that all SGRB jets have an isotropic kinetic energy $E_{\mathrm{K}}=10^{50}$ erg and a half-opening angle $\theta_{\mathrm{jet}}=0.2$ radians (11.5 deg), and that they are surrounded by a relatively tenuous interstellar medium with constant number density $n_{\mathrm{ISM}}=0.01$ cm*-3*. We fix the microphysical parameters333We refer here to the standard synchrotron afterglow model, and we set the microphysical parameters $p=2.5$, $\epsilon_{\mathrm{e}}=0.1$ and $\epsilon_{\mathrm{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 lightcurves available. so that the only remaining parameters needed to predict the afterglow lightcurve of the SGRB are the distance $d_{\mathrm{L}}$ and the viewing angle $\theta_{\mathrm{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 that

[TABLE]

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 (2016), and we use them as inputs to compute the lightcurve following Grossman et al. (2013) (using a constant grey opacity $\kappa=10$ cm2 g*-1*), assuming a blackbody 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_{\mathrm{ej}}$ and velocity vej to the merger through the fitting formulas of Dietrich & Ujevic (2016), we assume the H4 equation of state (Lackey et al., 2006; Glendenning & Moszkowski, 1991) which has a mid-range stiffness444This translates into mid-range values of the corresponding $M_{\mathrm{ej}}$ and vej. among those which are compatible with the observational constraints (Özel & Freire, 2016).

First, we need to derive the proper distributions of 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_{\mathrm{L,max}}$ out to which a NS-NS inspiral can be detected depends only on the binary plane inclination $\iota$ with respect to the line of sight, namely

[TABLE]

where $d_{\mathrm{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_{\mathrm{L}})\propto d_{\mathrm{L}}^{2}$ and $P(\iota)\propto\sin\iota$. By the above assumptions, the probability that a binary with luminosity distance $d_{\mathrm{L}}$ and inclination $\iota$ is detected is

[TABLE]

By Bayes’ theorem, the probability distribution of distance and inclination of a detected NS-NS inspiral is then $P(d_{\mathrm{L}},\iota\,|\,\mathrm{det})\propto P(\mathrm{det}\,|\,d_{\mathrm{L}},\iota)\times P(d_{\mathrm{L}})\times P(\iota)$, which gives

[TABLE]

The corresponding probability distribution of inclination for detected inspirals (which is obtained by marginalisation of Eq. 6 over $d_{\mathrm{L}}$) is then the well known

[TABLE]

For what concerns the distribution of 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_{\mathrm{lim}})$ we adopt the following Monte Carlo approach:

we construct our synthetic population of $N$ inspirals sampling distances and inclinations from $P(d_{\mathrm{L}},\iota\,|\,\mathrm{det})$ and masses $M_{1}$ and $M_{2}$ from the assumed normal distribution; 2. 2.

we compute the flux $F_{\mathrm{i}}=F_{\mathrm{A}}(t,d_{\mathrm{L,i}},\theta_{\mathrm{v}}(\iota_{\mathrm{i}}))$ of the jet afterglow or $F_{\mathrm{i}}=F_{\mathrm{M}}(t,d_{\mathrm{L,i}},M_{\mathrm{1,i}},M_{\mathrm{2,i}})$ of the macronova in the chosen band for each sample; 3. 3.

we estimate $P(F(t)>F_{\mathrm{lim}})$ as the fraction of the $F_{\mathrm{i}}$’s that exceed $F_{\mathrm{lim}}$.

Figure 2 shows the a priori detectability computed with the above method for $d_{\mathrm{L,max}}(0)=100$ Mpc (which corresponds roughly to the sky-position averaged aLIGO range for an optimally oriented NS-NS inspiral with the O1 sensitivity, see Abbott et al., 2016c) for Radio, Infrared and Optical observations (see the caption for details). The sensitivity of Optical observations has been chosen to match the limiting magnitude of the VST follow-up of GW150914, see Abbott et al. 2016a. 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_{\mathrm{K}}$, the ISM number density $n_{\mathrm{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 Fig. 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 about 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 $\mathcal{S}$ result in a posterior probability density $P(\xi\,|\,\mathcal{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\,|\,\mathcal{S})$, is essentially the $P(\xi\,|\,\mathcal{S})$ marginalised 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 §2.1.1 we already made use of two extrinsic parameters of the compact binary inspiral which are relevant for the SGRB afterglow and the macronova, namely the luminosity distance $d_{\mathrm{L}}$ and the binary inclination555the 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. $\iota$.

Recent works based on numerical simulations of NS-NS and BH-NS mergers (e. g.  Foucart, 2012; Giacomazzo et al., 2013; Hotokezaka & Piran, 2015; Kawaguchi et al., 2016; Dietrich & Ujevic, 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 equation of state 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_{\rm eff}$ of the black hole in the BH-NS case. Such information can be used to predict the observed lightcurve of the associated macronovas (e. g.  Kawaguchi et al., 2016; Dietrich & Ujevic, 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 lightcurve of the SGRB and/or of the dynamical ejecta macronova associated to the merger:

the luminosity distance $d_{\mathrm{L}}$ and the associated redshift $z$; 2. 2.

the orbital plane inclination $\iota$ with respect to the line of sight; 3. 3.

the component masses $M_{1}$ and $M_{2}$; 4. 4.

the effective spin $\chi_{\mathrm{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 spindown-powered transient described by Siegel & Ciolfi 2016) to those of the binary, whose distributions can be constrained by the GW signal.

In §4, we will use the fitting formulas provided in Dietrich & Ujevic (2016) to compute the posterior ejecta mass distribution associated to the example NS-NS inspiral treated in that section. We will refrain from deriving the SGRB jet energy from disk mass as suggested in Fig. 3, though, 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 §2.1 we introduced the idea of “a priori detectability” $P(F(t)>F_{\mathrm{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 $\mathcal{S}$ is observed, information it carries can be used to construct a posterior probability $P(F(t)>F_{\mathrm{lim}}\,|\,\mathcal{S})$ to better plan such observations. We will call it “a posteriori detectablity”. If the a priori detectability $P(F(t)>F_{\mathrm{lim}})$ is constructed using the prior distributions of the parameters of the EM transient model, the a posteriori detectability $P(F(t)>F_{\mathrm{lim}}\,|\,\mathcal{S})$ is obtained exactly the same way (as exemplified in §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 $\mathcal{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 which 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, distance, inclination, polarization angle, chirp mass and sky position of the binary share a certain degree of degeneracy: the same signal $\mathcal{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_{\mathrm{L}},\iota,M_{1},M_{2},...\,|\,\alpha,\mathcal{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_{\mathrm{lim}}\,|\,\alpha,\mathcal{S})$ which can be used as the basis of the EM follow-up strategy. We call this quantity “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:

• •

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 consist of scanning the localization uncertainty area, starting from the most probable sky location, until the source is found;

• •

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_{\mathrm{lim}})$ can be constructed, as shown in §2.1.1. This is the best follow-up timing information that can be constructed based on a priori knowledge only;

• •

After a GW signal $\mathcal{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_{\mathrm{lim}}\,|\,\mathcal{S})$ can be constructed. This exploits information contained in the GW signal, but it is still independent of the sky position;

• •

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 $\mathcal{S}$, compact binary inspiral parameters compatible with $\mathcal{S}$ and a particular sky position are in general different from those compatible with $\mathcal{S}$ and another sky postion. 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_{\mathrm{lim}}\,|\,\alpha,\mathcal{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 multi-dimensional kernel density estimation (KDE) technique, to be applied to the posterior samples obtained from a parameter estimation pipeline run on the gravitational wave 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 multi-dimensional 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\,|\,\alpha,\mathcal{S})$ of a quantity $q$ at sky position $\alpha$ is based on the concept of “inverse distance weighting” (Shepard, 1968): we assume that each posterior sample $\left\{\alpha_{\mathrm{i}},q_{\mathrm{i}},d_{\mathrm{L,i}},\iota_{\mathrm{i}},...\right\}$ contributes to the $P(q\,|\,\alpha,\mathcal{S})$ with a weight which is a decreasing function of the angular distance $\delta(\alpha,\alpha_{\mathrm{i}})$ between the posterior sample and the sky position $\alpha$. In particular, we assign a Gaussian weight to each posterior sample

[TABLE]

where the bandwidth $\sigma(\alpha)$ is taken as

[TABLE]

where $\left<\delta\right>$ is the arithmetic mean of the $\delta(\alpha,\alpha_{\mathrm{i}})$. The normalization of the weights is given by $\sum_{\mathrm{i=1}}^{N}w_{\mathrm{i}}=1$.

The ideas behind this method are simply that the closer the posterior sample is to sky-position $\alpha$, 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 (Eq. 9) is just “Silverman’s rule of thumb” (Silverman, 1982) for Gaussian KDE in a one dimensional parameter space (namely, the angular distance space).

The mean of $q$ at sky position $\alpha$ is thus computed as

[TABLE]

and similarly the variance

[TABLE]

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

[TABLE]

where $\mathrm{K}(x)$ is some kernel function, and $\sigma_{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

[TABLE]

where $F_{\mathrm{i}}(t)$ represents the flux (in the chosen band) at time $t$ of the lightcurve computed using the i-th posterior sample parameter values, $F_{\mathrm{i}}(t)=F(t,d_{\mathrm{L,i}},\iota_{\mathrm{i}},...)$. If we approximate the kernel functions with delta functions $K(x)\sim\delta(x)$, the expression becomes

[TABLE]

where $\mathrm{H}(x)$ is the Heaviside function.

3.3 The earliest, best and latest detection time maps

By the above method, 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 $\alpha$. By setting a minimum required detectability $\lambda$ one can define (for each sky position $\alpha$) a time interval during which $P(F(t)>F_{\mathrm{lim}}\,|\,\alpha,\mathcal{S})\geq\lambda$. If the detectability map $P(F(t)>F_{\mathrm{lim}}\,|\,\alpha,\mathcal{S})$ never reaches $\lambda$ 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

[TABLE]

where we set $p(t,\alpha)\equiv P(F(t)>F_{\mathrm{lim}}\,|\,\alpha,\mathcal{S})$ for ease of reading. Irrespectively of $\lambda$, the best detection time map can be defined as

[TABLE]

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 $\alpha$ is observed at a time as close as possible to $t_{\mathrm{B}}(\alpha)$, and in any case not earlier than $t_{\mathrm{E,\lambda}}(\alpha)$ or later than $t_{\mathrm{L,\lambda}}(\alpha)$. 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 (§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_{\mathrm{FoV}}$ at observing frequency $\nu_{\mathrm{obs}}$ and that the limiting flux $F_{\mathrm{lim}}$ for detection is independent of sky position and is always reached after an integration time $T_{\mathrm{int}}$. The outline of the algorithm is the following:

• •

we divide the skymap in patches, each representing a potential field to be observed;

• •

we define a list of available time slots (i. e. possible observing time windows) on the time axis, starting 24 hours 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 afterglow – but only that this method may not be applicable with very low latency666In the example of §4, actually, the only parameters needed to predict the possible SGRB afterglow light curves are $d_{\mathrm{L}}$ and $\iota$. Posterior distributions of these parameters (so called “extrinsic”) can be obtained with a low latency analysis tool such as BAYESTAR (Singer & Price, 2016).);

• •

starting from the patch with highest sky position probability, we check if the detectability map $P(F(t)>F_{\mathrm{lim}}\,|\,\alpha,\mathcal{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 not available anymore;

• •

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_{\mathrm{side}}$ of the HEALPix tessellation defines the number of pixels the sky is divided into, namely $N_{\mathrm{pixels}}=12N_{\mathrm{side}}^{2}$. The pixel area is then $A_{\mathrm{pixel}}=3438\,N_{\mathrm{side}}^{-2}$ deg2. Thus we replace the actual observations of duration $T_{\mathrm{int}}$ and field of view $A_{\mathrm{FoV}}$ with “pseudo observations” of area $A_{\mathrm{pixel}}$ and effective duration $T=T_{\mathrm{int}}\times A_{\mathrm{pixel}}/A_{\mathrm{FoV}}$, choosing $N_{\mathrm{side}}$ in order to minimize the difference between $A_{\mathrm{pixel}}$ and $A_{\mathrm{FoV}}$. The algorithm is thus implemented as follows:

we consider the posterior samples produced by a parameter estimation sampler (multiple sampling algorithms are implemented in LIGO’s LALINFERENCE parameter estimation tool777http://software.ligo.org/docs/lalsuite/lalinference/) applied to a simulated signal $\mathcal{S}$. These are points in the compact binary inspiral parameter space distributed according to the posterior probability density; 2. 2.

we find the 90% sky position confidence region of the source based on the posterior samples; 3. 3.

we divide the sky into pixels according to a HEALPix tessellation of order $N_{\mathrm{side}}$, and we consider only those pixels which fall inside the 90% sky position confidence region of the simulated signal; 4. 4.

we associate to each pixel $p$ the integral of the skymap probability density $P_{p}=\int_{A_{\mathrm{pixel}}(p)}P(\alpha\,|\,\mathcal{S})\,\mathrm{d}\Omega$ over the pixel area; 5. 5.

we associate to each pixel $p$ the earliest and latest detection times $t_{\mathrm{E,\lambda}}(p)$ and $t_{\mathrm{L,\lambda}}(p)$ averaged over that pixel; pixels for which the detectability never reaches $\lambda$ are excluded from the list of possible observations; 6. 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. 7.

we sort the pixels in order of decreasing $P_{p}$ and, starting from the first ($p=1$) pixel, we do the following:

1. (a)

we check that at least one available time slot is comprised between $t_{\mathrm{E,\lambda}}(p)$ and $t_{\mathrm{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; 2. (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_{\mathrm{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 to detect the transient. The inputs of the algorithm are:

the posterior sample list based on $\mathcal{S}$; 2. 2.

the detectability threshold $\lambda$; 3. 3.

the available time windows; 4. 4.

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

the HEALPix tessellation order $N_{\mathrm{side}}$, which should give a pixel area $A_{\mathrm{pixel}}$ close to $A_{\mathrm{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 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, a 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_{\mathrm{L}}=75$ Mpc, with orbital plane inclination $\iota=14^{\circ}$, at sky position $(\mathrm{RA},\mathrm{Dec})=(23^{\mathrm{h}}\,27^{\mathrm{m}}\,12^{\mathrm{s}}$, $-10^{\circ}\,30^{\prime}\,0^{\prime\prime})$, detected by the two-detector Advanced LIGO network on MJD 55483.27839 (i. e. at 06:40:53 of the 14th of October 2010 – 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_{\mathrm{K}}=10^{50}$ erg and half-opening angle $\theta_{\mathrm{jet}}=0.2$ rad (which is less than the viewing angle, thus the afterglow is orphan) is launched perpendicular to the orbital plane right after the merger, we computed its afterglow lightcurve assuming an ISM number density $n_{\mathrm{ISM}}=0.01\,\mathrm{cm}^{-3}$, adopting the same microphysical parameters as in §2.1.1, using BOXFIT v. 1.0 (van Eerten et al., 2011). The lightcurves at $\nu_{\mathrm{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_{\mathrm{K}}=10^{50}$ erg, in this example the only binary parameters which 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_{\mathrm{obs}}=1.4$ GHz produced using these samples (see the caption for additional information). The 90 percent sky position confidence area (represented by the larger red contours) covers approximately 1500 deg2.

4.1 Optical search

For our virtual optical EM follow-up, we adopt parameters inspired by the VST follow-up of GW150914 (Abbott et al., 2016a). We consider observations in the r band ($\nu_{\mathrm{obs}}=4.8\times 10^{14}$ Hz) with a detection limit $F_{\mathrm{lim}}=22.4$ AB magnitude, reached after $T_{\mathrm{int}}=100$ seconds (slew + integration). The field of view is $A_{\mathrm{FoV}}=1$ deg2. With the adopted flux limit, the optical lightcurve (Fig. 5) becomes too faint for detection after a few hours. Consistently, the detectability $P(F(t)>F_{\mathrm{lim}}\,|\,\alpha,\mathcal{S})$ is below 1 percent at all times $t>1$ d over the whole skymap except for a subregion of the skymap with a total area of 185 deg2. We run the follow-up strategy construction algorithm allowing 3 hours of available time per night from day 1 to day 15. We set $N_{\mathrm{side}}=64$ which gives $A_{\mathrm{pixel}}=0.84$ deg2 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 percent), 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 GW150914 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 to detect a relatively dim optical afterglow like that in Fig. 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_{\mathrm{int}}=1000$ s and $F_{\mathrm{lim}}=24.5$ mAB. 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 hours of telescope time. The $p=25$ observation (i. e. the 25th observation in descending order of $P_{p}$), scheduled $2\mathrm{d}\,2\mathrm{h}\,34\mathrm{m}\,45\mathrm{s}$ after the event, contains the EM counterpart. The flux at that time is 24.58 mAB, slightly dimmer than the required detection threshold, but detectable with a threshold S/N ratio of 5 under optimal observing conditions.

As explained in §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 hours 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 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 centres of the two large uncertainty regions (see Fig. 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 towards 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” image888Media release at http://www.ska.ac.za/media-releases/ in July 2016. We assume a field of view $A_{\mathrm{FoV}}=1.7$ deg2 at $\nu_{\mathrm{obs}}=1.4$ GHz. We conservatively estimate 50 $\mu$Jy rms noise for a $T_{\mathrm{int}}=1000$ s observation (slew + integration), assuming 16 working dishes. In a large area survey, usually a 10 sigma detection is required to avoid a large number of false alarms, thus we set $F_{\mathrm{lim}}=0.5$ mJy, i. e. we require the flux to be ten times the rms noise for the detection to be considered confident. We allow a maximum of 20 percent 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 hours 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$ d (382 deg2 in total). Setting $N_{\mathrm{side}}=64$, the follow-up construction algorithm outputs 337 “pseudo observations” (corresponding to 169 pointings, totalling 47 hours of telescope time), which are represented in Figure 8. The $p=35$ observation contains the counterpart. It is scheduled $5\mathrm{d}\,0\mathrm{h}\,9\mathrm{min}\,27\mathrm{s}$ after the event, when the flux of the EM counterpart is 0.52 mJy (see Fig. 5), which means that the afterglow is detected at better than 10 sigma. The first 35 (pseudo) observations make up 4.7 hours 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) lightcurve (see Figure 5) of the associated macronova with the same assumptions as in §2.1.1. Due to the rather large masses of the binary components and to their similar mass, the dynamically ejected mass is small ($M_{\mathrm{ej}}\approx 5.6\times 10^{-3}M_{\odot}$ according to the Dietrich & Ujevic (2016) fitting formula assuming the H4 equation of state). The lightcurve is thus quite dim. It peaks between the second and the third day, slightly brighter than 22.4 mAB in the J band. The effective temperature of the photosphere at peak is $T_{\mathrm{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 to VISTA: we choose a limiting flux $F_{\mathrm{lim}}=22.4$ mAB in the J band with $T_{\mathrm{int}}=1200$ s, and we set $A_{\mathrm{FoV}}=1.5$ deg2. Again, we assume that 3 hours 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 lightcurves exceed 22.4 mAB, thus we set $\lambda=0.5$ to limit the search to points of the skymap which 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 hours of telescope time), which are represented in Figure 9. The 79th observation in descending order of $P_{p}$ contains the counterpart. It is scheduled 6d 0h 18min 27s after the event, when the flux of the macronova in the J band is just below the 22.9 mAB and the effective temperature of the spectrum is around 2600 K. The flux is lower than $F_{\mathrm{lim}}$ (which was chosen to represent an indicative limit for achieving a S/N ratio $\sim$ 10), but according to ESO’s exposure time calculator999We queried the VIRCAM ETC at http://www.eso.org assuming a 1.2 airmass and a seeing of 0.8 arcsec., such emission would be detected through the VISTA telescope at ESO in the J band with a S/N ratio of 6, assuming a 1000 s integration, in optimal observing conditions. Our virtual Infrared follow-up would thus again result in a detection in a search with a threshold at S/N $\sim$ 5. To achieve it, 14.5 hours 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_{\mathrm{lim}}\,|\,\alpha,\mathcal{S})$ with the a priori detectability $P(F(t)>F_{\mathrm{lim}})$ computed in §2.1.1. This should simulate a search based on a priori information only. The results are the following:

• •

Optical search: the counterpart is in the field of view 5d 22h 53m 32s after the event, when the flux of the jet afterglow in the r band is as low as 26 mAB, which is definitely too faint for a detection;

• •

Radio search: the counterpart is in the field of view 20d 2h 18m 0s after the event, when the flux of the jet afterglow at 1.4 GHz is 120 $\mu$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 hours of telescope time for this single follow-up in order to include the observation that contains the counterpart;

• •

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 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 lightcurves. The possible lightcurves 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=(\mathrm{RA},\mathrm{Dec})$. In order to apply the idea to a practical example, we introduced an explicit method to construct the detectability maps (§3.1) and an algorithm which uses these maps to define an EM follow-up observing strategy (§3.4). We then applied the method to a synthetic example, showing that it improves significantly the effectiveness of the EM follow-up (§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 to mention some secondary details or to include 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 lightcurves are treated in a statistical sense, the models only need to represent correctly the peak flux and lightcurve general evolution. Fine details are lost in the processing of the lightcurves, and are thus unnecessary. Moreover, it is very straightforward to include our uncertainty on 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_{\mathrm{K}}$ and ISM density $n_{\mathrm{ISM}}$ of the SGRB afterglow, or the NS equation of state). 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 lightcurves 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 §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 can not 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: 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, 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 face again the same issues about 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 of 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 which 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_{\mathrm{L}}=75$ Mpc), the jet is slightly off-axis ($\iota=14^{\circ}$), the sky position uncertainty is large (more than $1500$ deg2) 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 centre 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 (§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:

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. 2.

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

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

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

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

Some recent works 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. (2016) developed and compared some ingenious algorithms that aim at maximizing 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. For what concerns 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 span 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 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 §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 a 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” lightcurves of an event must be computed (at least one per posterior sample, which means around $10^{4}$ lightcurves 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 lightcurves. Since the aim of the approach is to assess the detectability, rather than to fit the model to observational data, simple analytic models which capture the main features of the expected lightcurves are better than complex numerical models suited for parameter estimation. The macronova model by Grossman et al. (2013) is a good example: few minutes are sufficient to compute $10^{4}$ lightcurves on a laptop using this model. Parallelization is instead unavoidable101010In the specific case adopted in this paper, actually, the only two free parameters were the luminosity distance $d_{\mathrm{L}}$ and the viewing angle $\theta_{\mathrm{v}}$. Since the involved redshifts are very low, it was sufficient to compute a small number of base lightcurves for each observing frequency, at a fixed distance, with varying viewing angle. The actual lightcurves were then computed by interpolation of the base lightcurves 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 pre-computed lightcurves can be a general way to reduce the computational cost of the evaluation of the EM counterpart models. when such a large number of off-axis SGRB afterglow lightcurves 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 will have 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 to tune the models in order to better predict the detectability of subsequent events and may be an alternative way to get insights into the population properties of the EM counterparts.

6 Conclusions

The electromagnetic follow-up of a gravitational wave events is one of the major challenges that transient astronomy will face in the next 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 improve significantly the effectiveness of the search at least in some cases, as shown in the example in §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.

Acknowledgements

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. MB 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 A 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 10000 mock posterior samples whose underlying probability distribution in parameter space has a known analytic form. The parameter space is 4-dimensional, the parameters being right ascension $\mathrm{RA}$, declination $\mathrm{Dec}$, luminosity distance $d_{\mathrm{L}}$ and a fourth quantity $q$ with no physical meaning. The $\mathrm{RA}$, $\mathrm{Dec}$ and $d_{\mathrm{L}}$ parameters are independent. The right ascension is normally distributed with mean $12$h and sigma $30$min; the declination is uniformly distributed between $-7^{\circ}30^{\prime}$ and $+7^{\circ}30^{\prime}$; the luminosity distance is normally distributed with mean 100 Mpc and sigma 20 Mpc; the distribution of quantity $q$ depends on right ascension: its position-conditional distribution is normal, with mean and sigma given by

[TABLE]

The 90% position probability area of the posterior samples is about 435 deg2 wide and its shape is approximately a rectangle extending in right ascension from 11 to 13 hours and in declination from $-7^{\circ}5^{\prime}$ to $+7^{\circ}5^{\prime}$. Figure 10 shows the reconstructed position-conditional means and standard deviations computed in a set of points along the $\mathrm{Dec}=0^{\circ}$ 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 towards the use of posterior distributions of parameters other than the sky position to inform and improve the electromagnetic follow-up. In their approach, distance information encoded in the signal is used in conjunction with galaxy catalogues as the basis for a follow-up strategy based on pointing candidate host galaxies to maximise 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 posterior samples is small, i. e. at the borders and outside the 90% sky-position confidence region.