Searches after Gravitational Waves Using ARizona Observatories (SAGUARO): System Overview and First Results from Advanced LIGO/Virgo's Third Observing Run

The SAGUARO software suite ingests the VOEvents distributed by the NASA Gamma-Ray Coordinates Network (GCN)³⁰ system in real time. These alerts are employed by the LVC for disseminating GW event information, including the HEALPIX localization maps with distance constraints (e.g., Singer et al. 2016; Singer & Price 2016). The alert contents are described in the LIGO/Virgo Public Alerts User Guide³¹ and contain the classification probabilities for each GW event, which are split into five categories (see Kapadia et al. 2019, for more details): (1) terrestrial (not of astrophysical origin), (2) BNS (both components are neutron stars, 1 < M < 3 M_⊙), (3) MassGap (any component has a mass in the gap between neutron stars and black holes, 3 < M < 5 M_⊙), (4) NSBH (one component is a neutron star, 1 < M < 3 M_⊙ and one is a black hole M > 5 M_⊙), and (5) BBH (both components are black holes, M > 5 M_⊙). The alerts also contain two parameters, HasNS and HasRemnant, that indicate the likelihood that the event produces EM emission. HasNS indicates the probability that one of the components was a neutron star and HasRemnant indicates the probability that some material remained outside the final remnant compact object, as calculated by the Foucart et al. (2018) model.

Once an alert is received, our CSS search is automatically triggered if the classification parameters exceed our limits for triggering, the false alarm rate is better than our requirement, and the target is observable by CSS. We require that the classification for the GW event have a combined BNS, NSBH, and MassGap probability greater than 20% and a false alarm rate <12 yr⁻¹ in order to trigger. The GW event is determined to be observable if any of the CSS fields within the 90% probability region meet the following constraints: (1) airmass <2.5, (2) the projected distance on the sky to the Moon, d_moon, is governed by ${d}_{\mathrm{moon}}\gt (42\times {\theta }_{\mathrm{illum}}+3)$ deg, such that the d_moon limit increases with the moon illumination θ_illum (represented by a fractional number going from 0 for new to 1 for full), (3) Sun altitude <−12 deg. The exception to this was the first event of the LIGO/Virgo O3 run, S190408an, which was used as a full system test even though it was a clear BBH (see Section 5).

Once it is determined that a GW event meets our criteria for triggering, the SAGUARO software automatically inserts a selection of up to 12 fields (60 deg²) into the CSS observing queue with the option of manually triggering an additional 12 fields. These fields are selected to cover the highest-probability regions that are observable and for which we have template images. They are given a higher priority level that allows them to be observed immediately after the current CSS sequence is finished.

Once all four images of a field have been taken, basic processing including bias subtraction and flat-fielding are performed. Astrometric and photometric calibrations are then done with SCAMP (Bertin 2006, 2010a) using Gaia DR2, resulting in a typical standard deviation of ∼011 and ∼0.19 mag for the astrometry and zero-point, respectively.

A script watching for new data waits for all four images of a field to arrive before creating a median for that field using SWarp (Bertin 2010b). During the median creation, each image is background-subtracted to remove background variations before median combining. The background value is retained and used for scaling purposes on the median. The creation of medians allows for cleaner images on which we can search for transients, and removes artifacts such as cosmic rays that only appear on a single image. As an image set for a particular field can be separated by ∼30 minutes, there exists the possibility that not all four images for each field will be observed, due to weather, moon, or Sun constraints—we therefore have built in a time limit for our median creation algorithm to create medians with the available images if all four are not received. A real-time data-processing pipeline, described in more detail in Section 3.4, then processes these median images for transient detection.

