Observational constraints on the optical and near-infrared emission from the neutron star–black hole binary merger candidate S190814bv

Context. Gravitational wave (GW) astronomy has rapidly reached maturity, becoming a fundamental observing window for modern astrophysics. The coalescences of a few tens of black hole (BH) binaries have been detected, while the number of events possibly including a neutron star (NS) is still limited to a few. On 2019 August 14, the LIGO and Virgo interferometers detected a high-signiﬁcance event labelled S190814bv. A preliminary analysis of the GW data suggests that the event was likely due to the merger of a compact binary system formed by a BH and a NS. Aims. In this paper, we present our extensive search campaign aimed at uncovering the potential optical and near infrared electromagnetic counterpart of S190814bv. We found no convincing electromagnetic counterpart in our data. We therefore use our


Introduction
The discovery of the binary black hole (BH) merger event GW150914 (Abbott et al. 2016) was a major landmark in the history of physics. It was the first detection of gravitational waves (GWs) and the beginning of GW astronomy. The detection of the first confirmed binary neutron star (NS) merger, GW170817 (Abbott et al. 2017a), and the subsequent discovery of its electromagnetic (EM) counterparts -the short GRB 170817A (Abbott et al. 2017b;Goldstein et al. 2017;Savchenko et al. 2017) and the UV, optical, and IR transient AT2017gfo (Coulter et al. 2017a;Lipunov et al. 2017;Soares-Santos et al. 2017;Valenti et al. 2017) -was a second major breakthrough, and marked the beginning of multi-messenger astrophysics with GWs (Abbott et al. 2017c).
Following the success of the GW170817 follow-up campaign, considerable effort has been expended in mounting similar campaigns with the aim of discovering and characterising the counterparts of new GW events. To optimise the science return of the demanding observations of GW counterparts, a large fraction of the GW and EM community in member states of the European Southern Observatory (ESO) has gathered together to form the ElectromagNetic counterparts of GRAvitational wave sources at the VEry Large Telescope (ENGRAVE) consortium 1 . This paper introduces the collaboration and our first major campaign to search for an EM counterpart to a GW source, 1 http://www.engrave-eso.org/ the NS-BH event merger candidate S190814bv, reported during the O3 run of the LIGO Scientific Collaboration and the Virgo Collaboration (LVC) 2 . S190814bv was detected by the LVC on 2019 Aug 14 21:10:39.01 UT (MJD 58709.88240; The LIGO Scientific Collaboration and the Virgo Collaboration 2019a), and an alert was issued on 2019 Aug 14 21:31:40 UT (MJD 58709.89699), approximately 21 min after the merger. The source was localised to a 50% credible region of 133 deg 2 in the initial report (90% credible region of 772 deg 2 ), which was reduced to 5 deg 2 at 50% and 23 deg 2 at 90% half a day later (The LIGO Scientific Collaboration and the Virgo Collaboration 2019b) making this the best-localised candidate GW event so far. The source distance (as inferred directly from the GW observations) is 267 ± 52 Mpc. The estimated False-Alarm Rate (FAR) is extremely low, at 2.033 × 10 −33 Hz (1 per 1.559 × 10 25 years).
Preliminary parameter estimation (Veitch et al. 2015) indicated that the lighter object had a mass M 2 < 3 M , while the heavier object had a mass M 1 > 5 M , making this a NS-BH candidate according to the LVC classification criteria. The same preliminary analysis pointed to a negligible probability of any disrupted material remaining outside the final compact object (given by the parameter HasRemnant <1%), implying that in this case an EM counterpart was unlikely. We note that the classification of one of the components as a NS is based solely on the mass being <3 M and that a low mass BH is not ruled out by such low-latency classification.
No γ-ray, X-ray, or neutrino signal could be connected to the event (Molkov et al. 2019;Kocevksi 2019;Ohno et al. 2019;Pilia et al. 2019a,b;Sugizaki et al. 2019;Palmer et al. 2019;Evans et al. 2019;Cai et al. 2019;Svinkin et al. 2019; Ice-Cube Collaboration 2019; Ageron et al. 2019;Alvarez-Muniz et al. 2019). The relatively small localisation region led to a world-wide follow-up effort with optical and NIR telescopes (e.g. Gomez et al. 2019a;Andreoni et al. 2020;Dobie et al. 2019;Watson et al. 2020;Antier et al. 2020;Vieira et al. 2020). The ENGRAVE collaboration activated its search programmes to try to discover, or set limits on, an EM counterpart to S190814bv. ENGRAVE members ran wide-field searches for EM counterparts and the ATLAS, GOTO, GRAWITA-VST, Pan-STARRS, and VINROUGE projects also triggered their searches on this event. No promising EM counterpart was detected. In this paper we combine our ENGRAVE ESO/Very Large Telescope (VLT) data, a number of other narrow-field facilities, and the wide-field programmes to report the combined results of our search for a counterpart (Sects. 2 and 3). We place limits on the presence of a counterpart (Sect. 3.6), and discuss the implications of these limits for NS-BH mergers and future searches (Sect. 4). Unless otherwise specified, errors are given at 68% confidence level (1σ), upper limits are given at 3σ, and magnitudes are in the AB system. When needed, we assume a flat FLRW cosmology with H 0 = 70 km s −1 Mpc −1 and Ω m = 0.3.

Wide-field survey observations and results
We employed two different approaches to search for an optical or NIR counterpart to S190814bv. A number of wide-field facilities (with Field of View, FoV, of 1 deg 2 or more) were used to tile the LVC GW sky localization probability maps (skymap) with the aim of covering as much of the 2D probability of S190814bv localisation as possible. These telescopes (with apertures of 0.4 m-4 m) were situated in La Palma, Chile and Hawaii, giving a spread of latitude and temporal coverage. The second approach was to target specific galaxies in the 3D sky region with larger aperture (2 m-8 m) telescopes and smaller FoV cameras. In the next subsections we summarise the search for transients with the different facilities.

The search for transients with GOTO
The Gravitational wave Optical Transient Observer (GOTO 3 ) is a robotic array of wide-field optical telescopes sited at the Roque de los Muchachos Observatory in La Palma. It is operated by the University of Warwick on behalf of an international collaboration. The hardware is modular in design and optimised to autonomously respond to GW events, being able to cover large areas of sky quickly. At the time of S190814bv, GOTO was equipped with 4 active unit telescopes, each having an aperture of 40 cm at f /2.5 and featuring a 50-megapixel CCD detector. This corresponds to a plate-scale of 1 . 25 pix −1 and a FoV of 5.9 deg 2 per camera. Observations of the S190814bv error box were automatically scheduled on the basis of a ranked tiling pattern derived from the available skymaps (Dyer et al. 2018;Gompertz et al. 2020a). Each tile was observed using sequences of 60 s or 90 s exposures with the GOTO-L filter, which is a wide filter covering 400-700 nm (slightly wider than Sloan g+r combined). The bulk of the final localisation probability, as given by the LALInference skymap (LALInference.v1.fits), was covered in 8 observable tiles (Fig. 1), with 89.6% probability covered over the timespan MJD = 58710. 091-58710.230 (5.09-8.34 h after the GW event). Additional observations were obtained the following night, which covered 94.1%. The Moon was closer to the relevant tiles during this second night, affecting the zeropoints achieved in the exposures sets. Most tiles were observed multiple times (see Table B.2), though observing conditions were not optimal, given the presence of the Moon nearby and poor weather. The probability regions were also close to GOTO's lower declination limit, meaning some of the GW probability region could not be observed. Individual exposures are median combined in groups of 3-6 subsequent images and reached 5σ limiting magnitudes covering V = 17.4-19.1 mag. These are derived by photometrically calibrating our photometry to AAVSO Photometric All-Sky Survey (APASS 4 ) stars using the V band.
Raw GOTO images are transferred from La Palma to the Warwick data centre in real-time and processing begins minutes after acquisition using the GOTOphoto pipeline (GOTO collaboration, in prep.). Image-level processing includes detector corrections, astrometry tied to Gaia and photometric zeropoints using a large number of field stars. Difference imaging was performed on the median exposures using recent survey observations as reference. Source candidates were initially filtered using a trained classifier and cross-matched against a variety of catalogues, including the Minor Planet Center (MPC) 5 and Pan-STARRS1 3π survey . The classifier employs a random forest algorithm based on image features, largely following the procedure of Bloom et al. (2012). It was trained using an injected source data-set. Human vetting was performed on the resulting candidates using a web-based marshall interface. No viable optical counterpart candidates could be associated with S190814bv (Ackley et al. 2019).

