Continued radio observations of GW170817 3.5 years post-merger

We present new radio observations of the binary neutron star merger GW170817 carried out with the Karl G. Jansky Very large Array (VLA) more than 3\,yrs after the merger. Our combined dataset is derived by co-adding more than $\approx32$\,hours of VLA time on-source, and as such provides the deepest combined observation (rms sensitivity $\approx 0.99\,\mu$Jy) of the GW170817 field obtained to date at 3\,GHz. We find no evidence for a late-time radio re-brightening at a mean epoch of $t\approx 1200$\,d since merger, in contrast to a $\approx 2.1\,\sigma$ excess observed at X-ray wavelengths at the same mean epoch. Our measurements agree with expectations from the post-peak decay of the radio afterglow of the GW170817 structured jet. Using these results, we constrain the parameter space of models that predict a late-time radio re-brightening possibly arising from the high-velocity tail of the GW170817 kilonova ejecta, which would dominate the radio and X-ray emission years after the merger (once the structured jet afterglow fades below detection level). Our results point to a steep energy-speed distribution of the kilonova ejecta (with energy-velocity power law index $\alpha \gtrsim 5$). We suggest possible implications of our radio analysis, when combined with the recent tentative evidence for a late-time re-brightening in the X-rays, and highlight the need for continued radio-to-X-ray monitoring to test different scenarios.