One of the primary advantages of SAGUARO's transient discovery is the access to nearly three years of archival data with the current instrumentation, which enable the production of deep templates for image subtraction. We created templates that cover the entire CSS footprint of 25345 deg² (see Figure 1) with a median of ∼60 individual 30 s images contributing to each template after a series of quality cuts. First, we discard images that have <2000 detected point sources. This is done to ensure high astrometric precision for the templates. Second, we do not include images that have a sky brightness of <20 mag arcsec⁻². This cut was chosen based on the distribution of measured sky brightness values, and serves to exclude poor-quality data taken in adverse conditions (i.e., clouds or bright moon). Applying these cuts restricts the template creation to high-quality, deep images with a median 3σ limiting mag of 23.0 mag. Approximately 27% of the archival data were not included in the template creation. At the two extremes, 5% of the templates were created with more than 90 images per field, while 9% of the templates were only observed once and have ≤four images per field. For those fields with poor coverage (e.g., at very high decl. or close to the Galactic plane), no quality cuts were applied and all available images were used to create the template.

SAGUARO makes use of a data-processing pipeline written in Python (K. Paterson et al. 2019, in preparation will discuss the pipeline in detail) and an implementation of the image subtraction method called ZOGY³² (Zackay et al. 2016).

The data-processing pipeline controls the flow of the data, submitting images for image subtraction and creating detailed logs for each image. The pipeline also creates a mask for each image by identifying saturated stars and bad pixels. This allows regions with poor subtractions caused by bad pixels and saturated stars to be ignored during the transient detection.

Image subtraction is performed using the new median images and the templates discussed in Section 3.3. The difference image created by the image subtraction is then converted into a significance image, where the value of each pixel is represented by its significance (i.e., a pixel with value 5 will have a signal-to-noise ratio (S/N) of 5). Corrections, such as those for astrometric errors, are applied to the significance image to produce a corrected significance or "Scorr" image (see Zackay et al. 2016 for a detailed description on the products produced). As the significance or S/N is encoded directly in the Scorr image, it provides a direct way to find detections above a set threshold. For SAGUARO, we set the detection threshold on the subtracted images to 5σ. Thus, sources with a significance or S/N > 5 on the Scorr image will be flagged as transient candidates. A flux for each detection is obtained through point-spread function (PSF) photometry. The zero-point, calculated through SCAMP (Bertin 2006, 2010a) using Gaia DR2, is then used to convert this flux to magnitudes in Gaia G-band (roughly covering 330–1050 nm; see Weiler 2018). With the focus on transients, we cross-match detections against stellar sources in the templates to remove variable stars and poor subtractions associated with stars close to saturation. At present, SExtractor's CLASS_STAR parameter is used to filter stellar sources (defined as CLASS_STAR > 0.5). The remaining candidates are then loaded into the database for visualization and vetting.

For vetting in real-time, detailed information for each transient candidate is stored in a PostgreSQL database and postage stamp images are saved to disk. A Flask webserver allows visual inspection of candidates as they come in, and accommodates queries based on date, field ID and detection threshold. Candidates are sortable by S/N or machine-learning score (which gives the likelihood of the transients being real) to promote the most likely candidates.

Each candidate is automatically cross-matched against known moving objects from the Minor Planet Center (MPC³³ ) and known transients from the Transient Name Server (TNS³⁴ ). We also search for previous detections from the Zwicky Transient Facility (ZTF; Graham et al. 2019) and cross-match against galaxies from the Galaxy List for the Advanced Detector Era (GLADE; Dálya et al. 2018) catalog within the localization volume using a simple broker, the Steward Alerts for Science System (SASSy³⁵ ). This allows for lists of viable candidates to be quickly disseminated to the community through GCN notices within a few hours of taking the data.

As we show in Section 5, the EM counterpart search portion of SAGUARO is functioning and responding to LVC events in O3. We have also been ingesting all 1.5 m Mt. Lemmon telescope data during routine CSS operations in real time and searching for transients (e.g., Lundquist et al. 2019), independent of the main GW follow-up program. At the moment, the primary purpose of this wide-area search is to stress test our GW search pipeline, and gather a sufficient training set to improve our machine-learning algorithm for transient vetting.

