GOING THE DISTANCE: MAPPING HOST GALAXIES OF LIGO AND VIRGO SOURCES IN THREE DIMENSIONS USING LOCAL COSMOGRAPHY AND TARGETED FOLLOW-UP

The Advanced Laser Interferometer Gravitational-wave Observatory (LIGO; Aasi et al. 2015) began operations in 2015 (Abbott et al. 2016b) and almost immediately recorded the first-ever gravitational-wave (GW) signal from a binary black hole (BBH) merger, GW150914 (Abbott et al. 2016d). It should soon observe GWs from neutron star (NS) binary mergers (Abadie et al. 2010; Abbott et al. 2016f) too. These systems should present several kinds of electromagnetic (EM) transients that are detectable by existing and planned facilities (e.g., Metzger & Berger 2012). Joint broadband observations would tell the full story of these rare events and solve longstanding puzzles from the nature of short gamma-ray bursts (SGRBs; Paczynski 1986; Eichler et al. 1989; Narayan et al. 1992; Rezzolla et al. 2011) to the astrophysical sites of r-process nucleosynthesis (Rosswog et al. 2014; van de Voort et al. 2015; etc.) and enable these systems to be used as standard siren probes of the evolution history of the universe (Schutz 1986; Dalal et al. 2006; Nissanke et al. 2010). Metzger & Berger (2012) consider the radioactively powered kilonova (Li & Paczyński 1998) to be the most promising EM signature to find in coincidence with a LIGO event, though Barnes & Kasen (2013) have shown that the high optical opacities of the ejecta will cause this signature to be faint (M_i ≲ −13) and red (r − i ∼ 1) and to peak quickly, within days to weeks (e.g., Kasen et al. 2015). A consortium of partner gamma-ray, X-rays, optical, and radio facilities have embarked on an unprecedented campaign to search for EM counterparts of GW signals (Abadie et al. 2012b; Evans et al. 2012; Aasi et al. 2014; Abbott et al. 2016c).

GW localizations are currently ∼100–1000 deg² (Nissanke et al. 2013; Singer et al. 2014; Berry et al. 2015; Essick et al. 2015) and should shrink to ∼10–100 deg² (Abbott et al. 2016e) over years of detector upgrades and construction of additional detectors: Advanced Virgo (Acernese et al. 2015), KAGRA, and LIGO–India. Even the most accurate imaginable GW localizations of ≲10 deg² will be grossly larger than the ∼1'–10' fields of view (FOVs) of CCD cameras that are common on the world's largest optical and infrared telescopes. Consequently, robotic and low-overhead synoptic survey telescopes with primary mirror diameters of 1–8 m and FOVs of up to tens of deg² have been heralded as the most promising tools (Metzger & Berger 2012) for finding those fast and faint EM counterparts.

Kopparapu et al. (2008) and White et al. (2011) recognized that targeting individual galaxies can reduce the area to be searched while eliminating false positive candidates. A galaxy catalog was central to the follow-up strategy during Initial LIGO (Abadie et al. 2012a; Evans et al. 2012; Aasi et al. 2014). Although the increasing number of galaxies within the expanding GW range decreases the utility of this technique, Hanna et al. (2014) and Kasliwal & Nissanke (2014) argued that using both sky positions and distance estimates from GWs can still reduce the number of galaxies, especially in the early years when localizations are particularly coarse. Nissanke et al. (2013) showed that the region bracketed by the upper and lower limits of the 95% credible distance interval can reduce the volume, and hence the number of galaxies, by a factor of 60%.

While the Advanced LIGO commissioning plan calls for several steps in sensitivity and range, parallel construction of additional detectors will result in shrinking sky localization uncertainty (Veitch et al. 2012; Rodriguez et al. 2014). Intriguingly, Gehrels et al. (2016) point out that these effects roughly cancel each other, such that the typical GW error volume may scarcely vary over the next decade.

The obvious next step is to exploit the correlated three-dimensional (3D) structure of the reconstructed volumes. This immediately raises three questions. How accurately can distance be measured with GW detectors, particularly in the early two-detector configurations? What is the 3D shape of the reconstructed volumes? What is the minimal amount of information that is needed to faithfully describe these volumes?

