The afterglow and kilonova of the short GRB 160821B

ABSTRACT

GRB 160821B is a short duration gamma-ray burst (GRB) detected and localized by the Neil Gehrels Swift Observatory in the outskirts of a spiral galaxy at z = 0.1613, at a projected physical offset of 16 kpc from the galaxy’s center. We present X-ray, optical/nIR, and radio observations of its counterpart and model them with two distinct components of emission: a standard afterglow, arising from the interaction of the relativistic jet with the surrounding medium, and a kilonova, powered by the radioactive decay of the sub-relativistic ejecta. Broadband modelling of the afterglow data reveals a weak reverse shock propagating backward into the jet, and a likely jet-break at 3.5 d. This is consistent with a structured jet seen slightly off-axis (θ_view ∼ θ_core) while expanding into a low-density medium (n ≈ 10⁻³ cm⁻³). Analysis of the kilonova properties suggests a rapid evolution towards red colours, similar to AT2017gfo, and a low-nIR luminosity, possibly due to the presence of a long-lived neutron star. The global properties of the environment, the inferred low mass (M_ej ≲ 0.006 M_⊙) and velocities (v_ej ≳ 0.05c) of lanthanide-rich ejecta are consistent with a binary neutron star merger progenitor.

gravitational waves, nuclear reactions, nucleosynthesis, abundances, gamma-ray burst: general, stars: neutron

1 INTRODUCTION

Short duration GRBs were long suspected to be the product of compact binary mergers (Blinnikov et al. 1984; Goodman 1986; Paczynski 1986; Eichler et al. 1989; Narayan et al. 1992), involving either two neutron stars (NSs) or an NS and a solar-mass black hole (BH). The merger remnant, either a massive NS (Bucciantini et al. 2012; Giacomazzo & Perna 2013) or an accreting BH (Ruffert & Janka 1999; Baiotti, Giacomazzo & Rezzolla 2008; Kiuchi et al. 2009), launches a highly relativistic jet, which produces the GRB emission (Rezzolla et al. 2011; Paschalidis, Ruiz & Shapiro 2015; Ruiz et al. 2016). Later, the interaction of this jet with the surrounding medium produces the observed broadband (from radio to X-rays) afterglow spectrum via synchrotron radiation (e.g. Granot & Sari 2002). Within this framework, the study of the GRB afterglow probes the jet structure and geometry as well as the properties of the surrounding environment.

Another long-standing prediction of the NS merger model is the presence of a luminous, short-lived transient arising from the radioactive decay of freshly synthesized r-process elements (Li & Paczyński 1998; Metzger et al. 2010; Barnes & Kasen 2013; Tanaka & Hotokezaka 2013). Whereas the afterglow does not provide any direct link to the GRB progenitor, such radioactively powered transient, initially dubbed ‘mini-supernova’ and now more commonly referred to as ‘kilonova’, is the distinctive signature of compact binary mergers and clear signpost of heavy elements production in these systems (Lattimer & Schramm 1974, 1976; Freiburghaus et al. 1999; Goriely, Bauswein & Janka 2011; Roberts et al. 2011; Korobkin et al. 2012; Grossman et al. 2014; Rosswog et al. 2014). Evidence of kilonovae associated to short GRBs was only recently found (Tanvir et al. 2013; Yang et al. 2015; Jin et al. 2016). The time-scales and red colour of these first candidate kilonovae suggested that the merger ejecta were highly opaque, as predicted for the heaviest r-process elements (Barnes & Kasen 2013; Tanaka & Hotokezaka 2013). The high opacity causes any UV/optical emission to be significantly suppressed, thus naturally explaining the lack of kilonova detections in over a decade of Swift observations (Bloom et al. 2006; Kocevski et al. 2010; Kann et al. 2011; Pandey et al. 2019).

This slowly progressing field was revolutionized by the discovery of GW170817 and its electromagnetic counterparts (Abbott et al. 2017a, b). X-ray observations revealed the onset of an off-axis afterglow (Troja et al. 2017) powered by a relativistic structured jet (Mooley et al. 2018; Ghirlanda et al. 2019; Troja et al. 2019). Thanks to the proximity of the event, the associated kilonova AT2017gfo was characterized in great detail (e.g. Covino et al. 2017; Pian et al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Troja et al. 2017). In addition to the expected red emission, peaking in the nIR a few days after the merger, a distinctive feature of AT2017gfo was its luminous UV/optical light, peaking at early times and rapidly fading away (Evans et al. 2017). Although the origin of this blue component remains an open question, it immediately revealed a complex chemical composition and velocity structure of the merger ejecta, indicating that the phenomenology of a kilonova is likely determined by the interplay of multiple outflows emerging from the merger remnant (Kasen et al. 2017; Perego, Radice & Bernuzzi 2017; Radice et al. 2018; Wollaeger et al. 2018).

Fostered by the discovery of AT2017gfo and its luminous blue emission, several attempts were made to find similar cases in archival short GRB observations (Gompertz et al. 2018; Troja et al. 2018b; Rossi et al. 2019). Troja et al. (2018b) found that some nearby events have optical luminosities comparable to AT2017gfo. In particular, they showed that the short GRB150101B was a likely analogue to GW170817, characterized by a late-peaking afterglow and a luminous optical kilonova emission, dominating at early times. Other, although less clear, cases discussed in the literature include GRB060505 (Ofek et al. 2007) and GRB080503 (Perley et al. 2009). Overall, it seems plausible that kilonovae similar to AT2017gfo could have been detected in the optical, although not clearly identified prior to GW170817.

Whereas this observational evidence suggests that r-process nucleosynthesis is common in the aftermath of a short GRB, it does not inform on the production of the heaviest elements, i.e. those with atomic mass number A ≳ 140. NS merger ejecta with an electron fraction Y_e ≳ 0.25 do not have enough neutrons to push the nuclear chain past the second r-process peak at A  ≈ 130, thus producing only a blue and fast-fading kilonova (Li & Paczyński 1998; Metzger et al. 2010). An efficient production of lanthanides and actinides will increase the ejecta opacity, leading to a delayed and redder kilonova emission (Barnes & Kasen 2013; Tanaka & Hotokezaka 2013). Only the detection of such red component, peaking at IR wavelengths, is a clear signature of the production of heavy (A ≳ 140) r-process nuclei. Unfortunately, archival IR observations of short GRBs are sparse and mostly unconstraining, and the only possible nIR detection of a kilonova in a short GRB remains GRB130603B (Tanvir et al. 2013).