SAGUARO's wide-field CSS search reaches a limiting magnitude of G ≈ 21.3 mag. In Figure 2, we show the optical r-band and NIR K-band light curves for AT2017gfo, the optical counterpart of GW170817, (compiled in Villar et al. 2017; data from Andreoni et al. 2017; Arcavi et al. 2017a; Coulter et al. 2017; Cowperthwaite et al. 2017; Díaz et al. 2017; Drout et al. 2017; Kasliwal et al. 2017; Lipunov et al. 2017; Pian et al. 2017; Pozanenko et al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Troja et al. 2017; Utsumi et al. 2017; Valenti et al. 2017) for three distances: ≈40.7 Mpc, the distance of GW170817 (Cantiello et al. 2018), 100 Mpc, the typical expected distance for a BNS merger for LVC's O3, and 200 Mpc, the expected range for a face-on BNS merger at design sensitivity (Abbott et al. 2016a). Matched to SAGUARO's current discovery limits, our optical search maintains sensitivity to kilonovae of comparable luminosity to GW170817 out to 200 Mpc at δt ≲ 1 day. Once a promising kilonova candidate is identified, SAGUARO's follow-up efforts on larger aperture facilities (described in subsequent sections) will be employed to track the counterpart's temporal evolution. Indeed, it is the combined discovery and follow-up endeavor, the template coverage, and the extension to NIR wavelengths, that give SAGUARO an advantage in GW counterpart efforts.

Zoom In Zoom Out Reset image size

In the Northern Hemisphere, there are only a few optical >1 deg² wide-field efforts underway to search for optical counterparts to GW events. This includes efforts by Pan-STARRS (Chambers et al. 2016), ZTF (Graham et al. 2019), and ATLAS (Tonry et al. 2018). All three facilities are conducting optical counterpart searches following GW events in addition to their primary surveys. Among these searches, the SAGUARO discovery effort is most comparable to Pan-STARRs, which has a 1.8 m aperture and 7 deg² FOV. The discovery components of ZTF and ATLAS are ≈1–2 mag less sensitive, but both possess a FOV that is a factor of ≈6–7 times larger than SAGUARO/CSS.

When GW localizations cover ∼100–1000 s deg², potentially dozens of candidate counterparts will be uncovered, as was the case for the first two potential neutron star mergers of O3, S190425z, and S190426c (see, e.g., the summary of these events in Hosseinzadeh et al. 2019). While this large number of candidates can be winnowed down by making cuts on the transients' age, color, and association with galaxies at the appropriate distance, ultimately spectroscopic vetting must be done to uncover the true EM counterpart.

SAGUARO is committed to spectroscopic candidate vetting for very promising targets that are confirmed to have no detection before the GW event. Candidates from all programs, including our own Mt. Lemmon 1.5 m search, will be considered equally. When distance information is available, we prioritize targets that are localized near cataloged galaxies (e.g., Dálya et al. 2018) with distances consistent with that inferred from the GW signal. Candidate vetting is a dynamic process that benefits from real-time access to large aperture telescopes, as was the case for our team's vetting of PS19qp/SN2019ebq (Smith et al. 2019), which an initial spectrum suggested was consistent with a kilonova at the distance to S190425z (Nicholl et al. 2019a); we detail these observations further in Section 5.2.

Once a true GW counterpart is discovered, SAGUARO will spring into action to collect high-cadence optical and NIR light curves as well as spectral sequences while the transient is accessible to ground-based observatories.

