Introduction

Recent detections of gravitational waves (GWs) have revolutionized our understanding of the population of compact object binaries, impacting many areas of physics and astrophysics (Abbott et al. 2016(Abbott et al. , 2017a(Abbott et al. , 2019a(Abbott et al. , 2019b(Abbott et al. , 2020b(Abbott et al. , 2021a(Abbott et al. , 2021b)). While much can be learned from the GW signals alone, understanding the full astrophysical context of the merger event, including the association to its host galaxy, requires discovery of an electromagnetic (EM) counterpart. Mergers of two neutron stars are thus of particular interest, as they are predicted to produce radiation across the EM spectrum and have long been theorized to be the origin of short gamma-ray bursts (SGRBs; e.g., Eichler et al. 1989;Narayan et al. 1992;Fong & Berger 2013;Berger 2014). This was spectacularly confirmed by the discovery of the binary neutron star (BNS) merger GW170817 (Abbott et al. 2017a(Abbott et al. , 2017b)), which not only had associated gamma-ray emission (Abbott et al. 2017c;Goldstein et al. 2017;Savchenko et al. 2017), but also a bright kilonova (KN) detected in the ultraviolet, optical, and IR bands (Andreoni et al. 2017;Arcavi et al. 2017a;Chornock et al. 2017;Coulter et al. 2017;Cowperthwaite et al. 2017;Díaz et al. 2017;Drout et al. 2017;Kasliwal et al. 2017;Lipunov et al. 2017;Nicholl et al. 2017;Pian et al. 2017;Pozanenko et al. 2017;Smartt et al. 2017;Tanvir et al. 2017;Utsumi et al. 2017;Valenti et al. 2017;Villar et al. 2017Villar et al. , 2018) ) and a long-lasting synchrotron afterglow detected from radio through X-ray wavelengths (Alexander et al. 2017;Haggard et al. 2017;Hallinan et al. 2017;Margutti et al. 2017;Troja et al. 2017;Alexander et al. 2018;D'Avanzo et al. 2018;Dobie et al. 2018;Lyman et al. 2018;Margutti et al. 2018;Mooley et al. 2018aMooley et al. , 2018bMooley et al. , 2018c;;Nynka et al. 2018;Ruan et al. 2018;Troja et al. 2018;Fong et al. 2019;Ghirlanda et al. 2019;Hajela et al. 2019;Lamb et al. 2019;Piro et al. 2019;Troja et al. 2019Troja et al. , 2020)).

Mergers between a neutron star and a black hole are also predicted to result in detectable EM emission in some cases, in particular if the mass ratio of the binary is not too extreme (Kawaguchi et al. 2016;Metzger 2019). It is, however, an open question whether neutron star-black hole (NSBH) mergers also produce SGRBs (Murguia-Berthier et al. 2017;Gompertz et al. 2020), and how those SGRBs would compare to the cosmological SGRB population. The prompt gamma-ray emission from off-axis relativistic jets is highly suppressed due to relativistic beaming, and would likely be undetectable at the larger distances where most GW mergers will occur (Abbott et al. 2017c;Goldstein et al. 2017;Margutti & Chornock 2021). Thus, except for the small fraction of on-axis mergers, the best opportunity to determine whether NSBH mergers produce relativistic jets or outflows of sub-relativistic material is to search for synchrotron emission in the radio or X-ray, as was done for GW170817.

On 2019 August 14, Advanced LIGO/Virgo reported the detection of a new compact object merger candidate GW190814, with a preliminary false alarm rate of one in 10 25 yr (GCN 25324; LIGO Scientific Collaboration, Virgo Collaboration 2019a). It was initially classified as a MassGap event (meaning that the lighter member of the binary had a mass between 3-5M e ), but the classification was revised to an NSBH merger approximately 12 hr later (GCN 25333;LIGO Scientific Collaboration, Virgo Collaboration 2019b). This classification, together with the excellent localization (23 deg 2 with 90% confidence in the skymap provided by LIGO Scientific Collaboration, Virgo Collaboration (2019b) 13.5 hr post-merger), generated considerable interest and telescope investment from the astronomy community. Numerous follow-up efforts across the EM spectrum revealed no evidence for any counterpart (Dobie et al. 2019;Gomez et al. 2019;Ackley et al. 2020;Andreoni et al. 2020;Antier et al. 2020;Gompertz et al. 2020;Morgan et al. 2020;Page et al. 2020;Thakur et al. 2020;Vieira et al. 2020;Watson et al. 2020), broadly consistent with the highly unequal binary mass ratio revealed by the full GW analysis (Abbott et al. 2020b) and the small NS radius inferred from observations of GW170817 (Capano et al. 2020). The nature of the lighter M 2.59 0.09 0.08 -+  component-neutron star or black hole-remains unclear (Essick & Landry 2020;Tsokaros et al. 2020;Tews et al. 2021;Abbott et al. 2021b). Nevertheless, GW190814 provided an excellent test bed for various multiwavelength observing strategies, as its precise localization (tightened to 18.5 deg 2 in the final analysis by Abbott et al. 2020b) and large distance (241 45 41 -+ Mpc) will likely be typical of GW events discovered in O4 and beyond.

In particular, GW190814 prompted several independent searches for a radio counterpart. Unlike the optical sky, the variable radio sky is not well-characterized on timescales of months at the typical flux densities of plausible gravitational-wave transients (although the background rate of extragalactic radio transients is expected to be low; see, e.g., Metzger et al. 2015). This is largely a technological limitation: high-resolution radio imaging generally requires an interferometer, and most radio interferometers have very small fields of view. Thus, early widefield radio transient searches were often shallow (with sensitivity limits > few mJy) and had very limited temporal coverage, relying on just a few epochs to assess variability (e.g., Gregory & Taylor 1986;Levinson et al. 2002;Thyagarajan et al. 2011;Hodge et al. 2013). Alternate strategies were to obtain repeated deep observations of a small region of sky (e.g., Frail et al. 1994;Carilli et al. 2003;Ofek et al. 2011;Mooley et al. 2013;Hancock et al. 2016), or to repurpose data originally taken for other purposes (e.g., calibration fields, Bower et al. 2007;Bell et al. 2011;Frail et al. 2012). In all cases, the number of highly variable sources discovered was small, limiting the conclusions that could be drawn.

More recently, improvements in the mapping speed of existing radio facilities like NSF's Karl G. Jansky Very Large Array (the VLA) and the advent of new facilities with larger fields of view like ASKAP have enabled a new generation of sensitive wide-field radio surveys, refining estimates of the occurrence rates of radio transients on timescales of days to years (Mooley et al. 2016;Bhandari et al. 2018;Mooley et al. 2019). Comparison of new data to previous all-sky surveys has also begun to probe the population of radio transients and variables on timescales of decades (Mooley et al. 2016;Law et al. 2018;Nyland et al. 2021;Wołowska et al. 2021). These searches confirm earlier results that only a few percent of unresolved radio sources are highly variable on these timescales; most of the variables are active galactic nuclei (AGNs; Mooley et al. 2016;Bhandari et al. 2018;Radcliffe et al. 2019). The typical amplitude of AGN variability is small (few percent to factors of ∼few; see, e.g., Hovatta et al. 2008;Richards et al. 2011;Sarbadhicary et al. 2020), but extreme variability on decadal timescales is possible; for example, Nyland et al. (2021) recently discovered a population of quasars that have increased in brightness by 100% to >2500% on timescales of 20 yr.

While further work remains to fully characterize the variable radio sky on the timescales most relevant for GWs (months to years), the low background of transient and variable sources revealed by these searches suggests that GW counterpart searches in the radio sky may be promising. Several wide-field radio searches for GW merger counterparts have been previously employed even in cases when no bright radio counterpart is expected, to better understand the likely background rates of potential contaminating sources (e.g., Mooley et al. 2018d;Bhakta et al. 2021). However, GW190814 was the first event for which a significant fraction of the localization area could be covered to any significant depth by current radio facilities. A wide-field singlefrequency radio search covering 89% of the LIGO Scientific Collaboration, Virgo Collaboration (2019b) localization region was conducted with ASKAP at early times (2-33 days postmerger), ruling out the presence of an on-axis relativistic jet with isotropic-equivalent kinetic energy E iso = 10 51 erg within the observed region under standard assumptions about the jet microphysics (Dobie et al. 2019). This is consistent with the lack of bright X-ray or gamma-ray emission observed at early times, which would have been expected if such a jet were present (Palmer et al. 2019;Page et al. 2020;Watson et al. 2020).

Here, we present targeted late-time radio observations of 75 galaxies within GW190814ʼs localization volume, spanning 1-7 months post-discovery. All observations were taken with the NSF's Karl G. Jansky Very Large Array (the VLA). Our data are aimed at constraining the presence of highly off-axis initially relativistic jets, which might be expected in GW190814 given the measured high inclination of the system from the GW signal (46 deg;

Abbott et al. 2020b

). The timescale of our observations also allows for some limited constraints on the presence of slower-moving KN ejecta. A secondary goal is to characterize the background of variable and transient radio sources likely to be encountered in future radio searches for GW counterparts, with a focus on the implications for galaxy-targeted strategies. We present our observations in Section 2 and discuss our counterpart search and additional targeted observations of our most promising candidate in Section 3, ultimately concluding that this object is more likely to be an unrelated background source. In Section 4, we place limits on the existence of a radio counterpart to GW190814 and discuss the nature of the unrelated variable radio sources uncovered in our search. We conclude in Section 5 with some implications of our work for future radio searches for GW counterparts.

Observations

Galaxy Selection Criteria

To maximize the likelihood of observing the counterpart location with a minimal expenditure of telescope time, we selected a galaxy-targeted search strategy for GW190814. Galaxy-targeted searches are particularly appropriate for telescopes with smaller instantaneous fields of view (like the VLA). The galaxies we target include many of the most optically luminous galaxies within the localization volume, the same galaxies that have been prioritized for counterpart searches at optical and X-ray wavelengths. Such galaxies are attractive targets because they contain much of the stellar mass within the localization volume (and thus have the highest probability of actually containing the merger). However, they are also more likely to have detectable radio emission from unrelated sources (e.g., star formation, AGN activity), making searches for radio transients more challenging.

To select our targets, we generated a list of all galaxies from the Galaxy List for the Advanced Detector Era (GLADE) catalog (Dálya et al. 2018) with B-band luminosities L  0.1L * in the LALINFERENCE 90% localization volume that was circulated by the LIGO/Virgo collaboration 13.5 hr post-discovery (LIGO Scientific Collaboration, Virgo Collaboration 2019b). We then ranked this list based on a weighting of the galaxy's spatial position within the localization volume and the galaxy's B-band luminosity (a proxy for stellar mass), using the same procedure followed in Gomez et al. (2019) and Hosseinzadeh et al. (2019). We observed the top 75 galaxies on this list, out of a total of 723 galaxies. The observed galaxies contained ∼21% of the cataloged integrated stellar luminosity in the region. However, the GLADE catalog is known to be incomplete at the distance of GW190814 (Dálya et al. 2018;Abbott et al. 2020b). We therefore estimate the completeness of the catalog by integrating a Schechter B-band galaxy luminosity function down to 0.1L * to approximate the true total number of galaxies and the corresponding total integrated stellar luminosity contained in the region (see Gomez et al. 2019 for the exact function used). We find that the GLADE catalog is ∼50% complete down to 0.1L * at the distance of GW190814, in terms of the number of galaxies. We estimate that our 75 galaxies contain 14% of the total integrated stellar luminosity within the LALINFERENCE localization volume.

The final localization volume presented in Abbott et al. (2020b) shifted slightly compared to the LALINFERENCE localization volume (Figure 1, left panel) and shrank significantly, from 1.1 × 10 5 to 3.9 × 10 4 Mpc 3 . We therefore repeat the above calculations for the final localization volume. We find that 65 of our target galaxies remain within the final 90% localization volume presented in Abbott et al. (2020b), comprising 32% of the total integrated stellar luminosity within this region. Thus, assuming that the merger probability tracks stellar light, we have a roughly 32% chance that our observations covered the true position of the merger.

VLA Observations

We observed the top 75 galaxies in our ranked list in 2019 September (∼1 month post-merger; "epoch 1") and again in 2020 February/March (∼6-7 months post-merger; "epoch 2"). To facilitate scheduling, we split each epoch into three 3h 20 m scheduling blocks of 25 galaxies each, with the observing order chosen to minimize slew times. A complete list of the galaxies observed and the timing of each observation is given in We detect two significant point-like radio sources within our search region (dashed light blue line): "Candidate 1" (variable, not detected at 35 days post-merger), and "Source 2" (consistent with constant flux density). We also observe weak extended emission near the nucleus of ESO 474-035. At the distance of ESO 474-035 (d L = 271 Mpc), Candidate 1 would be at a projected offset of 87 kpc (blue line), at the upper end of the distribution measured from short GRBs. Unlike Source 2, which has a clear optical counterpart, Candidate 1 has no coincident optical emission to i  22.23 mag at 1.48 days after the merger, and no underlying, static host galaxy to r = 24.92 mag.

Table 1. The position of each galaxy and the size of the region searched for transients in each pointing is shown in comparison to the GW localization in Figure 1 (left panel). The search region from a sample pointing is also shown in Figure 1 (right panel).

To maximize the chance of detection, we utilized the sensitive C band receiver in 3-bit mode (4 GHz bandwidth, mean observing frequency 6 GHz) for all observations. We used 3C147 as our flux and bandpass calibrator for all of our observations and one of three phase calibrators, J0011-2612, J0118-2141, or J0120-2701, depending on which was closest to a given galaxy. We calibrated all data using the standard NRAO pipeline in CASA (McMullin et al. 2007) and imaged the data using the CASA task clean.

The first epoch of observations was taken when the VLA was in its most extended A configuration (beam size ∼0 33 at 6 GHz) and the second epoch when it was in its more compact C configuration (beam size ∼3 5 at 6 GHz). As any physically plausible radio afterglow at the distance implied by the GW signal (d 241

) is predicted to be an unresolved point source on the timescale of our observations, the differing resolution of the radio observations should not impact the flux density recovered for any (sufficiently isolated) bona fide counterpart. Nevertheless, the C-configuration data are more sensitive to diffuse emission from the merger host galaxy or other sources in the field, which caused additional challenges with the data imaging for a subset of our targets. This manifests as an elevated rms noise level in a small fraction of our images. In epoch 1, we achieved a median image rms of 12.2 μJy beam -1 , and in epoch 2, we achieved a median image rms of 17.6 μJy beam -1 . In both epochs, our typical time on source was ∼6 m 20 s for each galaxy.

Counterpart Search

Candidate Selection Criteria

We next searched each galaxy field for a radio counterpart to GW190814. To identify radio sources, we used Source Extractor (Bertin & Arnouts 1996) in combination with the distance from the GLADE catalogfoot_1 to locate all radio sources detected with >5σ significance in each image located within 100 kpc of each of our target galaxies. We chose this search radius because 95% of observed SGRBs are found within 100 kpc of their hosts (Fong & Berger 2013); thus, if GW190814 belongs to the population of mergers that produce SGRBs, then we expect that a real radio counterpart to GW190814 will likely be found within this area. This resulted in an initial list of 72 detected sources in epoch 1 and 102 detected sources in epoch 2. We then visually inspected all images to remove sources that were obvious imaging artifacts or clearly extended, leaving 75 sources detected in at least one epoch. We measured the flux densities of each detected source using the imtool fitsrc command within the pwkit package (Williams et al. 2017). All flux densities were extracted assuming a point-source fit. We found that, due to the differing rms noise level in the epoch 1 and 2 images, some sources that were only recovered by Source Extractor with >5σ significance in one of the two epochs may nevertheless be consistent with constant flux density. For each source detected in only one epoch, we therefore ran fitsrc at the source position in the other epoch to search for lower-significance emission. We then used the recovered flux density (or three times the image rms at the source position, if no emission was detected) to measure the variability of each object.

We assess the variability of our radio sources in several ways. Several previous wide-field searches for radio transients have characterized variability between two epochs in terms of the modulation index, m, and the variability t-statistic, V s (e.g., Mooley et al. 2016;Radcliffe et al. 2019;Bhakta et al. 2021). We define:

where S is the flux density of each source, as determined from a point-source fit with fitsrc, and

where σ is the total measurement error on the flux density. We include both the uncertainty derived from the point-source fit and an additional error term of 5% corresponding to the accuracy of the absolute flux calibration scale of the VLA (Perley & Butler 2017) in this quantity. We plot the distribution of our radio sources in V s and m in Figure 2.

Previous work has often focused on maximizing the purity of constructed samples of radio variables and transients, and thus has imposed fairly strict cutoffs for variability: Mooley et al. 16 Applying these criteria, we find 8 and 13 variable sources, respectively, in our sample, corresponding to respective variability fractions of 11% and 17%. This is significantly higher than previous blind untargeted searches, which have found that only a few percent of radio sources over large sky areas exhibit this level of variability (e.g., Carilli et al. 2003;Frail et al. 2012;Mooley et al. 2016;Bhandari et al. 2018;Dobie et al. 2019;Radcliffe et al. 2019;Sarbadhicary et al. 2020;Bhakta et al. 2021). We note that all 13 variables are coincident with bright galaxies in archival optical imaging, and 12 of them are coincident with the nuclei of our target galaxies to within astrometric errors. One natural explanation for the higher prevalence of variability in our sample compared to previous work is thus that the centers of (relatively) nearby galaxies selected as likely hosts for GW mergers are more likely than average regions of space to contain sources of detectable radio emission (e.g., weak AGN). We explore this further in Section 4.3.

One downside of optimizing for sample purity in our data is that GW radio counterparts are expected over a broad range of timescales: weeks to months for relativistic jet afterglows, or months to years or even decades for KN afterglows (e.g., Nakar & Piran 2011). Thus, with only two epochs of data, it is worth exploring additional, conservative methods to ensure that transients with variability timescales not well-aligned to our observing cadence are also discovered. For example, such techniques would be necessary to recover a GW170817-like radio transient in our 2, magenta star). Even an off-axis relativistic jet that peaks on the timescale of our second epoch may not satisfy the criterion for V s given the sensitivities achieved in our two epochs (e.g., sample model in Figure 2, blue star).

We are thus motivated to explore a less stringent criterion for variability, to emphasize completeness of our sample rather than purity. We therefore create a list of all radio sources inconsistent with constant flux density to within the >1σ measurement uncertainties calculated by fitsrc (using three times the image rms at the source position as an upper limit on the flux density for sources only detected in one epoch). We identified 34 potential radio counterparts using this criterion (Figure 2, black crosses). Of our 34 variable sources, 32 of them (including all 13 of the "highly Note.

a Galaxy selected for additional multifrequency follow-up, based on the presence of a variable radio source within 100 kpc discovered in our observations. (2021). However, these criteria may miss genuine GW counterparts: for example, GW170817ʼs radio counterpart is not obviously highly variable on the timescale of our observations (magenta star). Even jet emission that peaks at the time of our second epoch could be missed if observed with the cadence and sensitivity of our observations; one such model that we considered for GW190814 is shown by the blue star (a tophat jet with E K,iso = 5 × 10 51 erg, n = 0.09 cm -3 , ò e = 0.1, ò B = 0.01, p = 3.2, and viewing angle 45°). We utilize a less-restrictive variability definition in this work, resulting in a larger sample of possibly variable objects (black crosses). We select one of these variables for additional multifrequency follow-up, based on its large modulation index and its lack of an optical counterpart (Candidate 1, red star).

significant variables" satisfying the Bhakta et al. (2021) criteria discussed above) increased in brightness between epochs 1 and 2, while only two sources decreased in brightness. We do not expect this imbalance to result from differences in absolute flux calibration between the two epochs, as the variables were not preferentially observed in any single scheduling block (and as mentioned above, we already include a conservative additional 5% uncertainty term in the measurement errors used to compute V s , based on the known accuracy of the VLA flux calibration scale). Instead, the apparent increase in flux density of many of our sources may be partially explained by the reconfiguration of the VLA between epochs 1 and 2: the C-configuration epoch 2 data are more sensitive to diffuse low surface brightness emission (from, e.g., ongoing star formation) than the A-configuration epoch 1 data. (This is also consistent with the fact that a number of our radio sources appear point-like in epoch 1 and slightly extended in epoch 2, despite the epoch 2 data having ∼10× lower resolution).

The centers of galaxies may be particularly prone to such spurious detections of variability, as they may contain a superposition of extended and compact emission components from star formation and/or AGN activity. Indeed, the positions of 17 of our variable sources are consistent with the centers of their target GLADE galaxies to within astrometric uncertainties. While some models predict an enhanced rate of compact object mergers near the supermassive black holes at the centers of galaxies (Antonini & Perets 2012;McKernan et al. 2020;Perna et al. 2021;Zhu et al. 2021), this result suggests that GW counterparts near galactic nuclei will be particularly challenging to identify in the radio. We rule out these 17 sources as likely radio counterparts to GW190814.

Finally, we searched the VLA Sky Survey (VLASS; Lacy et al. 2020) Quick Look images and the PanSTARRS-1 data archive (PS1; Chambers et al. 2016) at the positions of the remaining candidates, to provide additional insight into their nature. Three sources are detected in VLASS epoch 1 data that predate GW190814, suggesting that they are unrelated to the merger. Furthermore, we found that all but six sources had a spatially coincident optical counterpart in the PS1 catalog, suggesting that these radio sources are also likely not plausible counterparts to GW190814. While the exact peak timescale of a hypothetical radio counterpart to GW190814 is uncertain, models of off-axis relativistic jets with θ obs ∼ 46°(GW190814ʼs binary inclination) and parameters typical of short GRBs are expected to show significant variability between the timescales of our two epochs (Figure 2, blue star). Of our remaining six candidates, the source that shows the largest |m| is located at R.A. = 00 h 52 h 45 418, decl. = -25°43′08 16 (J2000), within the field of the galaxy ESO 474-035, which was ranked sixth in our catalog of 75 galaxies (Figure 1). We discuss this source (hereafter "Candidate 1") and its evolution in more detail in the next section.

Modeling of Candidate 1 in the Field of ESO 474-035

We now focus on Candidate 1, to assess its viability as a counterpart to GW190814. Previous Magellan observations at 1.48 days after the GW trigger cover both the cataloged GLADE galaxy ESO 474-035 and the position of Candidate 1, and placed a limit of i > 22.23 mag on any transient optical emission at that time (Figure 1; see also Gomez et al. 2019). Moreover, deeper, pre-merger limits at the candidate position from the Legacy Survey (Dey et al. 2019) of r = 24.92 mag, derived from the 5σ depth of co-added DECam images (Brick ID 0131m257), place stringent constraints on an underlying source. For a satellite galaxy at the same distance as ESO 474-035, this translates to a constraint of L  4.1 × 10 6 L e . This is roughly four orders of magnitude below the luminosity of the Milky Way, ruling out all except the faintest dwarf galaxy regime (Simon 2019). Moreover, the lowest host galaxy luminosities derived for short GRBs are ≈10 9 L e (Berger 2014), well above the limit derived here. If instead there is a background galaxy at a luminosity similar to that of the Milky Way at this position, it would need to be at z  2 to be consistent with this limit, implying a luminosity of 3 × 10 41 erg s -1 for the radio transient, comparable to radioloud quasars (Kellermann et al. 2016). Thus, we find that, while a background quasar cannot be ruled out, if the origin of Candidate 1 is from a stellar progenitor, then it likely originated from ESO 474-035.

We triggered additional multifrequency radio observations of Candidate 1 at 1-12 GHz, which were carried out on 2020 May 5 (t = 266 days, "epoch 3"), with the VLA in C configuration. Our observations reveal a continued increase in the flux density at 6 GHz (Figure 3, left), confirming that the change in flux density is intrinsic to the source, rather than an artifact of the VLA configuration change between epochs 1 and 2. The broadband spectrum is optically thin, consistent with a single power law F ν ∝ ν β , where β = -1.1±0.06 (Figure 3, right). At the highest frequencies observed (8-12 GHz), we note that the emission appears partially resolved into two components, with centroids separated by ∼3″. At the distance of ESO 474-035 (271 Mpc), this would correspond to a physical separation of ∼4 kpc, orders of magnitude too large for a newly formed GRB jet. This may suggest that our candidate is instead a background double-lobed radio AGN jet undergoing a flaring event.

AGNs are one of the most common types of compact sources in the radio sky, and in the radio their flares have typical timescales of months to a few years (e.g., Hovatta et al. 2008). Thus, AGN flares are likely to be the largest source of contamination in searches for radio GW counterparts. We note that, prior to obtaining data with better spatial resolution in epoch 3, Candidate 1 initially appears very similar to expectations for a GW counterpart, and we may expect to find similar contaminants in future GW counterpart searches where high-resolution data are not available. We therefore invest some additional effort in fitting Candidate 1ʼs radio emission with models appropriate for GW counterparts, to see if we can distinguish Candidate 1 from a typical radio GW counterpart based on the physical properties required to fit the flux density alone. For this analysis, we use the combined flux density of the two resolved components when modeling the high-frequency epoch 3 data, and we assume that Candidate 1 is at the distance of ESO 474-035.

We consider two classes of radio GW counterpart models for Candidate 1: radio emission from collimated fast ejecta (i.e., an initially relativistic jet, possibly with some velocity structure) and from slower ejecta (the same material in which r-process nucleosynthesis occurs at early times, producing the KN optical transient; assumed to be quasi-spherical and moving at up to a few tens of percent of c). In both cases, the radio emission is synchrotron radiation arising from a population of electrons accelerated into a power-law distribution of energies, N(γ) ∝ γ -p for γ > γ m , as the merger ejecta shocks and interacts with the surrounding interstellar medium. The allowed parameter space is highly degenerate, as we observe only the rising portion of the light curve and a single power-law segment of the synchrotron spectral energy distribution. We find p = 3.2 ± 0.11 using the spectral slope computed from our multifrequency observations (β = -1.1±0.06), assuming the radio observations lie above ν m (the synchrotron frequency corresponding to γ m ) and below the cooling frequency (ν c ) (Granot & Sari 2002). This value is in line with some supernovae (see, e.g., Van Dyk et al. 1994;Chevalier 1998), although it disagrees with the precise value p 2.15 0.02 0.01 = -+ calculated for GW170817 (Hajela et al. 2019) and it is higher than previously observed in many cosmological SGRB afterglows, where typically 2 < p < 3 (Fong et al. 2015).

We first consider the possibility that Candidate 1ʼs radio emission is due to an off-axis relativistic jet. We fix p = 3.2 and θ obs = 46°(to put the jet in alignment with the best-fit binary inclination as derived from the GW signal; see Abbott et al. 2020b), and we assume that the fractions of energy carried by electrons and magnetic fields in the shock are ò e = 0.1 and ò B = 0.01, respectively. While this agrees with assumptions made in previous studies of cosmological SGRBs (e.g., Fong et al. 2015) and with previous theoretical and observational studies that find ò e ≈ 0.1 for relativistic shocks (e.g., Spitkovsky 2008; Beniamini & van der Horst 2017), we note that for most events at cosmological distances p, ò e and ò B are poorly constrained by the data directly and must be assumed. The full GW170817 data set suggests that ò B may be much lower than 0.01 in at least some jets (Margutti & Chornock 2021); if this is true for GW190814, then the constraints derived below are also affected. For example, models with lower values of ò e or ò B will have a larger total energy. We consider both tophat jets (in which all of the jet energy is contained within a narrow cone with opening angle 15°) and jets with a relativistic core surrounded by Gaussian wings of slowermoving material. We find that only the tophat jet models can reproduce the steep rise seen in Candidate 1ʼs 6 GHz light curve, and furthermore that we require a large isotropic-equivalent jet energy (E K,iso ∼ 8 × 10 53 erg) and a high density (n ∼ 0.5 cm -3 ); comparable to the largest energy values and in the top 40% of density values inferred for SGRBs (Fong et al. 2015). While we cannot entirely rule out the possibility that Candidate 1 is an offaxis relativistic jet launched by GW190814, we disfavor this possibility due to the large energy required and the high value of p. The high density would also be unexpected for such a highly offset transient, particularly due to the lack of any optical emission at the transient position to suggest a satellite galaxy or globular cluster environment for the transient.

We next compare the radio behavior to slow KN ejecta models, to determine whether they are consistent with our observations. Following the prescriptions of Schroeder et al. (2020), we modeled the 6 GHz lightcurve and multifrequency spectrum of Candidate 1 with a KN ejecta interaction model. We again fix p = 3.2, ò e = 0.1, and ò B = 0.01 in our modeling. Optical and near-infrared follow-up studies of GW190814 have placed constraints on the ejecta mass of M ej  0.04-0.1 M e (Gomez et al. 2019;Ackley et al. 2020;Andreoni et al. 2020;Kawaguchi et al. 2020;Morgan et al. 2020;Thakur et al. 2020;Vieira et al. 2020). We set the ejecta mass M ej to 0.01M e , as lower values of M e would cause the time of observation of epoch 3 to approach the deceleration time, t dec , of the ejecta. If t dec ≈ t obs,3 , the time of observation of epoch 3, the KN model light curve would start to decline, whereas we observe the 6 GHz light curve still rising through epoch 3.

Even after making these assumptions, several parameter degeneracies remain in our modeling. We therefore create a grid of 250 light-curve models exploring a range of values for the density n and ejecta energy E ej . We find a broad range of combinations that are consistent with our observations at the times of the two 6 GHz detections (Figure 4, solid line). Densities above n ∼ 1 cm -3 are ruled out, as these models would also require the light curve to begin declining by the time of our epoch 3 observation. As with the relativistic jet models, the ejecta energy required to match the observations is large, particularly for the lower densities expected given the large offset of Candidate 1. These models also require very high ejecta velocities, β 0 ∼ 0.80-0.94c for n ∼ 0.01-10 -4 cm -3 (Figure 5, red dashed line). While the distributions of ejecta velocities in some published KN models have tails to high velocities ∼0.9c (e.g., Figure 5, gray lines), the bulk of the energy must be carried by slower-moving ejecta (≈0.1-0.5c), to match the optical and infrared properties of the KN emission (Bauswein et al. 2013;Hotokezaka et al. 2013;Sekiguchi et al. 2016Sekiguchi et al. , 2016;;Ciolfi et al. 2017;Mooley et al. 2018a). The ejecta energy is only mildly sensitive to the assumed values of ò e and ò B ; if energy equipartition is assumed (ò e = ò B = 0.33), the required E ej decreases by a factor of 2-3 across the range of allowed densities. Conversely, models with lower values of ò e or ò B require a higher E ej . It is therefore difficult to explain the high ejecta energy required to fit our radio light curve with existing KN models.

In summary, both relativistic jet models and quasi-spherical KN ejecta models ultimately require very high energies and fast ejecta velocities to match the radio evolution of Candidate 1. A comparison to ejecta models consistent with GW170817 and cosmological SGRBs can be seen in Figure 5. It is clear that either the KN models used to explain Candidate 1 are probing a new regime of parameter space, one with higher energies and velocities than the other models that have been found to be consistent with the population of cosmological SGRBs and GW170817, or that some parameters (e.g., ò e and ò B ) differ from the standard values we assumed. Nevertheless, even for higher ò B values, it is difficult to construct a plausible physical scenario that accelerates sufficient ejecta to such high velocities, particularly for compact object mergers with highly unequal mass ratios like GW190814. Ultimately, we conclude that this analysis disfavors Candidate 1 as a radio counterpart to GW190814 on physical grounds. We suggest that similar analyses can be applied in future GW counterpart searches to discriminate between genuine radio GW counterparts and unrelated background AGN.

Discussion

Limits on a Relativistic Jet Launched by GW190814

Apart from Candidate 1, we detect no other convincing radio counterpart to GW190814 in our observations. We therefore next investigate the limits we can place on the existence of a relativistic jet (assuming that the merger occurred within one of the galaxies we targeted). We generate a grid of light curves at 6 GHz using afterglowpy (Ryan et al. 2020) for two relativistic jets: (1) a tophat jet with an opening angle θ jet = 15°, and (2) a Gaussian jet with a 15°jet core (θ jet ) and wings extending out to 6 θ jet . Similar jet opening angles are seen in many GRBs (Ryan et al. 2015). We compute the light curves for isotropic-equivalent kinetic energies, E K,iso = 2 × 10 51 , 5 × 10 51 , and 5 × 10 52 erg (representing a typical SGRB energy, as well as two more optimistic cases), over 0.1-1000 days. These light curves are computed on a fine density grid with varying circumburst density n = 10 -6 -10 3 cm -3 , and for viewing angles ranging from 0°(on-axis) to 90°. We fix p = 2.2, ò e = 0.1, ò B = 0.01, and d 241

Mpc. We then compare these light curves to the limits placed by our VLA observations at 38 and 208 days, using 5× the typical rms of our images in each epoch as an upper limit on the flux density of any counterpart at that time (Figure 6). These limits allow us to rule out the higher-density parameter space, as higher densities result in a brighter source, which would have been detected in our observations. Alternately, if the value of ò B is significantly lower than we assumed, then higher-density models may still be allowed by our data.

Compared to the best-fit viewing angle of θ obs = 46°calculated by LIGO/Virgo, we find our observations rule out densities n  0.4, 0.3, 0.02 cm -3 assuming a tophat jet model for E K,iso = 2 × 10 51 , 5 × 10 51 , and 5 × 10 52 erg, respectively. Compared to the circumburst densities found for the SGRB population (assuming ò e = 0.1 and ò B = 0.01, for consistency; see Fong et al. 2015), these constraints are comparable to the higher end of the population (at the ∼60%-85% level compared to SGRB densities). These limits are a strong function of the viewing angle; for the jet with E K,iso = 2 × 10 51 erg, the lower limit on the density ranges from 0.06-4 cm -3 for the full range of viewing angles 35°-64°a llowed by the GW analysis. For Gaussian jet models with the same energies, we find more constraining densities of n  0.1, 0.03, 9 × 10 -4 cm -3 (at the ∼40%-65% level of SGRBs) for θ obs = 46°(n  0.03-1 cm -3 for the 2 × 10 51 erg jet observed at θ obs = 35°-64°at d L = 241 Mpc). Thus, the assumed jet structure (tophat versus Gaussian) can affect the implied limits on circumburst density by a factor of a few. Our limits additionally become less constraining for more narrowly beamed jets (which peak earlier); we can rule out a tophat jet with E K,iso = 2 × 10 51 erg The allowed E ej vs. n phase space for Candidate 1, under the assumption that the radio emission is produced by the shock between 0.01M e of quasi-spherical ejecta and the ambient medium (with p = 3.2, ò e = 0.1, and ò B = 0.01). Densities higher than 1 cm -3 are ruled out by the 6 GHz light curve.

Figure 5. The E ej vs. Specific Momentum (Γβ) phase space for KN models consistent with Candidate 1ʼs radio evolution (red dashed line) in comparison to other ejecta models proposed for compact object mergers. Both quasispherical KN ejecta models and relativistic ejecta models struggle to reproduce the high energy required to match Candidate 1ʼs radio properties. The orange circles show the energy of the red, blue, and purple KN components associated with GW170817 from Villar et al. (2017), while the gray lines show two different models for the velocity distribution of quasi-spherical ejecta in this event (Mooley et al. 2018a). The purple shaded region is a representative range of maximum energies found for SGRB slow ejecta derived from late-time radio observations (Schroeder et al. 2020). The green line is the structured jet model for GW170817 from Margutti et al. (2018), while the blue shaded region is the beaming-corrected energy of the jet component in SGRBs (Fong et al. 2015). and θ jet = 10°viewed at θ obs = 46°for densities n > 1 cm -3 . In comparison, the search conducted by Dobie et al. (2019) can rule out comparable densities for a tophat jet with E K,iso = 10 51 erg, θ jet = 10°, ò e = 0.1, ò B = 0.01, and p = 2.2 only if the jet is viewed at 35°off-axis; for a jet that is 46°off-axis, the density is not constrained by their data. This emphasizes the importance of continuing radio transient searches to late times (6 months postmerger) to fully constrain the allowed parameter space for highly off-axis jets, matched to their later peak timescales.

Properties of Our Most Variable Sources: Extreme AGN

Flares?

The differing resolution and sensitivity of our data may impact the interpretation of sources with both small and large apparent flux density changes on the timescale of our observations, as discussed in Section 3.1. Nevertheless, five of of our 13 "highly variable sources" (selected using the definition of Bhakta et al. 2021) have strong detections in VLASS data predating 2019 August 18, and three have marginal detections, confirming that they are likely unrelated to GW190814. The brightest of these radio sources is coincident with ESO 474-026 (the most highly ranked galaxy on our prioritized list, at a distance of 244 Mpc), and appears strongly point-like in both our observations and the VLASS quicklook data. We therefore suggest that its variability is intrinsic, and that this radio source is a compact AGN. Interestingly, the flux density of this source increases dramatically, from 0.34 ± 0.03 mJy in epoch 1 to 1.31 ± 0.14 mJy in epoch 2. This factor of ∼4 increase in flux density over a timescale of ∼5 months is not unprecedented for AGN flares, although their typical fractional variability is much lower at 6 GHz (from a few percent to a factor of 2; see Hovatta et al. 2008). Additionally, Cox & Sparke (2004) reported a 6 GHz flux density of 2.5 ± 0.1 mJy for this source in VLA observations taken on 1994 October 2 at similar resolution (5″), suggesting that its large-amplitude radio variability is not unique to the period of our observations. ESO 474-026 has been classified as a polar-ring galaxy; its unusual optical morphology is most likely a short-lived state produced by a recent major galaxy merger (Reshetnikov et al. 2005;Spavone et al. 2012). We suggest that this merger event could have fueled the AGN responsible for the variable radio emission.

Some models predict an enhanced rate of BNS and NSBH mergers in the accretion disks surrounding AGN (McKernan et al. 2020;Perna et al. 2021;Zhu et al. 2021), but the resulting radio transient would be difficult to disentangle from the type of variability seen in ESO 474-026. Our results suggest that mergers in active galaxies will be particularly hard to discover and model in radio-only data sets, as radio AGN variability is still poorly constrained on the relevant timescales. This highlights the importance of continued monitoring of such sources, given the wide range of timescales expected for extragalactic radio transients. Ongoing and planned all-sky radio surveys are beginning to provide this vital long-term coverage; this particular source will next be observed in 2022 February during epoch 2.2 of VLASS.

Considerations for Future Galaxy-targeted Radio

Counterpart Searches

Finally, we briefly expand upon some implications of the high variability fraction of our sample of radio sources. Compared to untargeted searches, nearby galaxies are likely to have a higher surface density of detectable radio transients with unresolved, compact emission, including radio supernovae, tidal disruption events (TDEs), and AGN flares (e.g., Weiler et al. 1986Weiler et al. , 2002;;Berger et al. 2003;Hovatta et al. 2008;Romero-Cañizales et al. 2011;Alexander et al. 2015;Irwin et al. 2015;Metzger et al. 2015;Anderson et al. 2019;Figure 6. Portions of parameter space ruled out by the nondetection in our two epochs at ∼38 and ∼208 days for a 15°tophat (left) and Gaussian (right) jet, with E K,iso =(0.2, 0.5, and 5) × 10 52 erg, p = 2.2, ò e = 0.1, and ò B = 0.01. The shaded region shows the uncertainty from the distance. All space above and to the left of these limits is ruled out by our nondetections. Under these assumptions, we can rule out an energetic jet similar to those seen in cosmological short GRBs for a range of viewing angles and densities. The binary inclination angle and its associated uncertainty calculated from the full GW analysis (Abbott et al. 2020b) are shown by the dashed lines. Alexander et al. 2020), in addition to more diffuse emission associated with star formation. At 241 Mpc, our 5σ sensitivity corresponds to a radio luminosity of 4 × 10 27 erg s -1 Hz -1 in epoch 1 and 6 × 10 27 erg s -1 Hz -1 in epoch 2. This is sufficient to detect all known radio TDEs (Alexander et al. 2020) and the brightest ∼7% of core-collapse supernovae (Bietenholz et al. 2021). Given the small number of galaxies we observed, their fairly proximal distances, and the rates of transients per L * galaxy (Li et al. 2011;Metzger et al. 2015;van Velzen et al. 2020), we expect to find <1 serendipitous supernova or TDE in our search. However, our observations demonstrate that the contamination risk from AGN is much higher for observations with a cadence of several months, in agreement with previous work that shows radio AGN flares have a typical variability timescale of months to a few years (Hovatta et al. 2008;Richards et al. 2011).

Many galaxy-targeted searches for GW counterparts, including ours, use a prioritization scheme that ranks galaxies in part by their optical luminosity, with the intent of maximizing the fraction of stellar mass in the localization volume that can be observed with a finite amount of telescope time (e.g., Gehrels et al. 2016;Arcavi et al. 2017b;Ducoin et al. 2020). For radio interferometers in which slew times can be significant, observing strategies typically maximize a combination of probability covered and slew (Rana & Mooley 2019), but do not necessarily address the issue of heightened contamination in galaxies. By default, such galaxyranking algorithms also have the effect of prioritizing exactly the galaxies that are most likely to have detectable radio emission from unrelated processes such as star formation or AGN activity, even at distances of a few hundred Mpc. Indeed, we find that 23 of our target galaxies (31% of our sample) have detected radio emission coincident with the galaxy nucleus, and 17 of these nuclear radio sources are initially flagged as variable using our transient search criteria (i.e., we find that 74% of these nuclear sources exhibited significant variability on timescales of 5-6 months). We therefore conclude that identifying radio counterparts in or near the nuclei of the most massive galaxies will be challenging for low-cadence observations-particularly as most existing all-sky radio surveys are fairly shallow, so there are unlikely to be sufficiently deep premerger radio observations of the target galaxies to compare against.

We also find that, for radio variability searches conducted with the VLA or another reconfigurable radio interferometer, the array configuration needs to be carefully considered when planning GW follow-up and interpreting the data. As mentioned above, many of the most highly ranked galaxies in a GW localization volume may have detectable radio emission from other, unrelated sources; thus, lower-resolution observations taken in the VLA's C and D configurations (such as our epoch 2 observations) may suffer from issues of source confusion. In addition, while the VLA is equally sensitive to emission from point-like unresolved sources in all configurations, it is more sensitive to diffuse emission in its most compact configurations. Thus, sources identified as point-like in our high-resolution epoch 1 A-configuration data may in fact have more extended components that contribute to the fitted flux density measured in our lower-resolution epoch 2 C-configuration data (even when forcing a point-source fit), giving the appearance of variability. We see some evidence for this in our data, complicating our efforts to determine how much of the variability is intrinsic for these sources. This reinforces the importance of using variability selection criteria tailored for galaxy-targeted searches, and taking extra care near the nuclei of target galaxies where such extended emission is most likely to be detected. Deep template images of each galaxy in the relevant configuration(s) would be necessary to attempt to deconvolve a potential near-nuclear GW counterpart from the background variability of its host.

Conclusions

We carried out the first galaxy-targeted search for a radio counterpart to a GW merger event, GW190814. Although we detected several transient or variable sources, all are consistent with AGN variability or are artificial variables created by the different uv coverage of our two epochs of data, and are thus unlikely to be genuine radio counterparts to GW190814. Via additional monitoring of one initially promising candidate, we demonstrate that multifrequency radio observations can help distinguish background AGN flares from bona fide radio GW counterparts, as they may have different spectral indices and/or require unphysical parameter values to fit the shape of the radio light curve in the context of relativistic jet or KN afterglow models. For the 75 galaxies that we observed, comprising 32% of the stellar luminosity in the final localization volume, we can rule out a relativistic jet at the best-fit LIGO/Virgo viewing angle of ∼46°with isotropic-equivalent energies E K = 2 × 10 51 , 5 × 10 51 , and 5 × 10 52 erg propagating in an ISM-like constant-density medium of n  0.4, 0.3, 0.02 cm -3 for a tophat jet model, or n  0.1, 0.03, 9 × 10 -4 cm -3 for a Gaussian jet model (assuming θ jet = 15°, p = 2.2, ò e = 0.1, and ò B = 0.01). These limits will change if different values are assumed for θ jet , p, and the microphysical parameters; thus, care must be taken when cross-comparing the limits set by different studies in the literature.

Our results have a number of implications for future radio searches for GW counterparts and other radio transients. In particular, we find that searches specifically targeting nearby galaxies (where a compact object merger is most likely to occur) may encounter additional complexities. We identify a significantly higher fraction of variable radio sources in our galaxy-targeted search, in comparison to previous wide-field radio transient searches (e.g., Carilli et al. 2003;Frail et al. 2012;Mooley et al. 2016;Radcliffe et al. 2019;Sarbadhicary et al. 2020;Bhakta et al. 2021). This is likely due to a combination of several factors, some of which are relevant for any galaxy-targeted search (e.g., contamination from unrelated AGN variability), and some of which are specific to the VLA (e.g., the complications incurred by the regular reconfiguration of the VLA, which results in a nonuniform data set for any follow-up campaign lasting 3 months).

As Advanced LIGO, Virgo, Kagra, and future generations of GW detectors continue to improve their sensitivity, wide-field and galaxy-targeted radio counterpart searches will also need to improve their sophistication. In the future, sensitive, highresolution, wide-field radio facilities like ASKAP (Dobie et al. 2019), MeerKAT, and the Square Kilometer Array will likely be able to overcome some of the challenges encountered in this work (particularly if higher-frequency receivers are added). Nevertheless, our pilot study demonstrates that there is a niche for galaxy-targeted radio searches for nearby events with the VLA. Our observations were sufficient to rule out most plausible SGRB-like jet models for the galaxies we observed, with a modest investment of only 23 hr of telescope time. In the future, we would suggest adding additional epochs to enable at least two observations of each target galaxy per array configuration, to minimize the resolution effects encountered in this study and to better discriminate between AGN variability and the smooth rise and decline expected for a GW radio counterpart. To prevent the total time investment from being prohibitive, we also suggest that this strategy only be applied to events that are closer and/or better localized than GW190814, such that the total number of galaxies in the localization volume is 50. Such events will remain rare (we expect at most one such BNS or NSBH merger in O4; see Abbott et al. 2020a), but are nevertheless likely to be among the best-studied mergers. The increasing availability of deep template maps of the radio sky will also improve our ability to interpret the results of future radio counterpart searches, just as we found VLASS observations of our target galaxies to be useful in this work. Radio searches will remain the only way to discover EM counterparts to compact object mergers that occur in the daytime sky, whose optical and X-ray emission cannot be studied, and will thus remain an important tool for multimessenger studies, despite their challenges.