Searching for electromagnetic counterparts to gravitational-wave merger events with the prototype Gravitational-Wave Optical Transient Observer (GOTO-4)

ABSTRACT

We report the results of optical follow-up observations of 29 gravitational-wave (GW) triggers during the first half of the LIGO–Virgo Collaboration (LVC) O3 run with the Gravitational-wave Optical Transient Observer (GOTO) in its prototype 4-telescope configuration (GOTO-4). While no viable electromagnetic (EM) counterpart candidate was identified, we estimate our 3D (volumetric) coverage using test light curves of on- and off-axis gamma-ray bursts and kilonovae. In cases where the source region was observable immediately, GOTO-4 was able to respond to a GW alert in less than a minute. The average time of first observation was 8.79 h after receiving an alert (9.90 h after trigger). A mean of 732.3 square degrees were tiled per event, representing on average 45.3 per cent of the LVC probability map, or 70.3 per cent of the observable probability. This coverage will further improve as the facility scales up alongside the localization performance of the evolving GW detector network. Even in its 4-telescope prototype configuration, GOTO is capable of detecting AT2017gfo-like kilonovae beyond 200 Mpc in favourable observing conditions. We cannot currently place meaningful EM limits on the population of distant (⁠|$\hat{D}_L = 1.3$| Gpc) binary black hole mergers because our test models are too faint to recover at this distance. However, as GOTO is upgraded towards its full 32-telescope, 2 node (La Palma & Australia) configuration, it is expected to be sufficiently sensitive to cover the predicted O4 binary neutron star merger volume, and will be able to respond to both northern and southern triggers.

gravitational waves, (transients:) black hole mergers, (transients:) black hole - neutron star mergers, (transients:) neutron star mergers, (transients:) gamma-ray bursts

1 INTRODUCTION

Gravitational wave (GW) detections of compact object mergers are fast becoming common occurrences. During the first half of the advanced LIGO (aLIGO) – advanced Virgo (AdV) Collaboration (LVC) observing run 3 (O3a), which ran from 2019 April 1st to 2019 September 30th, 33 candidate compact object merger events were reported, with 21 of these classified as most likely due to binary black hole (BBH) mergers. The tantalizing prospect of finding an electromagnetic (EM) counterpart to GW events warrants follow-up, particularly in light of the rich scientific yield attained following the detection of both GW and EM signals from the binary neutron star (BNS) merger GW170817 (Abbott et al. 2017a). The dominant demographic of GW events, BBH mergers, is expected to be EM silent (Abbott et al. 2016), although numerous studies suggest that there is potential for an EM counterpart (Palenzuela, Lehner & Yoshida 2010; Moesta et al. 2012; Loeb 2016; Murase et al. 2016; Perna, Lazzati & Giacomazzo 2016; De Mink & King 2017; Janiuk et al. 2017; Chang & Murray 2018; McKernan et al. 2019), and a possible high-energy detection was claimed alongside GW150914 (Connaughton et al. 2016). Assessing observational constraints on the BBH merger population has value, because it helps to inform the follow-up strategy of ‘wide-fast’ survey telescopes chasing GW triggers, and because constraints on the EM signal will naturally become more stringent as the sample increases. Furthermore, there is not yet sufficient observational evidence to rule out the possibility of an EM counterpart in BBH mergers.

More promisingly, the O3a GW trigger population also contains candidate BNS and neutron star–black hole (NSBH) mergers, which are expected to exhibit an EM signature, at least in some cases. The prime candidate for such a counterpart is a kilonova (KN; Li & Paczyński 1998; Rosswog 2005; Metzger et al. 2010; Barnes & Kasen 2013; Metzger 2017), where heavy unstable elements are formed in the neutron-rich environment of the merger via rapid capture (r-process) nucleosynthesis (Lattimer & Schramm 1974; Eichler et al. 1989; Freiburghaus, Rosswog & Thielemann 1999), and subsequently produce thermal emission as they decay radioactively. We may also expect to observe a short duration (<2s; Kouveliotou et al. 1993) gamma-ray burst (sGRB; Blinnikov et al. 1984; Paczynski 1986; Eichler et al. 1989; Narayan, Paczynski & Piran 1992; Rosswog, Ramirez-Ruiz & Davies 2003; Belczynski et al. 2006; Fong & Berger 2013), wherein the merger launches a relativistic jet that can be detected at high energies, and subsequently produces a broad-band synchrotron afterglow as the ejecta decelerate and form shocks in the ambient environment (Blandford & McKee 1976). The KN signature from an NSBH merger is expected to be different to those produced by a BNS (Kawaguchi et al. 2016; Tanaka et al. 2018; Barbieri et al. 2019) – typically predicted to be brighter in the infrared (e.g. Metzger 2017; Kawaguchi, Shibata & Tanaka 2020), although see Foucart et al. (2019). There is also some evidence to suggest that NSBH mergers may be distinguishable in the observed sGRB population (e.g. Troja et al. 2008; Gompertz, Levan & Tanvir 2020).

KNe and sGRBs are already known to be linked to one another through coincident detections (Berger, Fong & Chornock 2013; Tanvir et al. 2013; Yang et al. 2015; Jin et al. 2016; Kasliwal et al. 2017; Gompertz et al. 2018; Jin et al. 2018, 2020; Troja et al. 2018a, 2019; Lamb et al. 2019b), with both BNS and NSBH KN models employed. KNe and sGRBs were also confirmed to be linked to BNS mergers by the detections of GRB 170817A (Abbott et al. 2017a; Goldstein et al. 2017; Hallinan et al. 2017; Margutti et al. 2017, 2018; Savchenko et al. 2017; Troja et al. 2017, 2018b; D’Avanzo et al. 2018; Lyman et al. 2018; Mooley et al. 2018) and the associated kilonova AT2017gfo (Blanchard et al. 2017; Chornock et al. 2017; Coulter et al. 2017; Covino et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Evans et al. 2017; Hjorth et al. 2017; Levan et al. 2017; Nicholl et al. 2017; Pian et al. 2017; Smartt et al. 2017; Soares-Santos et al. 2017; Tanvir et al. 2017; Villar et al. 2017) alongside GW170817 (Abbott et al. 2017b). In this instance, the sGRB afterglow was likely viewed somewhat away from the jet axis (Abbott et al. 2017b; Haggard et al. 2017; Kim et al. 2017; Lazzati et al. 2017; Lamb & Kobayashi 2018; Mandel 2018; Fong et al. 2019; Lamb et al. 2019a; Wu & MacFadyen 2019).

Locating EM counterparts to GW detections remains an extremely challenging task due to the very large localization uncertainties from the GW detectors. These uncertainties can span from hundreds to tens of thousands of square degrees, and while the localization performance will improve as the detector network expands (Abbott et al. 2020a), wide-field survey telescopes remain essential when searching for associated transients. The Gravitational-wave Optical Transient Observer (GOTO; see Steeghs et al., in preparation; Dyer et al. 2018) is one such facility. At design specifications, the project will include one northern node on La Palma (Spain), and one southern node (Australia), each of which will be equipped with 16 telescopes on two robotic mounts, with a field of view of ∼75 square degrees per node. The facility will provide rapid-response tiling of the large LVC error boxes backed by an ongoing sky patrol survey. By exploiting its high cadence and wide field, GOTO can quickly identify potential GW counterparts and flag them for further photometric and spectroscopic follow-up. During O3a, GOTO consisted of four wide field-of-view telescopes, each covering 4.8 square degrees, with a combined footprint of ≈19 square degrees. GOTO’s first installation is located at the Observatorio del Roque de los Muchachos on La Palma and is able to scan the visible Northern hemisphere sky once every ∼14 d. Hereafter, ‘GOTO-4’ specifically refers to the 4-telescope prototype configuration.

