A Case Study of On-the-fly Wide-field Radio Imaging Applied to the Gravitational Wave Event GW151226

The era of gravitational wave astronomy has begun. The Advanced Laser Interferometer Gravitational-wave Observatory's (aLIGO; Abbott et al. 2016a) first observing run (O1), which ran in the last quarter of 2015, reported two significant binary black hole (BBH) merger events, GW150914 and GW151226 (Abbott et al. 2016b, 2016c). The second observing run of aLIGO (O2) began in late 2016 and ended on 2017 August 25, with the Virgo detector joining on 2017 August 01 to form a three-detector network. Two more significant BBH mergers were detected, GW170104 and GW170814 (Abbott et al. 2017a; The LIGO Scientific Collaboration et al. 2017), and for the first time gravitational waves were detected from the coalesecence of two neutron stars, GW170817 (Abbott et al. 2017b). The notification of the discovery of BBHs triggered a world-wide, panchromatic search for their electromagnetic counterparts (EM; e.g., Abbott et al. 2016d, 2016e; Copperwheat et al. 2016; Cowperthwaite et al. 2016; Evans et al. 2016; Golkhou et al. 2017; Racusin et al. 2017). Thus far, no conclusive time-variable or quiescent emission has been found at any wavelength for BBHs. In contrast, the binary neutron star merger GW170817 was accompanied by EM signals detected at all wavelengths, including prompt gamma-ray emission (Goldstein et al. 2017), fast-fading optical/NIR (Coulter et al. 2017), and delayed X-ray (Troja et al. 2017) and radio emission (Hallinan et al. 2017; Mooley et al. 2018).

Prior to the detection of the EM counterpart to GW170817, radio emission was widely expected to arise on a wide range of timescales and luminosities from compact binary star mergers. Mergers involving neutron stars leave behind significant neutron-rich debris that settles into a disk. In the conventional picture, most of the debris disk is accreted by the newly formed black hole post-merger (leading to a short GRB in the case of binary neutron stars) and a small amount, about 0.01 M_⊙, is ejected (e.g., Rosswog et al. 1999, 2013; Bauswein et al. 2013; Hotokezaka et al. 2013; Radice et al. 2016). The forward shock into the ISM swept-up by the (sub-relativistic) ejecta is expected to generate broadband synchrotron emission. This gives rise to a milliJansky-level radio transient on timescales of months to years (e.g., Nakar & Piran 2011; Hotokezaka et al. 2016), peaking at frequencies around a few Gigahertz. In cases where relativistic jets are formed and beamed away from the observer, the deceleration of the jet through interaction with the ISM eventually opens up the emission cone into the observer's line of sight. Such orphan afterglows appear as radio transients on timescales of weeks to months (Hotokezaka & Piran 2015; Hotokezaka et al. 2016). This simple jet model likely needs to be modified as the X-ray and radio emission from GW170817 are best understood as the breakout of wide-angle, mildly relativistic outflow (consistent with a "cocoon") of material entrained by the jet (Gottlieb et al. 2017; Kasliwal et al. 2017; Lazzati et al. 2017; Mooley et al. 2018; Ruan et al. 2017).

For neutron star mergers, there is another possible channel for generating radio emission. A millisecond magnetar is a plausible merger remnant, where the magnetar wind drives a strong shock into the ejecta and the reverse shock results in a "plerion" (cf. the Crab Nebula). A strong plerionic radio emission, which is isotropic and independent of the ambient density, is expected on timescales of a few months (Piro & Kulkarni 2013). There is also a strong radio signal expected at late times (∼year timescale; even at small ISM densities) since the magnetar can drive the ejecta to relativistic velocities (e.g., Murase et al. 2017). Past searches for late-time radio emission from short GRBs have put constraints on the phase space of kinetic energy, ejecta mass, and ISM densities, and on the magnetar scenario (Metzger & Bower 2014; Horesh et al. 2016; Fong et al. 2016).

In the case of binary black holes, radio emission is not widely expected to arise from baryonic poor environments. However, if the BBH merger launches a relativistic jet into a dense, gas-rich environment, radio emission is expected at a level of the order of 10–100 μJy (at an assumed distance of ∼400 Mpc), on timescales of 10⁵ s and at frequencies around a few Gigahertz (Yamazaki et al. 2016; Kashiyama et al. 2017).