The search for transients with the VST
The VLT Survey Telescope (VST; Capaccioli & Schipani 2011) is a 2.6 m facility located at Cerro Paranal, Chile, and managed by ESO. The telescope is equipped with OmegaCam (Kuijken 2011), a 268-megapixel camera with a FoV of 1 × 1 deg 2 and a resolution of 0 . 21 pix −1 . The filter set includes Sloan ugriz filters. The telescope is operated in service mode. The GW trigger follow-up is performed using Guaranteed Time allocated to the Italian VST and OmegaCam Consortium.
The monitoring of the S190814bv sky area started on MJD 58710.36 (Grado et al. 2019a,b;Yang et al. 2019). For each pointing and epoch we obtained three dithered exposures for a total exposure time of 135 s. The pointings were visited up to five times during a period of two weeks (see Table B.2). All exposures were obtained using the r filter.
On the first night we imaged 15 deg 2 , covering 53.6% of the localisation probability of the preliminary BAYESTAR skymap (which was the only map available at that time), and 60.7% of the final LALInference skymap probability. Starting from the second epoch we revised the pointing list to optimise the sky coverage for the updated LALInference skymap. The survey area increased to 23 deg 2 covering a maximum of 87.7% of the localisation probability, as shown in Fig. 1.
Details of the image processing and candidate detection are given in Brocato et al. (2018) and Grado et al. (2020). The area identified by the skymap is not fully covered by VST archive observations, hence we used both Pan-STARRS1  and DECam (Abbott et al. 2018a) archive images as templates for the comparison. To select the candidates, we applied a random-forest machine learning algorithm trained on previous search instances (Yang 2018) and then visually inspected the candidates with the highest score. We detected a number of transients from which we removed objects detected only at one epoch and/or associated with stellar sources in the template images. The final list includes 27 transients (reported in Table 2). Out of these, 21 objects were already discovered by other surveys and had been registered on the on the Transient Name Server (TNS) 6 . In Table 2 we include also three transients reported in TNS by other groups that are detected on our images but are below the detection threshold of our search and were therefore not independently detected in our search.
All of these candidates show a slow evolution in the two weeks of the observing campaign (∆m < 1 mag between the first and last detection). Therefore we tentatively exclude that any of them are associated with the GW event (see Sect. 4.4 and Appendix C for more details on the rejected candidates). The limiting magnitudes of the stacked images were estimated by means of artificial star experiments. The limiting magnitude for each pointing, defined as the magnitude at which 50% of the artificial stars are recovered, is reported in Table B Table 1 with the probability contours of the initial skymap (BAYESTAR) and the refined skymap LALInference.

The search for transients with VISTA
The 4.1 m Visible and Infrared Survey Telescope for Astronomy (VISTA) is sited at Cerro Paranal and operated by the European Southern Observatory (ESO). The VISTA InfraRed CAMera (VIRCAM) has 16 detectors arranged in a sparse array, and conventionally six pointings are combined with offsets to form a contiguous "tile" of ∼1.6 deg 2 (Sutherland et al. 2015).
Observations were made with VISTA under the VINROUGE programme at three epochs, the first over several nights postmerger (beginning at MJD 58711.17), the second around a week later, and the final epoch roughly seven weeks post-merger which was used as our primary reference template. We only observed in the K s band (2.15 µm), to optimise our search for a red KN component. A large majority (>90%, see Table 1) of the LALInference localisation area (referred to here as "VISTAwide") was covered at all three epochs, as shown in Fig. 1. The single tile covering the highest likelihood region was re-imaged six times to provide deeper limits in that area (referred to as "VISTA-deep", enclosing ∼21% of the sky localisation probability). Full details of the area covered, timing and representative depth reached are given in Tables 1 and B.2. Initial processing of the data was performed using a pipeline based on the VISTA Data Flow System (VDFS; González-Fernández et al. 2018) modified for on-the-fly processing. Subsequently, the VINROUGE in-house pipeline for the Identification of GW counterparts through NIR Image A113, page 4 of 48 Notes. The start MJD refers to the start of observations on that night (for reference, the GW trigger occurred at MJD 58709.882). The given limiting magnitude is the median magnitude of the individual tiles that covered the probability listed. All times are in the observer frame.
Subtraction (IGNIS) was used to aid the search for transient sources. Using object lists generated from the VDFS pipeline for both the science and template images, positions were cross-checked to create a list in which the majority of objects visible across multiple epochs were removed, along with objects associated with error flags. Template and science images were paired based on area coverage, with templates resized to encompass the entire science area per image. Coordinates were aligned through a combination of the astrometry.net software solve-field, and SExtractor (Bertin & Arnouts 1996) and swarp to match template image positions directly to their corresponding science images. Template images were then subtracted from science images using the hotpants tool 7 (Becker 2015).
We searched for potential counterparts in the subtracted images using two approaches: first by eyeballing the regions around obvious galaxies, particularly those thought to be in the distance range of interest, and secondly through an automated search for sources.
Candidate transients from the automated search were culled based on various criteria, in particular, objects within low confidence regions (e.g. tile edges), with <5 σ detection significance, coincident with foreground stars or the bright cores of galaxies (for which the subtractions often left scars), or unusually sharp images suggestive of hot pixels rather than stellar sources. Moving sources were identified by reference to the MPC. A final check involved human vetting of remaining candidates (typically 1-10 per science image).
This process was repeated for all available science data (including the VISTA-deep field) across the three epochs, crossmatching them over with as many template files as could be attributed to each. Of the sources found in the automated procedure, all were deemed to be image artefacts. Similarly, no convincing new sources were found in the eyeball search, with the exception of the known transient AT2019noq, which was found to have AB magnitudes K = 20.12 ± 0.07 (at MJD 58711.23) and K = 20.06 ± 0.07 (at MJD 58719.25). This source was marginally below the adopted significance threshold in the subtracted image.

The search for transients with ATLAS
The Asteroid Terrestrial-impact Last Alert System (ATLAS) is a high cadence, near-Earth asteroid (NEA) survey with two telescopes located on two separate sites in Hawaii (Mauna Loa and Haleakala). The f /2 telescopes are 0.5 m aperture with 10.56k × 10.56k pixel CCD cameras (Tonry et al. 2018a). The plate scale is 1 . 86 pix −1 , giving each camera a FoV of 29.2 deg 2 . Both units scan the sky between −40 • < δ < +80 • with a cadence of approximately two days, weather permitting. ATLAS survey mode uses two composite filters -"cyan" and "orange" (c and o, respectively). Cyan covers the Sloan g and r filters and orange covers the Sloan r and i filters.
A typical NEA survey observing cycle is comprised of a sequence of 4 slightly dithered exposures (which we call quads), each lasting 30 s, with overheads and processing requiring an additional 10 s. The 4 exposures are typically separated by 15 minutes within a 1 h period to allow for detection and linking of fast-moving objects. ATLAS frequently adjusts this NEA optimised schedule to carry out similar sequences of quads over the sky area of a GW sky map (e.g. Stalder et al. 2017). Observations are processed by an automatic pipeline to produce de-trended, sky-flattened images. These are corrected astrometrically with respect to the ICRS using Gaia stellar positions, and corrected photometrically with respect to a custom built reference catalogue (Refcat2; Tonry et al. 2018b). Difference imaging is employed to identify transients in the survey data and source extraction and measurement are carried out as described in Tonry et al. (2018a). All detections with S /N ≥ 5 are read into a database at Queen's University Belfast and we require 3 or more detections at S /N ≥ 5 to form an object detection. After A&A 643, A113 (2020) such objects are defined, they are subject to various quality filters, machine learning algorithms and cross-matching to known minor planet, star and galaxy catalogues. ATLAS was serendipitously observing the S190814bv skymap region several hours before the GW detection during its normal survey mode. Hence any recent, young and bright transients would have been identified. Seven pointings of ATLAS covered the entirety of the LALInference skymap, and the first pre-discovery observation of the map started at MJD 58709.52 (8.7 h before S190814bv). The coverage continued until MJD 58709.635, or 2.8 h later (Fig. 1). Only the Haleakala telescope observed, in the cyan (c) filter. Some of the earliest exposures were affected by cloud cover and moonlight. No new transient objects which are not cross-matched with stars or known AGN were found in our images. ATLAS re-observed the field on the next two subsequent nights, in the o band (with the Mauna Loa unit). The second night of observations began at MJD 58710.602 (17.26 h after S190814bv) covering 99.8% of the localisation area within a 1 h period. The third night of observations began at MJD 58711.6, again covering 99.8% of the LALInference skymap probability. In none of the three post-event epochs did we find any new transients within the GW localisation area of ATLAS.