First, in the small to medium aperture range are the 1.5 m Kuiper, 1.8 m VATT and 2.3 m Bok telescopes, all of which are based in Southern Arizona. All three telescopes have imagers that can gather data on kilonovae, VATT and Bok have spectrographs, and Kuiper and Bok have access to SPOL, an imager/spectropolarimeter (Schmidt et al. 1992). For instance, in ∼1 hr exposure times the Bok B&C spectrograph can get high S/N spectra down to ∼19 mag; GW170817/AT2017gfo was brighter than 19 mag for ∼2 days. Similarly, the 90Prime 1 deg² camera (Williams et al. 2004) on the Bok telescope can image down to ∼24 mag in ∼30 minute exposures, facilitating kilonova optical light-curve follow-up for a week or more in the nearest events (Figure 2).

SAGUARO also has access to large aperture facilities: the 6.5 m MMT, the twin 6.5 m Magellan telescopes, the 2 × 8.4 m Large Binocular Telescope, and the two Keck 10 m telescopes. The spectral sequences that these facilities can provide will lend insight into the emission mechanisms of neutron star mergers and r-process element production. Additionally, all of these large aperture facilities will be used to obtain late-time light curves once any counterpart is too faint for spectroscopy, potentially down to ∼26.5 mag in the optical and ∼24.5 mag in the NIR. These same large aperture facilities will be used to study the host properties of the EM counterparts to GW events, as was briefly discussed in Section 2.

A candidate GW signal was identified using data from LIGO Livingston Observatory (L1), LIGO Hanford Observatory (L2), and Virgo Observatory (V1) on 2019 April 8 at 18:18:02.288 UTC (Ligo Scientific Collaboration & VIRGO Collaboration 2019a). This event, S190408an, had a >99% probability of being a BBH merger. The localization covered 387 deg² at a distance of 1473 ± 358 Mpc. S190408an is the first publicly announced astrophysical event in LVC's O3 and the entire localization region resides in the Northern Hemisphere (Figure 3), making it a prime target for SAGUARO's first trigger as an important test of the system.

The majority of the localization lies in an area not covered by the CSS survey due to its proximity to the Galactic plane. As a result, we only triggered three viable fields with template images. The localization was not immediately observable and the images were taken at a mid-time of δt ≈ 17.56 hr after the GW event. As shown in Figure 3, these fields covered 8.63 deg² of the 50% probability region and 13.45 deg² of the 90% probability region. These observations account for 13.8% of the 50% probability, 8.6% of the 90% probability, and 7.8% of the total probability. Poor weather conditions and the low elevation of the fields resulted in a 3σ limiting magnitude of 19.8 mag. The details of these observations are summarized in Table 1.

From these 3 fields, 5469 candidates above 5σ were detected. No known moving objects or transients were found in the median images after cross-matching these candidates against known moving objects from the MPC and known transients from TNS. Cross-matching to Toba et al. (2014) and Krawczyk et al. (2013) for known active galactic nuclei (AGNs) also found no matches. As generally expected for BBH mergers, no real astrophysical source associated with the GW event was found after candidate vetting.

A candidate GW signal was identified using data from L1 and V1 on 2019 April 25 at 08:18:05.017 UTC (Ligo Scientific Collaboration & VIRGO Collaboration 2019b). This event, S190425z, had a >99% probability of resulting from a BNS merger. The initial 90% localization covered 10183 deg² at a distance of 155 ± 45 Mpc. The 90% localization area was later improved to 7461 deg² at a distance of 156 ± 41 Mpc (Ligo Scientific Collaboration & VIRGO Collaboration 2019c). It is notable that S190425z was detected as a sub-threshold event in V1, and thus is considered a single detector event, explaining its large localization area with respect to the expected median value for O3.

Once the initial alert was received, the SAGUARO software automatically selected the 12 highest-probability fields to observe and placed these into the CSS queue. As Figure 4 shows, the localization has three distinct regions of high probability. Two of these regions are mostly visible from the southern hemisphere and one region is mostly visible in the north; the 12 selected fields were naturally in the northern region. These fields covered 58.5 deg² of the 50% probability region and 60.0 deg² of the 90% probability region. These fields accounted for 6.0% of the 50% probability, 3.4% of the 90% probability, and 3.0% of the total probability after the localization was updated (Table 1).