There are several challenges that must be overcome in any observational effort designed to detect EM counterparts. At radio wavelengths, the main barriers are achieving the necessary submillijasky sensitivity over the large aLIGO error regions of hundreds of deg² at Gigahertz frequencies needed to optimally test existing theoretical models. Equally important is having the ability to rapidly identify and reject any variable sources that could lead to misidentification, and to provide subarcsecond localization for any viable candidates (Metzger & Berger 2012; Metzger et al. 2015; Hotokezaka et al. 2016). In the face of such challenges, several alternative strategies have been adopted. Palliyaguru et al. (2016) and Bhalerao et al. (2017) have taken the approach of carrying out radio follow-up of optically selected candidates, identified from wide-field imaging surveys such as the Dark Energy Camera (DECam; Cowperthwaite et al. 2016) or the Palomar Transient Factory (PTF; Law et al. 2009). This approach avoids imaging large sky areas and it can make deep observations over a wide frequency range. Another approach has been to shift to lower frequencies in order to utilize existing wide-field instruments. Radio searches spanning 50–100 deg² were carried out with the Australian Square Kilometre Array Pathfinder (ASKAP; Bannister et al. 2016), the Low-Frequency Array (LOFAR; Broderick et al. 2016a, 2016b), and the Murchison Widefield Array (MWA; Kaplan et al. 2016) telescopes at frequencies of 863 MHz, 145 MHz, and 118/154 MHz, respectively. While these surveys did not have the sensitivity and angular resolution needed to detect and provide arcsecond-localization of any putative radio afterglows of the GW events, they were good proof-of-concept experiments for the follow-up of aLIGO sources using wide-field blind radio observations.

In this paper, we take a more direct approach by using the VLA On-the-Fly mosaicing capability to make sensitive, wide-field observations at arcsecond resolution and at Gigahertz frequencies. We use this new mode to carry out follow-up observations of GW151226. In Section 2, we describe the VLA observations and the data reduction carried over a 100 deg² region. The search for transients and the identification of variables is described in Section 3. This search is discussed in the context of future aLIGO-Virgo observing runs and radio follow-up programs in Section 4.

2.1. The Gravitational Wave Event GW151226

The GW151226 gravitational wave signal was initially identified by the GstLAL compact binary coalescence search (Messick et al. 2017) of the data from the LIGO Hanford and Livingston detectors at 2015 December 26.15, and localized by the BAYESTAR code (Singer & Price 2016), which is sensitive to compact binary star coalescence events. The false alarm rate for this event was reported as being lower than one per hundred years (Usman et al. 2016; Messick et al. 2017). Using Bayesian Markov-chain Monte Carlo and nested sampling to perform forward modeling of the full GW signal, including spin precession and regression of systematic calibration errors (LALInference; Veitch et al. 2015), the event was localized to within 1240 deg² (90% credible region), significantly improving over the BAYESTAR localization. The LALInference sky map of the gravitational wave event, together with our VLA survey region, is shown in the left panel of Figure 1.

Zoom In Zoom Out Reset image size

GW151226 marks the second direct detection of gravitational waves. Detailed offline analysis of the aLIGO data recovered the gravitational wave signal with a significance greater than 5σ. The initial (individual) and final black hole masses were estimated to lie between 5–22 M_⊙ and 19–27 M_⊙, respectively. The luminosity distance was estimated to lie between 250 and 620 Mpc, i.e., 0.05 < z < 0.12 (the ranges represent the 90% credible interval; Abbott et al. 2016f).

2.2. On-the-fly Mosaicing

2.3. Radio Observations

2.4. RFI Flagging, Calibration, Imaging, and Source Finding

Immediately after the observations for each epoch/region (each observing block) were complete, we downloaded the raw data from the VLA archive, onto the lustre file system at the NRAO AOC in Socorro. The data were then calibrated using the NRAO scripted pipeline (in CASA 4.5.0) with quacking¹¹ of the target fields removed. Due to substantial satellite-induced RFI, the data between 2.12 and 2.37 GHz (the second and third spectral windows) were fully flagged for all observations. Single scans (5500–5700 scans for each region per epoch) were then split from the calibrated data set into individual measurement sets. The calibration and splitting processes for each observing block took about 36 hr.