The search for transients with Pan-STARRS1
The Pan-STARRS system ) comprises 2 × 1.8 m telescopes on Haleakala, each with a 1.4-Gigapixel camera mounted at the Cassegrain f /4.4 focus of each unit. Here we describe observations with the Pan-STARRS1 telescope (PS1) and the camera GPC1. The GPC1 is composed of sixty Orthogonal Transfer Array devices (OTAs), each of which has a detector area of 4846 × 4868 pixels. The 10 micron pixels (0 . 26) give a focal plane of 418.88 mm in diameter or 3.0 degrees. This provides a FoV area of 7.06 deg 2 , and an active region of about 5 deg 2 (see Chambers et al. 2016, for a description of the focal plane gaps). The five filter system (generally denoted grizy P1 ) is described in Tonry et al. (2012) and Chambers et al. (2016). For filters in common, the PS1 filters have similar transmission profiles as those from the Sloan Digital Sky Survey (SDSS; Abazajian et al. 2009). Images from Pan-STARRS1 are processed immediately with the Image Processing Pipeline Waters et al. 2016). The existence of the PS1 3π Survey data ) provides a ready-made template image of the whole sky north of δ = −30 • , and we furthermore have proprietary i P1 data in a band between −40 • < δ < −30 • , giving reference sky images in the i P1 band down to this lower declination limit. All individual PS1 images have a reference sky subtracted from them  and sources with at least two detections with S /N ≥ 5 significance and spatially coincident to within 0 . 5 are detected and measured. PS1 typically observes in a quad sequence similar to ATLAS, with a set of 4 × 45 s exposures taken across a time span of 1 h to identify and link moving sources. The PS1 team can intervene at any moment and direct the telescope to observe a LVC GW sky map with a flexible choice of filter, exposure time, coverage, and dither and stack strategy. The difference images can be combined into deeper stacks or processed individually and the sources resulting from these are read into a large database at Queen's University Belfast. A series of quality control filters, machine learning algorithms and cross-matches against minor planet, stellar and galaxy catalogues are automatically run and human scanning occurs for all objects not coincident with known solar system objects, stars or catalogued AGN (see Smartt et al. 2016a,b, for more details). At the detection time of S190814bv (2019 Aug 14 21:10:39.01), Hawaii was in day time and PS1 began observing the field at 2019 Aug 15 12:40:37 UT, or 15.50 h after the LVC discovery time.
On the first night of observation the individual 45 s exposures (called "warps") were processed individually to search for any fading transient over the 2 h 33 min period of observation. Image sensitivities are estimated by injecting 500 point sources per skycell across a range of magnitudes and the limiting magnitude is defined when 50% of the sources are recovered (described in the content of database table DiffDetEffMeta in Flewelling et al. 2016). Each chip exposure is warped onto a pre-defined tesselation (called skycells, see Chambers et al. 2016), and the limits refer to these skycells. The 45 s exposures were combined into a nightly stack on the first three nights of observing. The stacks are made by median combining the warps of each skycell. On the first observing night 25 exposures were combined in i P1 and 31 exposures in z P1 in each skycell stack, giving a typical exposure time of 1125 s and 1395 s in i P1 and z P1 , respectively. We did not find any fading transient, but the true constraints are weak due to the dither strategy and fill factor.
For the two subsequent nights, we did not process the individual images, rather we combined all the z P1 band warps into a nightly stack. The effective exposure times were 12480 s and 13440 s on these respective nights. These were deeper than the 3π reference stacks in z P1 in this sky region, so over the next four weeks PS1 observed the region in z P1 to make a deeper and more uniform reference stack. The limiting magnitudes of the skycells on the three nights observing of S190814bv were calculated using the new, custom-made deeper reference-stack for template subtraction. The final sky coverage is plotted in Fig. 1.
All images were processed through the Image Processing Pipeline described above. Detections coincident with known stellar objects from the Gaia DR1 (Gaia Collaboration 2016), Guide Star Catalogue 8 , Two Micron All Sky Survey (2MASS, Skrutskie et al. 2006), SDSS DR12 (Alam et al. 2015) and PS1 ) catalogues were rejected. Additionally any objects coincident with known AGN were identified and excised from the transient search list. As discussed in Smartt et al. (2016a,b), the AGN identification is based mostly on the Véron-Cetty & Véron (2001) and MILLIQUAS 9 catalogues (Flesch 2015). The resultant objects are spatially cross-matched against known galaxies (mostly through the NASA Extragalactic Database, NED 10 ) and all are visually inspected. The objects discovered are listed in Table 3, along with their likely classification Srivastas et al. 2019;). None of these objects is a viable counterpart of S190814bv (see Sect. 4.4 for details on the candidate rejection). LVC skymap of S190814bv (LALInference.v1.fits) in greyscale, with galaxies selected through HOGWARTS over-plotted (black dots). Only the higher probability Northern region of the skymap is included, since the Southern region was not followed up by ENGRAVE. The size of the symbol of each galaxy (i.e. the black dots) is proportional to the probability of hosting the GW event (see Table B.1) given the skymap and a weighing scheme following Arcavi et al. (2017a). The various instruments are illustrated with different colours as in the figure legend, and the typical limiting magnitudes and filters used are given in Table B.1.

Galaxy targeted searches
In addition to the wide-field survey coverage, the unusually tight localisation map of S190814bv (5 deg 2 at 50%), and the distance estimate available from the GW signal (267 ± 52 Mpc), allowed us to define a coordinated programme of multi-wavelength observations of galaxies within the localisation region (e.g. Nissanke et al. 2013;Evans et al. 2016a;Gehrels et al. 2016). While these images cannot cover the whole 2D skymap, they can (often) significantly improve upon the depth of the wide-field surveys for a select number of high-luminosity galaxies (Fig. 2). To identify galaxies with the highest probability of hosting the GW event we utilised the HOGWARTS code 11 (Salmon et al. 2020), which ranks galaxies in the Galaxy List for the Advanced Detector Era (GLADE; Dálya et al. 2018) catalogue according to their probability of containing the corresponding merger given the 3D localization probability density (Singer et al. 2016), and based on the expectation that NS-BH merger rates follow the galaxy mass distribution (Arcavi et al. 2017a).
Since the expected colours and luminosities of the counterparts of NS-BH mergers (see Sect. 4.1) still have significant uncertainties (largely due to the lack of observational constraints), our goal was to obtain multi-colour (optical and NIR) imaging, which we prioritised over observing a greater number of galaxies. This strategy enabled our observations to be sensitive to counterparts that were either blue (e.g. disc-wind driven), or very red due to high lanthanide opacities in dynamical ejecta. While our observations targeted the most luminous galaxies weighted for Fig. 3. VLT/HAWK-I image of a galaxy targeted field. A number of catalogued galaxies (and at least one uncatalogued galaxy likely at the same redshift) are visible in the field. The insets show each individual galaxy as well as the resulting subtraction, demonstrating the absence of variable sources to the limits of the data in any of these possible host galaxies. Each galaxy inset is labelled by the HyperLEDA identifier from the GLADE catalogue, and the corresponding limiting magnitudes are listed in Table B.1. the localization probability, individual telescope pointings were refined in order to capture additional (lower luminosity) galaxies within the localisation volume of the LALInference skymap. We obtained a series of coordinated observations using the Gamma-Ray burst Optical and Near-IR Detector (GROND), the Liverpool Telescope (LT), the Nordic Optical Telescope (NOT), the Telescopio Nazionale Galileo (TNG), the Very Large Telescope (VLT) and the William Herschel Telescope (WHT). In total, over 400 multi-wavelength (grizJHK) images of the 67 most probable galaxies within the 3D volume were obtained in the ten days following the merger. When generating target lists for each telescope, we attempted to avoid unnecessary duplication of observations, while ensuring that the highest probability galaxies were observed to the greatest possible depth. In practice weather, seeing and other scheduling constraints meant that some duplication was unavoidable. Our global coverage is shown in Fig. 2 and a list of observed galaxies in order of decreasing probability (as defined in Sect. 3.7) is given in Table B.1. An example set of observations with the VLT High Acuity Wide field K s band Imager (HAWK-I; Pirard et al. 2004;Casali et al. 2006;Kissler-Patig et al. 2008;Siebenmorgen et al. 2011) is shown in Fig. 3.
Our techniques for searching for transient objects depended on the nature of the data available. All images obtained were manually, and rapidly, compared against existing optical survey data, in particular the PS1 3π survey (e.g. as was done in Coulter et al. 2017b, when AT2017gfo was first discovered). Given the brightness and proximity of the Moon at the time of the observations, only the VLT data exceeded the depth of the PS1 3π survey. For some observations these data remained the best comparison.
In most other cases, when reference images were subsequently obtained, we performed PSF-matched image subtraction using the hotpants code. The residual images were then manually inspected to identify any possible transient sources. We limited our search to the circle centreed on the nucleus of the investigated galaxy, with radius 1.5×R 25 , where R 25 is the galaxy isophotal radius at B = 25 mag arcsec −212 . In general, all galaxies are well subtracted except for some of the brighter nuclei that leave notable residuals and prevent the search for transients. We confirmed that no transients are identified in the difference images. To quantify the depth of these images we inserted artificial stars into the images in different positions within the galaxy search radius with a range of magnitudes and estimated their recovery in our difference images. The limiting magnitude is defined as the average magnitude of the faintest artificial stars that can be visually identified, where optical and NIR photometry is calibrated against the PS1 and 2MASS catalogues, respectively. In all cases these stars had S /N ∼ 3. We found that the limiting magnitude is fairly constant at different positions within the galaxy search radius with the exception of the nuclei.

VLT observations
The Very Large Telescope (VLT) is a facility operated by ESO on Cerro Paranal in the Atacama Desert of northern Chile (Arsenault et al. 2006), which consists of four individual 8.2 m telescopes (UT1-4). We obtained observations of 16 highpriority galaxies with the VLT using i and z imaging with the FOcal Reducer and low dispersion Spectrograph (FORS - Appenzeller et al. 1998, mounted at UT1), and NIR imaging for further 17 galaxies with HAWK-I (Kissler-Patig et al. 2008, mounted at UT4; see Fig. 3 for an example image) in the K s band (see Table B.1). All VLT data were reduced using the standard EsoReflex graphical environment (v2.9.1; Freudling et al. 2013). The observations were performed in one epoch of FORS imaging on 2019 Aug 16 and three epochs of HAWK-I imaging on 2019 Aug 16, 22-23 and 2019 Sept 23-24. Given the likely slow rise time of KNe in the NIR bands, the first two HAWK-I K s -band epochs were intended to be sensitive to the peak of the KN a few days after merger time. This complemented the FORS optical imaging within the first 24-48 h, which is more sensitive to early emission.
FORS observations consisted of 3 × 100 s observations in the i band, although one field was erroneously observed in the z band for the same exposure time. These images reached significantly deeper limiting magnitudes (i ∼ 23-24.5 mag) across the field than those obtained by smaller aperture telescope searches. The cores of some galaxies were, on occasion, saturated, removing our ability to detect transients close to the nucleus.

WHT observations
A series of optical and NIR observations of 17 galaxies in the sample were taken with the William Herschel Telescope (WHT, Boksenberg 1985) from 2019 Aug 14-22. Optical observations of 17 galaxies were obtained in the r band using the Auxiliaryport CAMera (ACAM, Benn et al. 2008) instrument on 2019 Aug 15, while NIR observations of 12 galaxies were taken in the K s band using the Long-slit Intermediate Resolution Infrared Spectrograph (LIRIS, Acosta-Pulido et al. 2002) over the following nights. Both the LIRIS and ACAM images were reduced using standard IRAF procedures 13 and the custom LIRIS package for LIRIS 14 .