Zoom In Zoom Out Reset image size

The observations of the fields started at δt ≈ 1.1 hr with respect to the GW event, ∼0.4 hr after it was added to the CSS observing queue, with a mid-time of δt ≈ 1.37 hr for all 12 fields. In favorable observing conditions, we reached a 3σ limiting magnitude of 21.3 mag. From the 12 fields observed, 2711 candidates above 5σ were detected. After cross-matching these candidates against known moving objects from the MPC, two known moving objects were found to have remained in the median-combined images. After cross-matching to the TNS, six previously discovered transients were found within the data: AT 2019bln (Tonry et al. 2019), SN 2017frl (Tonry et al. 2017; Xhakaj et al. 2017), AT 2019eaj (Fremling 2019), AT 2018cix (Fremling 2018), SN 2019aja (Stanek 2019), and SN 2019bzo (Brinnel et al. 2019; Nordin et al. 2019). Figure 5 shows an example of one such transient found. Further cross-matching to Toba et al. 2014 and Krawczyk et al. 2013 found 45 known AGNs associated with detections. Four candidates were found within the data: SN 2019eff (Chambers et al. 2019a), AT 2019efu (Chambers et al. 2019a), AT 2019ech (Chambers et al. 2019b), and AT 2019fgy (Chambers et al. 2019c). Only one candidate, SN 2019eff, was associated with a GLADE galaxy, but at a distance inconsistent with the reported range for the GW event. Spectral classification of SN 2019eff by Nicholl et al. (2019b) indicated a type IIb supernova. The remaining candidates were not followed up by any group.

Zoom In Zoom Out Reset image size

As the follow-up of S190425z developed, 69 GW counterpart candidates were reported in GCNs (Hosseinzadeh et al. 2019). One GW counterpart candidate reported by Pan-STARRS, PS19qp/SN2019ebq (Smith et al. 2019), appeared to be a promising kilonova candidate. An initial spectrum from the advanced Public ESO Spectroscopic Survey for Transient Objects (PESSTO) displayed a red, featureless continuum and narrow host lines at a redshift z = 0.037 (Nicholl et al. 2019a). This redshift corresponds to a distance consistent with the GW event (D ≈ 150 Mpc), and gave the transient an absolute magnitude of M_i ≈ −16.7 mag (corrected for Galactic extinction Schlafly & Finkbeiner 2011), roughly that of GW170817 at similar epochs. This candidate was not covered by the CSS fields triggered, but we were motivated by the distance, luminosity, and featureless spectrum seen by PESSTO, to trigger ToO imaging and spectroscopy on the LBT starting at δt ≈ 1.11 days after the GW trigger. The transient is well detected in the 60 s r-band acquisition image (Figure 6). Performing astrometry relative to Gaia DR2, we measure a location of R.A. = 17^h01^m1835, decl. = −07°00'10.5'' with a positional uncertainty of 0060 (1σ). Using standard tasks in IRAF, we measure a brightness for PS19qp/SN2019ebq of r = 20.44 ± 0.05 mag, translating to M_r ≈ −15.59 mag at z = 0.037 (or M_r ≈ −16.91 when accounting for the Galactic extinction in the direction of PS19qp/SN2019ebq of A_r = 1.313 mag; Schlafly & Finkbeiner 2011).

Zoom In Zoom Out Reset image size