Outside of GOTO, there has been a widespread, sustained effort by the observational community to identify EM counterparts to GW triggers throughout the LVC O3 observing run. Numerous wide-field follow-up missions have tiled GW error boxes searching for transients, including the All-Sky Automated Survey for Supernovae (ASASSN; Shappee et al. 2014), the Asteroid Terrestrial impact Last Alert System (ATLAS; Tonry et al. 2018), the Deca-Degree Optical Transient Imager (DDOTI; Watson et al. 2016), the Dark Energy Survey (DES; Dark Energy Survey Collaboration et al. 2016), the Global Rapid Advanced Network Devoted to the Multi-messenger Addicts (GRANDMA; Antier et al. 2020a), KMT-Net (Kim et al. 2016), the Mobile Astronomical System of TElescope Robots (MASTER; Lipunov et al. 2010), MeerLICHT (Bloemen et al. 2016), PanSTARRS (Kaiser et al. 2010), Searches After Gravitational waves Using ARizona Observatories (SAGUARO; Lundquist et al. 2019), the Télescope à Action Rapide pour les Objets Transitoires (TAROT; Boër 2001), the Visible and Infrared Survey Telescope for Astronomy (VISTA; Sutherland et al. 2015), the VLT Survey Telescope (VST; Capaccioli & Schipani 2011), and the Zwicky Transient Facility (ZTF; Bellm et al. 2019). No associated transients were identified (Anand et al. 2020; Antier et al. 2020a, b; Coughlin et al. 2020; Sagués Carracedo et al. 2020), but constraining limits were placed on a number of milestone events, including S190814bv, the first NSBH merger candidate identified in GW (Dobie et al. 2019; Gomez et al. 2019; LIGO Scientific Collaboration & Virgo Collaboration 2019; Ackley et al. 2020; Andreoni et al. 2020; Vieira et al. 2020; Watson et al. 2020), and several candidate BNS systems (Goldstein et al. 2019; Hosseinzadeh et al. 2019; Lundquist et al. 2019), including the unusually massive GW190425 (Coughlin et al. 2019; Hosseinzadeh et al. 2019; Lundquist et al. 2019; Abbott et al. 2020b).

In this paper, we investigate the 29 follow-up campaigns of GW triggers undertaken by GOTO-4 during the first half of the LVC O3 run (2019 April–October), before GOTO was upgraded to eight telescopes. The success of GOTO-4 in tiling the LVC error boxes is discussed. While no GW–EM counterpart candidate detections were made, we assess GOTO-4’s volumetric coverage in 3D by employing test sources to represent the expected EM signatures accompanying GW events, using the O3a events as a benchmark test sample. The observable horizons are compared to the distance distribution of O3a BNS events. The findings are used to inform future strategy, and highlight areas of focus for future upgrades.

In Section 2, we further discuss the data acquisition from the LVC archives and the GOTO pipeline. Section 3 introduces our test sources, and their application is detailed in Section 4. We show our results in Section 5, which are discussed in Section 6. Finally, we present our key conclusions in Section 7. We assume a cosmology of H₀ = 67.4 km s⁻¹ Mpc⁻¹, Ω_M = 0.315, and Ω_Λ = 0.685 (Planck Collaboration et al. 2018) throughout.

2 DATA SAMPLE

2.1 LVC superevents

