GALAXY STRATEGY FOR LIGO-VIRGO GRAVITATIONAL WAVE COUNTERPART SEARCHES

The advent of Advanced LIGO (aLIGO—Aasi et al. 2015) and Advanced Virgo (AdV—Acernese et al. 2015) heralds the dawn of a new age of discovery in which gravitational wave (GW) detections will supplement traditional electromagnetic (EM) detections, such as those by Swift (Gehrels et al. 2004). We refer to aLIGO/AdV in combination as LVC. These GW observatories will begin operating soon at a fraction of design capability, and within a few years should be able to detect neutron star–neutron star (NS–NS) mergers out to ∼450 Mpc and black hole–neutron star (BH–NS) mergers out to ∼900 Mpc (Aasi et al. 2013). These distances are referred to as the "horizon" for these events. Although black hole–black hole (BH–BH) mergers are detectable by aLIGO/AdV to even greater distances, since we are focused in this work on GW signals with accompanying EM counterparts, they will not be a part of our discussion.

It is important to distinguish between "horizon" and "range." The "inspiral range," the most commonly cited figure-of-merit regarding LVC sensitivity, is defined as the radius of a sphere whose volume equals the sensitivity volume within which a signal-to-noise ratio (S/N) $\geqslant \;8$ detection is achieved for a 1.4–1.4 ${M}_{\odot }$ NS–NS merger, averaged over all sky locations and binary inclinations. In this work we use μ to denote the NS–NS inspiral range. The antenna projection functions and associated averagings were discussed by Finn & Chernoff (1993) and Finn (1996). The inspiral range is not a hard upper limit, as face-on binary orbit mergers produce stronger signals (Dalal et al. 2006; Maggiore 2007; Schutz 2011; Nissanke et al. 2013). The maximum theoretical detection distance, the horizon, is 2.26 times the range (Finn & Chernoff 1993; see also Abadie et al. 2010). A simple calculation shows that in a general population with binaries of random inclinations and positions over the whole sky, the fraction of aLIGO detections one expects to pick up from beyond the nominal inspiral range is about half of the total (Nissanke et al. 2010, 2013; Singer 2015).

High S/N detections will permit important physical parameters to be measured that are not easily accessible through traditional means, such as the masses and spins of the merger components and the luminosity distance D_L to the merger. A joint EM/GW detection would provide powerful constraints on the merger physics. Various groups have begun to the examine the prospects of finding an EM counterpart to a GW trigger (Nuttall & Sutton 2010; Metzger & Berger 2012; Nissanke et al. 2013; Hanna et al. 2014; Chen & Holz 2015; Chan et al. 2015; Fan et al. 2015; Ghosh & Nelemans 2015). Evans et al. (2016) consider strategies for X-ray observations using Swift/X-ray Telescope (XRT).

In Section 2 we present an overview of short GRBs, including putative GW rates for NS–NS mergers, and a discussion of "kilonova" emission accompanying the decay of radionuclides. In Section 3 we describe a galaxy catalog—the Census of the Local Universe (CLU)—and quantify its completeness out to distances relevant for aLIGO. In Section 4 we estimate galaxy totals that are concomitant with putative aLIGO error volumes for the next 5 yr, if one restricts attention to bright galaxies. In Section 5 we discuss potential EM observations in the optical, X-ray, and radio which might follow GW detection, and in Section 6 we summarize our findings.

2.1. GW Radiation

2.2. GRBs and Beamed EM Radiation

2.3. Kilonova Emission

Any given galaxy catalog is generally not optimal for GW follow-up studies. Consider two extremes: The Two Micron All Sky Survey (2MASS) survey (Skrutskie et al. 2006) has good coverage in both the northern and southern skies, but does not go very deep (Huchra et al. 2012). The Millenium Galaxy Catalog, comprising spectroscopic redshifts of galaxies from 2dF or SDSS, is deep, but covers only a small slice along the celestial equator (Driver et al. 2005). An attempt to overcome these limitations led to the Gravitational Wave Galaxy Catalog = GWGC (White et al. 2011). Its only limitation for our current study is that it does not extend beyond 100 Mpc. One of us (M. M. Kasliwal et al. 2016, in preparation) has amassed a catalog based on the union of several existing catalogs—the CLU—which is suitable for GW+EM follow-up studies. As we shall show, for bright galaxies the CLU is complete out to the anticipated aLIGO GW inspiral range for NS–NS mergers up to 2020.