TNG observations
Optical and NIR images of a subset of 19 galaxies were carried out with the Italian 3.6 m Telescopio Nazionale Galileo (TNG, Poretti 2018), situated on La Palma, using the optical DOLoRes (Molinari et al. 1999) and near-infrared NICS (Oliva et al. 2001) instruments. Ten galaxies were observed in the r band with relatively short (120 s) exposures obtained on 2019 Aug 15 between 02:13 and 02:54 UT (D'Avanzo et al. 2019a). Image reduction was carried out using standard IRAF procedures. Astrometry was performed using the USNO-B1.0 15 catalogue. The typical upper limit is r ∼ 22.8 (3σ detection limit).
NIR observations of five galaxies were undertaken on 2019 Aug 16, Aug 20 and Sept 5 (usually between 03:00 and 05:00 UT). In addition, the galaxy ESO474-026 was observed on 2019 Aug 17. Each galaxy was observed for 1200 s total exposure time with the K s filter (D'Avanzo et al. 2019b). Image reduction was carried out using the jitter task of the ESO-eclipse package 16 . Astrometry and photometric zeropoints were calculated using the 2MASS 17 catalogue. The typical upper limit is K s ∼ 19.7 − 20.9 (3σ detection limit).

GROND observations
We observed 36 galaxies simultaneously in g , r , i , z , J, H, K s with GROND ), mounted at the 2.2 m MPG telescope at the ESO La Silla Observatory. For each galaxy we obtained an average exposure of 2.1 min in the optical bands and 3.9 min in the NIR bands, and the data were reduced using the GROND pipeline (Krühler et al. 2008), which applies bias and flat-field corrections, stacks images and provides an astrometric calibration. The observations reached typical 3σ detection limits of 20−22 mag in the r band and 17.5−19.5 mag in K s .

LT observations
The Liverpool Telescope (LT; Steele et al. 2004) is a 2 m fully robotic telescope on the Canary island of La Palma, Spain. A total of 19 galaxies were observed using the IO:O imaging camera. IO:O has a 10 × 10 arcmin 2 FoV and was operated with a 2×2 binning, providing a pixel scale of 0 . 3 pix −1 . The observations were made between 01:38 and 05:34 UT on 2019 Aug 15. For all fields, 2 × 150 s exposures in r band were obtained, and for some of the highest-probability candidates we also obtained 2 × 150 s exposures in i band. Reduced images were provided by the IO:O pipeline and stacked with SWarp 18 . Image subtraction of these data were performed using our own subtraction tools rather than hotpants, and detection limits were measured performing PSF photometry at fixed positions over a grid around the centre of each image to determine the median and standard deviation of the sky background as measured in each aperture. Using these measurements, with the calibration tied to the PS1 photometric standards in each field, we derived 3σ limiting magnitudes of 20.3 in both bands.

Galaxy catalogue incompleteness
While searches for EM counterparts targeting known galaxies within the localisation region of a GW are eminently feasible for 99% probability 95% probability 50% probability Fig. 4. 50, 95, and 99% probability regions for S190814bv. Boxes indicate the northern and southern regions of the map as discussed in the text. Galaxies which have a spectroscopic redshift in NED, and lie within the 95% contour at a distance ±3σ that of S190814bv are marked in black; galaxies with an inconsistent spectroscopic redshift are plotted in yellow. The inhomogeneous coverage of NED in the northern contour is clearly visible.
nearby events, such an approach becomes less effective as distance increases. There are two reasons for this: firstly, the density of galaxies per unit area on the sky increases such that tiling the GW map becomes more efficient; secondly, the completeness of galaxy catalogues drops off precipitously beyond 200−300 Mpc. Nonetheless, for S190814bv we used targeted deep optical and NIR observations of some of the most likely host galaxies in the probability map as reported by the HOGWARTS ranking tool. We thus needed to determine the completeness of the galaxy catalogues that HOGWARTS used for the position and distance of S190814bv.
To assess the completeness of GLADE, we queried NED for any galaxies within the 95% probability region of the LALInference skymap, and with a listed spectroscopic redshift. This resulted in 5209 galaxies, of which 1376 have a spectroscopic redshift within 3σ of the S190814bv distance luminosity marginalised over the whole sky (267 ± 52 Mpc). We plot the positions of these galaxies in Fig. 4. What is apparent from Fig. 4 is that the completeness of the NED database varies across the map, with a sharp drop in the number of galaxies above Dec = −25 • . This is almost certainly due to the lack of coverage of the 2dF galaxy redshift survey above Dec = −25 • (Colless et al. 2001).
We attempt to quantify this varying incompleteness in NED in order to determine what fraction of stellar luminosity and mass we have covered in the galaxy targeted search. To this end, we selected galaxy candidates from the PS1 3π catalogue, as this is the deepest, most homogeneous public imaging catalogue available over the whole skymap. PS1 reaches a limiting magnitude of 98% completeness for point sources of 22.5-23 mag in each of g P1 , r P1 , i P1 , with extended source completeness being about 0.5 mag brighter . For reference, an apparent magnitude limit of ∼22 mag corresponds to an absolute magnitude ∼−15.1 mag at a distance of 267 Mpc. This absolute magnitude is comparable to that of the Small Magellanic Cloud, and so at the distance of S190814bv PS1 is essentially complete to all galaxies of relevance. While we may miss some very low surface brightness dwarf galaxies, these contain so little stellar mass that they can be ignored for our purposes (see Sect. 3.7 for a discussion of this).
To create our galaxy candidate catalogue for S190814bv, we queried the PS1 database ) for all sources within the northern 95% localisation region. In order to select only extended objects, we require that g PSF − g Kron > 0.1 mag, and r PSF − r Kron > 0.1 mag, and in addition that the source has n detections > 10 within the PS1 catalogue. Finally, we limit ourselves to the brightest galaxies in the field, setting a threshold of r Kron < 20 mag, which is equivalent to an absolute magnitude of ∼−17 mag at the distance of S190814bv. We also mask out regions in our catalogue around the Sculptor Galaxy NGC 253 and globular cluster NGC 288, which both contain a large number of spurious detections in the PS1 catalogue.
Visual inspection of a random sample of sources from our extended source catalogue confirms that the majority ( 90%) are indeed galaxies (Fig. 5). The small number of sources brighter than r = 14 mag in the catalogue all appear to be saturated, bright stars rather than galaxies, so we impose a brightness cutoff at r = 14 mag. We are finally left with 23 466 candidate galaxies within the northern 95% localisation probability region of S190814bv. In order to better assess the issue of completeness, we cross-matched our PS1 galaxy catalogue against the NED and GLADE galaxy lists, requiring a matching radius of <1 . 5. We show the fraction of galaxies that have an associated NED or GLADE counterpart in Fig. 6. The GLADE completeness reaches a maximum of ∼80% for galaxies brighter than r ∼ 16.5 mag, but drops rapidly at fainter magnitudes, with a completeness of only ∼50% at r = 17.5 mag, and 20% for r < 18 mag; the NED completeness is substantially lower between r ∼ 15.5-17.5 mag.
The completeness of galaxy catalogues in the context of gravitational wave searches was also recently considered by Kulkarni et al. (2018). These authors employed a different approach to this work, by using nearby supernovae with known distances to generate a random sample of galaxies. From this, Kulkarni et al. then determine the fraction with extant spectroscopic redshifts. While Kulkarni et al. look at a somewhat closer distance (<200 Mpc), it is nonetheless encouraging that their "Relative Completeness Fraction" of 75% (which is implicitly weighted by host mass), is comparable to our completeness for the most luminous galaxies.