INTRODUCTION
GW170817 has been a milestone event for transient multi-messenger studies. It was the first binary neutron star (NS) merger observed by the LIGO and VIRGO detectors (Abbott et al. 2017), and so far it remains the only binary NS system from which gravitational waves (GWs) and a multi-wavelength (radio to gamma-ray) counterpart have been discovered (Abbott et al. 2020;Kasliwal et al. 2020;Paterson et al. 2021). The GW170817 NS-NS merger occurred at 12:41:04 on 2017 August 17 UTC, and its GW detection was followed by the detection of a γ-ray burst (GRB) by the Fermi and INTEGRAL satellites, ≈ 2 s after the merger. UV/optical/IR instruments subsequently identified the so-called kilonova counterpart (AT2017gfo), in the galaxy NGC 4993 at a distance of ≈ 40 Mpc, making GW170817/GRB170817a the closest short GRB with known redshift (e.g., Arcavi et al. 2017;Chornock et al. 2017;Coulter et al. 2017;Cowperthwaite et al. 2017;Drout et al. 2017;Kasliwal et al. 2017;Kasen et al. 2017;Kilpatrick et al. 2017;Pian et al. 2017;Shappee et al. 2017;Smartt et al. 2017;Tanvir et al. 2017;Valenti et al. 2017). Observations of the quasi-thermal UV/optical/IR emission from the GW170817 slow (∼ 0.1c − 0.3c), neutron-rich, kilonova ejecta were successful in verifying that mergers of NSs in binaries are production sites of heavy elements such as gold and platinum (e.g., Kasliwal et al. 2017;Kasen et al. 2017;Pian et al. 2017;Metzger 2017).
In addition to the quasi-thermal kilonova emission, a delayed non-thermal (synchrotron) afterglow from GW170817/GRB 170817a was first observed in the Xrays ≈ 9 days after the merger by the Chandra observatory (e.g., Troja et al. 2017;Haggard et al. 2017;Margutti et al. 2017). A radio afterglow detection with the Karl G. Jansky Very Large Array (VLA) followed, about two weeks after the merger (Hallinan et al. 2017). Further radio observations of the source proved decisive in narrowing down the morphology of the jet (Corsi et al. 2018;Dobie et al. 2018;Alexander et al. 2017;Margutti et al. 2018;Mooley et al. 2018b,a;Hajela et al. 2019), ruling out the simple uniform energy-velocity (top-hat) ejecta in favour of a structured jet, where the ejecta velocity varies with the angle from the jet axis (Mooley et al. 2018c;Lazzati et al. 2018;Ren et al. 2020;Ghirlanda et al. 2019). These observations, with the help of hydrodynamic simulations Nakar et al. 2018), set constraints on the opening angle of the jet core ( 5 deg), the observer's viewing angle (≈ 15 − 30 deg), the isotropic equivalent energy (∼ 10 52 erg) and the interstellar medium (ISM) density (∼ 10 −4 − 0.5 cm −3 ).
The extended radio follow-up of GW170817 up to 2.1 years after the merger had shown that the radio emission from the structured jet had faded below typical flux density sensitivities that can be reached with the VLA in a few hours of observing (Makhathini et al. 2020). Several theoretical scenarios, however, predict the possible emergence at late times of detectable electromag-arXiv:2103.04821v3 [astro-ph.HE] 18 Jun 2021 netic emission associated with the afterglow of the kilonova ejecta itself (e.g., Nakar & Piran 2011;Piran et al. 2013;Hotokezaka & Piran 2015;Kathirgamaraju et al. 2019;Hotokezaka et al. 2018;Margalit & Piran 2020). Indeed, numerical simulations show that during the NS-NS merger, a modest fraction of a solar mass is ejected from the system, and the total ejecta mass and velocity distribution of such ejecta depend on the total mass, mass ratio, and the nuclear equation of state (EoS) of the compact objects in the binary. While optical-UV observations are mostly sensitive to the low-end of the ejecta velocity distribution, radio (and X-rays) can probe the fastest moving ejecta tail, shedding light on whether NS-NS ejecta are broadly distributed in energy and velocity (as simulations seem to suggest) and providing indirect constraints on the nuclear EoS.
Tentative evidence for a very late-time re-brightening possibly associated with the kilonova afterglow of GW170817 has recently been reported in the X-rays (Hajela et al. 2020a(Hajela et al. , 2021Troja et al. 2020). On the other hand, relatively shallow observations (rms ≈ 4.3 µJy) with the VLA at 3 GHz had reported a lack of radio detection contemporaneous with the X-ray late-time rebrightening (Alexander et al. 2020). This radio nondetection was interpreted to be compatible with expectations from the simplest extrapolation of the X-ray excess to the radio band (Alexander et al. 2020). Here, we use much deeper VLA observations to show that the lack of a detection of radio emission in excess to that expected from the structured jet associated with GW170817 constrains kilonova ejecta models to a relatively steep energy-velocity ejecta profile.
Our work is organized as follows. We report our new observations in §2; in §3 we discuss our results within the kilonova afterglow model; finally, in §4 we conclude with a summary.

OBSERVATIONS AND DATA REDUCTION
As we describe in what follows, we have processed new and archival data of GW170817 obtained at radio and X-ray wavelengths between 2020 September and 2021 February. The full panchromatic afterglow light curve of GW170817 is available for download on the web 7 .

Radio Observations
Radio observations of the GW170817 field were carried out with the VLA on the dates listed in Table 1 at S band (2-4 GHz, nominal central frequency of 3 GHz) with the array in its B (September 2020) and A configurations (December 2020 -February 2021). Each observation was calibrated in CASA (McMullin et al. 2007) using the automated VLA calibration pipeline. The calibrated data were then manually inspected for further radio frequency interference (RFI) excision. We interactively imaged each observation using the CASA task tclean with one Taylor term (nterms=1) and robust weighting (robust=0.5), and derived the rms measurements using imstat. Table 1 lists the rms sensitivity reached in each observation, estimated within a region of 20 synthesized beams around the position of GW170817 (α = 13h09m48.069s, δ = −23d22m53.39s, J2000; Mooley et al. 2018c). We find no significant (> 3×rms) excess in a region of one synthesized beam around the position of GW170817 in any of the individual images. Next, we co-added and interactively imaged as above (nterms=1 and robust=0.5) all A-configuration 3 GHz VLA observations listed in Table 1. An rms of 1.3 µJy was reached at 2.8 GHz (Table 2) within a region of size approximately equal to 20 synthesized beams centered on the position of GW170817. Within one synthesized beam centered on the location of GW170817, we measure an imstat peak flux density value of ≈ 2.8 µJy at 2.8 GHz. With this procedure, several of the bright, extended sources present in the field left substantial deconvolution residuals. Thus, to mitigate the effects of deconvolution residuals, test the robustness of our measurement, and further improve our sensitivity, we imaged all 3 GHz data (both A and B configurations of the VLA) listed in Table 1 non-interactively with two Taylor terms (nterms=2), robust weighting (robust=0.5), single phase-only selfcal (solution interval of 4 minutes), and a cleaning threshold of 4 µJy. This yielded an image rms noise of 0.99 µJy (in a region of size equal to 20 synthesized beams around the position of GW170817; theoretical thermal noise ≈ 0.85 µJy) and a peak flux density value of 2.86 µJy (Table 2) within one synthesized beam centered on the location of GW170817.
A late-time VLA observation of the GW170817 field was also carried out in U band (nominal central frequency of 15 GHz) on 10 February 2021 (Table 1). We calibrated this dataset and interactively imaged the field (with nterms=1 and robust=0.5). No significant emis-sion is found at the location of GW170817 (Table 2). The rms measured in a region of size equal to 20 synthesized beams centered around the position of GW170817 is of ≈ 1.9 µJy at 15 GHz.

X-ray Data
We reprocessed and analyzed the Chandra ACIS-S observations of the GW170817 field obtained between 2020 December 10 and 2021 January 27 (obsIDs 22677, 24887, 24888, 24889, 23870, 24923, 24924;150.5 ks, PI Margutti) using the same procedure described in Makhathini et al. (2020). We find an unabsorbed flux density of (1.70±0.45)×10 −4 µJy at 2.4×10 17 Hz (1 keV; 1σ uncertainty) by combining the spectral products of all seven observations and fitting the data with an absorbed power-law model where the hydrogen column density N H has been fixed to the Galactic value, and where the photon index Γ (N Eν ∝ E −Γ ν ,where E ν = hν) has been fixed to 1.58 (Makhathini et al. 2020). To investigate if Γ is different from 1.58, we refitted the Chandra data leaving Γ as a free parameter. From the 2020 December-January 2021 data we find Γ = 2.71 +1.43 −1.08 . If we additionally combine the 96.6 ks of data obtained in 2020 March, we get Γ = 1.84 +0.80 −0.73 (90% uncertainties). Hence, in both cases, the value of Γ is consistent (well within the 90% confidence interval) with Γ = 1.58. Our results are also consistent with Hajela et al. (2020aHajela et al. ( , 2021; Troja et al. (2020Troja et al. ( , 2021.

DISCUSSION
As is evident from Figure 1, our late-time radio observations of GW170817 do not provide evidence for radio emission in excess to what is expected from the very late time tail of a structured jet afterglow model (black solid line). The X-ray observations, on the other hand, suggest a more pronounced statistical fluctuation or the possible emergence of a new component at higher frequencies (bottom panel in Figure 1).
The radio-to-X-ray spectral index as derived from the measured ratio of the late-time radio-to-X-ray flux densities (see Table 2) is Γ radio−X +1 = −0.535±0.024, within ≈ 2.1σ of the value adopted in Figure 1 (−0.584 ± 0.002) derived by Makhathini et al. (2020). The current measurements still carry too large uncertainties for claiming any clear evidence for a change in spectral behavior at late times. It may be possible however that the nondetection of a radio rebrightening associated with a kilonova afterglow, together with the tentative X-ray excess, suggest a flattening of the radio-to-X-ray spectrum at late times.
In general terms, it is not difficult to envision a scenario in which the electron index (p) for the ejecta responsible for the late-time X-ray excess differs from the one used to model the structured jet afterglow at earlier times (p = 2.07−2.14; Makhathini et al. 2020). The predictions of Fermi particle acceleration imply that the power-law index p expected at non-relativistic shock speeds is close to p ≈ 2, while at ultra-relativistic velocities one can have p ≈ 2.2 (e.g. Sironi et al. 2015, and references therein). Thus, a flattening of the radio-to-X-ray spectral index in GW170817, if confirmed by further follow-up, could support the idea that a non-relativistic ejecta component is starting to dominate the emission.
Differently from the above theoretical predictions for particle acceleration, non-relativistic ejecta observed in radio-emitting core-collapse supernovae typically have p = 2.5 − 3.2 (e.g., Chevalier 1998), pointing to steeper radio-to-X-ray spectra than that suggested by the X-ray excess observed in GW170817. Cases where a transition of the ejecta from the relativistic to the non-relativistic regime has been observed include GRB 030329 (e.g., Frail et al. 2005;van der Horst et al. 2008) and TDE Swift J1644+57 (Cendes et al. 2021). These two cases pointed to a slowly-increasing or constant value of p in the relativistic-to-non-relativistic transition, which again would correspond to a spectral steepening rather than a spectral flattening.
One key difference between GW170817 and other nonrelativistic flows is that we have seen a single powerlaw spectrum from radio, optical, and X-rays, i.e. synchrotron radiation from electrons with different Lorentz factors in many orders of magnitude. In this case, the spectral index should represent p. On the other hand, in other non-relativistic flows, we often determine p from the radio spectra closer to the minimum Lorentz factor γ m of the electron energy distribution, where the synchrotron emission may be dominated by thermal electrons around the typical Lorentz factor rather than accelerated electrons (e.g., Park et al. 2015;Maeda 2013). Thus, GW170817 offers an opportunity to test particle acceleration theory, and continued monitoring from radio-to-X-rays is key to this end. We also note that continued X-ray observations may offer an opportunity to probe the evolution of the cooling frequency in the Newtonian limit, when ν c ∝ β −3 t −2 , and thus constrain β (velocity of the ejecta in units of c) and the kinetic energy of the fast tail of the ejecta (see also Linial & Sari 2019).
An alternative explanation for the excess in X-rays (as compared to the radio) observed in GW170817 might be the possibility of a Compton echo of the X-rays from the prompt emission of GRB 170817A, scattering off surrounding dust (Beniamini et al. 2018). Given currently large uncertainties in the X-ray result, hereafter we focus on the constraints that the lack of a radio excess set on kilonova ejecta models.
Following Kathirgamaraju et al. (2019), the kilo-  Nynka et al. (2018), together with our latest measurement in the radio (3 GHz, latest yellow data point in the grey, shaded region) and X-rays (latest purple data point in the grey, shaded region) extrapolated to 3 GHz using the spectral index derived in Makhathini et al. (2020). The best fit structured jet model for GW170817 is also plotted (top panel, black line) along with the associated 1σ error region (blue shaded region). As evident from the lower panel, our radio measurement is compatible with the tail of the GW170817 jet within the large errors. On the other hand, the X-rays show a ∼ 2 σ excess and could indicate the onset of a new component (Hajela et al. 2020a(Hajela et al. , 2021Troja et al. 2020). nova blast wave drives a shock through the interstellar medium, resulting in synchrotron emission. Electrons are accelerated to a power-law distribution of Lorentz gamma factors γ e > γ e,m , with power-law index p. The energy in the kilonova blast wave is distributed as E(> βγ) ∝ (βγ) −α (with γ the Lorentz factor of the shocked fluid) and normalized to the total energy E at some minimum velocity β 0 such that E > (β 0 γ 0 ) = E. It is reasonable to assume that radio (GHz) observations are in between the minimum frequency, ν m (corresponding to γ m , see Nakar & Piran 2011), and the cooling frequency, ν c . In this case, the kilonova peak flux density reads (Nakar & Piran 2011): where Q x = Q/10 x is followed for all quantities (Q, all expressed in cgs units); B and e are the fractions of the total energy in the magnetic field and electrons respectively; n, the number density of the medium; d is the distance to the source; the normalization constant is calculated for p = 2.1. The time at which the kilonova afterglow emission peaks can be calculated as (Kathirga- (2) The blast wave can be approximated to be mildly relativistic before this peak, and therefore the rising part of the kilonova ejecta light curve can be easily modeled as (see Kathirgamaraju et al. 2019, and references therein): where: In Figure 2, we plot the rising portion of the 3 GHz kilonova light curves obtained following the above prescriptions, and setting β 0 = 0.3. This choice is motivated by the fact that observations in UV/optical/IR of the early kilonova, which only probe the slowest-moving material, point to speeds of ∼ 0.1c − 0.3c (see Section 1). Since we expect the radio to probe the fastest tail of the kilonova ejecta, we consider β 0 ∼ 0.3 a reasonable choice. We note however that smaller values of β 0 , though unlikely, would shift the radio light curve peak to later times, thus allowing for less steep values of α. We set E iso = 10 51 erg, d = 40 Mpc, while varying the power-law index α of the energy-speed distribution of the kilonova ejecta. The solid lines correspond to the choice p = 2.1, e = 7.8 × 10 −3 , B = 9.9 × 10 −4 n = 9.8 × 10 −3 cm −3 as derived from the modeling of the earlier-time panchromatic afterglow of the GW170817 structured jet (Makhathini et al. 2020). For comparison, the dashed lines show the case e = 10 −1 , B = 10 −3 n = 10 −2 cm −3 , p = 2.2, which corresponds to the generic case discussed in Kathirgamaraju et al. (2019). As evident from this Figure, to explain the absence of a kilonova detection in the radio one needs α 5 for the case where the density and micro-physical parameters are set equal to the ones measured for the structured jet afterglow. This constraint on α agrees with the predictions from numerical simulations described in Hotokezaka et al. (2018) and X-ray observations discussed in Hajela et al. (2020b). For the more generic parameters as in Kathirgamaraju et al. (2019), α 20.
For an equal mass ratio binary, a steeper energyvelocity distribution at a given βγ correlates with a stiffer NS EoS for a given cold, non-rotating maximum mass (compare e.g. SF Ho and LS220 in Figures 1 and 9 of Radice et al. 2018). For EoS with the same stiffness (i.e. with the same radii at NS masses of 1.4 M ), larger values of the cold, non-rotating maximum NS mass also correlate to steeper energy-velocity ejecta distribution (as long as a NS is formed even if for a short timescale of order ∼ 1 ms; compare e.g. BHBΛφ and DD2 in Fig. 9 of Radice et al. 2018). On the other hand, if a fast tail exists in the ejecta, one can robustly exclude stiff EoS and relatively high mass ratio scenarios due to the weak or absent core bounce in these scenarios (Nedora et al. 2021). Taken together, if future radio observations reveal a kilonova afterglow, these trends would favor moderate stiffness and mass ratio models. Given these considerations, the constraints we are setting here on α shed some light on the possible EoS, but cannot uniquely pinpoint it. An independent measurement of β via direct size imaging (once the ejecta becomes bright enough in the radio) together with constraints on α derived from light curve modeling, may help reducing degeneracies.

SUMMARY AND CONCLUSION
We have presented extensive late-time radio observations of the GW170817 field, carried out with the most extended configurations of the VLA. Combining the collected data we have built the deepest high-resolution image of the GW170817 field available so far. Our radio flux density measurements show that there is no evidence for emission in excess to the one expected from the afterglow of the GW170817 structured jet at 3 GHz and t ≈ 1200 days since merger. These results constrain the energy-speed distribution of the kilonova ejecta to be rather steep, with a power-law index of α 5 (for β 0 0.3). We finally commented on how the recent detection of a potential excess in the X-rays may hint to a flattening of the power-law index (albeit ≈ 2.1 σ in terms of significance) of the electron energy distribution of the kilonova ejecta compared to the value of this parameter as constrained by the earlier panchromatic afterglow observations. Further late-time monitoring of the GW170817 field with the VLA is likely to reveal whether a kilonova afterglow is emerging.
We are grateful to Dale A. Frail for contributing to shaping this work with many insightful discussions. A.B. and A.C. acknowledge support from the National Science Foundation via Award #1907975. K.P.M. and G.H. acknowledge support from the National Science Foundation Grant AST-1911199. D.L.K. was supported by NSF grant AST-1816492. D.L. acknowledges support from the National Science Foundation via Award #1907955. The National Radio Astronomy Observatory is a facil-ity of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.