The short duration GRB160821B, thanks to its low redshift and a rich multiwavelength data set, is an excellent testbed for kilonova searches and afterglow studies. The presence of a kilonova was discussed by Kasliwal et al. (2017) and Jin et al. (2018), who did not find conclusive evidence for it due to the sparse data set used in both studies. Here, we present a comprehensive broadband analysis of this event, including data from Swift and XMM–Newton in the X-rays, the Gran Telescopio de Canarias (GTC), the William Herschel Telescope (WHT), the Hubble Space Telescope (HST), the Keck I telescope, and the Discovery Channel Telescope (DCT) in the optical/nIR, and the Jansky Very Large Array (VLA) in the radio. Thanks to the good temporal and spectral coverage of our data set, we resolve two emission components, which we identify as the GRB afterglow and its associated kilonova. We can exclude dust as the origin of the observed red colour, and interpret it as evidence for a lanthanide-rich kilonova emission.

Throughout the paper, we adopt a standard ΛCDM cosmology (Planck Collaboration VI 2018). Unless otherwise stated, the quoted errors are at the 68 per cent confidence level, and upper limits are at the 3σ confidence level.

2 OBSERVATIONS AND DATA ANALYSIS

GRB 160821B triggered the Swift Burst Alert Telescope (BAT; Barthelmy et al. 2005) at 22:29:13 UT on 2016 August 21, hereafter referred to as T₀. Based on its duration of T₉₀  = 0.48 ± 0.07 s (Lien et al. 2016), it is classified as a short burst. According to the BAT GRB catalogue¹, its spectrum is relatively soft, and best described (χ² = 59 for 59 dof) by a power law with a high-energy exponential cut-off at 50 keV. The time-averaged fluence is (1.10 ± 0.10) × 10⁻⁷ erg cm⁻² s⁻¹ in the 15–150 keV energy band. The Fermi GBM Burst catalogue² reports a fluence of (2.0 ± 0.2) × 10⁻⁷ erg cm⁻² s⁻¹ over the broader the 10–1000 keV energy band. Using this value and the redshift z ∼ 0.1613 from the candidate host galaxy (see Section 2.2.6), we derive an isotropic-equivalent energy E_γ,iso =  (1.3 ± 0.3) × 10⁴⁹ erg in the 10–1000 keV energy band, at the lower end of the observed distribution of short GRB energetics.

Its short duration and possible low redshift made it a prime target for kilonova searches, triggering an intense multiwavelength campaign, as we describe in detail below. A log of the observations is reported in Table 1.

Table 1.

Observations of GRB160821B.

Date	T–T0	Telescope	Instrument	Exposure	Band	AB mag	Flux densitya
MJD	(d)			(s)			(μJy)
Optical/nIR
57621.940	0.002	Swift	UVOT	147	wh	>21.9	<8
57621.942	0.004	Swift	UVOT	209	u	>21.4	<12
57622.013	0.076	GTC	OSIRIS	270	r	22.67 ± 0.10	3.4 ± 0.3
57622.017	0.080	GTC	OSIRIS	270	i	22.39 ± 0.07	4.3 ± 0.3
57622.020	0.083	GTC	OSIRIS	180	z	22.28 ± 0.06	4.7 ± 0.3
57623.006	1.07	WHT	ACAM	1440	r	23.83 ± 0.25	1.2 ± 0.3
57623.027	1.09	WHT	ACAM	1680	z	23.6 ± 0.3	1.4 ± 0.4
57623.878	1.941	GTC	CIRCE	540	H	>23.8	<1.1
57623.895	1.958	GTC	CIRCE	1800	J	>24.0	<0.9
57623.921	1.984	GTC	CIRCE	600	Ks	>23.3	<1.7
57623.958	2.021	GTC	OSIRIS	800	g	25.67 ± 0.15	0.22 ± 0.03
57623.965	2.028	GTC	OSIRIS	720	r	25.12 ± 0.12	0.36 ± 0.04
57623.973	2.036	GTC	OSIRIS	450	i	24.56 ± 0.12	0.58 ± 0.06
57623.980	2.043	GTC	OSIRIS	420	z	24.31 ± 0.17	0.72 ± 0.11
57625.564	3.627	HST	WFC3	2484	F606W	26.02 ± 0.06	0.157 ± 0.009
57625.631	3.694	HST	WFC3	2397	F160W	24.53 ± 0.08	0.57 ± 0.04
57625.697	3.760	HST	WFC3	2397	F110W	24.82 ± 0.05	0.44 ± 0.02
57625.929	3.992	GTC	OSIRIS	450	g	>25.6	<0.24
57625.934	3.997	GTC	OSIRIS	420	i	>25.8	<0.19
57626.234	4.297	Keck I	MOSFIRE	160	Ks	24.0 ± 0.4	0.9 ± 0.3
57626.922	4.985	GTC	OSIRIS	800	r	26.49 ± 0.20	0.101 ± 0.019
57629.402	7.465	Keck I	MOSFIRE	145	Ks	>23.9	<0.9
57630.321	8.383	Keck I	MOSFIRE	110	Ks	>23.7	<1.2
57631.924	9.987	GTC	OSIRIS	720	i	>26.0	<0.15
57631.935	9.998	GTC	OSIRIS	1200	g	>25.8	<0.20
57631.950	10.013	GTC	OSIRIS	960	r	>26.2	<0.13
57632.325	10.388	HST	WFC3	1863	F606W	27.9 ± 0.3	0.028 ± 0.008
57632.383	10.446	HST	WFC3	2397	F110W	26.9 ± 0.4	0.07 ± 0.02
57632.449	10.512	HST	WFC3	2397	F160W	26.6 ± 0.3	0.08 ± 0.02
57645.088	23.151	HST	WFC3	1350	F606W	>27.2	<0.05
57645.108	23.171	HST	WFC3	1497	F110W	>26.6	<0.09
57645.154	23.217	HST	WFC3	2097	F160W	>25.7	<0.19
57721.2	99.2	HST	WFC3	5395	F110W	reference	–
57725.3	103.3	HST	WFC3	2484	F606W	reference	–
58333.7	711.7	HST	WFC3	2796	F160W	reference	–
Radio
57622.11	0.17	VLA	–	3600	C	–	26 ± 5
57623.06	1.13	VLA	–	3600	C	–	<15
57632.01	10.07	VLA	–	6480	X	–	<11
57639.04	17.10	VLA	–	6460	X	–	<33
X-ray
57721.995	0.057	Swift	XRT	185	0.3–10 keV	–	0.15|$^{+0.08}_{-0.06}$|
57722.001	0.063	Swift	XRT	566	0.3–10 keV	–	0.05|$^{+0.03}_{-0.02}$|
57722.008	0.070	Swift	XRT	784	0.3–10 keV	–	0.05 ± 0.02
57722.064	0.126	Swift	XRT	363	0.3–10 keV	–	0.08|$^{+0.04}_{-0.03}$|
57722.068	0.130	Swift	XRT	396	0.3–10 keV	–	0.07|$^{+0.04}_{-0.03}$|
57722.074	0.136	Swift	XRT	865	0.3–10 keV	–	0.037|$^{+0.018}_{-0.014}$|
57722.132	0.195	Swift	XRT	862	0.3–10 keV	–	0.035 ± 0.019
57722.222	0.285	Swift	XRT	1273	0.3–10 keV	–	0.025|$^{+0.014}_{-0.010}$|
57722.264	0.327	Swift	XRT	584	0.3–10 keV	–	0.05 ± 0.02
57722.278	0.340	Swift	XRT	1584	0.3–10 keV	–	0.029 ± 0.011
57722.356	0.419	Swift	XRT	3908	0.3–10 keV	–	0.014 ± 0.05
57722.962	1.024	Swift	XRT	9008	0.3–10 keV	–	(3.6 ± 0.2)× 10−3
57724.725	2.327	Swift	XRT	8777	0.3–10 keV	–	<4.6 × 10−3
57625.879	3.942	XMM-Newton	EPIC/PN	10880	0.3–10 keV	–	(2.3 ± 0.3)× 10−3
57631.913	9.976	XMM-Newton	EPIC/PN	22665	0.3–10 keV	–	(2.9 ± 1.5)× 10−4
57737.570	15.172	Swift	XRT	26000	0.3–10 keV	–	<10−3