GW triggers from individual LVC analysis pipelines are aggregated into ‘superevents’, which are announced on the Gamma-ray Co-ordinates Network (GCN) and presented on the Gravitational-wave Candidate Event Database (GraceDB¹). Initial position reconstruction of the GW source is performed by the BAYESTAR algorithm (Singer et al. 2014; Singer 2015). The algorithm outputs a Hierarchical Equal Area isoLatitude Pixelization (HEALPix(add footnote: http://healpix.sf.net); Górski et al. 2005) all-sky map of the posterior probability, as well as the location, scale, and normalization of the conditional distance distribution for each pixel across a grid of millions. At a later time, the BAYESTAR reconstruction may be superseded by a volume reconstruction using LALInference (Aasi et al. 2013; Veitch et al. 2015), the Advanced LIGO Bayesian parameter estimation library. Skymaps are made available for download in the form of a Flexible Image Transport System (FITS; Wells, Greisen & Harten 1981) file.

Using the HEALPix maps, it is possible to calculate the probability contained in a particular region of sky, or to construct a probability density distribution with distance along a given line of sight (for comprehensive recipes, see Singer et al. 2016). Table 1 presents the 29 LVC superevents that were followed up by GOTO-4 during the first half of the LVC O3a run (excluding a further three triggers that were followed up but later retracted by the LVC). The mean and standard deviation of the all-sky posterior probability distance distribution (Singer et al. 2016) is shown, along with the classification of each event and the associated false alarm rate (FAR) from GraceDB. We also include the LVC astrophysical classification, which can be BBH (both binary components are constrained to >5 M_⊙), BNS (both binary components are constrained to <3 M_⊙), NSBH (one object is constrained to <3 M_⊙, and the other to >5 M_⊙), or MassGap (at least one of the binary components is constrained to be between 3 and 5 M_⊙). The terrestrial (noise) classification category is neglected. In all cases, we perform our analysis with the most recently released probability map, and the given distances, classifications, and FARs correspond to this map.

Table 1.

The sample of LVC superevents that was followed up by GOTO. The distance and σ_dist columns represent the posterior mean and standard deviation of the distance to the source, marginalized over the whole sky (Singer et al. 2016). The classification probabilities and FAR are taken from GraceDB. Note that ‘missing’ probability (i.e. cases where the given probabilities do not sum to 100 per cent) will have been assigned to the ‘terrestrial’ (noise) category.

			Classification probability
Event	Distance	σdist	pBBH	pNSBH	pBNS	pMassGap	FAR	Announcement
	Mpc	±Mpc	per cent	per cent	per cent	per cent	(yr−1)	GCN
S190408an	1473	358	>99	0	0	0	8.86 × 10−11	24096
S190412m	812	194	100	0	0	0	5.30 × 10−20	24098
S190421ar	1628	535	97	0	0	0	0.47	24141
S190425z	156	41	0	0	>99	0	1.43 × 10−5	24168
S190426ca	377	100	0	6	24	12	0.61	24237
S190510g	227	92	0	0	42	0	0.28	24442
S190512at	1388	322	99	0	0	0	0.06	24503
S190513bm	1987	501	94	<1	0	5	1.18 × 10−5	24522
S190517h	2950	1038	98	<1	<1	2	0.07	24570
S190519bj	3154	791	96	0	0	0	0.18	24598
S190521g	3931	953	97	0	0	0	0.12	24621
S190521r	1136	279	>99	0	0	0	0.01	24632
S190630ag	1059	307	94	<1	0	5	4.54 × 10−6	24922
S190706ai	5263	1402	99	0	0	0	0.06	24998
S190707q	781	211	>99	0	0	0	1.66 × 10−4	25012
S190718y	227	165	0	0	2	0	1.15	25087
S190720a	869	283	99	0	0	0	0.12	25115
S190727h	2839	655	92	<1	0	3	4.35 × 10−3	25164
S190728q	874	171	34	14	0	52	7.98 × 10−16	25187
S190814bv	267	52	0	>99	0	>1	6.40 × 10−26	25324
S190828j	1803	423	>99	0	0	0	2.67 × 10−14	25497
S190828l	1609	426	99	0	0	0	1.46 × 10−3	25503
S190901ap	242	81	0	0	86	14	0.22	25606
S190910d	632	186	0	98	0	0	0.12	25695
S190915ak	1584	381	>99	0	0	0	0.03	25753
S190923y	438	113	0	68	0	0	1.51	25814
S190924h	548	112	0	0	0	99	2.82 × 10−11	25829
S190930s	709	191	0	0	0	95	0.09	25871
S190930t	108	38	0	74	0	0	0.49	25876

a

Under the assumption that this source is astrophysical in origin, the classification probability becomes NSBH 12 per cent: MassGap 5 per cent: BNS 3 per cent (Ligo Scientific Collaboration & VIRGO Collaboration 2019z).

2.2 Data collection with GOTO

Although we here present an offline analysis of the events, the data collection is driven by a low-latency realtime response. The LVC all-sky probability map is automatically ingested by the GOTO sentinel (Dyer et al. 2018), which analyses the map and produces a schedule of exposures across the observable probability region. These exposures are taken on a fixed grid of tiles on the sky to compare GOTO’s new exposure to a reference image of that same patch of the sky. Observations are scheduled as soon as the initial map is available, and updated whenever a map update is released.

GOTO has the capacity to do multiband photometry having three colour filters. However, for initial follow-up GOTO uses its wide L-band filter. This collects photons between 3750 and 7000 Å, roughly equivalent to a combination of GOTO’s B,G, and R filters and comparable to the combined passband of the SDSS g and r filters.

Each tile is comprised of images from all four telescopes and has a combined field of view of ≈19 square degrees. The default strategy is to visit each tile at least twice per event, where each visit contains a set of 3 × 60 s exposures, which are conducted back-to-back in the wide L-band filter. The three exposures in each set are reduced and then median combined to create the science image. A reference image of the corresponding tile is then subtracted from the science image to identify any new transients. For BNS events, tiles are repeated in order to provide many passes over the skymap over the course of several nights.

2.3 Image processing and data mining

Data are processed on a dedicated cluster of machines located in a Warwick server room. The cluster includes high CPU core count processing nodes, storage arrays, and data base nodes all connected via 10 Gb ethernet. An automated process flow ensures new image frames enter the process queue automatically as they are downloaded from the observatory. For the prototype, a single high-end processing node is sufficient to keep up with processing in realtime (about four 50 megapixel images each minute). The software stack has been developed by the consortium, and performs image level processing, astrometry and photometry calibration, image alignment and subtraction, and source/feature detection. The results are ingested into a PostgreSQL data base in realtime to allow vetting of candidates with short latencies (Steeghs et al., in preparation).

The initial stages of GOTO’s image reduction pipeline cover per-image bias subtraction, dark subtraction, flat-field correction, overscan correction and trimming. Following a source detection pass, an astrometric solution is found using astrometry.net (Lang et al. 2010) and photometric zero-points are derived through comparison with either APASS V filter (Henden et al. 2015) or PS1 g-filter magnitudes for a large number of field stars. The systematic uncertainty in the zero-points varies between Unit Telescopes (UTs) and where a source falls on a given image, but is typically better than 0.15 mag. A set of exposures (usually 3) are then median-combined, correcting for any astrometric offsets between exposures. For each median stack, a reference image is identified (if available), which then triggers the difference imaging stage, using the HOTPANTS tool (Becker 2015).

Features are detected on the subtracted images and passed to a random forest classifier. This makes use of a number of source attributes such as flux, full-width half maximum, local noise level, etc. in an attempt to filter out artefacts and cosmics. During the prototype phase, GOTO employs a temporary classifier that is similar in structure to that of Bloom et al. (2012) and it was trained using an injected source data set. However, this is an area of active development within the collaboration and will be replaced in the next iteration of the pipeline (Mong et al., in preparation; Killestein et al., in preparation). Low scoring sources are marked bogus leaving human judgement to get any remaining high-confidence sources. Features with a reasonable classifier score are then ingested into the GOTO ‘Marshall’, which presents source and contextual information via a browser for human inspection. This entire process is completed approximately 10–20 min after the images are taken in the current prototype pipeline. Human candidate vetting takes place in real-time alongside the follow-up campaign, which typically lasts for several days following a trigger. Any announcements are disseminated via a GCN circulars and/or TNS submissions. For the purpose of this paper, the data are mined after a campaign has been completed using a script which pulls all observations linked to each event. The observations are analysed and their metadata is taken to assess follow-up performance.

3 TEST SOURCES

In order to assess our coverage along the radial distance of the LVC skymap, we must define a test source to recover. For this analysis, we define three physically motivated light curves, as well as a reference source with a constant magnitude. Our physically motivated sources are:

• A gamma-ray burst afterglow, viewed along the jet axis (see Section 3.1).

• A gamma-ray burst afterglow viewed off-axis.

• A Bazin function (Bazin et al. 2011) representing a kilonova-like evolution (see Section 3.2).

These three phenomena are associated with compact object mergers, though none are necessarily expected to accompany the merger of a BBH. However, we note the possible Fermi-GBM detection of GW150914 (Connaughton et al. 2016), and the subsequent theoretical works that attempt to link (weak) GRBs with BBH mergers, as well as some additional optical phenomena (Murase et al. 2016; Perna et al. 2016; De Mink & King 2017; Janiuk et al. 2017; Chang & Murray 2018). A kilonova-like emission profile remains the only EM accompaniment to have been detected alongside a GW signal within the wavelength range and time-scales covered by GOTO (GW 170817/AT2017gfo; Abbott et al. 2017a), albeit from a neutron star origin. A constant source of m_L = 19 is also included in our analysis as a comparison case. This allows us to measure how well GOTO-4 would have performed in retrieving a persistent and reasonably bright new object in a given search field.

Our model light curves are all constructed in the g and r filters. The mean of these two models provides a close match to the GOTO’s L band, which is used during GW follow-up. Light curves at a representative 100 Mpc distance are shown in Fig. 1.

Figure 1.

Light curves of our test sources at an assumed distance of 100 Mpc.

3.1 Gamma-ray burst models

Our GRB afterglow model light curves are constructed following Sari, Piran & Narayan (1998). We assume fairly typical physical parameters for sGRBs (e.g. Fong et al. 2015; Gompertz et al. 2015). The isotropic equivalent energy in gamma rays, |$E_{\gamma , \rm iso} = 10^{52}$| erg, the fraction of energy contained in the emitting electrons, ϵ_e = 0.1, and the fraction of energy contained in the magnetic fields, ϵ_B = 10⁻². The circumburst environment is assumed to have a constant density with radial distance, with a particle density of n = 10⁻³ cm⁻³. Electrons in the forward shock are assumed to be accelerated into a power-law distribution of Lorentz factors with an index of p = 2.2. The half-opening angle of the jet is set to θ_o = 0.087 rad (5°). These parameters result in a jet break at t ≈ 1.82 d (cf. Granot et al. 2018).

In the off-axis case, we use the analytical solution from Granot et al. (2018) for the jet break time

$$\begin{eqnarray*} t_{\rm jb} = 0.7(1+z)\bigg (\frac{E_{51}}{n}\bigg)^{1/3} \bigg (\frac{\theta _{\rm o}}{0.1}\bigg)^2 \hbox{d} \end{eqnarray*}$$

(1)

and peak flux time

$$\begin{eqnarray*} t_{\rm peak}(\theta _{\rm obs}) = \bigg (\frac{\theta _{\rm obs}}{\theta _{\rm o}}\bigg)^2 t_{\rm jb} \hbox{ d,} \end{eqnarray*}$$

(2)

where |$E_{51} = \theta _{\rm o}^2 E_{\rm iso}/2$| is the beaming-corrected energy in units of 10⁵¹ erg. The flux is assumed to be zero at t < t_jb then rises smoothly until its maximum at t = t_peak. At t > t_peak, the flux evolves following the standard on-axis evolution. All physical parameters are set to be the same as the on-axis case, including the jet half-opening angle. We assume that the observer is located at an angle of θ_obs = 0.174 rad (10°) from the jet axis. This choice reflects a more optimistic case that provides a (comparatively) bright signal in the absence of a GRB prompt trigger. As the observer moves further from the jet axis, the expected signal becomes fainter, and the light-curve peaks later.

3.2 Kilonova model

Our Bazin function model is based on the kilonova AT2017gfo. Following the method of Gompertz et al. (2018), we fit Bazin or exponential functions to the full data set² (Andreoni et al. 2017; Arcavi et al. 2017; Cowperthwaite et al. 2017; Díaz et al. 2017; Drout et al. 2017; Evans et al. 2017; Lipunov et al. 2017; Pian et al. 2017; Shappee et al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Valenti et al. 2017; Villar et al. 2017; Pozanenko et al. 2018) to obtain the phenomenological evolution of the g (including g, F475W and V filters) and r (including r, F606W, F625W, and R filters) bands. The rise and peak of the g and r bands are unconstrained for AT2017gfo; in both filters the light curve is best fit with an exponential profile that has a decay time of τ_{f, g} = 0.94 and τ_{f, r} = 1.47 d, and a normalization of A_g = 679.5 and A_r = 631.5 μJy. To avoid overestimating the flux at early times, we modify both functions from exponential to Bazin by adding a rise time of τ_r = 0.1 d, and t₀ = t_max − τ_rln (τ_f/τ_r − 1) = 0.3 d, where t_max is the peak emission time. The resulting light curves are consistent with the data, although the rise and peak parameters are not strictly constrained. Our L-band magnitudes are derived by assigning an equal weighting to the g- and r-band models.

4 METHOD

We first assess the combined GOTO-4 sky coverage for each event. There are two sources of overlap in the observations: first, between individual pointings (tiles) on the sky, and secondly, each tile is the combined footprint of 4 UTs, which themselves overlap one another (see Fig. 2). While this setup is advantageous because it eliminates the chip gap problem seen in other survey missions (and aids intra-telescope calibration) it makes calculating sky coverage more complex, because simply summing the on-sky footprint of the each image (1/telescope) results in a lot of double (or more) counting. Fortunately, a natural solution to double counting coverage is available in the form of the pixel grid in the LVC healpix maps. The healpy³query_polygon routine returns the indices of pixels contained within a user-defined polygon, meaning that we can generate a pixel coverage list by entering the coordinates of the corners of each individual UT snapshot. GOTO-4 sky coverage is then calculated by summing the area of all unique pixel instances in the list. Similarly, the total probability covered is calculated by summing the probability contained within every unique pixel covered. We calculate coverage with the ‘inclusive’ keyword set to false when running query_polygon, meaning that only pixels with centres that lie within our observed area are counted. The typical pixel size of the healpix maps is ≈10 square arcmin; much smaller than the total areas counted.

Figure 2.

Four example GOTO-4 tiles (blue, black, green, and red), each of which is comprised of 4 UT snapshots. There is significant overlap in observation both from overlapping tiles and from overlapping UT fields within a given tile. These effects must be accounted for to avoid ‘double-counting’ the covered probability.

4.1 Volumetric coverage

To assess our volumetric coverage, we compare our test source models (Section 3) at each observing epoch to the limiting magnitude of the science image in order to assess from how far away they could be detected. Model magnitudes are corrected for Galactic extinction using an all-sky foreground reddening, (E(B − V)), healpix map based on the results of Schlegel, Finkbeiner & Davis (1998), which is available for download from the Goddard Space Flight Center website.⁴ Extinction corrections are applied based on which pixel of the reddening map each LVC probability map pixel falls into. The two maps typically have comparable resolution (within a factor of two), but variability in the dust maps on scales of less than ≈10 arcmin² is lost due to being averaged over the pixel. Schlegel et al. (1998) values of E(B − V) are converted to the more recent findings of Schlafly & Finkbeiner (2011) using E(B − V)_SF11 = 0.86 × E(B − V)_SFD98 (Schlafly et al. 2010; Schlafly & Finkbeiner 2011).

Next, we calculate the horizon out to which our extinction-corrected test-source magnitude could be detected by GOTO-4 for each observation. The horizon is defined as the distance from which the extinction-corrected test-source magnitude would be equal to the 5σ upper limit of a given image, m_lim. This limit is based off the science frame. This determines the optimal depth the transient search can achieve from image subtraction (assuming the reference image can see deeper). m_lim is stored in the FITS header of the GOTO-4 observation, and is defined using the pipeline zero-point derived limiting magnitude. Where there are multiple visits to the same LVC probability map pixel, we take the visit with the farthest horizon, and discard the duplicates. Our definition of the horizon implicitly treats m_lim as a ‘hard’ cutoff, where everything brighter is detected and everything fainter is missed. In practice, m_lim represents the average magnitude limit for a given image, which means that in some cases brighter objects would be missed within the horizon, and fainter objects could be detected beyond it (e.g. due to background noise fluctuations, proximity to bright stars, etc). However, the recovery curve, while not a box function, is found to be very steep (Fig. 3), making our hard limit assumption a useful working approximation. In brief, the recovery curve was built by injecting a large number of mock transients into the images, using a large sample of images from the O3a observations and covering a wide range of magnitudes. These fields are then image subtracted, using two different methods, and an automatic source extraction test is done using sextractor. If a source is found within 2 pixels of its injected location it is counted as recovered. Full details of the process are available in Cutter et al. (in preparation). The output shown here is expected of a typical GOTO exposure. For our analysis we show both the HOTPANTS curve that is close to the prototype pipeline performance, as well as the ZOGY in Parallel implementation that is under development.

Figure 3.

The fraction of artificially injected test sources recovered versus magnitude, down to the average m_lim (vertical blue line). In the optimal case nearly 80 per cent of sources are recovered at this limit (horizontal line). The steep turnover indicates that a hard cutoff at m_lim is a suitable working approximation for the observable horizon.

This again supports our working assumption in terms of using our O3 campaigns as guidance towards future strategy, and performance will be further improved. We also note that variations in m_lim across the images are significantly larger than this particular effect, and is the dominant contribution in terms of the overall coverage achieved.

Finally, LVC probability map pixels are sorted into groups of equal observable horizon, where their probability density functions (PDFs) are summed (cf. Singer et al. 2016), and the combined PDF of each group is integrated out to their shared horizon. Our full volumetric probability coverage is then the sum of all of the groups.

4.2 Coverage of an AT2017gfo-like event

While it is important to assess our coverage of all existing GW events, the main targets of EM follow-up facilities like GOTO are BNS mergers that are (relatively) nearby, since these events are the most likely to produce a detectable signature in the form of a KN. AT2017gfo remains the only KN detected alongside a BNS merger that was confirmed by GW detections, and hence is the gold standard. We therefore also assess our coverage and follow-up strategy by calculating the maximum observable horizon of an AT2017gfo-like event in each LVC probability map, and how much of the 2D probability coverage GOTO obtained can be recovered when increasing the assumed distance to the test source. Note that this test is done in addition to the test described in Section 4.1, which was performed on all four test models including our KN analogue.

The method is largely the same as in Section 4.1, except that the test source is always the composite g- and r-band Bazin function fit to AT2017gfo. Instead of constructing and integrating the PDFs from the LVC skymap, we calculate the expected extinction-corrected magnitude of AT2017gfo in each healpix pixel for distance increments out to 500 Mpc. These magnitudes are compared to m_lim, and the column in each line of sight is considered to be covered for the given distance where the magnitude is brighter.

5 RESULTS

Results of the 29 GOTO-4 follow-up campaigns from the first half of the LVC O3 run are shown in Table 2. On average, GOTO-4 began observations 8.79 h after receiving an LVC preliminary alert (9.90 h after trigger), and tiled 732.3 square degrees of the sky per LVC superevent, encompassing 45.3 per cent of the localization probability. When observations were unconstrained at the time of the alert (i.e. night time, with the field above the altitude limit), GOTO-4 began observations less than a minute after receiving the LVC preliminary notice. As an example, Fig. 4 illustrates our follow-up of S190923y, which represents the ‘mean’ campaign given that we covered 723.7 square degrees and 39.4 per cent of the probability in this case (see Table 2). We also calculate the total probability available to GOTO-4 from its site on La Palma for each event, i.e. the probability that lies within the region of the sky above the telescope’s minimum altitude constraint (30°) and away from the Sun. These are shown as pA_vis in Table 2. On average, 64.4 per cent of the total probability was observable for a given event over a 2 to 3 night period, meaning that GOTO-4 tiled ≈70.3 per cent of the available probability per event on average.

Figure 4.

GOTO-4 follow-up of S190923y (25814). Each blue box represents one GOTO-4 tile, with the observing strategy prioritising the highest probability tiles first (darker shaded regions are higher probability). The area covered and pA values in Table 2 comprise the sum of the physical area and probability contained within these tiles, respectively, after accounting for overlap. The shaded area indicates unobservable regions, either from altitude limits or Sun constraints. pA_vis (Table 2) indicates the total fraction of probability that is observable (i.e. that lies within the non-grey region). Axes are in RA/Dec.

Table 2.

GOTO-4 coverage of the LVC probability maps. δt_trig is the time between the GW trigger and the first GOTO-4 observation. δt_alert is the time between receiving the LVC preliminary notification and the first GOTO-4 observation. pA is the percentage of the total probability that was tiled by GOTO-4. pA_vis represents the percentage of the total probability that was visible to GOTO-4 from its site in La Palma, accounting for Sun constraints and altitude limits. The 3D Coverage columns indicate the volumetric coverage for each of the test sources defined in Section 3. The KN Range columns indicate the horizon out to which 90, 50, and zero per cent of the 2D probability coverage is retained in a search for an AT2017gfo-like event. † denotes that a BAYESTAR map was used; no LALInference map was available.

	Response time		2D coverage			3D coverage				KN range
Event	δttrig	δtalert	Area	pA	pAvis	pVbazin	pVGRB	pVoff-axis	pVc19	D90	D50	D0
	(h)	(h)	(deg2)	(per cent)	(per cent)	(per cent)	(per cent)	(per cent)	(per cent)	(Mpc)	(Mpc)	(Mpc)
S190408an†	11.4	10.8	156.1	20.2	23.8	1.20 × 10−5	1.47 × 10−2	2.82 × 10−7	3.22 × 10−5	31	70	135
S190412m†	15.0	14.0	295.2	94.4	94.7	8.68 × 10−3	3.48	0	1.07 × 10−2	107	117	151
S190421ar	48.3	29.1	114.3	8.88	36.6	4.92 × 10−5	3.97 × 10−3	3.89 × 10−7	3.49 × 10−4	57	61	66
S190425z	12.4	9.50	2667.1	22.0	38.1	5.90	20.6	2.57 × 10−3	8.10	46	134	227
S190426c	5.30	5.00	772.7	54.1	70.2	1.10 × 10−2	8.98	0	1.42 × 10−2	4	44	136
S190510g	1.42	0.40	116.1	0.21	0.55	2.06 × 10−3	0.21	0	3.60 × 10−2	48	55	57
S190512at	2.78	2.50	315.1	87.1	92.4	8.52 × 10−5	0.37	0	1.26 × 10−4	22	60	154
S190513bm†	0.55	0.05	116.2	28.5	76.3	1.35 × 10−5	0.59	0	2.51 × 10−5	56	83	120
S190517h†	15.9	15.2	112.7	14.8	51.6	1.40 × 10−6	1.25 × 10−4	0	1.62 × 10−6	49	67	84
S190519bj†	5.35	4.35	664.8	84.7	85.3	2.41 × 10−6	9.55 × 10−4	0	3.64 × 10−6	43	69	161
S190521g	0.13	0.05	393.2	43.7	86.7	8.30 × 10−6	7.57 × 10−2	0	1.11 × 10−5	94	107	126
S190521r†	15.2	15.1	720.7	91.9	92.9	3.85 × 10−6	1.17 × 10−3	0	7.32 × 10−6	9	51	93
S190630ag	2.40	2.40	1170.3	60.9	79.5	1.33 × 10−3	19.0	1.66 × 10−7	3.09 × 10−3	71	112	150
S190706ai	0.33	0.03	543.9	36.7	48.5	8.03 × 10−6	1.07	1.67 × 10−8	2.86 × 10−5	55	94	168
S190707q	12.4	11.7	722.9	34.4	59.3	2.06 × 10−5	2.77 × 10−2	0	2.54 × 10−5	18	53	122
S190718y†	6.58	6.10	242.5	61.2	72.9	1.12	28.9	1.54 × 10−2	2.45	10	27	90
S190720a	0.08	0.04	1358.3	62.1	73.3	1.89 × 10−4	9.51	7.67 × 10−7	5.45 × 10−4	42	54	163
S190727h	15.0	14.9	714.7	42.3	93.5	5.72 × 10−7	6.03 × 10−5	0	1.43 × 10−6	52	66	140
S190728q	14.8	14.5	146.9	89.5	94.0	5.55 × 10−4	1.03	0	8.62 × 10−4	114	124	139
S190814bv	1.83	1.50	717.9	94.1	99.1	1.23 × 10−2	89.6	2.33 × 10−6	2.12 × 10−2	55	61	81
S190828j	16.1	15.8	442.2	9.11	81.6	1.01 × 10−5	2.30 × 10−3	6.45 × 10−8	1.27 × 10−5	34	105	149
S190828l	16.9	16.5	453.6	1.94	50.5	5.60 × 10−5	9.20 × 10−3	4.66 × 10−7	7.34 × 10−5	127	138	154
S190901ap	0.12	0.04	2523.5	38.3	45.3	0.34	30.2	8.40 × 10−4	1.16	62	88	144
S190910d	0.13	0.03	1675.0	41.2	85.1	5.43 × 10−3	17.6	0	1.87 × 10−2	28	69	148
S190915ak	29.9	29.8	18.2	0.08	0.08	3.63 × 10−11	2.39 × 10−9	0	8.42 × 10−11	10	10	15
S190923y†	13.8	13.7	723.7	39.4	59.7	1.91 × 10−2	8.95	0	2.29 × 10−2	46	95	120
S190924h	2.97	2.90	281.3	70.2	73.1	4.52 × 10−5	26.4	5.05 × 10−8	3.59 × 10−4	61	75	101
S190930s	6.28	6.20	2139.9	92.2	92.2	2.20 × 10−3	14.2	1.06 × 10−6	4.48 × 10−3	13	89	142
S190930t†	12.8	12.7	918.2	6.84	9.91	1.24	6.55	1.06 × 10−3	2.01	48	109	130
Mean	9.90	8.79	732.3	45.3	64.4	0.30	9.91	6.87 × 10−4	0.48	48	79	126
Median	6.58	6.20	543.9	41.2	73.1	8.52 × 10−5	1.03	0	3.59 × 10−4	48	70	136

Notably, GOTO-4 tiled more than 500 square degrees on more than half of its campaigns (15), and more than 1000 square degrees on 6 different occasions. This achievement clearly highlights the value of rapid-response, wide-field instruments when searching for poorly localized transient events. GOTO-4’s volumetric coverage is also presented in Table 2, expressed as the percentage of the three-dimensional probability map probed for each test source. The recovery performance is predictably diminished; EM follow-up facilities are not designed to search for predominantly BBH mergers at a mean distance of 1.3 Gpc. However, GOTO-4 was able to cover almost 10 per cent of the total probability volume on average for an on-axis GRB model. We can rule out our test case on-axis GRB in the observed area in three events (S190510g, S190814bv, and S190930t) to better than p < 0.05. In order to place an approximate limit on emission from a BBH merger, we take the mean of the least constraining observation within one day after trigger across every BBH follow-up campaign. This is a 5σ upper limit of m_L ≳ 19 on emission from a BBH merger inside the GOTO-4 field of view and during the GOTO observing epochs at times of less than 1 d after GW trigger.

Such wide searches do uncover a significant number of candidate detections. Following the automated classifier mentioned previously, and an automatic check against known MPC objects (which otherwise would dominate), we are typically presented with on average one high confidence transient candidate per 15 deg² during the O3a searches. Contextual information from external catalogues is used to sort these and flag interesting candidates in the context of the GW search. False positives include un-catalogued variable stars and flare stars, which were ruled out as counterparts via comparison with PS1 images. Many other candidates were disregarded as potential counterparts using the on-sky proximity and redshift of any nearby galaxies through the GLADE catalogue (Dálya et al. 2018), or via a pre-trigger reported detection either in our own detections tables, or already reported by others as a known transient in the Transient Name Server (TNS).

Some of the GOTO photometry of unassociated and uncatalogued transients were reported to the TNS where they were high confidence. However, since O3a was very much a test run for the project and its transient classification stack, we relied on several layers of human confirmation and also required multiple detections. Furthermore, automated TNS submission was not in place and generally we did not report photometry if other projects had already done so, unless it turned out to be a true object of interest. We intend to speed-up this process and improve our classifier robustness going forward so that we can increase our chances of being the first to report.

The expected number of events of each type in our sample are presented in Table 3. We present the results for the localization map as a whole, for the area tiled by GOTO-4, and for those events closer than 250 Mpc, which are our main targets of interest with GOTO. For off-axis GRBs, KNe, and our constant 19th magnitude test sources, the expected rates are essentially unchanged when limiting our results to those events within 250 Mpc. This fact indicates that distant (>250 Mpc) events do not contribute to the probability of detecting sources of this type. In contrast, the chances of detecting an on-axis GRB are greatly diminished when limiting the search volume; their much greater luminosities mean that they can be detected from much further away.

Table 3.

The number of detections of each type of signal we would expect to have seen in our sample and highlighted sub-samples if all LVC superevents contained our test sources. ‘Whole map’ assumes that the test source could be anywhere in the map, whether the region was observed or not. ‘Obs’ assumes that the test source occurred within the region covered by GOTO-4. ‘Nearby’ refers to events that occurred within 250 Mpc. ‘BBH’, ‘BNS’, and ‘All non-BBH’ are given for the ‘Whole map’ scenario, but are limited to events where the given classification was assigned the highest probability by the LVC, excluding terrestrial.

	Events	GRB	Off-axis GRB	KN	19th mag
Whole map	29	2.87	10−4	0.09	0.14
Obs	29	7.04	3 × 10−4	0.49	0.90
Nearby	5	0.86	10−4	0.09	0.14
Nearby, obs	5	4.15	3 × 10−4	0.49	0.90
BBH	17	0.34	10−8	10−4	2 × 10−4
BNS	5	0.89	10−4	0.07	0.12
All non-BBH	12	2.53	10−4	0.09	0.14

The GW population is dominated by distant BBH mergers, but for nearby (<250 Mpc) triggers where the majority of the probability map can be observed and tiled in a two night campaign, GOTO-4 can expect to detect a KN similar to AT2017gfo in fewer than 5 campaigns if one is present in each skymap. This is due to the still relatively poor localization performance of the current GW detector network, coupled with the relatively limited coverage of the GOTO-4 prototype. This will improve substantially in the future (Abbott et al. 2020a), as further discussed below.

Table 2 includes the horizons out to which we retain 90, 50, and 0 per cent of our 2D coverage for an AT2017gfo analogue, following the method in Section 4.2. For example, if GOTO covered 100 per cent of the probability in the LVC skymap in 2D, D₉₀ represents the distance (in Mpc) out to which we probe deep enough to still cover 90 per cent of the probability – where the lost 10 per cent lies along sightlines that were not probed with sufficient depth to recover the transient. D₅₀ indicates that 50 per cent of the observed probability is retained, and D₀ represents the distance at which no GOTO-4 observation was sufficiently deep to recover an AT2017gfo-like KN. This ‘completeness’ versus distance for each event is also shown in Fig. 5, along with the mean of all events. These results show that with reasonable observing conditions, GOTO-4 can detect an AT2017gfo-like event out beyond 100 Mpc, and is capable of achieving 200 Mpc in a favourable line of sight. Of further encouragement is the duplication factor; on average, GOTO-4 observed a given LVC skymap pixel 4.8 times during a campaign (due to a combination of repeat visits and tile/UT overlap), meaning that there is a great deal of scope to improve the observable horizon through image stacking. In some cases, however, the observable horizon is held back by poor observing conditions. These conditions include clouds, high airmass, a bright moon, or even an unfavourable Galactic pointing resulting in high extinction along the line of sight.

Figure 5.

The fraction of the 2D probability coverage that is retained with increasing distance. Each dotted line represents an LVC follow-up campaign with GOTO-4. The thick black line is the mean of the 29 events, and the grey area is the 1σ standard deviation.

6 DISCUSSION

6.1 General constraints

We find that limits on EM emission from the O3a population of BBH mergers using GOTO-4 are currently unconstraining. Under the assumption that every event does contain one of our test sources (an assumption that is extremely unlikely to be true in practice), we would require more than 10⁵ BBH events to place any meaningful limits with our O3a configuration. However, with information on the inclination of the system, we may be able to constrain BBH mergers as GRB progenitors after fewer than 50 events. GOTO-4 covered almost 10 per cent of the LVC probability volume on average for an on-axis GRB model during O3a, and hence we would expect to find an on-axis GRB in ≈one in 10 events if one is present in each case. These constraints will tighten as coverage is significantly improved going forward, thanks to the scaling up of GOTO and the improved localization precision from the LVC. More generally, GOTO-4 places a mean 5σ limit of m_L ≳ 19 in the observed area within one observer frame day from trigger for all BBH events (where observations were taken inside this window). This figure comes from taking the mean value of the least constraining tile from each campaign. However, it does not include foreground extinction, and the observations do not cover the full 24 h period.

Our analysis suggests that the current GOTO-4 follow-up strategy for GW triggers is not suited for placing constraints on associated off-axis GRBs. This limitation arises for two reasons: the first is that they are faint compared to the other expected transients (see Fig. 1) and the second is that the GOTO follow-up strategy focuses on the first two or three nights after the trigger (as is appropriate for KNe), whereas off-axis GRBs will peak several days later than this. The time of observation is the more impactful of the two factors, and we still expect to find untriggered off-axis GRBs in the survey data, because their intrinsic brightness lies well within GOTO’s capability. In the future, GOTO will be able to go deeper faster as more telescopes are deployed, and the LVC localizations will improve, leaving smaller error regions to search. Our findings here suggest that the GOTO GW follow-up strategy may benefit from returning to a candidate field after ∼5–10 d post-trigger, in order to better constrain the possibility of an off-axis GRB. Such discoveries were indeed made during O3 (Kool et al. 2019; Perley et al. 2019), highlighting the scientific potential for revising the strategy.

Our KN model is based on AT2017gfo, a confirmed BNS merger (Abbott et al. 2017b). NSBH mergers are also expected to produce KNe, and their emission profiles may be quite different. If the NS plunges directly into the BH without material being disrupted (the fate of the NS depends primarily on the binary mass ratio and the BH spin; Shibata & Taniguchi 2011; Kawaguchi et al. 2015; Kyutoku et al. 2015; Foucart, Hinderer & Nissanke 2018), then no EM transient is expected at all. The LVC analysis suggested that this was the case for all of the NSBH merger candidates during O3a (S190814bv, S190910d, S190923y, and S190930t; Ligo Scientific Collaboration & VIRGO Collaboration 2019t, w, x, y). If tidal stripping occurs during inspiral, then material may remain outside of the BH event horizon, and a KN can be produced. In the optical frequencies at which GOTO-4 observes, this KN may in fact be fainter than an AT2017gfo analogue, unless a significant (≥0.02 M_⊙) amount of matter is ejected (e.g. Kawaguchi et al. 2020; Zhu et al. 2020). However, for a larger ejecta mass, this signal may in fact be brighter instead. A KN from an NSBH merger may therefore be brighter than AT2017gfo, fainter, or absent entirely. Like with BNS mergers, we emphasize that a significant diversity of light curves are likely to be produced in nature, and that our use of AT2017gfo (as the best-studied and most concrete example) is as a known representative.

6.2 The binary neutron star population

When focusing on nearby (<250 Mpc) events, or those classified as a BNS merger (where BNS is the most probable classification, excluding terrestrial), the prospects of a detection are naturally greatly improved. Table 2 shows that the average GOTO-4 follow-up campaign probes sufficiently deeply to detect AT2017gfo (D_L = 41 ± 3.1 Mpc; Hjorth et al. 2017) because the mean D₉₀ = 48 Mpc. This means that 90 per cent of observed sightlines, by probability contained, are probed deeply enough to detect AT2017gfo at a distance of 48 Mpc. Similarly, per Table 2, the average response time of 9.90 h is fast enough to catch the peak of the emission, and typically 64.4 per cent of the available probability is covered. In favourable conditions the telescope can detect a similar event out to ∼200 Mpc. Notably, Table 3 shows that the chances of detecting a KN is a strong function of the 2D coverage, given the large jump in probability when we calculate our expectations for the surveyed area. This fact will be further exploited by the expansion of the GOTO network, both in terms of an increased number of telescopes as well as a second node in the Southern hemisphere. When fully deployed, this amounts to 2 × 16 scopes, a factor of 8 more than the GOTO-4 prototype presented here. This will allow faster tiling and more comprehensive and deeper coverage.

Table 3 indicates that in its current configuration, with no stacking, GOTO-4 can expect to find one AT2017gfo-like KN per few tens of events observed within 250 Mpc. This relatively low estimated discovery rate is primarily due to the O3a BNS sample being more distant on average than pre-O3 expectations, where the ‘most likely’ estimate of Abadie et al. (2010) suggested ∼4 BNS mergers per year within 100 Mpc (with expectations ranging from 0.04 to 42 events). More recent rate estimations following the LVC O1 and O2 runs predicted 0.5–16 BNS mergers per year (Abbott et al. 2019), and the most recent estimate from the LVC, which includes GW190425 (here S190425z, The LIGO Scientific Collaboration et al. 2020), constrained the rate of BNS mergers within 100 Mpc to be between 1 and 12 events per year.

A naïve estimate, based on the 5 candidate BNS mergers during the first half of O3 (∼6 months) and using the most distant detection (377 Mpc; S190426c) to define the radius of the observed volume, implies a rate of ∼4.5 × 10⁻⁸ Mpc⁻³ yr⁻¹. This indicates that based on the observed population of candidate BNS mergers during the first half of O3, only 0.1 BNS merger events within 100 Mpc are expected; towards the lower end of the Abadie et al. (2010), Abbott et al. (2019), and The LIGO Scientific Collaboration (2020) rate estimates. However, our estimate assumes that the GW detector network was always on during O3a, and does not account for the variation in the measured candidate masses. It also implicitly assumes that all of the candidate BNS events are indeed astrophysical. Such an assumption has further implications for their classifications – for example, under the assumption that S190426c is of astrophysical origin, it is found to more likely be an NSBH merger (Ligo Scientific Collaboration & VIRGO Collaboration 2019z). Our naïve estimate therefore provides a rough guide as to the volumetric BNS merger rate, but should be treated with caution.

Fig. 6 shows the distance distribution of events classified as BNS by the LVC alongside the distribution of observable horizons (assuming a Bazin model) achieved by GOTO-4 during LVC follow-up. We make 100 000 draws from the GOTO-4 horizon PDF, and compare it to 100 000 draws from both the weighted event distance PDF (the weighting for each BNS is |$1/\sigma _{\rm dist}^2$|⁠, see Table 1), and the unweighted event distance PDF including GW170817. Assuming each draw is a unique event with unique follow-up, we find that GOTO-4 reaches the depth required for a detection in 18.3 (22.5) per cent of cases for the weighted (unweighted  + GW170817) distributions.

Figure 6.

The PDF of the horizons achieved for a Bazin model during the 29 GOTO-4 LVC follow-up campaigns (red). These are compared to the weighted PDF (weights |$= 1/\sigma _{\rm dist}^2$|⁠) of the 5 BNS merger events during this time (green, solid), their unweighted PDF (dashed), and the unweighted PDF when the distance to GW170817 (41 Mpc; Hjorth et al. 2017) is included (dotted). The weighted PDF including GW170817 is not shown because the distance errors for this event are far smaller than the other 5 due to its EM detection.

6.3 Factors that influence performance

We further investigate events with unusually high or low observable horizons. The average maximum horizon achieved was 126 ± 40 Mpc (mean ± 1σ standard deviation in the sample). The LVC superevents with the worst observable horizons are S190915ak, S190510g, and S190421ar. At 15 Mpc, S190915ak is the worst by far, and the explanation is an unfortunate combination of high airmass (2) and a bright moon (91.5 per cent illumination). These poor conditions were coupled with a single observing epoch taken almost a full day after peak, when our test model had decayed almost a full magnitude from maximum brightness. For S190510g, we attained a maximum horizon of 57 Mpc. The airmass was 1.78 though the illumination was low at just 30 per cent. In this case, the largest detrimental factor was actually observing too early; all epochs were within 0.1 d of the trigger, preceding the rapid rise to peak of the test model (Fig. 1). The contemporaneous model was two magnitudes below peak. S190421ar had a maximum observable horizon of 66 Mpc, an airmass of 1.88 a fairly high moon illumination of 77 per cent, and also endured high winds, which were on average 21 km h⁻¹ during exposure. This event also suffered from higher extinction along the line of sight compared to the other campaigns, with a median value of A_V ≈ 0.2 mag; more than double the typical median. Furthermore, this event was only observed two days after trigger; two magnitudes below the model peak.

At the high end of the observable horizon distribution, only S190425z is more than 1σ above the mean. The maximum observable horizon achieved in this case was 227 Mpc. The moon illumination was 51 per cent, and the airmass was low, at 1.27. While these conditions are clearly better than the three aforementioned cases, perhaps the biggest contributor to the distant horizon is the fact that the bulk of the observations were taken right at the peak of our model flux. The next highest horizon in the sample was 168 Mpc for S190706ai, which had an airmass of 1.79 and a moon illumination of 34 per cent. These conditions are very similar to S190510g, with the key difference being that S190706ai was observed much closer to the model peak (though not as close as 190425z).

The depth of m_lim is clearly affected by the weather, wind shake, airmass, moon brightness, and the telescope optics, among other things. However, the indication from our analysis of the O3a campaign is that while observing conditions do play a role in determining from how far away we can expect to detect a KN (and can render a campaign entirely unconstraining to our models, as in S190915ak), the timing of the observation with respect to peak flux is in fact the dominant variable. This is reflected in the mean m_lim values; while S190915ak clearly experienced very poor observing conditions, with m_lim = 15.0, the other four events discussed have mean m_lim ranging between 19.3 (S190421ar) and 20.0 (S190425z). Despite this, their observing horizons vary greatly due to the differing proximity of the observations to the model peak. For our Bazin model based on AT2017gfo, observing as close to 0.5 d after trigger as possible is therefore highly desirable. One major caveat to this is that KN evolution at early times (∼a few hours) is largely unknown due to a lack of observations. Our model at these times is therefore not well constrained (see Section 3.2), and hence the poor recovery performance of events with only very early observations may be pessimistic. Additionally, AT2017gfo is only one (well studied) event. The indication from cosmological GRBs is that KNe show considerable diversity in their emission (Fong et al. 2017; Gompertz et al. 2018; Ascenzi et al. 2019; Rossi et al. 2020), meaning that the model employed here may be too optimistic/pessimistic in its peak magnitude, or evolve faster/slower than any given future event. AT2017gfo is fainter than all but one sGRB KN candidate, but in several cases non-detections that probe deeper than the AT2017gfo models imply that there is room for a fainter, undetected population (Gompertz et al. 2018; Pandey et al. 2019).

6.4 Future prospects

In addition to the factors discussed in Section 6.3, improvements in the limiting magnitudes can be made with a simple increase in exposure time. There are two obvious ways that the horizons presented in this paper can be improved: stacking the existing exposures, and/or increasing the exposure time for future follow-up campaigns. The former method could potentially yield depth increases of ∼0.85 mag, since our stacking gains scale as |$\approx 2.5 \log \sqrt{N}$| where N is the number of stacks and our mean duplication factor is 4.8. The second method to increase depth will happen naturally as GOTO approaches design specifications; adding more telescopes means that the footprint of an individual tile/pointing becomes larger, and hence the GW probability regions can be tiled faster. The result is more time for repeat visits and/or longer exposure times, as well as more recent reference tiles from a higher survey cadence. Additionally, as the GW detector network expands to design sensitivity, the localization precision of GW triggers (in particular nearby BNS mergers) is expected to improve (Abbott et al. 2020a), meaning that it may no longer be necessary to tile many thousands of square degrees. The most recent estimates (Abbott et al. 2020a) for the fourth LIGO/Virgo run (O4), which begins in January 2022 and will include the Kamioka Gravitational-wave Detector (KAGRA) suggests that the BNS (BBH) localization precision will improve from a median |$270^{+34}_{-30}$| (⁠|$280^{+30}_{-23}$|⁠) square degree 90 per cent credible interval in O3 to a median |$33^{+5}_{-5}$| (⁠|$41^{+7}_{-6}$|⁠) square degree 90 per cent credible interval during O4. These factors will improve the depth of the observations on average. In particular, the flexible design of GOTO enables it to point both mounts at a single tile to improve depth once the localizations have improved to the point where breadth is no longer an issue.

With regards to the sensitivity required, the expected BNS range during O4 (Abbott et al. 2020a) is 160–190 (aLIGO), 90–120 (AdV), and 25–130 Mpc (KAGRA). For comparison, Fig. 7 shows GOTO’s expected probability recovery fraction with distance in increments of m_lim. This incorporates the real observing conditions and follow-up schedule for each of the 29 events in O3a; it differs from Fig. 5 only in that m_lim is held constant for every tile. We find that attaining m_lim of magnitude 22 (Fig. 7, green line) will provide the range necessary for the predicted O4 BNS distribution for all but the most distant events. This is attainable for GOTO with an increase in exposure time or further image stacking (see e.g. Steeghs et al., in preparation).

Figure 7.

The mean fraction of 2D probability coverage that is retained with distance for our sample using an assumed m_lim for all tiles. ‘CALLIM5’ is the GOTO fits header keyword for m_lim. In addition to greater depth overall, we expect more visits to produce a more ‘box-like’ plot, with flatter sections out to greater distances, followed by faster drop offs. This is because more repeats along a given line of sight will provide deeper alternatives to poor observing epochs with shallow horizons, which create the sloping curvature seen here.

The mean response time of 8.79 h from the time that the LVC preliminary notice is received is impeded by notices that are received during the La Palma day. The forthcoming upgrade to add a southern node in Australia will further improve GOTO’s response time by increasing the window in which one of the facilities can respond (for overlapping latitudes), as well as greatly increasing the range of observable declinations.

7 CONCLUSIONS

We find that even in its 4-telescope prototype configuration, GOTO covered 70.3 per cent of the available probability region on average when following up LVC superevents. The mean area covered was 732 square degrees, but up to 2667 square degrees were covered in a single campaign. Over 500 square degrees were observed during 15 of the campaigns, including more than 1000 square degrees in 6 of them. In cases of well-timed events that fell in unconstrained tiles, GOTO-4 began observations less than a minute after receiving the LVC alert. On average, GOTO began observations 8.79 h after receiving an alert (9.90 h after the GW trigger). Despite no detections of GW–EM counterpart candidates in O3a, the telescope therefore comfortably fulfilled its role as a wide-field GW follow-up facility.

During the first half of O3, a typical GOTO observing campaign covered 45.3 per cent of the LVC probability map, and could have unearthed an AT2017gfo-like KN up to 126 Mpc away. GOTO-4 achieves a maximum depth sufficient to recover our KN test source in one-in-five follow-up attempts for a distance drawn from the BNS distance distribution in O3a. However, we find that due to their distance, it is not possible to place model-constraining limits on EM emission from the distant (>250 Mpc) population of BBH mergers detected by the LVC unless they house on-axis gamma-ray bursts; GOTO-4 was able to cover almost 10 per cent of the total probability volume on average for an on-axis GRB model, and rule one out on three occasions. We also note that future GOTO GW follow-up strategy would benefit from returning to a candidate field at ∼5–10 d post-trigger in order to better constrain the possibility of an off-axis GRB.

Based on our findings, the primary focus for GOTO should now be improving depth. The average duplication factor of 4.8 visits per LVC probability pixel and the recent upgrade to 8 telescopes will help to achieve this goal. We find that reaching a 5σ limit of 22 mag will provide KN coverage of almost the entire BNS volume probed by GW interferometers during O4 (though this ‘volume’ will be anisotropic, and many mergers will be detected at or near the horizon, where most of the volume is). This is expected to be comfortably within reach once GOTO attains its full 2 node, 32 telescope configuration. Efforts to improve depth will be aided by improved localizations from the expanding network of GW detectors, which will allow longer GOTO exposures and/or more repeated visits, and the option to point both mounts at a single target. The second GOTO node in Australia will also allow very rapid response in a greater fraction of cases, as well as enable follow-up of triggers whose fields lie predominantly at southern latitudes. The versatile design of GOTO means that it can evolve as the GW detectors do.

ACKNOWLEDGEMENTS

We thank the anonymous referee for helpful comments that improved the manuscript.

BPG, AJL, and KW have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 725246, TEDE, PI Levan).

DS, KU, and JL acknowledge support from the Science and Technology Facilities Council (STFC) via grants ST/T007184/1,ST/T003103/1, and ST/P000495/1.

Armagh Observatory & Planetarium is core funded by the Northern Ireland Government through the Dept. for Communities.

EP acknowledges funding from the Spanish Ministry of Economics and Competitiveness through project PGC2018-098153-B-C31.

RPB, MRK, and DMS acknowledge support from the ERC under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 715051; Spiders).

RLCS, POB, and RAJEF acknowledge funding from STFC.

The Gravitational-wave Optical Transient Observer (GOTO) project acknowledges the support of the Monash-Warwick Alliance; Warwick University; Monash University; Sheffield University; the University of Leicester; Armagh Observatory & Planetarium; the National Astronomical Research Institute of Thailand (NARIT); the University of Turku; Portsmouth University; and the Instituto de Astrofísica de Canarias (IAC).

Some of the results in this paper have been derived using the healpy and HEALPix packages.

DATA AVAILABILITY

Data products will be available as part of planned GOTO public data releases.

Footnotes

1

https://gracedb.ligo.org/

2

Curated on kilonova.space (Guillochon et al. 2017)

3

https://healpy.readthedocs.io/en/stable/index.html

4

lambda.gsfc.nasa.gov/product/foreground/fg_sfd_get.cfm