In Singer et al. (2014), we elucidated the localization areas and shapes that we expect in early Advanced LIGO and Virgo. In a similar spirit and using the same catalog of simulated events, we now reveal the shape, scale, and overall character of the 3D reconstructed volumes—enabled by a new and extremely efficient encoding of the 3D probability distributions. We provide a representative sample of GW volume reconstructions in a format that could be made available beginning with the second Advanced LIGO observing run. Complementing existing technologies to detect (Cannon et al. 2012; Messick et al. 2016) and localize (Singer & Price 2016) GW mergers within minutes of data acquisition, our approach can provide distance-resolved GW sky maps in near real time.

All classes of instruments can benefit from the new 3D localizations, but they are especially powerful for conventional, large-aperture telescopes with narrow-FOV instruments, and particularly at near-infrared wavelengths, where large-FOV cameras are scarce (Emerson et al. 2006). Our galaxy-targeted strategy of monitoring the most probable ∼100 galaxies for several nights following a GW trigger could be implemented on large infrared facilities. However, even a pilot program on robotic 2 m telescopes could be surprisingly powerful for the first few Advanced LIGO–Virgo observing runs.

The amplitude, or signal-to-noise ratio (S/N), of a GW signal is determined by a degenerate combination of inclination and distance.¹⁴ The Malmquist bias leads to a broad universal distribution of binary inclination angles, peaking at 30°(Schutz 2011). The distance of a GW source can generally be estimated with ∼30% fractional uncertainty (Cutler & Flanagan 1994; Berry et al. 2015).