Note. ^aOptical/nIR fluxes were corrected for Galactic extinction due to the reddening E(B − V) =0.04 along the sightline (Green et al. 2018). X-ray fluxes were corrected for Galactic absorption N_H = 5.7 × 10²⁰ cm⁻² (Willingale et al. 2013), and converted into flux densities at 1 keV using the best fit photon index Γ = 1.88.

2.1 X-rays

2.1.1 Swift/XRT

Observations with the X-Ray Telescope (XRT; Burrows et al. 2005) on-board Swift started 57 s after the trigger, and monitored the afterglow for the following 23 d for a total net exposure of 198 s in Windowed Timing (WT) mode and 39 ks in Photon Counting (PC) mode. Swift data were retrieved from the public on-line repository³ (Evans et al. 2009) by using custom options for the light-curve binning and standard settings for the spectral extraction. Spectra were modelled within xspec v.12.10.1 (Arnaud 1996) by minimizing the Cash statistics (Cash 1979). The Galactic absorption column was modelled with the function tbabs and fixed to the value N_H = 5.7 × 10²⁰ cm⁻² (Willingale et al. 2013). The time-averaged spectrum, from 4 to 74 ks, is well described by an absorbed power-law function with photon index Γ = 1.8 ± 0.2. Based on this model, we derive a counts-to-unabsorbed flux conversion factor of 4.4 × 10⁻¹¹ er g cm⁻² ct⁻¹. The derived X-ray fluxes are reported in Table 1. No evidence for intrinsic absorption was found at the galaxy’s redshift of z ∼ 0.1613, and we place a 3σ upper limit |${N_{{\rm H},z} \lt 5.5\times 10^{21}{\text{ cm}^{-2}}}$| using the model ztbabs.

2.1.2 XMM–Newton

The X-ray afterglow was also observed with XMM–Newton at 2 epochs, on 2016 August 25 (obsID: 0784460301) and on 2016 August 31 (obsID: 0784460401; PI: Tanvir). The PN (Strüder et al. 2001) and two MOS (Turner et al. 2001) CCD cameras operated in Full Frame mode and with the thin optical blocking filter. The PN, MOS1, and MOS2 exposure times were, respectively: 24.6, 27.1, 27.1 ks for the first observation and 34.1, 36.0, 35.9 ks for the second observation.

Data were retrieved from the public archive, and processed using the XMM–Newton Science Analysis Software (SAS; Gabriel et al. 2004) version 16.1.0. After removing intervals of high particle background, the effective exposures were reduced by 55 per cent for the first observation and 30 per cent for the second observation. For our analysis, we selected only events with FLAG = 0 and PATTERN≤4 and PATTERN≤12 for PN and MOS, respectively.

The afterglow is detected with high significance during the first XMM epoch. Source spectra were extracted from a circular region with radius of 20  arcsec, and the background was estimated from two independent source-free boxes around the afterglow position. We extracted 57 source counts from the PN exposure, and modelled the spectrum as described in Section 2.1.1 with an absorbed power-law function of photon index Γ = 1.88 ± 0.24. During the second epoch, the source is still visible with a detection likelihood (DET_ML = 7.9) slightly above the standard detection threshold adopted for the XMM–Newton source catalogues (DET_ML = 6). By summing PN and MOS data, we extract a total of 47 source counts. The corresponding X-ray flux was derived by assuming the same spectral parameters of the first epoch. Our results are reported in Table 1.

2.2 Optical/nIR

2.2.1 Swift/UVOT

Observations with the Ultra-Violet Optical Telescope (UVOT; Roming et al. 2006) on-board Swift started 76 s after the trigger. The GRB position was initially imaged with the white and u filters, but no counterpart was detected. We used the zero-points provided by Breeveld et al. (2011) to convert UVOT count rates into the AB magnitude system (Oke 1974). The corresponding 3σ upper limits are reported in Table 1. Subsequent observations used all the optical and UV filters, and are reported in Breeveld & Siegel (2016).

2.2.2 Gran Telescopio Canarias (GTC)

We observed the GRB afterglow (PI: A. Castro-Tirado) with the 10.4 m Gran Telescopio de Canarias (GTC), located at the observatory of Roque de los Muchachos in La Palma (Canary Islands, Spain), equipped with the Optical System for Imaging and low intermediate Resolution Integrated Spectroscopy (OSIRIS) and the Canarias InfraRed Camera Experiment (CIRCE) instruments. Deep optical images in the g, r, i, and z filters were taken over five different nights, starting as early as 1.8 h after the trigger, in order to characterize any spectral evolution of the GRB counterpart (Xu et al. 2016). A single epoch of nIR imaging in the J, H, and K_s filters was carried out 2 d after the trigger.

Data were reduced and aligned in a standard fashion using a custom pipeline, based mainly on Astropy and photutils python libraries. We combined the images by weighting each individual frame based on its depth, and applying a 3σ clipping algorithm. Photometric zero points and astrometric calibration were computed using the Pan-STARRS catalogue (Chambers et al. 2016). We then performed point spread function (PSF) matching photometry of the afterglow on the combined images. To construct the empirical PSF, we selected bright and isolated stars close to the afterglow, and combined them weighting by their flux. Background and host light contamination were left as free parameters in the fit. The resulting values are reported in Table 1.