In order to show completeness, we must adopt a model for galaxy number density in the local universe. The Schechter luminosity function (Schechter 1976) provides a useful description of the space density of galaxies as a function of their luminosity, ${\rho }_{{\rm{gal}}}(x){dx}={\phi }^{*}{x}^{a}{e}^{-x}{dx}$, where $x=L/{L}^{*}$ and L* is a characteristic galaxy luminosity where the power-law form of the function truncates. It has proven to be applicable over up to 10 mag in deep surveys (e.g., Bonne et al. 2015). The CLU catalog is amassed from many different surveys. One of the CLU data columns, ${b}_{{\rm{tc}}}$, consists of apparent B-magnitudes m_B for entries where they are available, and pseudo-m_B values for sources from other bands, for instance 2MASS. Hence for this work (${\phi }^{*}$, L, ${L}^{*})\to ({{\phi }_{B}}^{*}$, L_B, ${{L}_{B}}^{*}$), Physically, ${L}_{B}^{*}$ represents the turn-over in the distribution between a power-law for low x and an exponential for high x.

We adopt the following values derived from B-band measurements of nearby field galaxies: ${\phi }^{*}=(1.6\pm 0.3)\times {10}^{-2}{h}^{3}$ Mpc⁻³, $a=-1.07\pm 0.07$, ${L}_{B}^{*}=(1.2\pm 0.1)\times {10}^{10}{h}^{-2}{L}_{B,\odot }$, with a corresponding ${M}_{B}^{*}=-19.7\pm 0.1+5{\mathrm{log}}_{{\rm{10}}}h=-20.47$ (e.g., Norberg et al. 2002; Liske et al. 2003; González et al. 2006, and references therein). By comparison, for the Milky Way galaxy ${M}_{B}=-20.42$. This is based on a Milky Way B-band luminosity $2.3\times {10}^{10}{L}_{B,\odot }$ (Carroll & Ostlie 1996), where ${L}_{B,\odot }$ is the solar B-band luminosity, and an absolute solar B magnitude of 5.48 (Allen 1973). We take h = 0.7 based on the latest weighted overlap between the Planck results and the rest of astronomy (Ade et al. 2014).

Integrating over luminosity gives integrated number density ${\rho }_{{\rm{0,gal}}}={\int }_{{x}_{1}}^{\infty }{\rho }_{{\rm{gal}}}(x){dx}$. Although ${\rho }_{{\rm{0,gal}}}\to \infty$ as ${x}_{1}\to 0$ for $a\lt -1$, the integrated luminosity density diverges only for $a\lt -2$. One has

$$\begin{eqnarray}&&{\displaystyle \int }_{{x}_{1}}^{\infty }{\phi }^{*}{L}^{*}{x}^{a+1}\mathrm{exp}(-x){dx}={\phi }^{*}{L}^{*}{\rm{\Gamma }}(a+2,{x}_{1}),\end{eqnarray} \tag{ 1 }$$

where Γ is the incomplete gamma function. For $a=-1$ the total luminosity density is ${\phi }^{*}{L}^{*}{\rm{\Gamma }}(2+a)={\phi }^{*}{L}^{*}=1.9\times {10}^{8}\;{{hL}}_{B,\odot }$ Mpc⁻³. Dividing by a Milky Way B-band luminosity yields a density $\sim 6\times {10}^{-3}$ MWE Mpc⁻³, where MWE = Milky Way equivalent galaxy. For $a=-1$, half of the luminosity density is contributed by galaxies with ${x}_{1/2}\gt 0.693$. For the power law of interest in this study, $a=-1.07$, the cutoff lies at ${x}_{1/2}\gt 0.626$, or ${M}_{B1/2}=-19.97$. This corresponds to ∼0.66 of the Milky Way luminosity. To arrive at this ${x}_{1/2}$ value we used the fact that ${\int }_{{x}_{1}}^{\infty }{x}^{a+1}{e}^{-x}{dx}=1.04559$ for ${x}_{1}=0$ and $a=-1.07$, and half this value $1.04559/2$ is achieved for ${x}_{1}=0.626$.