We then ran a distributed imaging process on 70 cores across five computer nodes at the NRAO AOC using CASA task clean, such that each core was imaging a single scan at any given time. To enable linear mosaicing, we convolved the images with fixed circular synthesized beams¹² of 8 arcsec (for the C config observations) or 3 arcsec (for the B config observations). The imaging of each scan was itself done over two iterations. In the first iteration, we cleaned using Briggs weighting and a robust parameter of zero in order to minimize imaging artifacts. Images with sources brighter than 20 mJy were then used to derive and apply a single phase self-calibration procedure. Clean boxes were then derived based on the sources found using the ia.findsources tool in CASA. This was followed by clipping of the visibilities in amplitude (in all spectral windows; in order to minimize the effect of RFI) at a level equal to the sum of the mean and three times the standard deviation of the visibility amplitudes in the spectral window with the least RFI. The second imaging iteration used the self-calibrated visibilities to image using natural weighting and clean boxes. During each cleaning step, we used 700 clean iterations, two Taylor terms, a cyclefactor parameter of 7, and a clean stopping threshold of 0.3 mJy. The pixel size was chosen so as to sample the synthesized beam across at least four pixels. The images had a center frequency of 3.0 GHz, except for the observation on 2017 June 27 (where the data had substantial RFI between 2 and 3 GHz) for which the center frequency is 3.3 GHz.¹³ A small fraction of the single-scan images had strong spike-shaped artifacts (there is a known bug in CASA clean when using two or more Taylor terms, which is currently being fixed), and for these scans we chose the Briggs-weighted images (output during the first imaging iteration described above) instead of the ones with two Taylor terms. Based on our inspection of the single-scan images, we expect that the Briggs weighting reduces the integrated flux density of extended sources in the survey region by a fraction less than or equal to ∼10%. For unresolved sources, there will be no change in flux density due to the introduction of the Briggs-weighted images.

Linear mosaicing of the single-scan images was then carried out using FLATN in AIPS, followed by cropping of the mosaic into 4096 × 4096 pix² subimages. We used the primary beam parameters from Perley et al. (EVLA Memo 195) during the linear mosaicing step. The imaging and linear mosaicing processes together required about 6 hr per region/epoch. We then made a background noise map for each subimage using RMSD in AIPS, which was then supplied to SAD for the cataloging of sources down to¹⁴ 6σ. The cumulative noise plots made from the background rms noise maps is shown in Figure 2.

Zoom In Zoom Out Reset image size

As we are interested in only unresolved and partially resolved sources, we chose only those sources from each epoch that had an integrated-to-peak flux density ratio of <1.5 and prepared a merged catalog, consisting of 5307 sources. For sources in the merged catalog that were not detected in some of the epochs, we set the peak flux density to be equal to the peak pixel value at the location of the source. For those sources that have integrated flux density lower than the peak flux density, we set the integrated flux density to be equal to the latter. This cataloging process required about 1 hr. The total time required from downloading the raw data for each observing block and producing a final source catalog was 43 hr.

We used the catalog of 5307 sources from Section 2.4 to carry out a search for transient sources that appeared or disappeared in one or more of the three epochs. No transients were found to a 6σ limit of 670 μJy (50% completeness threshold for the merged catalog over three epochs and 100 deg²).

The same catalog was used to investigate variability following a process described in more detail in Mooley et al. (2016). In short, for each source, we calculate the variability statistic, V_s = ΔS/σ for each epoch and a modulation index $m={\rm{\Delta }}S/\bar{S}$, where S is the flux density, $\bar{S}$ is the average flux density, ΔS is the flux density difference, and σ is the rms noise. Significant variables were identified as those sources that have V_s larger than 4σ, and the absolute value of the modulation index, $| m|$, larger than 0.26 (i.e., a fractional variability, f_var > 0.3). The constraint on V_s ensures fewer than one false positive will be detected as a variable source in our search, assuming Gaussian statistics. Considering that our flux scale is accurate to only 3%–5% (Thyagarajan et al. 2011; Mooley et al. 2013), the artificial variability induced on account of our imaging parameters (∼10%; see Section 2.4 and also Mooley et al. 2016), and our usage of the true primary beam instead of the smeared OTFM beam, the constraint on the fractional variability is needed to minimize the number of false positives.