On August 23 (≈T₀+2 d), a GTC (+ OSIRIS) spectrum (3 × 1500s) with the R1000B grism and a 1|${^{\prime\prime}_{.}}$| slit covering the 3700–7500 Å range was gathered with the slit being placed in order to cover the GRB location. Data were reduced and calibrated using standard routines. No absorption or emission lines can be detected superimposed on the continuum emission.

2.2.3 William Herschel Telescope (WHT)

Optical imaging with the Auxiliary port CAMera (ACAM) started on 2016 August 22 at 23:55 (≈T₀ + 1.1 d; Levan et al. 2016). Observations were carried out in the r filter for a total integration of 24 min and in z filter for a total exposure of 28 min. Seeing during the observations was around 0.8 arcsec. The data were reduced within IRAF⁴ following standard procedures (e.g. bias subtraction, flat-fielding, etc.). Aperture photometry was performed using the photutils⁵ package, and calibrated to nearby Pan-STARRS sources (Chambers et al. 2016). A 20 per cent systematic uncertainty was added to our measurements to account for contaminating light from the nearby candidate host galaxy. Our measurements are reported in Table 1.

2.2.4 Hubble Space Telescope (HST)

We activated our program to search for kilonovae (GO14087,GO14607; PI: Troja) and, starting on August 22, obtained several epochs of deep imaging with the IR and UVIS channels of the Wide Field Camera 3 (WFC3). A complete log of the observations is reported in Table 1. Data were processed through the STSCI pipeline, and standard tools within the stsci_python package on AstroConda⁶ were used to align, drizzle and combine the exposures into the final images. The resulting plate scales were 0.067 arcsec pixel⁻¹ for the WFC3/IR images and 0.033 arcsec pixel⁻¹ for the WFC3 F606W images. Late-time observations of the host galaxy were used as reference templates. Aperture photometry was performed on the subtracted images, and the tabulated zero points were used to to determine the source brightness.

The GRB counterpart is detected in all filters in the earlier 2 epochs, whereas is no longer visible in subsequent visits. Its position, determined from the optical frames, is RA (J2000) = 18:39:54.56, Dec (J2000) = + 62:23:30.35 with a 1σ uncertainty of 0.04 arcsec, and lies 5.7 arcsec from the centre of a bright face-on spiral (Fig. 1), which is likely the GRB host galaxy. At a redshift z ∼ 0.1613, this corresponds to a projected physical offset of 16.40 ± 0.12 kpc.

Figure 1.

Top panel: The field of GRB160821B observed with HST using the F606W (blue), F110W (green), and F160W (red) filters. The white circle marks the position of the GRB counterpart, located 16 kpc from the centre of its galaxy. Bottom panel: Optical spectrum of the GRB host galaxy taken with the DeVeny spectrograph on the 4.3m DCT. The brightest lines are identified. The crossed circles mark the position of strong telluric features.

2.2.5 Keck

NIR observations with the MOSFIRE instrument on the Keck I telescope were taken at three different epochs, as previously reported by Kasliwal et al. (2017). A possible detection was found during the first epoch at T₀ + 4.3 d, whereas the source was undetected at later times. We retrieved the archival data, and independently analysed them using standard procedures for CCD data reduction. For the first epoch, four frames considered particularly noisy were removed, whereas the remaining ones were aligned using SCAMP (Bertin 2006) and stacked with SWarp (Bertin et al. 2002) into a final image with 160 s of total exposure. Our analysis confirms the presence of a weak signal (3.5σ) at the GRB position. A magnitude of K_s = 22.12 ± 0.38 was derived by performing aperture photometry calibrated to nearby point sources from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006). Our result is in agreement with the value quoted by Kasliwal et al. (2017). We used the offsets from Blanton & Roweis (2007) to convert the 2MASS Vega magnitudes to the AB system, as quoted in Table 1.

2.2.6 Discovery Channel Telescope (DCT)

The DeVeny spectrograph on the 4.3-m DCT was used on 2017 March 19 with the 300 g mm⁻¹ grating in the first order and a 1|${^{\prime\prime}_{.}}$|5 slit to obtain a spectrum for the GRB host galaxy, covering the 3600–8000 Å range at a dispersion of 2.2 Å per pixel. Standard iraf procedures were used for reduction and calibration.

The final spectrum is shown in Fig. 1 (bottom panel). Several nebular emission lines are visible, including those at λ_obs ≈ 4329, 5814, 5759, and 5647 Å, associated with [O ii], H β and [O iii] transitions, respectively, and identify the galaxy to be at z = 0.1613 ± 0.0004, consistent with the preliminary estimate of Levan et al. (2016). The bright line at 7622 Å, associated with the Balmer H α line, falls within the telluric A band. Lines properties were derived by modelling them with Gaussian functions using the splot task in iraf.

2.3 Radio

Radio observations were carried out with the Karl J. Jansky Very Large Array (JVLA) at four different epochs, on 2016 August 22 and 23 (project code: 15A-235; PI: Fong), and on 2016 September 1 and 8 (project code: 16B-386; PI: Gompertz) in B configuration at the centre frequencies of 6 GHz with a bandwidth of 2 GHz (the former 2 epochs) and 10 GHz with a bandwidth of 4 GHz (the latter 2 epochs). The primary calibrator was 3C 286 and the phase calibrator was J1849 + 6705 for all the four epochs. The data were downloaded from the VLA Archive and calibrated using the JVLA CASA pipeline v1.3.11 running in CASA v4.7.2. The data were then split, imaged and cleaned using CASA in interactive mode: Briggs weighting, robustness = 0.5 and 1000 clean iterations were performed. The results are presented in Table 1. The afterglow was detected during the first epoch at a flux of 26 μJy, consistent within the uncertainties with the value reported in Fong, Alexander & Laskar (2016). By performing a 2D Gaussian fit to the source region we obtain a position of RA (J2000) = 18:39:54.56, Dec (J2000) = 62:23:30.32 with an error of 0.07 and 0.08 arcsec, respectively. This is consistent with the optical position (Section 2.2.4). For the other epochs, a 3σ flux density upper limit is provided.

3 DISENTANGLING THE AFTERGLOW AND KILONOVA EMISSION