Probability covered by our targeted search
As demonstrated in the previous section, while the GLADE catalogue is somewhat incomplete in terms of galaxy number in the localisation volume, it contains the majority of the most luminous (and hence most massive) galaxies. In order to quantitatively estimate the efficiency of our galaxy-targeted search, in terms of covering the GW localisation probability, we proceed here to assign a definite probability of being the actual host of S190814bv to all catalogued galaxies in the volume, accounting for the mentioned incompleteness. The full list of targeted galaxies and the corresponding observations are reported in Table B.1. Let us consider galaxies as point-like objects, and let i be an index that runs on all galaxies that are located within the volume V α that contains a given fraction α of the GW 3D posterior localisation probability density P 3D (RA, Dec, d L ) (i.e. the 3D skymap). We assume the probability P gal,i that the merger has taken place within galaxy i to be proportional to the product between P 3D (x i ), namely the GW localisation probability density per unit volume at the galaxy position x i = (RA i , Dec i , d L,i ), and R NS−BH,i , that is the NS-BH merger rate in galaxy i (which in principle depends on its present properties and on its history). Artale et al. (2019) have shown, combining state-of-the art compact binary population synthesis models and cosmological simulations, that the NS-BH rate in galaxies at low redshift correlates almost linearly with the galaxy total stellar mass (R NS−BH ∝ M 0.8 ), with some scatter driven by differences in galaxy merger histories, specific star formation rate and metallicity evolution. Based on these results, for simplicity we assume R NS−BH,i ∝ M i and we use the galaxy K s -band luminosity L K as a proxy for galaxy mass, so that R NS−BH,i ∝ L K,i . This leads to which is similar to the galaxy ranking score used by Arcavi et al. (2017b), but with B-band replaced by K s -band luminosity (which is a better tracer of galaxy mass). In order to compute the normalisation constant A, we impose the condition where N is the total number of galaxies in the volume V α (equality (2) is strictly valid only when α = 1, but in practice it remains correct to an excellent approximation for α close to one). Since our catalogue only contains a fraction of the actual galaxies in the volume, we need to split the sum on the RHS of Eq.
(2) into two terms where N cat is the number of GLADE galaxies within V α . Assuming the remaining K s -band luminosity (present in the volume, but missing from the catalogue) to be uniformly distributed within the volume, we can approximate the last term as where L cat is the total K s -band luminosity in GLADE galaxies, L TOT is the total K s -band luminosity in the localisation volume, and P 3D V α is the GW 3D localisation probability density averaged over the volume, that is P 3D V α = α/ |V α |, where |V α | is the extent of the localisation volume (e.g. in comoving Mpc 3 ). This finally gives the normalisation constant as The total K s -band luminosity in the localisation volume can be estimated as is the K s -band local luminosity density (Hill et al. 2010).
In order to construct the volume that contains 95% of the localisation probability from the latest public 3D GW localisation probability density based on GW parameter estimation (the LALInference skymap; The LIGO Scientific Collaboration and the Virgo Collaboration 2019b), we employ a 3D greedy binning approach. The 3D bins are defined by dividing the sky into tiles using a healpix grid (with N side = 1024), and further dividing the distance coordinate into 3000 linearly-spaced bins between 0 and 700 Mpc. The probability contained in each bin is assigned based on the 3D skymap (following Singer et al. 2016). The bins are then summed in order of decreasing probability density until the enclosed probability equals 95%, which defines the desired localisation volume. The extent of the obtained volume is |V 95% | ≈ 1.4×10 5 Mpc 3 , which gives L TOT ∼ (6.9±1.5)×10 13 L . 1061 GLADE galaxies fall within this volume. Only 45% of these have a K s -band measurement reported in the catalogue, due to the 2MASS magnitude limit. To circumvent this problem we utilise our own VINROUGE observations to obtain K s -band magnitudes for a large fraction of the galaxies. These data were processed through the VISTA Data Flow System (González- Fernández et al. 2018), which provides outputs in the same form as for other VISTA public surveys, including catalogue counts and photometric calibration (per tile) for each observation. We therefore determine and extract K s -band magnitudes for sources on each tile and cross match the resulting catalogues with our GLADE output. This results in 876 matches, providing an 82% completeness. For the brightest galaxies we use the K s −J colour from 2MASS and the redshifts reported in GLADE to k-correct the VISTA magnitudes, in order to compute the corresponding luminosities. For the remaining galaxies, we use the median k-correction (which amounts to 0.10 mag). We finally compute the K s -band luminosity as where 3.27 is the absolute K s -band magnitude of the Sun (Willmer 2018) (in the Vega system),K is the k-corrected Cumulative K s -band luminosity distribution of GLADE galaxies in the S190814bv localisation region, and of our targeted galaxies, compared to the expected distribution in the 95% localisation volume assuming the K s -band luminosity density and distribution from Hill et al. (2010). The grey area in the latter distribution represents the 1σ uncertainty.
K s -band magnitude of the galaxy, and d 0 = 10 pc. Summing over all galaxies with a K s -band magnitude measurement, we obtain L cat ∼ 7.8 × 10 13 L , that is, the GLADE catalogue for S190814bv is essentially complete in terms of K s -band luminosity. This can be seen in Fig. 7, which compares the cumulative K s -band luminosity distribution of GLADE galaxies in the S190814bv 95% localisation volume (red line) and that of our targeted galaxies (blue line) with the expected distribution in the same volume, based on a Schechter fit to the local galaxy luminosity distribution (Hill et al. 2010). This comparison indicates that, despite the incompleteness of current catalogues, galaxytargeting-based searches are still viable out to these distances, as already suggested, for example, by Hanna et al. (2014), Evans et al. (2016b) and Gehrels et al. (2016). All the galaxies in our targeted search (see Table B.1) apart from five have a measured K s -band magnitude reported in GLADE (from 2MASS). In other words, even though the HOGWARTS code selects the galaxies based on their B-band luminosity, the resulting sample is generally bright in K s -band as well. We compute L K (following Eq. (6)) and therefore P gal,i (Eq. (1)) using our VISTA magnitudes, as explained above. The resulting distribution of covered probability as a function of limiting magnitude in different bands is shown in Fig. 8. The sum of the probabilities over the targeted galaxies in our search amounts to ∼50%. This does not enable us to place stringent limits on the properties of the putative EM counterpart of S190814bv using the galaxy-targeted search alone, but it nevertheless shows that targeted searches still have a reasonable chance of detecting a counterpart at ∼250 Mpc.