We used epoch E1 as the reference epoch and compared E2 and E3 independently with this epoch using the V_s and m for finding significant variable sources. Based on our inspection of the ratios of source flux densities between E2 and E1, we had to apply a multiplicative correction factor of 1.07 to the source flux densities in E2. Similarly, our comparison of E3 versus E1 instructed us to apply a multiplicative correction factor¹⁵ of 1.30 to the source flux densities in E3. This correction factor made the distribution of the variability statistic V_s symmetrical, but also broadened it to a Gaussian-like distribution (σ = 2.5, hence our 4σ variability selection criterion in this case is V_s > 10) rather than a Student-t distribution (see Figure 9 of Mooley et al. 2016). During the variability search between E3 and E1, we also discovered that a significant fraction of the sources in the final source catalog were resolved out in E3, and an abnormally large number of variables appeared below ∼1 mJy (this can again be attributed largely to the angular resolution differences). Hence, for the comparison between E3 and E1 we restricted our search to only those sources (total of 2782) having an integrated-to-peak flux density ratio¹⁶ of <1.2.

The plot of the variability statistic versus the modulation index is shown in Figure 3. We found 72 significant variable sources between E1 and E2 and 42 variables between E1 and E3, having fractional variability larger than 30%. This indicates a variability fraction of 1.5 ± 0.2% for timescales between a few weeks and a few months.

Zoom In Zoom Out Reset image size

We have used the VLA to image a 100 deg² error region of the O1 aLIGO BBH merger event GW151226. While a bright EM counterpart was not expected in this case, we have demonstrated a near real time ability to conduct blind searches in a phase space where the predicted afterglow signature is strongest. As we noted in Section 1, a late-time radio afterglow peaking on timescales of hundreds of days is one of the more robustly predicted afterglow signatures from neutron star mergers (Nakar & Piran 2011). Such searches are best carried out at Gigahertz frequencies, since synchrotron self-absorption suppresses the signal below 1 GHz for up to several years post-merger. Estimates of the radio afterglow signal currently suffer from uncertainties in the circum-merger gas density, and a full search must be sensitive to the low densities (10⁻³ cm⁻³), as seen for some short gamma-ray bursts (Fong et al. 2014), and the density within galactic disks, ∼1 cm⁻³ (Draine 2011), as expected for the Galactic binary pulsar population.

Metzger & Berger (2012) have listed the "cardinal virtues" that must be met for an experiment to be considered a competitive follow-up effort for detecting the predicted signatures of compact binary star mergers. These are (1) that they are detectable with current observing facilities with a reasonable time allocation, (2) that they accompany a significant fraction of GW events, (3) that they can be distinguished from other astrophysical transients, and (4) that they can provide arcsecond sky localization. Below we demonstrate how this experiment addresses the challenges and meets almost all of these requirements.

Since the field of view of the VLA at 3 GHz is small (0.06 deg²), we take advantage of the superb instantaneous sensitivity of the VLA and use the OTFM method, rapidly slewing the antennas (16/minute) in order to image a 100 deg² error region with only 7.5 hr of on-sky integration time. These observations were made at three epochs approximately 50, 101, and 185 days post-merger (see Table 1). These are an appropriate range of timescales for searching for the rise and fall of the expected late-time radio afterglow. The typical rms noise achieved for each epoch was approximately 120 μJy/beam. This is an astrophysically interesting sensitivity level for detecting afterglows circum-merger gas densities >0.1 cm⁻³ and distances of 100–200 Mpc (Hotokezaka et al. 2016).

We searched for any transient sources that appeared in the 100 deg² region over the three epochs. None were found. This implies an upper limit to the areal density of transients brighter than 700 μJy of 10⁻² deg⁻². Our limit is similar to that obtained from the 3 GHz Stripe 82 Pilot project (Mooley et al. 2016). The types and timescales of the different radio transients at frequencies of a few Gigahertz are shown in Figure 4. On the timescales that were sampled in this experiment, the expected transients include tidal disruption events (TDE), orphan GRB afterglows, and Type II core-collapse supernovae. If the gravitational wave event is contained within our survey region, we can place an upper limit of about 10²⁹ erg s⁻¹ Hz⁻¹ to any late-time radio afterglow from the merger event at an assumed mean distance of 440 Mpc. At an areal density of 10⁻² deg⁻², none of these transient source classes are expected to be a significant contaminant in a blind radio survey for the late-time afterglow from a binary compact remnant merger Mooley et al. (2016). The quietness of the radio sky stands in contrast to optical counterpart searches (e.g., Copperwheat et al. 2016; Cowperthwaite et al. 2016; Palliyaguru et al. 2016), for which there were a number of transients that required spectroscopic follow-up before being eliminated as candidates. Radio observations of GW sources also have an advantage over the optical and X-ray in cases where solar/lunar observing constraints are present, and dust-obscured environments in the host galaxy itself, or along lines of sight through our own Galaxy. For example, continued radio monitoring of GW170817 was carried out during a crucial period in the evolution of the afterglow light curve, while optical and X-ray telescopes were constrained by the Sun, thus allowing us to distinguish between the off-axis jet versus cocoon models (Mooley et al. 2018).