The kilonova AT2017gfo was characterized by a quasi-thermal spectrum peaking in the optical/UV and rapidly evolving towards redder wavelengths (e.g. Drout et al. 2017; Pian et al. 2017), in overall agreement with kilonova models (e.g. Kasen, Fernández & Metzger 2015; Kasen et al. 2017). The colour evolution of GRB 160821B (Fig. 2) is consistent with a similar behaviour. By using a simple power-law model, F_ν ∝ν^−β to describe the data, we derive β = 0.70 ± 0.20 at T₀ +2 h, rather standard for an afterglow spectrum with ν_m < ν < ν_c (Granot & Sari 2002). Data from the first night of observations show a flatter spectrum, β  ≈ 0.30, although with larger uncertainties. The optical/nIR counterpart displays a much redder colour, β ≈ 1.4–1.8, between 2 and 4 d, and then return to a more typical value of β = 1.1 ± 0.30. Evidence of a possible spectral break at λ  ≈ 10 000 Å is seen in the GTC OSIRIS and CIRCE data at T₀ + 2 d, but no longer visible at later times. Whereas spectral breaks are observed in many GRB afterglows, the change in spectral slope Δβ ≳ 1.0 implied by our observation is hard to reconcile with standard afterglow theory. Substantial colour variation is often observed during the early (≲10³ s) afterglow phases (e.g. Melandri et al. 2017), but is uncommon on longer time-scales (Li et al. 2018). The observed colour evolution in GRB 160821B appears atypical for a GRB afterglow, and instead consistent with the onset of a kilonova peaking in the optical at ≲1 d, then rapidly shifting to longer wavelengths. In our last observation, the intermediate value of the spectral index suggests a significant contribution from the underlying afterglow. This implies that the IR emission started to fade a few days after the burst, as observed in AT2017gfo, and is therefore shorter lived than the candidate kilonova component in GRB130603B (Tanvir et al. 2013).

Figure 2.

Colour evolution of the optical/nIR counterpart, compiled including data from GTC (circles), WHT (diamonds), Keck (star), and HST (squares).

No significant emission from the kilonova AT2017gfo was detected at X-ray or radio wavelengths (Hallinan et al. 2017; Troja et al. 2017). Based on this evidence, we use the X-ray and radio emission to trace the underlying afterglow component, and compare it to the optical/nIR data set in order to detect any excess from the kilonova. In order to model the afterglow, we consider a standard scenario (e.g. Zhang & Mészáros 2004, and references therein) in which the interaction between the jet and the environment generates two shocks: a highly relativistic forward shock (FS) propagating into the outer medium, and a mildly relativistic reverse shock (RS) travelling backward into the ejecta. The shocked electrons are accelerated into a power-law distribution, N(E) ∝ E^−p, and emit their energy via synchrotron radiation. The resulting broadband spectrum is characterized by four quantities: the self-absorption frequency ν_a, the synchrotron frequency ν_m, the cooling frequency ν_c, and the peak flux F_pk, where we use the subscript FS and RS to distinguish the two spectral components.

Template light curves for the kilonova emission were synthesized by interpolating over the rest-frame spectra of AT2017gfo, and stretching times by a factor of (1 + z). We used the spectroscopic data from VLT/X-Shooter (Pian et al. 2017; Smartt et al. 2017) and HST (Troja et al. 2017) for t > 1.5 d, and photometric measurements from Drout et al. (2017), Evans et al. (2017), Tanvir et al. (2017), and Troja et al. (2017) at earlier times.

3.1 Basic constraints to the afterglow

The early X-ray afterglow displays a bright and rapidly fading light curve (Fig. 3). The measured X-ray flux decreases from 2 × 10⁻¹⁰ erg cm⁻² s⁻¹ at 250 s post-burst to 2 × 10⁻¹² erg cm⁻² s⁻¹ at 400 s, implying a temporal slope α  ≈ 9, where F_X ∝ t^−α. This sharp decay cannot be reproduced by standard FS models and is generally attributed to the sudden cessation of central engine activity (Troja et al. 2007; Rowlinson et al. 2010). Other scenarios, invoking a bright RS emission (e.g. Uhm & Beloborodov 2007; van Eerten 2014), are not consistent with the simultaneous optical upper limits, nor with the constraints from the radio afterglow data (see below). We therefore attribute this first phase of the afterglow to the continuous activity of the central engine (see also Lü et al. 2017), and exclude it from subsequent modelling.

Figure 3.

Temporal evolution of the early X-ray afterglow of GRB160821B (red circles) and, for comparison, GRB090515 (grey circles; Rowlinson et al. 2010). The sharp drop in flux (∝t⁻⁹, dashed line) and the deep upper limit from Swift/UVOT rule out an external shock origin for the observed X-ray emission.

After the first orbit, the X-ray afterglow follows a simple power-law decay with slope α₁ = 0.84 ± 0.08. The condition β_X ≳ β_opt ≈ 0.7 is satisfied for an electron spectral index p ≈ 2.3 when the cooling frequency |$\nu _{\rm c}^{\rm FS}$| lies within or above the X-ray band. Given its slow temporal evolution, ν_c ∝ t^−0.5, its passage does not affect the optical/nIR data over the time span of our observations. This additional evidence adds support to the idea that the observed colour evolution is not related to the GRB afterglow. A comparison to the basic closure relations for FS emission (Zhang & Mészáros 2004) shows that the temporal slope α₁ is shallower than the predicted value 3β/2 ≈ 1.05. We interpret this flattening of the light curve as a viewing angle effect (e.g. Ryan et al. 2015) due to the contribution from the jet lateral structure. For this reason, we adopt a Gaussian jet profile, as described in Troja et al. (2017, 2018a), to better model the afterglow evolution.

XMM–Newton data rule out the presence of an early (≲1 d) temporal break, proposed by Jin et al. (2018) based on a smaller data set. They are instead consistent with an uninterrupted power-law decay up to 4 d (3.5 d rest frame) and no spectral evolution. The measured spectral index is β = Γ − 1 = 0.88 ± 0.24, consistent with the earlier XRT spectrum. In the second epoch of XMM observations, the measured flux is three times lower than the predictions based on the simple power-law temporal decay. The lower flux is consistent with a steeper power-law decline of slope α₂ ∼ 2.3, and may indicate that a temporal break, likely a ‘jet-break’ (Rhoads 1999; Sari, Piran & Halpern 1999), indeed occurred after 3.5 d.