The effective distance r_eff of a GW signal is the maximum distance at which it could have produced the observed S/N. The horizon distance r_H of a detector is the farthest distance at which the most favorably oriented source (at the detector's zenith, and with a binary inclination of ι = 0°) would register a threshold S/N (generally defined as S/N = 8). We give approximate formulae for the horizon distance in Equations (1) and (2) of Singer et al. (2014). The range is the direction- and orientation-averaged distance of sources detectable at a threshold S/N, r_R ≈ r_H/2.26 (Finn & Chernoff 1993; Schutz 2011).

The effective distance description has two major limitations. First, the source may sometimes lie beyond the effective distance because of measurement noise. Worse, there is no obvious way to describe the probability enclosed within, say, the two-dimensional (2D) 90% credible region on the sky and the effective distance; this number is always ≤90%. Second, notwithstanding the large fractional distance uncertainty, there is nontrivial structure to the full 3D reconstructed volumes that can be exploited to reduce the volume under consideration.

During Advanced LIGO's first observing run (O1), the network (Abbott et al. 2016e) consisting of LIGO Hanford Observatory (LHO) and LIGO Livingston Observatory (LLO) tends to produce probability sky maps consisting of one to two long, thin sections of a great circle (Kasliwal & Nissanke 2014; Singer et al. 2014). We provide this as an illustration of the main features for a two-detector network. We assume, as in (Abbott et al. 2016e), a BNS¸ range of 54 Mpc, though the range during O1 was about 40% better. The corresponding 3D geometry is shown in Figures 1 and 2. The degenerate arcs correspond to either one or two thin, rounded, slightly oblique petals, about 1°–5° wide, 10°–100° broad, and 10–100 Mpc deep. The "forked tongue" sky localization features due to the degeneracy of the sign of the binary inclination angle (Singer et al. 2014) are evident as narrow crevices running along the outside edges of the petals. The shape irresistibly suggests a tree ear fungus or a seed of the jacaranda tree.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The O2 configuration (Abbott et al. 2016e), which may include LHO and LLO with improved sensitivity at a BNS range of ∼100 Mpc, as well as Advanced Virgo, leads to more compact and elaborate combinations of petal-shaped regions. In the most favorable three-detector cases where the area on the sky is localized to a single compact region, the reconstructed volume is a spindle a few degrees in radius and ∼100 Mpc long.

Although the reconstructed regions are highly structured, the posterior probability distribution along a given line of sight is simple and generally unimodal; once again, a consequence of the Malmquist bias and the universal distribution of binary inclination angles.

This intuition leads us to suggest that the conditional distribution of distance is well fit by an ansatz whose location parameter $\hat{\mu }({\boldsymbol{n}})$, scale $\hat{\sigma }({\boldsymbol{n}})$, and normalization $\hat{N}({\boldsymbol{n}})$ vary with sky location ${\boldsymbol{n}}$:

$$\begin{eqnarray}p(r| {\boldsymbol{n}}) & = & \displaystyle \frac{\hat{N}({\boldsymbol{n}})}{\sqrt{2\pi }\hat{\sigma }({\boldsymbol{n}})}\exp \left[-\displaystyle \frac{{(r-\hat{\mu }({\boldsymbol{n}}))}^{2}}{2\hat{\sigma }{({\boldsymbol{n}})}^{2}}\right]{r}^{2}\\ \mathrm{for}\ r & \geqslant & 0.\end{eqnarray} \tag{ 1 }$$

This form is equivalent to the product of a Gaussian likelihood and a uniform-in-volume prior. We show that this is a good fit in Section 6 of the Supplement.

The outputs of the LIGO–Virgo localization pipelines are Hierarchical Equal Area isoLatitude Pixelization (HEALPix) all-sky images whose N_pix pixels give the posterior probability ρ_i that the source is contained inside pixel i. We add three additional layers: ${\hat{\mu }}_{i}=\hat{\mu }({{\boldsymbol{n}}}_{i})$, ${\hat{\sigma }}_{i}=\hat{\sigma }({{\boldsymbol{n}}}_{i})$, and (for convenience) ${\hat{N}}_{i}=\hat{N}({{\boldsymbol{n}}}_{i})$. The first layer, ρ_i, is unchanged and still represents the 2D probability sky map.

The probability that a source is within pixel i and at a distance between r and r + $d$r is ρ_i times Equation (1):

$$\begin{eqnarray}&&P(r,{{\boldsymbol{n}}}_{i})\,{dr}={\rho }_{i}\displaystyle \frac{{\hat{N}}_{i}}{\sqrt{2\pi }{\hat{\sigma }}_{i}}\exp \left[-\displaystyle \frac{{(r-{\hat{\mu }}_{i})}^{2}}{2{{\hat{\sigma }}_{i}}^{2}}\right]{r}^{2}\,{dr}.\end{eqnarray} \tag{ 2 }$$

The sky map is normalized¹⁵ such that

$$\begin{eqnarray}&&\displaystyle \sum _{i=0}^{N-1}{\rho }_{i}=1,\end{eqnarray} \tag{ 3 }$$

so Equation (2) is also normalized such that

$$\begin{eqnarray}&&\displaystyle \sum _{i=0}^{N-1}{\int }_{0}^{\infty }P(r,{{\boldsymbol{n}}}_{i})\,{dr}=1.\end{eqnarray} \tag{ 4 }$$

The r² term is necessary in Equations (1), (2) so that the probability density per unit volume vanishes at the origin. Equation (2) should be thought of as the probability distribution in spherical polar coordinates. If, however, one needs to perform a calculation in Cartesian coordinates, one converts using volume element, given by

$$\begin{eqnarray}&&{dV}={r}^{2}\,{dr}\,{\rm{\Delta }}{\rm{\Omega }}=\displaystyle \frac{4\pi }{{N}_{\mathrm{pix}}}{r}^{2}\,{dr}.\end{eqnarray} \tag{ 5 }$$

The r² cancels in the resulting probability density per unit volume:

$$\begin{eqnarray}&&\displaystyle \frac{{dP}}{{dV}}={\rho }_{i}\displaystyle \frac{{N}_{\mathrm{pix}}}{4\pi }\displaystyle \frac{{\hat{N}}_{i}}{\sqrt{2\pi }{\hat{\sigma }}_{i}}\exp \left[-\displaystyle \frac{{(r-{\hat{\mu }}_{i})}^{2}}{2{{\hat{\sigma }}_{i}}^{2}}\right].\end{eqnarray} \tag{ 6 }$$

Sky maps for compact binary merger candidates are produced by two codes with complementary sophistication and speed. The first is BAYESTAR, which rapidly triangulates matched-filter estimates of the times, amplitudes, and phases on arrival at the GW sites (Singer & Price 2016). The second is LALInference, which stochastically samples from sky location, distance, and component masses and spins (Veitch et al. 2015). Both methods directly sample the full 3D posterior probability distribution. The ansatz parameters are extracted using the method of moments as elaborated upon in Section 5 of the Supplement.

We use this encoding to demonstrate the utility of the 3D structure of the GW posteriors. LIGO provided 2D localizations during O1 but did not calculate or distribute low-latency GW distance estimates. (A directional distance estimate for GW151226 produced two weeks after the event (LSC & Virgo 2016; Smartt et al. 2016) did help to rule out the redshift of a Pan-STARRS optical transient candidate.) Without a 3D sky map, one could have provided the horizon distance r_H calculated from the detectors' sensitivity or the effective distance r_eff based on the signal's S/N. The new 3D sky maps allow us to find 90% credible volumes that are 2–30 times smaller (10th to 90th percentile) than the volume within the 2D 90% credible area and the horizon distance, or 1–7 times smaller than the volume within the 2D 90% credible area and the effective distance.

However, the 3D localizations truly shine when we minimize not the volume, but rather the total exposure time required to observe every galaxy within the 90% credible volume to a given flux limit. We neglect intrinsic scatter in absolute magnitude; the resulting conservative figure of merit allows us to focus on the effect of the distance posterior itself.

If BNS mergers have hosts that are similar to short gamma-ray bursts (SGRBs), then their local rates are likely traced by a combination of recent star formation (measured by B-band luminosity) and stellar mass (measured by K-band luminosity; e.g., Leibler & Berger 2010; Fong et al. 2013). As in Gehrels et al. (2016), we attempt to mitigate these concerns as the limited completeness of galaxy catalogs by considering only the brightest galaxies. If we assume a B-band Schechter function with α = −1.25, ϕ* = 1.2 × 10⁻² h³ Mpc⁻³, and L* = 1.2 × 10¹⁰ h⁻² L_⊙ (Longair 2008) as a proxy for the BNS merger rate, then we find that 2.8 × 10⁻³ galaxies per Mpc³ comprise half of the total luminosity (see Figure 3). The areal density of galaxies out to a distance r and with a luminosity greater than L is

$$\begin{eqnarray}{N}_{\mathrm{gal}} & \approx & 0.28\ {\deg }^{-2}\\ & & \times {\left(\displaystyle \frac{r}{100\mathrm{Mpc}}\right)}^{3}\left(\displaystyle \frac{{\phi }^{* }}{4\times {10}^{-3}\,{\mathrm{Mpc}}^{-3}}\right){\rm{\Gamma }}\left(\alpha +1,\displaystyle \frac{L}{{L}^{* }}\right),\end{eqnarray} \tag{ 7 }$$

at most a handful per deg² within distances that are relevant for BNS mergers.

Zoom In Zoom Out Reset image size

As a proof of concept, in the right panel of Figure 2, we show the potential host galaxies that are consistent with a simulated early Advanced LIGO BNS localization. We use the galaxy group catalog of Tully (2015). This may not be an ideal galaxy catalog for GW follow-up; among other reasons, being derived from 2MASS, its magnitudes trace mass rather than recent star formation. Better alternatives that should be available in the near future include the Census of the Local universe (D. Cook et al. 2016, in preparation; Gehrels et al. 2016) and the Galaxy List for the Advanced Detector Era.¹⁶ However, pending the availability of larger compilations of galaxy catalogs, it serves to illustrate our idea because it is ∼50% complete out to ∼150 Mpc. Furthermore, incompleteness is encoded self-consistently by placing a variable bandwidth weighted kernel at the position of each galaxy.

Our control observing strategy takes all galaxies within a given 2D credible region, out to an optimal direction-independent distance inferred from the estimated mass of the source, the loudness of the signal, and the sensitivity of the detectors. It employs the same exposure time for every galaxy. This naive 2 + 1-dimensional (2+1D) construction is already far more sophisticated than the effective distance or horizon distance that could have been available in O1.

Different facilities demand different figures of merit. For FOVs ≪ 0.5 deg², most observations pick out single galaxies. In this case, the optimal selection of galaxies exploits both the upper and lower limits of the reconstructed volume, and the total exposure time is dependent on the number density of galaxies. For FOVs ≫ 0.5 deg², most observations pick out many galaxies; an exposure tuned to reach a fixed luminosity limit for a distant galaxy also captures all of the galaxies at intervening distances. In this case, we tune the exposure time to reach the most distant galaxy in each field, and the total exposure time is independent of the number density of galaxies. The instrument sensitivity is also important. In the source limited regime, the exposure time to reach a fixed limiting luminosity for a source at a distance r varies as r². For sky background limited imaging, the timescales as r⁴. The three cases of greatest interest are summarized below:

$$\begin{eqnarray}T=\mathop{\min }\limits_{D:{\int }_{D}P\,{dV}=0.9}\left\{\begin{array}{lr}\displaystyle \frac{T* }{\omega }{\displaystyle \int }_{D}{\left(\displaystyle \frac{r}{r* }\right)}^{4}d{\rm{\Omega }} & I:\mathrm{large}\,\mathrm{FOV},\\ & \mathrm{sky}\,\mathrm{limited}\\ T* \phi {\displaystyle \int }_{D}{\left(\displaystyle \frac{r}{r* }\right)}^{2}{dV} & {II}:\mathrm{small}\,\mathrm{FOV},\\ & \mathrm{source}\,\mathrm{limited}\\ T* \phi {\displaystyle \int }_{D}{\left(\displaystyle \frac{r}{r* }\right)}^{4}{dV} & {III}:\mathrm{small}\,\mathrm{FOV},\\ & \mathrm{sky}\ \mathrm{limited}.\end{array}\right.\end{eqnarray} \tag{ 8 }$$

Here, T* is the time required to reach a fiducial flux limit for a source at a distance r*, ϕ is the average number density of galaxies, and ω is the FOV of the instrument in steradians. The integral is over solid angle (case I) or volume (cases II and III) and is minimized over all possible regions D that contain 90% posterior probability. Table 1 presents the results of this exercise: volumes, numbers of galaxies, and exposure times. Each case typifies a search for a different EM signature with a different kind of instrument:

Case I describes a kilonova search with a synoptic optical survey instrument. We adopt the Zwicky Transient Facility (ZTF) our example and assume a FOV of 47 deg². Adopting M_i = −13 as the absolute magnitude of a kilonova, we derive from Table 1 of Kasliwal & Nissanke (2014) a fiducial ZTF exposure time of T* ≈ 6600 s for a limiting magnitude of i ≈ 23.5 at a distance of r* = 200 Mpc. In O1, a ZTF-like survey would have taken ∼0.4 hr to tile the entire region to an appropriate depth in exposures that average ∼2 minutes in duration. The number of exposures decreases to a handful in the HLV configuration in Advanced LIGO's second observing run (O2), but the exposure duration increases to 35 minutes per field.

Case II models an afterglow search with the Swift X-ray Telescope (XRT) (considered in detail by Evans et al. 2016). Following Kanner et al. (2012), we adopt a fiducial exposure time of T* = 100 s at r* = 200 Mpc. In O1, we find an average exposure time of 7 s per field, increasing to 25 s per field in O2. These are unrealistically short due to Swift's slew rate (3° s⁻¹) and overhead (25 s; Evans et al. 2016). However, even adding an overhead of ∼30 s per galaxy, we suggest that a galaxy-targeted Swift campaign could be accomplished within a few 90 minute orbits.

Case III describes a kilonova search with a traditional small FOV optical telescope. Since we need to examine tens to hundreds of galaxies, we restrict ourselves to fully robotic telescopes dedicated to time-domain science and select the Las Cumbres Observatory Global Telescope (LCOGT) and Liverpool 2 m telescopes as our models. Using the LCOGT exposure time calculator,¹⁷ we find an exposure time of T* ≈ 2500 s to reach a depth of i = 23.5 mag at S/N = 5 during half moon phase on the Spectral camera. In O1 we find an average exposure time of about 0.2 minutes, for a total exposure time of 0.4 hr. Overhead would dominate. The average exposure time increases to ∼2.5 minutes in O2. The total exposure time increases dramatically from 0.4 to over 13.4 hr due to the rapid increase in exposure time with distance and the modest increase in the number of galaxies. However, the total exposure time is still short enough that two 2 m telescopes working in coordination could conceivably monitor all of the selected galaxies at a 1–2 night cadence.

To focus on the utility of distance and structure information, we neglected several details. Chen & Holz (2014) address optimal selection of which GW events to follow up. We set aside the feasibility of imaging tens to hundreds of targets in rapid succession with Swift, though this is being implemented (Evans et al. 2016). We ignored the question of whether the XRT exposure time can be finely tuned from one target to the next, though our results justify doing so if possible. We ignored variation in observability conditions such as Sun and Earth avoidance, weather, South Atlantic Anomaly passages, and moon phase, details which are better left to a facility-specific paper. We did not consider optimizing field selection given limited time resources, for which we refer the reader to Chan et al. (2015).

We did not address the significant issues of galaxy catalog completeness, construction of a galaxy prior from an incomplete galaxy catalog (Fan et al. 2014), or the suitability of any particular galaxy catalog, though targeting the brightest and most massive galaxies (those with L ≳ L*) should simplify these concerns (Gehrels et al. 2016).

Antolini & Heyl (2016) proposed using photometric redshifts to optimize GW follow-up, improving completeness at the expense of accurate distances. Improved completeness mitigates the concern that BBH mergers like GW150914, which are probably formed in low-metallicity environments (Abbott et al. 2016a), may display host environment preferences similar to long gamma-ray bursts and probably do not track with massive galaxies. However, our present work focuses on NS binary mergers, which are expected to be found in fairly eclectic host environments due to the long time delay between formation and evolution of the binary and its GW-driven inspiral into the LIGO band. Therefore, the completeness of existing spectroscopic redshift surveys seems tolerable for our approach and for BNSs.

One possible concern for small-FOV telescopes is the offsets that mergers may have from their host galaxies due to supernova (SN) kicks. Fong & Berger (2013) find that SGRBs have a median offsets of 4.5 kpc from their hosts. For even an improbably nearby merger at z = 0.005 or a luminosity distance r = 22 Mpc, the projected radius of the 24' XRT FOV is 74 kpc, and of the 10' LCOGT imager, 31 kpc, encapsulating well over 90% of SGRB offsets.

In principal, our exposure time estimates should account not only for sky background, but additional image subtraction background due to the light of the host galaxy. Offsets should help here too because SGRBs are typically found at separations >1.5 times the effective radii of their host galaxies.

One particular advantage of galaxy-targeted searches is the reduction in false positives. In a magnitude-limited snapshot, SNs of Types Ia, Ibc, and II are found in proportions of 68.6%:4.3%:27.1% (Li et al. 2011). Assuming a volumetric rate of 3 × 10⁻⁵ Mpc⁻³ yr⁻¹ (Li et al. 2011), an average absolute magnitude of −19 over a duration of 1 week, and a limiting magnitude of 23.5, we calculate an areal rate of about 2 Type Ia SNs per deg². A typical wide-field GW follow-up campaign searching an area of ∼100 deg² will be contaminated by hundreds of SNs. However, if we consider only 10' × 10' patches around 100 nearby galaxies, the on-sky footprint of <3 deg² translates to a background of merely 6 SNs Ia and 3 core-collapse SNs.

This points toward a new role in GW follow-up for small-FOV, large-aperture telescopes, in addition to and beyond their role of vetting candidates identified by synoptic surveys. Setting aside scheduling and proposal processes for the time being, suitable facilities for our strategy should have primary mirror diameters of 4–10 m and optical or infrared imagers with ∼10' FOVs that are either permanently installed at one of the foci or are rapidly deployable (e.g., mounted on Nasmyth platforms). Candidates include ALFOSC on the Nordic Optical Telescope, LMI on the Discovery Channel Telescope, WIRC on the Hale Telescope at Palomar, FourStar on Magellan, GMOS on Gemini North and South, FLAMINGOS-2 on Gemini South, FORS2 on VLT, LRIS at Keck, or GTC equipped with OSIRIS. As a pathfinder, we encourage deploying existing 2 m class robotic optical telescopes in this manner during the early Advanced LIGO and Virgo observing runs.

We thank the Aspen Center for Physics and NSF grant #1066293 for hospitality during the conception, writing, and editing of this paper. We thank P. Shawhan and F. Tombesi for detailed feedback on the manuscript. Supplementary material, including a sample of reconstructed GW volume FITS files, will be made available at https://dcc.ligo.org/P1500071/public/html. See the Supplement (Singer et al. 2016) in the journal for more details. This is LIGO document P1500071-v7.