Discussion
NS-BHs are hybrid merger events that offer insights into a range of behaviours that are not accessible through other mergers. They have both a larger total mass and a larger chirp mass than NS-NS systems. Thus they should produce a stronger GW signal that can be observed out to greater distances. No extant NS-BH systems are known and their range of masses, and astrophysical rates still therefore have very few observational constraints. Population synthesis models show NS-BH systems may be some- what rarer than NS-NS (e.g., Abadie et al. 2010;Dominik et al. 2015;Eldridge & Stanway 2016;Kruckow et al. 2018;Giacobbo & Mapelli 2018;Neijssel et al. 2019), but with significant uncertainties (Belczynski et al. 2016). There is also tentative evidence that they may contribute to the known population of cosmological short GRBs (Gompertz et al. 2020b).
During the last phase of the NS-BH coalescence, the NS can be partially or totally disrupted by the BH tidal field or swallowed directly by the BH without any significant mass left outside the merger remnant (Lattimer & Schramm 1976;. EM emission is expected when the disruption occurs before the NS reaches the innermost stable circular orbit (ISCO) of the BH. Tidal disruption depends on the mass ratio of the two compact objects, on the BH spin, and on the NS Equation of State (EoS, Kyutoku et al. 2011;Foucart 2012;Foucart et al. 2018). Simulations in Newtonian gravity show that the NS can also be disrupted over several orbits (Rosswog 2005;Davies et al. 2005). The properties of the progenitors dictate the mass ejected in tidal tails, and the potential formation of a disc wind. In turn, the properties of the ejecta (mass, electron fraction, entropy, and expansion velocity) determine the nucleosynthetic outcome, and hence the contribution that such binaries may make to the heavy element budget of the Universe (see e.g. Rosswog et al. 2017;Just et al. 2015;Roberts et al. 2017). The presence and properties of ejecta and disc determine also the possible formation of a relativistic jet and hence electromagnetic emission as a short GRB. The EM counterparts could also be much more varied than in the NS-NS case (Rosswog 2017).
Finally, the combination of a GW-detected NS-BH binary with an EM counterpart would enable a standard siren measurement of the Hubble constant and other cosmological parameters (e.g., Schutz 1986;Nissanke et al. 2010) out to larger distances than attainable via NS-NS binaries. In what follows, using our limits, we place constraints on the properties of the putative KN and GRB jet that might have been associated with S190814bv.  Fig. 9. 3σ or 50% completeness upper limits from the wide-field instrument follow-up campaign. The data used are referenced in Sect. 4.1. The kilonova models, representing an AT2017gfo-like evolution, are shifted to the luminosity distance measure from the LVC skymap (267 Mpc; LALInference.v1.fits), and the shaded regions represent the 1σ confidence interval (±52 Mpc). Absolute magnitudes assume a distance of 267 Mpc. Foreground extinction is not included.
merger (Abbott et al. 2017a), it is nonetheless prudent to compare it to our limits for S190814bv because it is the only high-confidence KN to date. Foreground 19 and host galaxy extinction is assumed to be negligible in this analysis.  2017). We find that some of the early VST observations were deep enough to detect a KN of similar brightness to AT2017gfo if one occurred within the 1σ distance confidence interval. The first VISTA-deep observation also constrains an AT2017gfo-like KN down to the S190814bv distance, and several PS1 frames constrain the near end of the distance distribution. However the large distance to this event precludes a strong statement on whether an AT2017gfo-like event would have been detected by PS1, VISTA or VST. Our deepest limits do exclude KNe (within the relevant frames) similar to those which have been claimed to accompany GRB 130603B (Tanvir et al. 2013;Berger et al. 2013), GRB 050709 (Jin et al. 2016), GRB 060614 (Yang et al. 2015, and GRB 150101B (Gompertz et al. 2018;Troja et al. 2018b). These were all brighter than AT2017gfo at similar epochs to our sampling (Gompertz et al. 2018). Note, though, that some claimed KNe are fainter than AT2017gfo, such as those accompanying GRB 160821B (Lamb et al. 2019b;Troja et al. 2019) and GRB 070809 (Jin et al. 2020).
Our galaxy-targeted observations are able to place significantly deeper limits over a fraction of the error box. In partic- 19 The typical value of E(B−V) over the skymap is <0.1. 20 As described on the ENGRAVE webpage, http://www. engrave-eso.org, these Bazin et al. (2011) model fits are purely phenomenological (cf. Gompertz et al. 2018), and describe the temporal evolution of AT2017gfo when shifted to the luminosity distance (267 Mpc) of S190814bv. Our galaxy targeted limits for S190814bv, alongside the equivalent AT2017gfo KN models. Due to the different distances of the observed galaxies, the data and models are presented in absolute magnitudes. For limits below the model lines, our observations would have uncovered a transient comparable to AT2017gfo, had it been present. We also show the apparent magnitude of the data and models when shifted to the luminosity distance of S190814bv (267 Mpc). Foreground extinction is not included. ular, early observations, which exceeded r > 22 mag, are well below the expected brightness of an AT2017gfo-like KN. Data taken in the IR on timescales of 5−10 days reaching K > 21 mag are also competitive. This is illustrated in Fig. 10, which shows the comparison (in absolute magnitudes) between the galaxytargeted limits and AT2017gfo. Solid lines in the figure show our phenomenological fits to the AT2017gfo light curve in the listed bands, converted to absolute magnitudes at the distances of the targeted galaxies. Downward pointing triangles show the absolute limiting magnitudes of our galaxy-targeted observations. As can be seen, early observations from the WHT and TNG in the r-band (Sects. 3.2 and 3.3), and later observations in the IR from HAWK-I are the most constraining (Sect. 3.1), and we are confident in these cases that our observations would have uncovered a KN similar to AT2017gfo if it had been present in the targeted galaxies.
In order to assess more quantitatively the ability of our observations to uncover a putative AT2017gfo-like transient, we can combine our wide-field and galaxy-targeted observations as described in Appendix A. Assuming as our EM counterpart model an AT2017gfo-like event whose flux is scaled by a constant factor, we can derive the covered probability as a function of the ratio between our limiting flux in the most constraining observations and that of AT2017gfo, which is shown in Fig. 11, where we show the covered probability in the r, i, z and K bands (red, purple, blue and cyan lines, respectively), and the combined probability (black line) that corresponds to having a constraining observation in at least one band. For the few, highest probability galaxies the most constraining limit arises from our galaxy targeted programme, but for the majority of the localisation volume the most constraining limits are through wide field observations, in particular from VST (r band), PS1 (z-band) and VISTA (K s -band). Our search is therefore sensitive to an AT2017gfolike KN over ∼40% of the localisation probability (as defined in Sect. 3.7), and over ∼80% to a transient with the same temporal behaviour, but brighter by a factor of 2. The most constraining observations (due to both depth and coverage) are those in the r and K s bands. . Covered probability (defined as described in Appendix A) at which we are sensitive to KN of a given brightness relative to AT2017gfo, based on both our galaxy targeted observations and widefield limits. Coloured lines show the covered probabilities in four different bands, listed in the legend. The black line shows the combined covered probability.

Constraints on the ejecta and on the binary properties
There are fundamental differences between the merger of two NSs (e.g. Ciolfi et al. 2017;Radice et al. 2018;Wollaeger et al. 2018;Shibata & Hotokezaka 2019) and that of a NS and a BH (e.g. Foucart et al. 2013;Kyutoku et al. 2013;Fernández et al. 2017;Tanaka et al. 2014). In the latter case less extreme values for the NS-BH mass ratio, larger BH spin and stiffer NS EoS favour the disruption of the NS before the ISCO, enabling the formation of the accretion disc, tidal tails, and unbound ejecta. This material forms different components from which EM signals can originate (Rantsiou et al. 2008;Pannarale & Ohme 2014;Foucart 2012;Hinderer et al. 2019;Foucart et al. 2019;Barbieri et al. 2019). The KN emission for a given merger is a function of the mass deposited in the various components of the KN, including low-electron-fraction tidal tails (Foucart et al. 2014;Kiuchi et al. 2015;Roberts et al. 2017;Kyutoku et al. 2018) and the neutrino-and viscosity-driven less neutronrich winds (Fernández & Metzger 2013;Just et al. 2015). It is therefore relevant to compare the observational limits on any KN emission from S190814bv with the expectations of NS-BH models.
We use, for that purpose, the multi-component, anisotropic, NS-BH-specific KN model presented in Barbieri et al. (2019), which builds on the NS-NS KN model of Perego et al. (2017). In this model, three outflow components produce KN emission: (1) the tidal ejecta (which are concentrated close to the orbital plane and have the shape of a crescent); two disc-related winds, namely (2) the neutrino-driven wind from the inner part of the accretion disc and (3) the viscosity-driven wind that results from small scale turbulence of magnetic origin inside the disc. For simplicity, since the neutrino-driven wind (2) is expected to unbind only a small fraction of the disc mass in NS-BH remnants, we neglect that component. As a further simplification, we fix the average (root mean square) velocity v and the (grey) opacity κ of the remaining two components to plausible values, namely v t = 0.3c and κ t = 15 cm 2 g −1 for the tidal ejecta and v w = 0.1c and κ w = 5 cm 2 g −1 for the viscous disc wind, based on their expected velocity and composition: the tidal ejecta are typically expected to retain a very low electron frac- tion Y e < 0.2 (e.g. Fernández et al. 2017) leading to efficient r-process nucleosynthesis of Lanthanides and hence a high opacity κ t > 10 cm 2 g −1 (Tanaka et al. 2020); disc wind outflows feature a wider range of Y e , due to viscous heating and neutrino irradiation from the inner part of the disc. Differently from NS-NS mergers, though, the absence of shocks and of intense neutrino production by a meta-stable neutron star remnant are likely to cause the disc wind to remain significantly neutron-rich (e.g. Fernández & Metzger 2013;Just et al. 2015). This justifies our choice of κ w = 5 cm 2 g −1 , which is appropriate for outflows with intermediate Y e ∼ 0.25 -0.35 (Tanaka et al. 2020). We assume θ v = 30 • as the viewing angle 21 -measured with respect to the total angular momentum axis -which is the most likely value for a GW-detected inspiral (see Schutz 2011). We are left with the total masses of the two components, M t and M w , as free parameters. By requiring the resulting light curves to be compatible with our upper limits (both from the galaxytargeted and from the wide-field searches), following the method outlined in Appendix A, we obtain the constraints on M t and M w shown in Fig. 12. The colour map in the figure shows the confidence level at which we can exclude each pair (M w , M t ), denoted as 1 − P(M w , M t ). The region to the lower left of the white dashed line is constrained only by galaxy-targeted observations, while the outer region is constrained mostly by widefield observations. Large tidal ejecta masses M t > 0.05 M are excluded with high confidence >95%, and we can exclude the region M t > 0.01 M and M w > 0.1 M at approximately one sigma confidence. By employing numerical-relativity-based fitting formulae that link the properties of the outflows to those of the progenitor binary (Foucart et al. 2018;Kawaguchi et al. 2016), the limits can also be translated into constraints on the NS-BH binary intrinsic properties, again following Barbieri et al. (2019, see also Barbieri et al. 2020). By assuming the disc wind mass to be 30% of the total disc mass (e.g. Fernández & Metzger 2013;Just et al. 2015;Fernández et al. 2019), we take our representative limits on the disc and tidal ejecta masses (corresponding to the region excluded at 1σ confidence in Fig. 12) to be M disc < 0.3 M and M t < 10 −2 M , respectively. Figure 13 shows the NS-BH parameter space allowed by our limits, for three  Fig. 13. Constraints on the BH spin (a BH ) and the BH mass (M BH ) of the NS-BH binary, assuming remnant disc and tidal ejecta mass limits of M disc < 0.3 M and M t < 10 −2 M , which correspond to approximately 1 sigma exclusion confidence (see text). The SFHo EoS has been adopted to compute the NS tidal deformability.  Figure 14 shows the corresponding limits assuming the DD2 EoS (Typel et al. 2010;Hempel et al. 2012), which is stiffer than SFHo. These two EoSs are representative of the uncertainties in the NS EoS obtained from present nuclear and astrophysical constraints (e.g. Oertel et al. 2017), as well as from constraints derived from GW170817 (Abbott et al. 2018b).
Based on these results, we can therefore exclude that the progenitor NS-BH binary produced a large amount of ejecta. This is consistent with the negligible probability for remnant material left after the merger as reported by the LIGO/Virgo Collaboration in low latency (The LIGO Scientific Collaboration and the Virgo Collaboration 2019a). This indicates that most likely the NS was not disrupted by tidal forces during the final part of the inspiral towards its BH companion, which disfavours high (aligned) BH spins and small mass ratio (or both), as shown quantitatively in Figs. 12-14. At very low mass ratios q = M BH /M NS , the maximum allowed BH spin actually increases with decreasing BH mass. Although the NS tidal disruption is more likely for lower BH masses, the frac-tion of unbound mass (on which our constraint is tighter) is much smaller . As a caveat, we note that the dynamical ejecta fitting formula by Kawaguchi et al. (2016) was calibrated only for mass ratios 3 ≤ M BH /M NS ≤ 7 (but see Foucart et al. 2019). We note that a small or absent amount of mass left outside the remnant BH was found to be the most likely outcome of BH-NS mergers in the population synthesis simulations described in Zappa et al. (2019).

Constraints of GRB afterglow-like emission
By assuming that S190814bv launched a short GRB jet, and that all short GRBs have a similar jet structure to that seen in GW170817, we can use the upper limits on any prompt γ-ray emission to constrain the inclination of the system (e.g. Salafia et al. 2019;Saleem et al. 2020;Song et al. 2019). We employ the two jet structures of Lamb et al. (2019a) (see also Resmi et al. 2018;Salafia et al. 2019), both of which are compatible with the afterglow of GRB 170817A, namely a Gaussian and a two-component structure. For both structures, the central-core isotropic-equivalent kinetic energy is E K,iso = 10 52 erg and the Lorentz factor is Γ = 100. The two-component structure has "wings" with 10% of the core kinetic isotropic equivalent energy and Γ = 5. The core half-opening angles are θ c = 0.09 rad for the Gaussian structure, and θ c = 0.07 rad for the two-component structure. Figure 15 shows the reported Fermi/GBM upper-limit (Kocevksi 2019 -pink line), assuming a soft, ∼1 s duration burst (see Goldstein et al. 2016) at 267 Mpc. Using the two jet structure models described above, the observed isotropic equivalent γ-ray energy for an off-axis observer can be found using the method in Ioka & Nakamura 2019 (which is equivalent to that described in Salafia et al. 2015). We assume a 10% efficiency for energy dissipated as γ-rays by the jet and include opacity due to pair-production where the Lorentz factor is Γ 20−30 following the method in Lamb & Kobayashi (2016, 2017, see also Matsumoto et al. 2019). The top panel of Fig. 15 shows the resulting isotropic-equivalent emitted gamma-ray energy, as a function of the viewing angle, for the Gaussian model (orange dotted line) and the two-component model (blue dash-dotted line) respectively. If there had been a successful GRB 170817Alike jet, the figure shows that the system should be inclined at >10 • , or with a (θ v − θ c ) 5 • where θ c is the jet's core opening angle and θ v is the viewing angle from the central axis. The cosmological population of short GRBs typically have an isotropic γ-ray energy in the range 10 49 E γ,iso 10 52 erg (Fong et al. 2015), thus a successful-GRB producing jet (if any) may have had a lower efficiency or core energy than those assumed here and the lower-limit on the off-axis angle could be smaller.
The afterglow for each structure is shown in the bottom panel of Fig. 15 at an inclination of 10 • and 20 • (thick and thin lines, respectively). The r band VST upper limits are shown as red triangles and the i-and z-band PS1 upper limits as purple and blue triangles, respectively. For the afterglow light curves we assume an ambient density n = 10 −3 cm −3 , microphysical parameters ε B = ε e 2 = 0.01, and an electron distribution index of p = 2.15. In Fig. 15 we only show the r band light-curve, noting that the difference in magnitude for r−z is δm r−z ∼ 0.2 for our model parameters. Our model afterglow light curve is too faint to be constrained by the upper limits. However, for these parameters we can rule out an environment with an ambient density

Comparison to other studies
S190814bv has also been the target of further follow-up reported by other groups. Gomez et al. (2019a) present a study of 96 GLADE galaxies in the 50% error region, which represent 70% of the integrated luminosity in the covered region. They estimate that they cover 25% of all galaxies in the overall localisation area, and are complete down to 0.75 L * (Schechter 1976) in the region they cover. The typical limiting magnitude they obtain after image subtraction is i = 22.2 mag, which is equivalent to M i = −14.9 mag and therefore fainter than AT2017gfo at the distance of S190814bv, at an observing time of ≈ 36 hrs after the GRB. These limits are comparable to the ones we reach with VST in r at a similar time. They rule out KNe with M ej > 0.01 M in their observations, but the incomplete coverage prevents their observations from being constraining at a high confidence. Andreoni et al. (2020) present the results of the EM counterpart search by the GROWTH collaboration, using both DECam wide-field tiling observations and targeted spectroscopic and photometric observations of detected candidate transients. All candidates are found to be unrelated SNe. The wide-field coverage is very complete (>98%), with the most constraining limit reached at 3.4 days, z > 22.3 mag. This limit is comparable to our PS1 z band limits at a similar time. Comparing their limits to different models, they constrain, depending on model and distance, the ejecta mass to M ej < 0.03 · · · 0.1 M , which is consistent with our results.
Watson et al. (2020) present shallow wide-field observations with the DDOTI imager, covering the entire main probability region (90% of the total probability) down to an unfiltered limit of w > 17.9 mag half a day after trigger. They find no candidate transients. They are able to rule out typical on-axis sGRBs but would not have detected a KN similar to AT2017gfo at the distance of S190814bv. Dobie et al. (2019) present their search for radio transients using ASKAP. They find a single significant transient, AT2019osy, which they suggest is likely associated with a low luminosity AGN. ENGRAVE observations of AT2019osy will be presented in a companion paper in preparation. Antier et al. (2020) report on rapid early follow-up by the Global Rapid Advanced Network Devoted to the Multi-messenger Addicts (GRANDMA), reaching limits of 17−18.5 mag within the first two hours, earlier but shallower than the GOTO limits presented in this work. These place no constraints on the potential KN emission associated with S190814bv.
Finally, during the revision of this manuscript a preprint was circulated by Vieira et al. (2020), describing the wide-field optical search by the Canada-France-Hawaii Telescope (CFHT). The search reaches limits comparable to ours, despite covering a lower total localisation probability. The corresponding constraints on the putative KN ejecta masses are similar to ours, even though we caution that they are obtained using non BHNSspecific KN models.

Ruling out identified transients as counterparts
The worldwide intensive efforts to identify a counterpart to S190814bv had led to the identification of multiple transients within the error localisation, even though this remains one of the smallest regions available. In part this is due to the deep observations that were capable of identifying transients sources fainter than 22 mag. In addition to the transients identified here through our searches, additional counterparts have been found by other groups (Andreoni et al. 2020;Gomez et al. 2019a;Dobie et al. 2019;Vieira et al. 2020). In total approximately 75 unique optical transients were identified. In principle, each of these should be considered a potential counterpart unless it can be ruled out through follow-up observations. There are various routes that such an approach can take. Firm reasons for rejection include: 1. The identification of the transient source in imaging taken prior to the detection of the GW event (Pre.Det) 22 . 2. A spectrum of the source or host galaxy that places the source outside of the plausible 3D-GW volume (i.e. too distant, or too close, (Spec.Host.z)). 3. A spectrum which identifies the source as a different kind of transient event, for which the progenitors are known (SN). 4. The source is actually moving, normally because it is an asteroid, but in one case a high proper motion star (Ast, HPM). In addition there are further indications which can be used to disfavour sources, but offer a less secure rejection of their association with S190814bv such as: 5. A photometric redshift which is inconsistent with the 3D-GW volume (Phot.Host.z). 6. A lightcurve which does not match the expectations for the counterparts of NS-NS or NS-BH mergers, but is in keeping with a supernova (SN?) 7. No obvious underlying host galaxy (No.Host) 8. A source which is nuclear in its host galaxy and therefore likely to be related to AGN activity (AGN?).
These latter scenarios (5-8) are not as robust as 1-4 since each has potential pitfalls. Photometric redshifts often have significant associated uncertainty and are prone to catastrophic failure. It is possible that a photometric redshift which formally places the event outside the GW horizon could be in error. The use of lightcurves requires some assumptions as to the nature of the electromagnetic emission from the GW event. Since we have only a single well sampled kilonova, and a handful of events identified superposed to short-GRB afterglows, this provides a limited observational comparison. Furthermore, we have no kilonova clearly associated to black-hole neutron star mergers. Nonetheless there are strong reasons to expect low ejecta masses and hence strong limits on the associated luminosities and timescales, hence the photometric evolution can provide a constraint. The lack of an obvious host galaxy would at first sight suggest a distant object, or a large kick to the progenitor. In most cases we would expect to be able to identify a host within ∼100 kpc of the transient location. The absence of such a host would disfavour an association with the GW event. However, it should also be noted that some short GRBs arise from an apparently "hostless" population (e.g. Berger 2010;Tunnicliffe et al. 2014), and such an event could be missed. Finally, while most nuclear activity is due to either AGN activity or nuclear starbursts, there are suggestions that mergers could be driven at much higher rates within accretion discs around supermassive black holes (Bartos et al. 2017;Stone et al. 2017). Hence, while nuclear events would apparently be disfavoured as counterparts this should not rule them out. Indeed, the transient AT2019osy, identified as a nuclear radio transient, was an object of interest in the error region of S190814bv. For this reason we separate events which are ruled out for one of the reasons 5-8 from those firmly ruled out via 1-4. In Appendix C we present a summary of all transient sources identified in the error region of S190814bv by our searches and those of others. We also indicate the reasons that each of these can be rejected (or not). Of the 73 sources presented there 36 are ruled out robustly, 29 are unlikely based on weaker constraints, the remaining 8 have little information to make such distinctions. There is some value in ascertaining if any of these eight events could be plausible counterparts through future observations, for example to obtain host redshifts. However, we also note that these events are not ruled out due to a paucity of observational constraints, rather than any particular diagnostics which would indicate they are likely related to S190814bv. Indeed, given the small error localisation of S190814bv and its relatively high distance, this list of transients provides some indication of the challenge that will remain in identifying robust EM counterparts to GW sources even in the 4-detector era.

Conclusions
S190814bv was unique amongst the GW detections to date in having an exceptionally small error box. This in turn made it plausible to search for EM emission via targeting of known galaxies. However, such searches were hampered by the large distance to the event (267 ± 52 Mpc compared to ≈40 Mpc for GW170817). Although unprecedented in previous observations, systems like S190814bv may well be more common in the future thanks to both the increasing sensitivity (hence range) of the detectors, and the addition of further GW observatories such as KAGRA in the Kamiokande underground site, Japan, and LIGO-India (Abbott et al. 2018c). Hence, it is relevant to consider what may be the most effective route to the identification of counterparts in this era (corresponding to O4 and beyond). It is striking that at these distances relatively sensitive wide-field searches such as PS1, VST and VISTA do not, in general, reach sufficiently deep limits to constrain a KN comparable in luminosity to AT2017gfo. This suggests that the current generation of widefield facilities may not be especially well suited to the majority of candidates in the future, where observations may need to reach r > 23 mag to probe a reasonable fraction of KN parameter space (see Sagués Carracedo et al. 2020 andCoughlin et al. 2020 who reach similar conclusions). While some wide-field facilities may be able to attain sufficient depth over a significant fraction of future events (e.g. DECam, BlackGEM, Rubin Observatory LSST) it may well be the case that events at ∼300 Mpc may only be detectable by 8 m class observatories. The requirement to observe such events will depend sensitively on where the true event rate of NS-NS and NS-BH lies. At the higher end, optical/IR observers can focus on more nearby events. However, should the event rate lie at the lower end it is necessary to consider if 8 m telescope resources may be needed to identify counterparts, and ELT-like resources required for their follow-up. It therefore continues to be case that the effort should be expended on extending the GW detector network such that the 3D probability volumes for the GW events can become tractable for such observations to be plausible.   For an EM counterpart model defined by an intrinsic light curve dL/dν(ν, t) (specific luminosity at a given rest-frame frequency ν, as a function of rest-frame post-merger time t), one can in principle use the galaxy-targeted and wide-field observations to exclude the presence of such emission in a given galaxy. We define here a framework that allows us to combine the results of different searches with a heterogeneous range of telescopes.
We work under the simplifying assumption that each observation has a well-defined limiting flux, above which we can exclude a detection with high confidence. A galaxy-targeted search consists of a set of observations of N gal galaxies, each observed N obs,i times (i here runs on the galaxies). Each observation takes place at a time t i, j post-merger, and reaches a limiting flux F lim,i, j in a band whose central frequency is ν i, j . We can exclude that galaxy i hosted the putative EM counterpart as long as for any index j running over the N obs,i observations of that galaxy. In that case, the observations contribute a total of P gal,i to the confidence at which the particular EM counterpart model can be excluded. Formally, calling ξ the set of parameters and assumptions that define a particular EM counterpart model, we can exclude ξ with a confidence defined by where H is the Heaviside step function. The quantity E i, j (ξ) is 1 if at least one observation of galaxy i is constraining (i.e. it satisfies inequality reported in Eq. (A.1)), and 0 otherwise. In the absence of a detection, the quantity defined by Eq. (A.2) is most commonly referred to as the "covered probability" with respect to a particular source model. In the case of a wide-field search, one can define the corresponding exclusion confidence as 23 For wide-field observations we assume the localisation probability density to follow the GW 3D skymap, that is, we do not distribute such probability to galaxies as in the galaxy-targeted search case. This relies on the assumption that the galaxy density averages out on scales as large as those probed by wide-field observations. Figure B.1 shows that this is a good approximation in our case, and this removes the uncertainty on the catalogue completeness.
Here the sum runs over a number N tiles of non-overlapping sky tiles (e.g. a healpix tessellation), each observed N obs,i times, at epochs t i, j , with limiting fluxes F lim,i, j in bands whose central frequencies are ν i, j . P tile,i is the (2D) skymap probability density integrated over the tile, while dP/dr defines how this probability density is distributed over luminosity distance r ≡ d L (i.e. skyposition-conditional distance probability density of the tile in the 3D skymap - Singer et al. 2016). Here E i (r, ξ) equals 1 up to the distance beyond which the putative EM counterpart becomes too faint to be detected by the wide-field observations of tile i, and 0 for longer distances. This general framework allows for combining constraints from different wide-field searches.
The results from wide-field and galaxy-targeted searches can be combined conservatively by taking the most constraining between the two for each particular EM counterpart model. GW probability versus galaxy-weighted probability in widefield observations. Comparison between the LALInference GW sky localisation probability density (P 2D ) integrated over the tiles of our wide-field observations and the sum of individual galaxy probabilities P gal (Sect. 3.7) over galaxies that fall in the same tile. The data are reported in Table B.2. Due to the wide field of view of these facilities, the galaxy density averages out, and the probabilities computed in the two ways are very similar. We note that the larger scatter for Pan-STARRS observations is due to the fact that observations are divided into smaller "sky cells". Grouping sky cells in the same observation would reduce the scatter.

Appendix B: Data tables
In this appendix, we provide tables which list exhaustively both our galaxy-targeted (Table B.1) and wide-field (Table B.2) observations. Figure

Notes.
For each galaxy we report: HyperLEDA name (except for those followed by 2MX in parentheses, which are from 2MASSX), right ascension, declination and luminosity distance as listed in the GLADE catalogue * , K s -band magnitude as measured in our own VISTA/VINROUGE data, K s -band luminosity, 3D skymap probability density per unit volume at the galaxy position, and localisation probability (Eq. (1)), followed by a list of observations of that galaxy, for which we specify the telescope, the post-merger time at which the observation was performed, the band, and the limiting magnitude.          Notes. For each facility we report a list of unique tiles, identified by a Tile ID and by the RA and Dec of the tile centre, along with the LALInference sky location probability density (P 2D ) integrated over the solid angle subtended by the tile, and the sum of the galaxy probabilities P gal over galaxies that fall within the tile (see Fig. B.1 for a comparison plot). For each tile, we then report the MJD, Filter and Limiting magnitude of the observations carried out of these tiles.
A113, page 28 of 48 K. Ackley et al.: Optical/near-infrared constraints on a NS-BH merger candidate

Appendix C: Candidate counterparts
The list of publicly reported candidate counterparts for S190814bv is given in Spectroscopic redshift is consistent with GW, while DECAM photometry on TNS appears to show a rapid decline (2.6 mag over one day). However, if associated with S190814bv, then the absolute magnitude three days after the GW event would be i = −18.1.   tifier CFHT0054-2345zau) as a potential counterpart after the first submission of this manuscript, and subsequently added to Table C.1 upon revision. On 08/21 we observed a similar footprint with the Pan-STARRS2 telescope in the i P1 filter. We did not carry out an independent transient search in these data, since the facility was still in science commissioning mode, but we stored the data to provide additional photometry for any interesting source. We cross-checked these data to confirm the proposed i-band transient found by the CFHT search: we have a nearby 5σ detection on MJD 58716.573 (5 days after the CFHT discovery), but it is 0 . 6 from CFHT0054-2345zau, closer to the galaxy's core (0 . 48 separation) and, if real, its measured i P1 magnitude is 22.2 ± 0.2. However, this is clearly offset from the CFHT object, and no excess flux is visible at the position of the CFHT object to i P1 < 22.3. Vieira et al. find i = 21.5 just 1 day later, which is therefore incompatible with our PS2 images. There is no detection in the z P1 data of PS1 either, hence we consider both the CFHT object and the excess flux in the PS2 i P1 -band to be bogus artefacts from image subtraction. Notes.
We only list candidates that were first reported to the TNS after 7 Aug (one week prior to the GW event), and that were discovered up to two weeks after the GW event. The "Comments" column lists the reason why candidates were discounted: Non.Det.
= no subsequent detection after first report; Ast   Bogus?* contributed to governance as a GC member provided comments on the manuscript. GP served on the on-call operations team. SPi is a member of the spectroscopy team and coordinator of the EPO working group. SPo is a core member of the GOTO collaboration. AP provided comments on the manuscript. MAPT contributed WHT observations. MLP provided comments on the manuscript. JQV provided comments on the manuscript. FR contributed obtaining limiting magnitudes of difference images. GR contributed to the GOTO observations. ARa is the PI of the GROND ToO time project. ARe provided Pan-STARRS and ATLAS science analysis. TMR contributed to the reduction and analysis of the HAWK-I data. SSR created and ran the subtraction pipeline for the VISTA imaging. ARo contributed to follow-up with TNG and VST. SR contributed to editing the introduction and Sect. 4.1. NBS contributed to FORS and HAWK-I observations and data reduction and served on the oncall operations team triggering VLT observations. ASC served on the on-call operations team triggering VLT observations. OSS developed and implemented the framework for the comparison of observational limits and theoretical models, led the writing of Sects. 3.7, 4, and Appendix A, produced Figs. 7, 8, 11-14 and B.1, and compiled Tables B.1 and B.2. LSa contributed the ranked galaxy lists from the HOGWARTs algorithm and web application. RS discussed the results and commented on the manuscript. SS provided comments on the manuscript. LSb contributed to the design of X-shooter observations. PScha was a member of the writing team on-call for this event and contributed to editing the manuscript. PSchi is the VST PI and provided guaranteed time VST observations. ASBS carried out observations with Pan-STARRS. TS is a GROND team member and reduced GROND images. SJS contributed to governance as the chair of the GC, co-leads the Pan-STARRS and ATLAS surveys, and contributed manuscript text. MS served on the oncall operations team. KWS developed and ran the Pan-STARRS and ATLAS transient science servers at QUB. JS contributed to governance as a GC member and provided comments on the manuscript. SS helped compile the table of PS1 candidates. ERS contributed as a member of the GOTO team. RLCS was a member of the writing team. DS contributed to governance as a GC member, is the GOTO PI & liason, was the WHT UK PI, contributed data analysis and provided content. GS served on the on-call writing team. CWS provided Pan-STARRS and ATLAS science analysis. NRT is PI of VINROUGE and contributed to governance as a GC member, contributed to VISTA analysis and paper writing. VT contributed to data analysis. JLT co-leads the ATLAS science surveys. MT served on the on-call operations triggering VLT observations and writing teams. KU processed GOTO data -image stacking and photometric calibration. AJvdH provided comments on the manuscript. SDV contributed to management as an EC member, served on the on-call operations team, and coordinated part of the VLT observations. NAW provided comments on the manuscript. DW served on the writing team on-call for this event, on the on-call operations team and contributed to the interpretation of the data. KWie served on the on-call operations team. KWii developed and maintained the computing environment and software for the imaging working group. LW served on the on-call operations team. SY contributed to transient detection in VST data. SXY provided comments on the manuscript. DRY built and helps manage much of the ENGRAVE communication infrastructure.