Radio observations show a fading afterglow since early times. For standard FS emission, this may indicate that ν_m is already below the radio range and that radio, optical, and X-rays belong to the same spectral segment, as observed for GW170817 (D’Avanzo et al. 2018; Lyman et al. 2018; Troja et al. 2019). However, for GRB160821B the flat radio-to-optical spectral index β_OR ≈ 0.3 at 4 h rules out such regime, and shows that |$\nu _{\rm m}^{\rm FS} \gg$| 6 GHz. In this case, the fading radio light curve could be explained by an early jet-break (e.g. GRB140903A; Troja et al. 2016), a flaring episode (e.g. GRB050724A; Berger et al. 2005), or an RS component from the shock-heated relativistic ejecta (e.g. GRB051221A; Soderberg et al. 2006). Simultaneous X-ray data allow us to exclude the first two options, favouring the RS scenario. For a short GRB with T₉₀ ≈ 0.5 s, the more likely scenario is the thin-shell case leading to a Newtonian RS (Lloyd-Ronning 2018; Becerra et al. 2019). The early radio observations constrain the RS peak flux |$F_{\rm pk}^{\rm RS}\gtrsim$| 26 μJy and frequency |$\nu _{\rm pk}^{\rm RS}\lesssim$|6 GHz at 3.5 h. Given the predicted evolution of the RS characteristic frequency |$\nu _{\rm m}^{\rm RS} \propto t^{-1.5}$| (Kobayashi 2000), the RS peak was well below the optical range at 100 s, and therefore did not contribute to the early X-ray emission. Optical upper limits, down to wh ≳ 21 mag at 100 s (Fig. 3), also exclude the presence of a bright RS component at early times.

Radio observations at later times help constrain the FS peak flux, |$F_{\rm pk}^{\rm FS}\lesssim$|50 μJy, and synchrotron frequency, |$\nu _{\rm pk}^{\rm RS}\gtrsim$| 100 GHz at 1 d. For a homogeneous circumburst medium, the radio flux should rise as t^1/3, and eventually become visible at late times. The lack of detection at 10 and 17 d therefore supports the presence of an earlier jet-break.

When viewed within the effective core opening angle θ_core, that is the viewing angle θ_view ≲ θ_core, the afterglow of an FS with lateral structure is similar to that of a top hat jet viewed on-axis and the jet break time-scales roughly as |$\theta _{\rm core}^{8/3}$|⁠, with deviations as θ_view approaches θ_core (Rossi, Lazzati & Rees 2002). Numerical simulations of initially top-hat jets show that the observed jet-break time-scales with these angles as |$t_j = 3.5 {\rm ~d}\times (1+z)E_{50}^{1/3} n_{-3}^{-1/3} ((\theta _{\rm core} + \theta _{\rm view})/0.2)^{8/3}$|⁠, where E₅₀ is the on-axis isotropic equivalent energy in units of 10⁵⁰ erg and n₋₃ is the circumburst medium number density in units of 10⁻³ cm⁻³ (van Eerten, Zhang & MacFadyen 2010). Then, using fiducial values E₅₀ = 1 and n₋₃ = 1, the observed break in the light curve is consistent with a jet break if θ_core + θ_view ≈ 0.2 rad.

3.2 Broadband afterglow modelling

Based on the preliminary constraints discussed in Section 3.1, we described the afterglow emission using a standard FS model and a structured jet with a Gaussian profile of width θ_core. We included in the fit the X-ray data, the late radio upper limits, and the first epoch of GTC optical/nIR observations (≈T₀ + 0.08 d). The rest of the optical/nIR data (from T₀ + 1 d to T₀ + 10.5 d; Table 1) were treated as upper limits to the afterglow flux. The early radio detection, dominated by RS emission, was also not included.

Prior to the fit all the data were corrected for Galactic extinction. Given the evidence for negligible intrinsic absorption in the X-ray spectra, and also considered the GRB location in the outskirts of its host galaxy, we assumed |$A^{\rm host}_{\rm V}\approx$| 0. A small amount of extinction may affect our estimates of the spectral slope, but not the evident colour evolution (Fig. 2).

We followed the same procedure used in Troja et al. (2018a, 2019), and performed a Bayesian Markov Chain Monte Carlo (MCMC) model fit to data from synthetic detections generated from our jet model. Despite the large number of free parameters, the FS properties are reasonably well constrained:

$$\begin{eqnarray*} {\rm log~} E_{\rm 0,iso} ~({\rm erg}) & = & 50.4 ~(-0.3, +0.7) \\ {\rm log~} n ~({\rm cm}^{-3}) & = & -2.8 ~(-1.4, +1.1) \\ {\rm log~} \epsilon _e & = & -0.17 ~(-0.22, +0.12) \\ {\rm log~} \epsilon _B & = & -2.1 ~(-1.0, +1.0) \\ p & = & 2.31 ~(-0.05, +0.08) \\ \theta _{\rm core} ~({\rm rad}) & = & 0.14 ~(-0.06, +0.36) \\ \theta _{\rm view} ~({\rm rad}) & = & 0.17 ~(-0.08, +0.25), \end{eqnarray*}$$

where E_{0, iso} is the isotropic equivalent blastwave energy, n the ambient density, ϵ_e and ϵ_B the shock microphysical parameters, and θ_view the observer’s angle with respect to the jet-axis. The resulting prompt radiative efficiency is |$\eta _{\gamma } = E_{\gamma ,\rm iso} / \left(E_{\gamma ,\rm iso} + E_{0,\rm iso} \right)$| = 0.05|$^{+0.05}_{-0.04}$|⁠. Viewing angle effects are taken into account into our modelling and the derived |$E_{0,\rm iso}$|⁠, but not in the observed |$E_{\gamma ,\rm iso}$|⁠. Accounting for them would increase its value by a factor 3, that is η_γ ≈ 0.15.

Our best-fitting model and its uncertainty are shown in Fig. 4. In this model, early emission at T₀ + 0.1 d is dominated by a small patch around the line of sight with Lorentz factor ∼ 8. The patch widens as the jet decelerates, through the jet break, until at T₀ + 10d whole jet is in sight and has a Lorentz factor of roughly 1.5. Whereas at early times (<T₀ + 0.1 d) it provides a good representation of the broadband data set, at later times it does not naturally account for the drastic colour change of the optical/nIR data. It therefore underpredicts the observed emission at nIR and, to a less extent, optical wavelengths. In Fig. 4, we also overplot the template light curves for AT2017gfo in the F606W and F110W filters, rescaled to match the observed fluxes. The close resemblance between the temporal evolution of the optical /nIR excesses and the kilonova light curves suggests that they share a common origin.

Figure 4.

Multicolour light curves of GRB160821B compared to the standard FS (dashed line) + RS (dot–dashed line) scenario. The FS model is the best fit to the broadband data set, whereas the RS is described as a fast fading power law of slope ≈ 2 (Kobayashi 2000). The jet-break time t_jet is shown by the thick vertical line. The shaded areas show the 68 per cent uncertainty in the model. Excess emission at optical and nIR wavelengths is compared with the template kilonova light curves of AT2017gfo (solid line). The redshifted optical light curve of SN1998bw (dotted line; Galama et al. 1998) is also shown for comparison. Errors are 1σ, downward triangles are 3σ upper limits. For plotting purposes, r, z, and K_s data were rescaled using the observed colours (Fig. 2) in order to match the F606W and F110W filters, respectively.