Figure 1 presents sky maps of the CLU catalog, showing all the galaxies, and also those for which $x\gt 1$, where $x={L}_{B}/{{L}_{B}}^{*}$. The dark strips evident in the top panel indicate individual deep surveys which make up the CLU. Restricting the sample to intrinsically bright galaxies cleans it up considerably and reveals large scale structure.

Zoom In Zoom Out Reset image size

Figure 2 provides a measure of the completeness as gauged by the Schechter function. We compare the CLU catalog with the GWGC and 2MASS redshift survey out to 200 Mpc using 12 radial slices. Within each slice we bin the data in $x=L/{L}^{*}$ and compare with the Schechter function weighted by the volume of the given slice. For CLU and GWGC $x={L}_{B}/{{L}_{B}}^{*}$, whereas for 2MASS $x={L}_{K}/{{L}_{K}}^{*}$. We adopt M_K^* = −23.1 + 5 log₁₀h = −23.87. Since GWGC ends at 100 Mpc there are no data in the more distant bins. For all three catalogs there is a progressive loss of fainter galaxies with distance relative to Schechter. For all radii the CLU follows the Schechter function for $x\gtrsim {x}_{1/2}$. In some of the panels in Figure 2 one sees spikes at $x\simeq 10$. These may be due to misidentified local objects spuriously placed at larger distances.

Zoom In Zoom Out Reset image size

In this paper we propose to utilize the brighter galaxies in each bin, which tends to mitigate the effects of incompleteness brought about by a progressive loss of fainter galaxies with increasing distance. The CLU catalog is fairly complete relative to Schechter above ${x}_{1/2}$. This is shown in Figure 3, where we compare with GWGC and 2MASS. In this work we will argue for including only the brighter galaxies in GW+EM follow-up studies. As noted earlier, this implies cutting galaxies below ${x}_{1/2}$. For the three fiducial distances of interest in the next section, 60 Mpc, 120 Mpc, and 180 Mpc, the CLU catalog is complete above this cut line at a level of ∼100%, ∼80%, and ∼40%, respectively. We stress that this completeness only includes galaxies for which $x\gt {x}_{1/2}$. Since these three distances are based on the NS–NS inspiral ranges for 2015, 2017, and 2020, an additional factor which must be included for the third year stems from the fact that the CLU catalog ends at 200 Mpc, namely the fact that about half the detections would be expected to arise from between the range and the horizon. Therefore for the third year the completeness is ∼20%. Both Figures 2 and 3 show some regions of "overcompleteness" below 100 Mpc. This is probably due to local over-densities such as the Virgo Supercluster.

Zoom In Zoom Out Reset image size

Our choice of ${x}_{{\rm{cut}}}={x}_{1/2}$ is motivated by a trade-off between including enough galaxies so that we encompass a reasonable fraction of putative sGRB host galaxies (e.g., Berger 2014; see his Figure 8) on the one hand, and not having incompleteness become too great an issue on the other hand. Taking ${x}_{{\rm{cut}}}={x}_{1/2}$ picks out about half of the sGRB host galaxies shown in Berger (2014), therefore the true efficiency factors for the three target years would be ∼50%, ∼40%, and ∼10%, respectively, after convolving the $x\gt {x}_{1/2}$ completeness with the fact that ${x}_{1/2}$ lies near the median in the observed sGRB x distribution.