We obtained 2 × 600 s exposures of PS19qp/SN2019ebq at 2019 April 26 11:00 UT using the Multi-Object Double Spectrographs (MODS; Pogge et al. 2010) and a 12 slit. Here we focus on the red-side spectrum, which spanned a usable range of ∼5650–9200 Å, and can be seen in Figure 6; note that we have corrected for Milky Way extinction, with a color excess of E(B − V) = 0.5 mag (Schlafly & Finkbeiner 2011). Two broad absorption features are apparent, at ≈7500 Å and 8200 Å, likely due to O i λ7774 and the Ca ii IR triplet, respectively—these spectroscopic signatures are a clear sign that PS19qp/SN2019ebq is a normal supernova. Spectral classification with the Supernova Identification (SNID; Blondin & Tonry 2007) software package indicates a type Ib/c supernova after maximum light; we plot a +15 day spectrum of SN2004aw (Taubenberger et al. 2006), a type Ic, in Figure 6 to show that this classification is reasonable. We also show a spectrum of AT2017gfo, the optical counterpart of GW170817, at +1.4 days after the GW170817 merger (Smartt et al. 2017; Valenti et al. 2017), clearly showing that PS19qp/SN2019ebq is an unrelated supernova (see also Carini et al. 2019; Dimitriadis et al. 2019; Jencson et al. 2019; McCully et al. 2019; Morokuma et al. 2019).

A candidate GW signal was identified using data from L1, L2, and V1 on 2019 April 26 at 15:21:55.337 UTC. This candidate, S190426c, initially had a 49.3% probability of resulting from a BNS merger, a 23.7% probability of being a MassGap event, a 12.9% probability of resulting from a NSBH merger, a 0.0% probability of being a BBH merger, and a 14.0% probability of being of terrestrial origin (Ligo Scientific Collaboration & VIRGO Collaboration 2019d). The source classification was later significantly revised to a 51.6% probability of being a NSBH merger, 21.5% probability of being a MassGap event, and a 12.9% probability of being a BNS merger, with the BBH and terrestrial probabilities unchanged (Ligo Scientific Collaboration & VIRGO Collaboration 2019e). The distance of the event was also updated to 377 ± 100 Mpc with the 50% probability region covering 214 deg² and the 90% probability region covering 1131 deg².

Although a GW170817-like source would not be detectable by our SAGUARO search at this distance, little is known about the expected emission from NSBH mergers. While Setzer et al. (2019) showed that kilonovae from NSBH mergers should be much dimmer than those from BNS mergers when viewed face-on, the mass ejected is expected to be highly anisotropic, making the viewing angle for these mergers important (Kyutoku et al. 2013). Simulations from Rosswog et al. (2017) found that NSBH kilonovae reach similar optical peak brightness to BNS kilonovae, while Tanaka et al. (2014) found that higher ejecta masses could make the kilonovae from NSBH events more luminous than those with BNS mergers. With such a wide range of predictions for the emission from NSBH mergers, follow-up of these sources is crucial to place the first observational constraints on models.

Once the GW alert was received, the 12 fields shown in Figure 7 were selected to be observed. The highest-probability regions of the localization were either north of our +60 deg decl. limit or fell in an area near the Galactic plane that is not covered by the CSS survey due to crowded fields. As a result, we canceled the initial trigger and triggered the telescope (selecting the same fields from the initial localization) after the localization map was updated (Ligo Scientific Collaboration & VIRGO Collaboration 2019f). These fields were observed at a mid-time of δt ≈ 1.74 days, covering 13.4 deg² of the 50% probability region and 58.9 deg² of the 90% probability region. These fields account for 4.3% of the 50% probability, 5.6% of the 90% probability, and 5.1% of the total probability after the localization was updated (e.g., Ligo Scientific Collaboration & VIRGO Collaboration 2019f). A slightly lower elevation and brighter sky resulted in a 3σ limiting magnitude of 20.8 mag.

Zoom In Zoom Out Reset image size

From the triggered fields, 11307 candidates above 5σ were detected. No known moving objects, transients, or AGNs were found in the median images after cross-matching these candidates against known moving objects from the MPC, known transients from TNS, and known AGNs (Krawczyk et al. 2013; Toba et al. 2014). One candidate was found within the data (AT 2019eij, Delgado et al. 2019) and was not associated with a GLADE galaxy.