Fig. 5 reports the spectral energy distribution at three different epochs. Whereas at early times (T₀ + 2 h) the broadband spectrum can be described by an afterglow model, at later times emission at optical/nIR wavelengths is brighter and redder than the model predictions (dashed line). This excess can be modelled with a blackbody spectrum of decreasing temperature, from 4500 K at 2 d to 3300 K at 4 d, and radius R ≈ 10¹⁵ cm.

Figure 5.

Broadband spectral energy distribution of GRB160821B. Downward triangles are 3σ upper limits. A standard non-thermal afterglow, including a forward shock (dashed line) and a reverse shock (dotted line) component, can reproduce the early data (T₀ + 2 h). At later times, an additional emission component is visible in the optical/nIR range. It can be described by a blackbody spectrum (solid line) with decreasing temperature.

4 KILONOVA PROPERTIES

4.1 Comparison to AT2017gfo and GRB130603B

We interpret the red excess detected in GRB160821B as kilonova emission from fast-moving lanthanide-rich ejecta. This is only the second case of a short GRB with a kilonova detection in the nIR, where the emission is determined by the heaviest elements (A ≳ 140). In Fig. 6, we compare its properties to the kilonova AT2017gfo and to the other candidate kilonova in GRB130603B. The IR emission of GRB130603B is brighter and longer lived, and its properties are not a good match, as reported in our preliminary analysis of this event (Troja et al. 2016). The observed emission resembles more closely the colour and temporal evolution of AT2017gfo, although less luminous by 1 mag. Our results are consistent with the conclusions of Jin et al. (2018), and provide a better temporal and spectral coverage of the candidate kilonova.

Figure 6.

Rest-frame nIR light curves for AT2017gfo (red solid line) and the two candidate kilonovae in GRB 130603B (star) and GRB 160821B (circles). Error bars are 1σ statistical uncertainties. The hatched area shows the uncertainty in the afterglow subtraction.

The complexity of these systems, in particular the poorly known nuclear physics involved, lead to large uncertainties in the modelling of kilonovae light curves and spectra (e.g. Rosswog et al. 2017). We do not attempt here a systematic comparison to the various models presented in the literature, but estimate the basic explosion parameters based on the scaling relations of Grossman et al. (2014). The observed nIR luminosity, L_nIR ≈ 2 × 10³⁹ erg s⁻¹, and time-scales, t_pk ≲ 3 d, imply a low ejecta mass M_ej ≲0.006 M_⊙ and high velocity v_ej ≳ 0.05c for an opacity κ ≈10 g cm⁻³. This result agrees well with the constraints v ≈ R_bb/t ≈ 0.1–0.2 c, derived from simple blackbody fits (Fig. 5). The ejecta mass is substantially lower than the values inferred for other GRB kilonovae (Tanvir et al. 2013; Yang et al. 2015; Jin et al. 2016), and comfortably within the range of dynamical ejecta from double NS mergers.

Our analysis also finds evidence for an early blue excess, although with larger uncertainties. It is suggestive that the luminosity and time-scale of this blue component are consistent with the early optical emission in AT2017gfo (Fig. 4, top right panel). The blue colour and early onset require a larger mass (M_ej ≈ 0.01 M_⊙) of lanthanide-poor material, produced, for example, by the merger remnant.

4.2 Effects of a long-lived NS

The merger of two NSs can lead either to a stellar-mass BH or to a hypermassive highly magnetized NS (Giacomazzo & Perna 2013; Piro, Giacomazzo & Perna 2017). The latter is thought to significantly affect both the kilonova colours and afterglow evolution through its continuous energy injection and strong neutrino irradiation (Gao et al. 2015; Kasen et al. 2015; Lippuner et al. 2017; Radice et al. 2018). Indeed, the red colours of GW170817/AT2017gfo and its smooth afterglow light curve, mostly consistent with a standard FS emission, were used to exclude a wide range of possible NS configurations (Margalit & Metzger 2017; Ai et al. 2018). Only a short-lived proto-magnetar or a long-lived NS with a weak poloidal field could be consistent with the electromagnetic properties of GW170817 (Li et al. 2018; Metzger, Thompson & Quataert 2018; Yu, Liu & Dai 2018; Piro et al. 2019).

GRB160821B allows us to study the link between the GRB central engine and the kilonova properties. Its X-ray afterglow shows evidence of a long-lasting central engine. In particular, the early X-ray emission shows a phase of nearly constant flux followed by a sharp decay (Fig. 3) often interpreted as a signature of a magnetar engine (Fan & Xu 2006; Liang, Zhang & Zhang 2007; Troja et al. 2007; Lyons et al. 2010; Rowlinson et al. 2010; Dall’Osso et al. 2011). The sudden cessation of X-ray emission at t ≈ 200 s may be caused by the NS collapse into a BH (Troja et al. 2007).

Within this framework, the observed plateau luminosity, L_X ≈ 5 × 10⁴⁷ erg s⁻¹, and its lifetime, T  ≈ 200 s, can be used to infer the magnetar properties (Zhang & Mészáros 2001). For an isotropic emission, Lü et al. (2017) derived an initial spin period P₀ ≈ 9 ms and a dipolar magnetic field B₀ ≈ 3 × 10¹⁶ G, in overall agreement with the assumption of a magnetar which suffered from dominant energy loss via gravitational wave radiation (Fan, Wu & Wei 2013). Beaming would substantially affect these results, and yield unphysical values of magnetic field and period. This indicates that, if the power source of the early X-ray emission is a newborn spindown NS, then the magnetar-driven outflow is nearly isotropic, in agreement with the recent claim of magnetar-driven fast X-ray transients (Xue et al. 2019).

Prolonged irradiation from the NS remnant will affect the ejecta composition and velocity profile, resulting in a bluer and short-lived kilonova emission. However, the colour evolution of GRB160821B supports the presence of a red kilonova evolving on time-scales similar to AT2017gfo, although less luminous by a factor of 2–5. While the lower IR luminosity could be an effect of the long-lived NS, the red excess shows that, despite it, a good amount of lanthanide-rich material was formed and released into the ambient medium, e.g. by tidally stripped matter ejected before the merger (Korobkin et al. 2012).

5 ENVIRONMENT

Late-time HST observations (Table 1) place deep upper limits on any galaxy underlying the GRB position. Such faint host galaxies are rarely observed in short GRBs (e.g. Piranomonte et al. 2008) and imply, at least, a moderately high redshift z ≳ 1. Similar distance scales are disfavoured by the presence of a kilonova, for which luminosities >>10⁴³ erg would lead to implausible properties of the merger ejecta. We consider more natural to link GRB160821B with the nearby bright galaxy (Fig. 1), whose spectrum and disturbed morphology are consistent with a star-forming spiral galaxy, possibly undergoing a merger phase. At the measured offset of 5.7 arcsec, the probability of a chance alignment between the GRB and the galaxy is 2 per cent. However, we find that the overall galaxy properties suggest an environment typical of short GRBs, thus strengthening the case for a physical association.