Another caveat arises from the fact that by utilizing B-band data we are implicitly taking B-band luminosity as a proxy for the compact object merger rate. Although this assumption is generally supported by sGRB observations (Berger 2014), it may in fact not be a universally good proxy (e.g., de Freitas Pacheco et al. 2006; O'Shaughnessy et al. 2010; Hanna et al. 2014).

As regards to identifying a candidate host galaxy for a GW event, a galaxy catalog is only half of the picture. Of course the starting point is the GW localization. Initial localizations in the ∼2015 time frame are expected to be ∼500 deg², decreasing to ∼20 deg² by ∼2020 (Aasi et al. 2013; Singer et al. 2014; Berry et al. 2015). In terms of the depth, the NS–NS inspiral range should increase from ∼60 Mpc to ∼180 Mpc over the same time frame. A projected timeline of NS–NS inspiral range versus date is given by Aasi et al. (2013). One may obtain a good estimate using simple considerations. In the CLU catalog there are ${N}^{*}=27559$ galaxies brighter than ${{M}_{B}}^{*}$, and ${N}_{1/2}=47438$ brighter than ${M}_{B1/2}$. If one restricts ${N}_{1/2}$ based on the limited sky areas and volumes relevant for the three target years, i.e., 500 deg² × 60 Mpc for 2015, 100 deg² × 120 Mpc for 2017, and 20 deg² × 180 Mpc for 2020, one obtains ${N}_{{\rm{gal}}}\;=$ 26.0, 36.4, and 18.8, respectively. As noted previously, since these distances represent putative NS–NS inspiral ranges, for a localization at μ we must roughly double the totals to take into account going from range to horizon if the distance error is large (i.e., plus or minus a factor of two).

Realistic idealizations for aLIGO localizations have been undertaken by several recent groups that attempt to quantify the "volume reductions" which might be realized. Nissanke et al. (2010) use a Markov Chain Monte Carlo (MCMC) technique to calculate distance measurement errors associated with aLIGO localizations of astrophysical populations of NS–NS and BH–NS binaries, considering both isotropically oriented as well as beamed events. They take as their starting point a precise sky localization based on a coincident EM detection of the same GW event. They present Fisher-matrix-derived linear scalings for $[{\rm{\Delta }}{D}_{L}/{\rm{True}}\ {D}_{L}]$ for the two populations, assuming four GW detectors. If the EM emission from the NS–NS merger providing the coincident EM signal is isotropic, they find that, in combination with the precise sky position, the distance to NS–NS binaries can be measured with a fractional error of ∼20%–60%, with most events clustered near ∼20%–30%. If the EM emission from the NS–NS merger is beamed, with a ∼25° opening angle, then the error on the distance is reduced by a factor of ∼2 and much of the high error tail is eliminated. BH–NS events are measured more accurately: the distribution of fractional distance errors lies in the range ∼15%–50%, with most events clustered near ∼15%–25%. Nissanke et al. (2013) carry out extensive end-to-end simulations, looking at GW sky localization, distance errors, and volume errors using NS–NS and BH–NS mergers. They compare MCMC-derived distance measures with 3D volume measures. They outline optimal strategies to prepare for identifying EM counterparts of a GW merger.

L. P. Singer et al. (2016, in preparation) have created full 3D position reconstructions for a large population of simulated early aLIGO NS–NS events. They provide a simple approximation for the 3D distance distribution and qualitatively describe the shapes that emerge. In the present work, as we are more interested in the impact of the galaxy catalog we use an even simpler description: we assume that the localization subtends a given solid angle, and is a shell between two constant radii.

L. P. Singer et al. (2016, in preparation) calculate full 3D aLIGO reconstructions based on BAYESTAR, a rapid reconstruction method for BNS mergers, and LALInference, a full Bayesian parameter estimation code (Singer et al. 2014; Singer 2015; Veitch et al. 2015), for the near-future of aLIGO—2015 and 2016. In this work we consider a longer time frame, i.e., up to the full achievement of aLIGO design sensitivity. In order to calculate galaxy sky counts in error boxes of given areas on the sky, we take a simplified approach compared to Nissanke et al. (2010, 2013), and Singer (2015). Namely, rather than using a realistic 3D reconstruction, for a given line of sight (LOS) we take a simple top-hat windowing function. For our three putative aLIGO target years—2015, 2017, and 2020—we adopt $\mu =60$, 120, and 180 Mpc, respectively, as fiducial LVC NS–NS inspiral ranges.

We consider 1000 randomly selected LOSs over the sky, and then search the CLU catalog to find galaxies within an angular separation that would place them inside an error box on the sky of (i) 500 deg², (ii) 100 deg², or (iii) 20 deg² for the three cases. Although the actual sky-projected aLIGO localizations will be complicated, the more important factor is simply the total sky area involved. Only galaxies are considered for which their luminosity places them above the 50th percentile mark ${x}_{1/2}$. The CLU galaxy count in a radial bin is incremented only if $| (r-\mu )/\mu | \lt 0.3$ (Hanna et al. 2014), i.e., if the galaxy lies within a thick spherical shell with ${\rm{\Delta }}\mu /\mu =0.6$, where the NS–NS inspiral range $\mu =60$, 120, and 180 Mpc, respectively. This has the effect of essentially doubling the galaxy counts derived by truncating the integration at the range μ, thereby mimicking the effect of the range-to-horizon mapping (Hanna et al. 2014).

Figure 4 shows the results of this experiment. Each point is the mean of the 1000 individual LOS values, and the error bar is the standard deviation. Two main points are worth noting. The most obvious is simply that larger error boxes yield more galaxies. The second is that galaxy counts are only added out to $1.3\mu$, by construction. In reality, one would continue to pick up sources out to the horizon at $2.26\mu$. The total number of sources for $d\lt 2.26\mu$ is twice that for $d\lt \mu$ despite the factor 11.5 increase in volume (as many binaries are not at favorable inclination). To fully model this would require a variable efficiency factor with distance between the range and horizon, hence the motivation for our simple top-hat method. The CLU totals out to $r=\mu$ (for a localization in each year at μ) are ${N}_{{\rm{gal}}}=17.7\pm 5$, $22.5\pm 4$, and $10.5\pm 2$, respectively. As noted previously, the CLU is ∼100% complete above ${x}_{1/2}$ at 60 Mpc, ∼80% complete at 120 Mpc, and ∼40% complete at 180 Mpc. In addition, the NS–NS range-to-horizon mapping roughly doubles the totals for localizations $r\simeq \mu$ with a large radial uncertainty.

Zoom In Zoom Out Reset image size

For consistency we may consider galaxy counts derived directly from the Schechter function. The galaxy density above ${x}_{1/2}$ is ${\int }_{{x}_{1/2}}^{\infty }{\phi }^{*}{x}^{a}{e}^{-x}{dx}=2.35\times {10}^{-3}$ Mpc⁻³. If we then multiply this by a volume ${\rm{\Delta }}V=(4/3)\pi {\mu }^{3}({\rm{\Delta }}{\rm{\Omega }}/4\pi )$, where, for the three target years (μ, ${\rm{\Delta }}{\rm{\Omega }}$) = (60 Mpc, 500 deg²), (120 Mpc, 100 deg²), and (180 Mpc, 20 deg²), respectively, we obtain ${N}_{{\rm{gal}}}=25.8$, 41.3, and 27.8.

The results given in Figure 4 reveal the dramatic reduction in tiling requirements brought about by restricting a search methodology to the brighter galaxies. Were one simply to tile the entire 2015, 2017, and 2020 putative error boxes A=500, 100, and 20 deg² using an EM detector with a field of view (FOV) of $\delta A=0.1$ deg², i.e., $\sim 19\times 19$ arcmin, the number of tilings required would be $A/\delta A=5\times {10}^{3}$, 10³, and $2\times {10}^{2}$, respectively.

Concentrating on selected bright galaxies reduces the tiling effort considerably. The optimal size of a "tile" is dictated by observations of short GRBs—in particular their locations within their host galaxies. Fong et al. (2010) use precise HST localizations to study the cumulative distribution of projected physical offsets for short GRBs with sub-arcsecond positions and find them to lie within ∼100 kpc of their host galaxy centers. For a typical distance of interest in this study, 100 Mpc, this corresponds to ∼10⁻³ radian or ∼0.057 deg—a projected area on the sky of ∼0.010 deg² or ∼37 arcmin². This can be covered by an EM detector with a rectangular FOV of at least 7 × 7 arcmin or 49 arcmin². For such detectors the number of tilings will reduce to the number of galaxies as shown in Figure 4, i.e., ∼20, ∼20, and ∼10, respectively. This represents a reduction in requisite tilings by ∼1–2 orders of magnitude over a brute force methodology, i.e., simply covering the entire GW error box.

Furthermore, the underlying strategy of focusing attention only on bright galaxies is strengthened by the observation that short GRBs tend to lie in the larger, brighter galaxies (in contrast to long GRBs, which lie preferentially in dwarf irregulars). Fong et al. (2010) compare the projected physical offset distribution between short and long GRBs and find that, when normalized to the sizes of their host galaxies, the distributions are indistinguishable. However, the absolute length scales differ by a factor of ∼5 (see also Fong & Berger 2013). Berger (2014, see his Figure 8) plots all known sGRB host galaxy $x={L}_{B}/{{L}_{B}}^{*}$ values; they span a range $0.1\lesssim x\lesssim 2$. Our adopted cut value in this study ${x}_{1/2}\simeq 0.6$ is roughly at the median of the observed distribution, therefore our bright galaxy strategy would pick up about half the sGRBs shown in Berger's sample. Therefore the total effective completeness would be reduced by another factor ∼2 below that given in the previous section, which only considered completeness for $x\gt {x}_{1/2}$. The trade-off against lowering our cut value so as to include a greater fraction of putative sGRB host galaxies, i.e., using, for example, ${x}_{1/3}$ or ${x}_{1/4}$, is two-fold: (i) the number of galaxies would increase rapidly and (ii) incompleteness would become a more severe issue. A further caveat involves the issue of large SN kicks imparted to NS–NS binaries potentially leading to hostless events (Kelley et al. 2010; Behroozi et al. 2014).

5.1. Optical

5.2. X-Ray

5.3. Radio

In this work we have argued for a strategy in which only the brighter galaxies are considered in GW-EM follow-ups. We show that the CLU catalog is fairly complete for $x\gt {x}_{1/2}$ out to 200 Mpc. By weighting the galaxy counts within projected sky areas with a simple model for the aLIGO radial localization, we find that only about 18 ± 5 (for 2015), 23 ± 4 (for 2017), or 11 ± 2 (for 2020) galaxies need to be considered, assuming error boxes of 500, 100 and 20 deg², respectively. This results, if one restricts attention to galaxies in the upper 50th percentile, in integrated luminosity density. Furthermore, there are numerous EM detectors with the property that one tile (i.e., one FOV) would encompass the region of interest. These facilities—optical, X-ray, and radio detectors—have carried out GRB follow-up previously and could participate in future tiling observations. Having one tile per galaxy reduces the number of requisite tiles down to the number of galaxies, a reduction in tiling effort by ∼1–2 orders of magnitude.

The CLU efficiencies for ${\rm{x}}\gt {x}_{1/2}$ vary from ∼100% at the NS–NS inspiral range μ = 60 Mpc to ∼40% at μ = 180 Mpc. For our third target year, 2020, the efficiency is cut by an additional factor of two due to the fact that ∼half of the detections at any given time come from beyond the range, and the CLU is truncated beyond 200 Mpc. In addition to the CLU incompleteness, our adopted cut ${x}_{{\rm{cut}}}={x}_{1/2}$ lies at roughly the median in the observed sGRB host galaxy ${L}_{B}/{{L}_{B}}^{*}$ distribution (Berger 2014), therefore we lose another factor of two in potential sGRB hosts. Lowering our ${x}_{{\rm{cut}}}$ value further would allow us to cover more of the expected sGRB host galaxy ${L}_{B}/{{L}_{B}}^{*}$ distribution but would exacerbate the CLU incompleteness issue. An additional caveat arises from the hostless sGRBs with optical afterglows, and the (majority) of sGRBs for which no afterglows have been observed. If these lie further from the host galaxy centers than the ∼100 kpc maximal offset indicated by Fong and Berger et al. (2013), then our strategy could miss the kilonova emission.

M.M.K. acknowledges generous support from the Carnegie-Princeton Fellowship. S.N. and L.P.S. thank the Aspen Center for Physics and the NSF Grant #1066293 for hospitality during the editing of this paper. S.N. acknowledges generous support from the Radboud University Excellence Initiative. We thank internal LIGO reviewer Ilya Mandel for excellent feedback on all aspects of the paper.