Zoom In Zoom Out Reset image size

While no transients were found, about 1.5% of the compact, persistent radio sources in the field exhibited significant variability. Most radio sources vary to some degree or another. Below ∼1 GHz, the variations are dominated by plasma effects such as interstellar scintillation, while at high frequencies the variations are intrinsic, likely originating from changes in the black hole–disk–jet environments (Altschuler 1989). The frequency and level of variability is in agreement with previous Gigahertz studies (Mooley et al. 2013, 2016; Hancock et al. 2016) and is consistent with normal AGN activity¹⁷ on a timescale of a few months. In principle, a transient occurring in a star-forming galaxy could be mistakenly identified as an unresolved or partially resolved variable source. In practice, with the 8 arcsecond synthesized beam and the rms noise of ∼120 μJy/beam, the experiment does not have sufficient surface brightness sensitivity to detect star-forming galaxies within the LIGO volume (Condon 2015) for GW151226. However, host-galaxy contamination may be severe in the cases where merger afterglows are located within ∼100 Mpc and observations are carried out in VLA D config (see Figure 5). Hotokezaka et al. (2016) looked at contamination as a function of galaxy distance for VLA B config and L band (1.4 GHz), and they find that contamination due to starburst emission and bright AGNs will be ∼5%. Nevertheless, an interferometer such as the VLA not only provides a subarcsecond position for a potential merger event, but it also acts as a spatial filter of extended host-galaxy emission.

Zoom In Zoom Out Reset image size

While we have shown that Gigahertz radio searches for late afterglows are a promising method to search for the EM signature of a GW event, the experiment could still be improved. In Section 2.4, we noted that 1.8 days were required to process each epoch. Lower latency would be possible if the data calibration stage was further optimized. The benefits of such improvements are not immediately obvious. The radio afterglow reaches its peak months after the gravitational event, well after any putative optical signature has faded away. The follow-up of any radio candidate would be with a goal of identifying a host galaxy and measuring its redshift, and as such could be done at a more leisurely pace.

A weakness of radio searches is that the large uncertainty in the expected circum-merger densities predicts a wide range of peak flux densities at Gigahertz frequencies (see Section 1). Figure 6 shows the cumulative distributions of peak flux densities of neutron star merger (simulated population of gravitational wave sources; Hotokezaka et al. 2016). It is possible to achieve noise values of 50–75 μJy/beam and thus broaden the parameter space of the late-time afterglow models. Such deeper searches will be able to access up to 40% of the merger afterglows. For this blind survey approach, and with finite VLA telescope time, the trade off is to image smaller fields of view. Such decreases in the size of the GW error regions are expected in the upcoming runs jointly with LIGO and Virgo, and eventually runs including KAGRA and LIGO-India (Nissanke et al. 2013; Kasliwal & Nissanke 2014; Abbott et al. 2016g). Another approach for increasing sensitivity is to significantly reduce the total area searched by doing targeted searches. This approach involves deep radio follow-up of transient candidates identified at other wavelengths (Palliyaguru et al. 2016; Bhalerao et al. 2017) or using a master galaxy catalog in the localized volume (Gehrels et al. 2016; Singer et al. 2016). Such a strategy has been successful in the case of GW170817. However, as the distance range for the gravitational wave detectors increases, blind surveys or even hybrid approaches between targeted and blind searches will become important.

Zoom In Zoom Out Reset image size

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. K.P.M. acknowledges support from the Oxford Centre for Astrophysical Surveys, which is funded through the Hintze Family Charitable Foundation. K.P.M. is currently a Jansky Fellow of the National Radio Astronomy Observatory. We thank the anonymous referee for valuable comments that helped improve this manuscript.