The galaxy’s absolute B-band magnitude is M_B ≈ −19.9 AB, approximately 0.7|$L^{*}_B$| for a late-type galaxy (Type 3 + 4; Zucca et al. 2009), and its colours are consistent with those of an irregular galaxy (Fukugita, Shimasaku & Ichikawa 1995). It is detected in the WISE bands with W1(3.4 μm) = 20.52 ± 0.10 and W2(4.6 μm) = 21.0 ± 0.3 AB mag, from which we estimate a stellar mass log (M/M_⊙) ≈ 8.5 (Wen et al. 2013), at the lower end of the short GRB distribution (Berger 2014).

We used standard emission line diagnostics to infer the average properties of the putative host. The prominent nebular emission lines are indicative of on-going star formation. The [O ii] line luminosity gives a star formation rate SFR([O ii]) ≳  1.5 M_⊙ yr⁻¹ (Kennicutt 1998), consistent with the estimate SFR(FUV) ≈  1.2 M_⊙ yr⁻¹ derived from the UV luminosity M_w2 ≈ −18.5 AB mag (Hao et al. 2011). Line ratios, log ([O iii]/H β) ≈ 0.4 and log ([O iii]/[O ii]) ≈ −0.1, are substantially lower than in long GRB host galaxies, and within the range of short GRB hosts (Levesque & Kewley 2007). The solution based on the R₂₃ (Pagel et al. 1979) and O₃₂ indicators is degenerate, and we can only constrain the upper metallicity branch (Kobulnicky & Kewley 2004), for which we find 12 + log(O/H) ≈ 8.7, consistent with a solar metallicity (Asplund et al. 2009).

The projected physical offset is 16.4 kpc from the galaxy’s centre, at the higher end of the offset distribution for short GRBs (Troja et al. 2008; Fong & Berger 2013; Tunnicliffe et al. 2014). As only a small fraction of stars is expelled during galaxy interactions (e.g. Behroozi, Loeb & Wechsler 2013), the observed offset is likely the result of an intrinsic kick imparted to the progenitor upon birth. For a stellar age of 300 Myr and a velocity dispersion v_disp ≈ 120 km s⁻¹, a minimum kick velocity v_kick ≳ 80 km s⁻¹ would be required to reproduce the observed distance. However, a more detailed treatment, accounting for the galaxy’s potential as well as the possibility of multiple orbits around the galaxy, finds that natal kicks larger than 150 km s⁻¹ are needed to explain the large offsets of some short GRBs (Behroozi, Ramirez-Ruiz & Fryer 2014).

An alternative explanation for short GRBs with large offsets was discussed by Salvaterra et al. (2010), who investigated compact binaries dynamically formed in globular clusters. At the distance of the candidate host, HST upper limits can only exclude the presence of luminous cluster systems, such as those found in M87, but are otherwise unconstraining. The range of circumburst densities from afterglow modelling (Section 3.2) is also compatible with an origin in a globular cluster (Freire et al. 2001).

6 CONCLUSIONS

GRB160821B is a nearby (z ∼ 0.1613) short GRB with bright X-ray, radio, and optical/nIR counterparts. The X-ray emission shows evidence for continued energy injection from a long-lived central engine, active up to 200 s after the burst. At later times, the X-ray afterglow is consistent with standard forward shock emission, whereas the radio signal was likely dominated by a weak reverse shock. Optical/nIR observations show a clear evolution towards red colours, consistent with the onset of a lanthanide-rich kilonova similar to AT2017gfo. The identification of this additional component is challenging due to the contamination from the underlying bright afterglow, and required an uncommonly rich data set to be disentangled. Within the sample of short GRBs, this is only the second kilonova detected in the nIR. Its low luminosity implies a low mass of lanthanide-rich ejecta, possibly as an effect of a long-lived NS remnant.

ACKNOWLEDGEMENTS

This work is dedicated to the memory of John K. Cannizzo, a friend and colleague with whom we shared many stimulating conversations about kilonovae and short GRBs.

We thank the anonymous referee for his/her insightful and constructive comments. We gratefully thank R. Clavero and S. Dichiara for assistance, B. Metzger, and J. Barnes for help in the preliminary modelling. We acknowledge the use of public data from the Swift data archive. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. These results also made use of Lowell Observatory’s Discovery Channel Telescope. Lowell operates the DCT in partnership with Boston University, Northern Arizona University, the University of Maryland, and the University of Toledo. Partial support of the DCT was provided by Discovery Communications. LMI was built by Lowell Observatory using funds from the National Science Foundation (AST-1005313). The work is partly based on the observations made with the Gran Telescopio Canarias (GTC), installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias, in the island of La Palma.

ET acknowledges financial support provided by the National Aeronautics and Space Administration through HST-GO14357, HST-GO14087, HST-GO14607, and HST-GO14850 grants from the Space Telescope Science Institute, operated by the Association of Universities for Research in Astronomy, Incorporated. AJCT, YDH e IM acknowledge financial support from the State Agency for Research of the Spanish MCIU through the ‘Center of Excellence Severo Ochoa’ award for the Instituto de Astrofísica de Andalucía (SEV-2017-0709). The development of CIRCE at GTC was supported by the University of Florida and the National Science Foundation (grant AST-0352664), in collaboration with IUCAA. RSR acknowledges support by Italian Space Agency (ASI) through the Contract no. 2015-046-R.0 and by AHEAD the European Union Horizon 2020 Programme under the AHEAD project (grant agreement no. 654215). YDH also acknowledges the support by the program of China Scholarships Council (CSC) under the grant no. 201406660015. JBG acknowledges the support of the Viera y Clavijo program funded by ACIISI and ULL. GN and AT acknowledge funding in the framework of the project ULTraS (ASI–INAF contract no. 2017-14-H.0). GR acknowledges the support of the University of Maryland through the Joint Space Science Institute Prize Postdoctoral Fellowship.

Footnotes

1

https://swift.gsfc.nasa.gov/results/batgrbcat/

2

https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermigbrst.html

3

http://www.swift.ac.uk/xrt_products/

4

iraf is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation.

5

https://photutils.readthedocs.io/en/stable/

6

AstroConda is a free Conda channel maintained by the Space Telescope Science Institute (STScI). It provides tools and utilities required to process and analyse data from the Hubble Space Telescope (HST), James Webb Space Telescope (JWST), and others: http://AstroConda.readthedocs.io/.