Short GRB 160821B: a reverse shock, a refreshed shock, and a well-sampled kilonova

TL;DR

This paper presents detailed multi-wavelength observations of short GRB 160821B, revealing complex afterglow features, a kilonova with specific ejecta masses, and evidence for a refreshed shock and reverse shock, advancing understanding of neutron star mergers.

Contribution

It provides the most well-sampled kilonova light-curve without a gravitational wave trigger and identifies a refreshed shock and reverse shock in the afterglow of a short GRB.

Findings

Detection of a kilonova with specific ejecta masses.

Evidence for a refreshed shock at >1 day post-burst.

Observation of a reverse shock in the radio afterglow.

Abstract

We report our identification of the optical afterglow and host galaxy of the short-duration gamma-ray burst GRB 160821B. The spectroscopic redshift of the host is $z=0.162$, making it one of the lowest redshift sGRBs identified by Swift. Our intensive follow-up campaign using a range of ground-based facilities as well as HST, XMM and Swift, shows evidence for a late-time excess of optical and near-infrared emission in addition to a complex afterglow. The afterglow light-curve at X-ray frequencies reveals a narrow jet, $\theta_j\sim1.9^{+0.10}_{-0.03}$ deg, that is refreshed at $>1$ day post-burst by a slower outflow with significantly more energy than the initial outflow that produced the main GRB. Observations of the 5 GHz radio afterglow shows a reverse shock into a mildly magnetised shell. The optical and near-infrared excess is fainter than AT2017gfo associated with GW170817, and is well explained by a kilonova with dynamic ejecta mass $M_{\rm dyn}=(1.0\pm0.6)\times10^{-3}$ M$_{\odot}$ and a secular (postmerger) ejecta mass with $M_{\rm pm}=(1.0\pm0.6)\times10^{-2}$ M$_\odot$, consistent with a binary neutron star merger resulting in a short-lived massive neutron star. This optical and near-infrared dataset provides the best-sampled kilonova light-curve without a gravitational wave trigger to date.

Tables

Table 1. Table 1: Optical and near-IR photometry of the sGRB 160821B afterglow

$\Delta t$ (d)	$t_{\exp }$ (s)	Telescope/Camera	Filter	$A{B}_0$	Source of photometry
0.95	14×300	TNG/DOLoRes	$g$	24.02±0.16	This work
2.02	7×120	GTC/OSIRIS	$g$	25.56±0.16	This work
3.98	10×120	GTC/OSIRIS	$g$	25.98±0.15	This work
6.98	21×120	GTC/OSIRIS	$g$	26.90±0.18	This work
0.05	6×300	NOT/AlFOSC	$r$	22.58±0.09	This work
0.07	6×300	NOT/AlFOSC	$r$	22.52±0.06	This work
0.08	3×90	GTC/OSIRIS	$r$	22.53±0.03	This work
1.06	6×240	WHT/ACAM	$r$	23.82±0.07	This work
1.95	9×300	NOT/AlFOSC	$r$	24.81±0.07	This work
2.03	5×120	GTC/OSIRIS	$r$	24.80±0.06	This work
3.64	4×621	HST/WFC3/UVIS	F606W	25.90±0.06	This work
4.99	27×120	GTC/OSIRIS	$r$	26.12±0.25	This work
10.40	4×621	HST/WFC3/UVIS	F606W	27.55±0.11	This work
23.20	1350	HST/WFC3/UVIS	F606W	¿27.34	This work
0.08	3×90	GTC/OSIRIS	$i$	22.37±0.03	This work
2.04	5×90	GTC/OSIRIS	$i$	24.44±0.10	This work
4.00	3×90	GTC/OSIRIS	$i$	25.70±0.38	This work
9.97	18×90	GTC/OSIRIS	$i$	¿25.59	This work
0.08	3×60	GTC/OSIRIS	$z$	22.39±0.02	This work
1.08	6×240	WHT/ACAM	$z$	23.60±0.15	This work
1.99	9×300	NOT/AlFOSC	$z$	23.90±0.23	This work
2.04	7×60	GTC/OSIRIS	$z$	24.34±0.24	This work
3.76	2397	HST/WFC3/IR	F110W	24.69±0.02	This work
10.53	2397	HST/WFC3/IR	F110W	26.69±0.15	This work
23.18	1498	HST/WFC3/IR	F110W	¿27.34	This work
0.96	33×20	GTC/CIRCE	$H$	23.83±0.35	This work
3.71	2397	HST/WFC3/IR	F160W	24.43±0.03	This work
10.46	2397	HST/WFC3/IR	F160W	26.55±0.23	This work
23.23	2098	HST/WFC3/IR	F160W	¿27.21	This work
4.3	45×30.8	Keck/MOSFIRE	$K$	24.04^+0.44_-0.31	Kasliwal et al. (2017b)

Table 2. Table 2: Swift (top) and XMM-Newton (bottom) X-ray observations in the 0.3–10 keV band, of the sGRB 160821B afterglow after the first hour.

$t$	0.3-10 keV flux
(d)	(${10}^{-14}$ erg cm-2 s-1)
${0.06}_{-0.01}^{+0.01}$	${59.6}_{-10.8}^{+10.8}$
${0.14}_{-0.02}^{+0.06}$	${45.8}_{-7.50}^{+7.50}$
${0.30}_{-0.03}^{+0.03}$	${32.1}_{-7.45}^{+9.46}$
${0.34}_{-0.01}^{+0.01}$	${28.0}_{-6.00}^{+7.42}$
${0.42}_{-0.02}^{+0.13}$	${13.1}_{-2.99}^{+3.74}$
${1.02}_{-0.30}^{+0.39}$	${3.44}_{-1.10}^{+1.49}$
${2.33}_{-0.67}^{+2.11}$	$\le 2.53$
${3.91}_{-0.12}^{+0.12}$	${1.70}_{-0.21}^{+0.21}$
${9.95}_{-0.17}^{+0.17}$	${0.51}_{-0.20}^{+0.20}$

Table 3. Table 3: Radio data used in the analysis

$t$	$\nu$	Flux density	Source
(d)	(GHz)	(mJy)
$0.15$	$5.0$	$0.035$	Fong et al. (2016)
$1.10$	$5.0$		Fong et al. (2016)
$10.06$	$9.8$	$0.016\pm 0.004$	This work
$17.09$	$9.8$		This work

Equations

$${\Gamma }_2\gtrsim 47.4{\left(\frac{1+z}{t_c}\right)}^{1/2}{\left(\frac{E_{\mathrm{k},\mathrm{iso}}}{10^{51}\mathrm{erg}}\right)}^{1/6}{\left(\frac{n}{10^{-4}{\mathrm{cm}}^{-3}}\right)}^{-1/6}{\theta }_j^{1/3},$$

$$\rho(r,t)\propto\left\{\begin{array}[]{lr}r^{-3}~{}t^{-3}&0.025c\leq r/t\leq 0.15c,\\ r^{-6}~{}\zeta(\theta)~{}t^{-3}&0.15c\leq r/t\leq 0.9c.\end{array}\right.$$

$$\zeta (\theta )=0.01+\frac{0.99}{1+e^{-20(\theta -\pi /4)}},$$

Full text

(stars:) gamma-ray burst: individual GRB 160821B, stars: neutron

††journal: ApJ††facilities: GTC(OSIRIS/CIRCE), HST(WFC3), NOT(ALFOSC), Swift(XRT and BAT), TNG(DOLORES), VLA, WHT(ACAM)

1 Introduction

Short-duration gamma-ray bursts (sGRBs) are widely thought to result from the merger of a binary neutron star (BNS) or a neutron star and a stellar mass black hole system. A fraction of the neutron star matter disrupted during the inspiral or collision will undergo rapid accretion onto the remnant object and launch an ultra-relativistic jet (e.g. Nakar, 2007; Gehrels et al., 2009). Energy dissipation within such a jet produces a GRB, and, as this outflow decelerates, an external shock forms producing broad-band afterglow emission. This progenitor model is supported by the fact that well-localised sGRBs (mainly the sample discovered by the Neil Gehrels Swift Observatory, hereafter referred to as Swift) appear to be produced in a wide range of stellar populations, including those with no recent star formation, and on occasions at large distances (10s of kpc in projection) from their putative host galaxies (e.g. Fong et al., 2013b; Tunnicliffe et al., 2014).

A further signature of compact binary mergers involving neutron stars is via the observation of a slower transient, variously called a ‘macronova’ (Kulkarni, 2005), ‘kilonova’ (Metzger et al., 2010), or ‘merger-nova’ (Gao et al., 2015) (in this paper we shall use the term kilonova). A kilonova is powered by the radioactive decay of heavy, unstable, neutron-rich species created from decompressed neutron star material which is ejected during the merger (e.g. Li & Paczyński, 1998).

The first compelling observational evidence for such a kilonova was the case of sGRB 130603B, for which excess near-infrared emission was detected in Hubble Space Telescope (HST) imaging at about one week in the rest frame after the event (Tanvir et al., 2013; Berger et al., 2013). That this excess appeared in the near-IR tallied with predictions that the same heavy r-process elements created in the kilonova should produce dense line-blanketing in the optical, leading to emission appearing in the near-IR in the days to weeks following the merger (Barnes & Kasen, 2013; Kasen et al., 2013; Tanaka & Hotokezaka, 2013). A further interest in these events comes from the fact that this process of radioactive decay naturally leads to stable r-process elements, thus potentially explaining the abundances of more than half the elements in the universe heavier than iron (e.g. Lattimer & Schramm, 1974; Freiburghaus et al., 1999; Rosswog et al., 2018). Mapping the diversity and evolution of kilonova events over cosmic time is therefore an essential ingredient to quantifying their global contribution to nucleosynthesis.

At a redshift $z=0.36$ (de Ugarte Postigo et al., 2014), identifying the kilonova emission in the afterglow to sGRBs 130603B was challenging and would not currently be feasible at higher redshifts, where the bulk of well-localised sGRBs have been found. Indeed, state-of-the-art modelling of neutron-star binary mergers suggests that ejection of sufficient material to create a kilonova as bright as this is unlikely to happen in most mergers, and may require special circumstances such as a high mass-ratio for the components of the binary (e.g. Hotokezaka et al., 2013; Just et al., 2015; Sekiguchi et al., 2016). Nonetheless, following this discovery, and based on archival data, possible kilonova signatures were identified via a late-time $I$-band excess emission in two earlier GRBs; namely sGRB 050709 at $z=0.16$ (Jin et al., 2016), and GRB 060614 at $z=0.125$ (Yang et al., 2015). More recently, it has been proposed that the optical counterparts identified for sGRB 070809 at $z=0.22$ (Jin et al., 2019, although note that the host identification, and therefore redshift, in this case is rather uncertain) and sGRB 150101B at $z=0.13$ (Troja et al., 2018) may have been dominated by kilonova emission. For GRB 060614 the claim is particularly controversial in that its prompt duration, $T_{90}\sim 100$ s, is much longer than the canonical $T_{90}\leq 2$ s for a sGRB. However, the absence of an accompanying bright supernova combined with it exhibiting an initial spike of gamma-rays with duration of only a few seconds has led to speculation that it could have been produced by a compact binary merger (Gal-Yam et al., 2006; Gehrels et al., 2006; Perley et al., 2009; Kann et al., 2011).

The recent multi-messenger observation of the BNS merger GW170817, discovered via gravitational waves and associated with a burst of $\gamma$-rays, GRB 170817A, detected by Fermi and INTEGRAL (Abbott et al., 2017a, b; Goldstein et al., 2017; Savchenko et al., 2017), provided an opportunity to test directly the merger progenitor model. GRB 170817A appeared faint when compared to the cosmological sample of sGRBs and by considering the compactness problem and lack of an early afterglow indicates that the burst of $\gamma$-rays is unlikely to be a typical sGRB seen off-axis (e.g. Lamb & Kobayashi, 2018; Ziaeepour, 2018; Matsumoto et al., 2019), however, Ioka & Nakamura (2019) show that the observed GRB emission likely originates from a ‘mid’-region of a structured outflow. The rapid decline and super-luminal motion of the late-time afterglow to GW170817 offer strong support for the sGRB - BNS association (Ghirlanda et al., 2018; Mooley et al., 2018; van Eerten et al., 2018; Lamb et al., 2019). Additionally, a kilonova was seen to follow GW170817, and monitored intensively at UV, optical and near-infrared wavelengths (e.g. Coulter et al., 2017; Tanvir et al., 2017; Evans et al., 2017; Smartt et al., 2017; Pian et al., 2017; Andreoni et al., 2017; Kasliwal et al., 2017a; Cowperthwaite et al., 2017). By scaling the well-sampled GW170817 kilonova lightcurve to the distance of sGRBs with afterglows, attempts have been made to investigate the diversity of the kilonova population (Gompertz et al., 2018; Ascenzi et al., 2018; Rossi et al., 2019).

Here we report a search with HST, XMM-Newton, and ground-based telescopes including the Gran Telescopio Canarias (GTC), the Nordic Optical Telescope (NOT), the Telescopio Nazionale Galileo (TNG), the William Herschel Telescope (WHT), and the Karl G. Jansky Very Large Array (VLA) for afterglow and kilonova emission accompanying sGRB 160821B, associated with a morphologically disturbed host galaxy at $z=0.162$. We supplement these data with publicly available and/or published in other sources Swift, VLA, and Keck data. Throughout we assume a flat universe with $\Omega_{\rm m}=0.308$ and $H_{\rm 0}=67.8$ km s*-1* Mpc*-1* (Planck Collaboration et al., 2016). Optical and near-IR magnitudes are reported on the AB system. In §2 we report the observations at X-ray, optical, near-IR, and radio frequencies plus the identification of the afterglow and the host. The results, interpretation and afterglow and kilonova modelling are shown in §3. We discuss these results in §4 and give concluding remarks in §5.

2 Observations

2.1 Discovery of sGRB 160821B

The Burst Alert Telescope (BAT) onboard Swift triggered on sGRB 160821B on 2016 Aug 21 at 22:29 UT. The reported duration of the burst was $T_{90}(15-350\,{\rm keV})=0.48\pm 0.07$ s (Palmer et al., 2016). The burst was also detected by Fermi/GBM, from which a somewhat longer duration of $\approx 1$ s was found (Stanbro & Meegan, 2016). Lü et al. (2017) performed a joint fit to the Swift/BAT and Fermi/GBM data, finding the total fluence in the 8-10000 keV band of $(2.52\pm 0.19)\times 10^{-6}$ erg cm*-2*. This corresponds to an isotropic energy, assuming the redshift of $z=0.162$, of $E_{\gamma,\rm iso}=(2.1\pm 0.2)\times 10^{50}$ erg, fairly typical of the population of short GRBs with measured redshifts (Berger, 2014).

2.2 Afterglow identification

After slewing, the X-ray Telescope (XRT) on Swift detected a fading afterglow which provided a refined localisation, and from the X-ray spectrum found no evidence for significant absorption beyond that expected due to foreground gas in our Galaxy (Sbarufatti et al., 2016). As described below, our early optical imaging identified the afterglow of the burst and a prominent nearby galaxy at a separation of about 5.7 arcsec (Xu et al., 2016).

With a magnitude of $r\approx 19.4$ (Section 2.3), the probability of the chance alignment of an unrelated galaxy of this brightness or brighter this close to the line of sight is $P_{\rm chance}\approx 1.5$% (using the formalism of Bloom et al., 2002) and although low, is not entirely negligible. However, the absence of any faint underlying quiescent emission in our final HST epochs (see Section 2.2), which might otherwise suggest a higher redshift host, adds support to our working hypothesis that this is the host galaxy of sGRB 160821B.

The Nordic Optical Telescope, located in the Canary Islands (Spain), began optical observations at 23:02 UT, only 33 minutes post-burst. These revealed an uncatalogued point source within the X-ray error region, presumed to be the optical afterglow (Xu et al., 2016). The best astrometry came from our HST images, and gave a position of RA(J2000) = 18:39:54.550, Dec(J2000) = +62:23:30.35 with an uncertainty of $\approx 0.03$ arcsec in each coordinate, registered on the GAIA DR2 astrometric reference frame (Gaia Collaboration, 2016; Gaia Collaboration et al., 2018). Fong et al. (2016) reported a detection of the radio afterglow at 5 GHz with the VLA, which provided a burst location of RA(J2000) = 18:39:54.56, Dec(J2000) = +62:23:30.3 (reported error 0.3 arcsec), consistent with our HST localisation.

2.3 Host galaxy and redshift

The position of the proposed host galaxy measured from our HST images is RA(J2000) = 18:39:53.968, Dec(J2000) = +62:23:34.35. We obtained spectroscopy of this galaxy with the WHT using the Auxiliary Port Camera (ACAM), in observations beginning on 2016 Aug 22 at 22:57 UT (Levan et al., 2016). The data were reduced using standard IRAF routines. The resulting 2D and 1D extracted spectra are shown in Figure 1, with emission lines of H$\alpha$, H$\beta$, [S II] and [O III] providing a redshift of $z=0.1616\pm 0.0002$. The slit was aligned to cross both the nucleus of the main galaxy and a fainter blob of emission to the north, labelled ‘B’ and ‘C’ respectively on Figure 2. The latter turned out to be a higher redshift galaxy111For completeness, we note that the impact parameter of the GRB from this background galaxy is $\approx 50$ kpc, and it has a $P_{\rm chance}\approx 40\%$, confirming that it is not a good alternative host candidate. at $z=0.4985\pm 0.0002$, the spectrum of which is also shown in Figure 1.

At a redshift $z=0.162$ the separation between afterglow and host corresponds to 16.4 kpc in projection, which is consistent with the offset distribution found for other sGRBs (Fong et al., 2013a; Tunnicliffe et al., 2014).

Morphologically, the host appears to be a face-on, disturbed spiral galaxy (Figure 2). The extended, warped appearance of the central bulge suggests an ongoing merger, and the nebular emission lines are consistent with active star formation. It is interesting to note, although most likely coincidental, that the hosts of both sGRB 130603B and GRB 170817A were also notably disturbed (Tanvir et al., 2013; Levan et al., 2017).

The foreground extinction corrected magnitude of the host from the HST imaging (with the flux from the $z=0.5$ background galaxy subtracted) is $r_{606,0}=19.4$. This corresponds to an absolute magnitude of $M_{r}=-20.0$, which is $\sim L^{*}/3$ with respect to the Loveday et al. (2015) “blue” (star-forming) galaxy population.

The $r$-band 25 mag arcsec*-2* isophote has a radius of $\approx 3.5$ arcsec, corresponding to a linear scale of $\approx 10$ kpc. However it is possible to trace lower surface brightness emission from the galaxy out to the GRB location, albeit at a faint surface brightness level of $\approx 27$ $r$ mag arcsec*-2*.

2.4 Further optical and near-infrared monitoring

sGRB 160821B is amongst the lowest redshift sGRBs found by Swift to date. This, combined with its comparatively low foreground Galactic extinction of $A_{V}=0.118$ mag (Schlafly & Finkbeiner, 2011), motivated an intensive follow-up monitoring campaign.

Further optical and near-IR imaging was obtained with the NOT, the GTC, and the WHT over the next several nights. These data were reduced using standard procedures, and calibrated photometrically using Pan-STARRS (optical) and 2MASS (near-IR) stars in the field.

Observations with the HST using the Wide Field Camera 3 (WFC3), were obtained in the F606W filter (a wide filter spanning approximately the $V$ and $r$ bands), the F110W filter (a wide $YJ$ band) and the F160W filter ($H$ band) from several days to several weeks post-burst (Troja et al., 2016). We adopted the standard photometric calibration for these bands222http://www.stsci.edu/hst/wfc3/phot_zp_lbn, and aperture corrections were determined using bright point sources on the frames.

In all cases, interactive aperture photometry was performed using the Gaia software333 http://astro.dur.ac.uk/~pdraper/gaia/gaia.html. Care was taken to obtain sky estimates close to the position of the transient, since the background was not entirely free of light from the host galaxy.

These observations revealed the counterpart to be initially steady in brightness during the observations made on the first night, but thereafter it faded monotonically in all bands. In the third HST visit, at $\sim 23$ days, no emission is detected at the burst location, which was confirmed by a final visit at $\approx$100 days. A summary of the results of all our optical and near-IR photometry for the sGRB 160821B afterglow, together with selected magnitudes reported elsewhere, is presented in Table 1.

2.5 X-ray monitoring

Swift/XRT monitoring continued for 2.5 days, showing evidence for a significant break to a steeper rate of fading around 0.4 days. Our XMM-Newton observations comprised two visits at approximately 4 and 10 days post-burst. The first visit produced a very significant detection, and was above a simple extrapolation between the last Swift visits. This is discussed further in §3.

A summary of the X-ray observations is presented in Table 2.

2.6 Radio monitoring

The 5 GHz radio detection in 1 hour of observations at 3.6 hours after the burst had a reported flux density of $\sim 35$ $\mu$Jy; an additional observation with the same telescope at 26.5 hours post-burst returned a 3$\sigma$ upper limit of 18 $\mu$Jy (Fong et al., 2016).

Late-time radio observations of the GRB 160821B field were carried out with the VLA, at a central frequency of about 10 GHz and nominal bandwidth of 4 GHz. The first observation started on 2016 September 01 at 23:24:16 UT; the second observation started on 2016 September 08 at 00:10:33 UT. Data were calibrated using the automated VLA calibration pipeline available in the Common Astronomy Software Applications (CASA). After calibration, data were inspected for flagging, and then imaged using the CLEAN algorithm available in CASA. For each of the observations, we estimated the maximum flux density measured within a circular region centered around the position of GRB 160821B and with a radius of 0.6 arcsec (comparable to the nominal FWHM of the VLA synthesized beam in its B configuration at 10 GHz). If the maximum peak density found within this region is above $3\times$ the image rms, then we report the measured flux density value and assign to it an error obtained by adding in quadrature the image rms and a 5% absolute flux calibration error. On the other hand, if the maximum flux density within the selected circular region does not exceed the $3\times$ rms, we report an upper-limit with value equal to $3\times$ the image rms. Radio data444We note that the measured radio flux at $\sim 17$ days is $\sim 31~{}\mu$Jy and only just below $3\times$ the image rms. The presented upper-limit at this time, $<33~{}\mu$Jy, is likely an underestimate, where the flux at the GRB location plus $2\sigma$ would give a limit of $<53~{}\mu$Jy. are listed in Table 3.

3 Light-curve Behaviour, Interpretation, and Modelling

In this section we describe the behaviour of the light-curve at the various observed frequencies. Additionally, we give our interpretation of this behaviour before estimating the light-curve with physically motivated models. These models provide parameter estimates for the various contributing emission components.

3.1 X-ray frequency light-curve behaviour

A period of extended emission555Due to the lack of a clear or consistent definition for extended emission in GRBs, we follow Kisaka, & Ioka (2015) who define extended emission as X-ray emission with a duration $\sim 10^{2}$ s and indicative of a long-lasting central engine. We additionally note that sGRB 160821B is included in the sample of sGRBs with EE by Kisaka et al. (2017) and Kagawa et al. (2019). (EE) follows the sGRB 160821B prompt emission for a duration of $\sim 200$–300 s. Following the rapid decline of the EE, Swift/XRT and XMM-Newton observations show a shallower decline between $\sim 0.01$ and 10 days; as expected from an afterglow. However, this late-time X-ray flux deviates from the expected power-law decline of a simple afterglow model. The flux level drops below that expected from a power-law decay between $\sim 0.3$ to 4 days. Re-binning the Swift/XRT data into photon bins with a lower minimum count, the behaviour of the X-ray light-curve is more clearly revealed; see Figure 3 where the grey markers show the data using the typical minimum photon count per bin and the black markers show the re-binned flux levels (a triangle indicates an upper-limit). A photon index $\Gamma=1.7$ is assumed, which is consistent with both Swift/XRT ($\Gamma=2.0^{+0.7}_{-0.6}$) and XMM-Newton ($\Gamma=1.4^{+0.5}_{-0.4}$). Horizontal error-bars indicate the duration of the observations at each point. The re-binned data reveal a break in the X-ray light-curve at $\sim 0.35$ days, where the flux drops significantly for all the following data, and the flux level at 2–3 days is comparable to the XMM-Newton observed flux level at $\sim 4$ days.

3.2 Behaviour at optical and near-infrared frequencies

Figure 4 shows the spectral energy distribution of all the optical data from Table 1, where we have averaged together points taken in the same filter at close to the same time. The colour evolution of the transient exhibits a trend from blue in observations taken roughly one day after the burst to a much redder colour in all subsequent detections. This is immediately indicative of an emerging kilonova component which itself is evolving from blue to red on time-scales of days; (see e.g. Perego et al., 2014; Tanaka et al., 2018; Wollaeger et al., 2018). $r$- and $H$-band data are shown in Figure 3 for comparison with a typical power-law decline extrapolated from the power-law used to show the behaviour at X-ray frequencies (see §3.1). The deviation from a power-law with an excess in blue and then red is evident; the behaviour at optical and near-IR is distinct from that at 1 keV.

We note that while treating the F606W magnitudes as $r$-band in principle introduces a systematic error, the measured $g$-F606W colour is flat (consistent with our interpretation below that the optical light is afterglow dominated at these times), indicating that colour corrections would be smaller than the photometric errors. (Furthermore, even for our kilonova models, at the time of those epochs, the predicted difference between F606W and the $r$-band is $\lesssim 0.2$ AB mag.)

3.3 At radio wavelengths

Radio observations show a fading source between $\sim 0.1$ and 1 day, but a detection at $\sim 10$ days indicates continued radio afterglow emission as shown by the red contours in Fig. 2 (Note: the small apparent offset between the radio and optical positions is consistent with the effects of noise in the radio map, given the low $S/N$). The late afterglow is limited by a non-detection at $\sim 17$ days.

3.4 Interpretation

A kilonova component is likely to peak in the optical within one to two days post-merger, leading us to expect the $r$-band flux to be dominated by afterglow at the early ($\sim 0.1$ days) and late ($\sim 10$ days) epochs. Inspection of the spectral energy distribution at $\sim 0.1$ days between the X-ray (1 keV) and the $r$-band optical data reveals $\beta=0.66\pm 0.03$, where $F_{\nu}\propto\nu^{-\beta}$, and is consistent with $\beta=0.68\pm 0.07$ at $\sim 10$ days in agreement with this expectation (see Figure 3). Using the broader spectral index limits at $\sim 10$ days, and assuming a temporal decline as $F_{\nu}\propto t^{-\alpha}$, where $\alpha=3(p-1)/4$, the power-law behaviour for the limits on $p$ from $p=2\beta+1$ is shown. A break in the light-curve at $t_{j}\sim 7$ days is required, where $\alpha=-p$ at $t>t_{j}$; this break will be achromatic. The X-ray light-curve drops significantly below the lower-limit ($p=2.23$) power-law extrapolated to earlier times from $\sim 4$ days.

The X-ray light-curve exhibits an earlier break at $t\sim 0.35$ days, and a late-time excess. Afterglow variability is discussed in Ioka et al. (2005), and such an excess is expected from either a refreshed shock where a slower shell catches up with the initial decelerating outflow (e.g. Panaitescu et al., 1998; Zhang & Mészáros, 2002), or a structured jet with an angle-dependent energy and Lorentz factor distribution (e.g. Lamb & Kobayashi, 2017). By assuming the jet structures used to model the afterglow to GRB 170817A in Lamb et al. (2019), where on-axis the resultant GRB would have been consistent with the short GRB population (e.g. Salafia et al., 2019), then from the observed $\gamma$-ray energy of GRB 160821B we can estimate the system inclination following Ioka & Nakamura (2019). For a Gaussian structure with GRB 170817A-like core energy $[\log_{10}(E_{c})=52.4^{+0.4}_{-0.5}]$, then to reproduce the prompt $\gamma$-ray energy of GRB 160821B, the system should be inclined at $\sim\theta_{c}+(3\pm 2)^{\circ}$ (see also Troja et al., 2019); for a two-component jet $[\log_{10}(E_{c})=52.0^{+0.6}_{-0.9}]$ then the opacity of the low-$\Gamma$ second component must be considered (e.g. Lamb & Kobayashi, 2016) and the expected inclination would be $\sim\theta_{c}+(1.5\pm 1.5)^{\circ}$. For a structured jet, however, a late-time re-brightening in the afterglow is only expected for some structure profiles and at higher inclinations, $\sim(3-5)\times\theta_{c}$666A late excess/re-brightening is not expected from a Gaussian profile structure. where bright $\gamma$-ray emission is not expected (see Lamb & Kobayashi, 2017, 2018; Gill & Granot, 2018; Beniamini & Nakar, 2019; Matsumoto et al., 2019). Considering the bright GRB we assume that GRB 160821B is on-axis or very close to on-axis, where the resultant afterglow would behave similarly to the on-axis case regardless of the jet structure (see Lamb & Kobayashi, 2017). For our working model we favour a refreshed shock scenario with two shells where $\Gamma_{1}>\Gamma_{2}$, here the subscript indicates the shell order. If the jet breaks at $t\sim 0.35$ days, then the apparent break at $t>4$ days is indicative of a turnover in the light-curve following a significant energy injection episode.

The extended emission at X-ray frequencies lasting until $\sim 200$–300 s post sGRB 160821B supports continued engine activity beyond the timescale of the GRB. This X-ray emission is consistent with an outflow episode driven by fallback accretion onto a spinning black hole (Rosswog, 2007; Metzger et al., 2008; Nakamura et al., 2014; Kisaka, & Ioka, 2015; Yu et al., 2015; Kisaka et al., 2017). A peak or break time of $\sim 4$ days for the refreshed shock indicates that the bulk Lorentz factor of the outflow when the second shell catches the first should be low, with $\Gamma(t)\sim 10$ and the second shell will have a Lorentz factor much lower than the value typically expected for a successful GRB, $\Gamma_{2}\ll 100$. Energy dissipated within a low-$\Gamma$ outflow is not expected to be emitted at $\gamma$-ray energies; $\gamma$-rays injected into the outflow will be coupled to the plasma and these photons will adiabatically cool and thermalise due to scattering. The effect of these processes is to suppress any resulting emission which will have a spectral peak at $\sim$X-ray frequencies. Photons that fail to escape from a low-$\Gamma$ jet will be reabsorbed by the outflow and contribute to the jet kinetic energy driving the afterglow (Kobayashi, & Sari, 2001; Kobayashi et al., 2002; Lamb & Kobayashi, 2016). The energy-loss by the photon distribution and re-absorption by the outflow will result in a very low value for the emission efficiency, $\eta$. This low-$\Gamma$ X-ray extended emission producing shell follows the initial, high-$\Gamma$, GRB producing shell which will decelerate as $\Gamma_{1}(t)\propto t^{-3/8}$ as it sweeps-up the ambient medium. However, the second shell encounters very little material and will catch up with the forward shell when $\Gamma_{1}(t)\sim\Gamma_{2}/2$ (Kumar & Piran, 2000). The energy of the second shell refreshes the forward shock resulting in a re-brightening of the afterglow (e.g. Granot et al., 2003).

Although limited, the observations at radio frequencies place tight constraints on any possible afterglow, and the afterglow parameters will be constrained by the detection and upper-limits at 1–10 days. The early radio detection at $\sim 0.1$ days, brighter than the following upper-limits and flux at $\sim 10$ days, is likely the result of a reverse shock (e.g. Mészáros & Rees, 1997; Sari & Piran, 1999; Kobayashi, 2000; Kobayashi, & Sari, 2001; Resmi & Zhang, 2016; Lamb & Kobayashi, 2019). Given the X-ray to optical spectral index $\beta\sim 0.66$, the 5 GHz radio emission at $\sim 0.1$ days is below the characteristic synchrotron frequency $\nu_{m}$; if the $\sim 0.1$ days radio emission at 5 GHz belongs to the forward shock, then as $F_{\rm 5GHz}=F_{\nu,{\rm max}}(\nu_{\rm R}/\nu_{m})^{1/3}$ and considering the flux at X-ray frequencies is $F_{\rm X}=F_{\nu,{\rm max}}(\nu_{\rm X}/\nu_{m})^{-\beta}$, then $\nu_{m}\sim 6.4\times 10^{14}(F_{\rm X}/F_{\rm 5GHz})^{\sim 1}$ Hz giving $\nu_{m}\sim 10^{12}$ Hz. As $\nu_{m}\propto t^{-3/2}$ and $t^{-2}$ for the afterglow before and after the jet break, the 5 GHz radio emission will brighten until a peak when $\nu_{m}=5$ GHz or the jet breaks; in either case, the upper-limit of $18~{}\mu$Jy at $\sim 1$ day post sGRB 160821B  rules out the earlier detection being due to the forward shock. This is the first successfully modelled candidate of a reverse shock in an unambiguous sGRB afterglow and indicates that, in some cases, emission from the reverse shock can be bright despite previous non-detections (Lloyd-Ronning 2018; however see Becerra et al. 2019 where a reverse shock was recently claimed for the candidate short GRB 180418A). Any afterglow model that can explain the behaviour at X-ray frequencies and the early and late optical and near-IR should also be consistent with the detection and limits at radio frequencies.

The afterglow at both radio and X-ray frequencies can constrain the behaviour at optical and near-IR. These observations indicate an excess in blue at early times followed by a reddening; this behaviour is indicative of a kilonova. Previous studies of sGRB 160821B have been restricted to much smaller photometric data-sets and consequently have only drawn weak conclusions about the possibility of a kilonova component and the nature of the afterglow (Jin et al., 2018; Kasliwal et al., 2017b; Gompertz et al., 2018). Here, we use the X-ray, early optical and radio constraints on the afterglow emission to interpret the kilonova contribution at optical and near-IR frequencies. We use the latest kilonova light-curve models based on numerical-relativity simulations to constrain the dynamical and post-merger ejecta masses (e.g. Kawaguchi et al., 2018).

3.5 Afterglow Modelling

We use the analytic solution for a relativistic blast-wave from Pe’er (2012), and the method for generating afterglow light-curves from Lamb et al. (2018) to estimate the broadband afterglow for a given set of parameters. We use the observed data to constrain several of the GRB afterglow parameters. As the optical flux at $\sim 10$ days could still have some kilonova contribution, we use the 1 keV to $r$-band spectral slope at $\sim 0.05$ days to estimate $p$, where $\beta\sim 0.66$ giving $p=2.3$. If we assume a prompt efficiency of $\eta\sim 0.1$–0.15 (Fong et al., 2015), then the isotropic equivalent kinetic energy in the initial outflow is $E_{{\rm k,iso}}\sim(1$– $2)\times 10^{51}$ erg. Throughout, we fix $\varepsilon_{B}=0.01$ for the forward shock, consistent with the range for short GRBs (Fong et al., 2015).

The optical flux is approximately flat between 0.05 and 0.07 days; this flatness combined with a likely reverse shock in the radio at the same time indicates that these points coincide with the deceleration timescale for the outflow. By fixing the ambient density to $n=10^{-4}$ cm*-3*, consistent with the location in the outskirts of the host galaxy (see Figure 2), the Lorentz factor of the GRB outflow can be estimated; $\Gamma_{0}\sim 18~{}[t_{d}/(1+z)]^{-3/8}(E_{{\rm k,iso}}/10^{51}~{}\mathrm{erg})^{1/8}~{}(n/10^{-4}~{}\mathrm{cm}^{-3})^{-1/8}\sim 55$–60, where $t_{d}\sim 0.06$ days is the deceleration time. Similarly, the break at $t_{j}\sim 0.35$ days can be used to estimate the jet half-opening angle, $\theta_{j}\sim 0.05~{}[t_{j}/(1+z)]^{3/8}~{}(E_{{\rm k,iso}}/10^{51}~{}\mathrm{erg})^{-1/8}~{}(n/10^{-4}~{}\mathrm{cm}^{-3})^{1/8}\sim 0.033$ rad, or $\sim 1.9$ deg. As the break time dominates the opening angle estimation, we can put weak limits on this value of $1.9^{+0.10}_{-0.03}$ degrees (these small errors are only the formal fit uncertainty given this choice of jet model and decomposition of the light curve; the systematic errors from uncertainties in the model assumptions are much greater, and poorly quantifiable), this narrow jet is consistent with the opening angle range for short GRBs (Jin et al., 2018).

The forward shock is refreshed at $\sim 1$ day, peaking at $\sim 3$ days and then declining as $\sim t^{-p}$. We assume that the second shell has the same half-opening angle as the first. As the jet has broken, side-ways expansion could widen the initial blast-wave and the second shell will only refresh the blast-wave with an opening angle $\leq\theta_{j}$. By assuming that the radius of the blast-wave is roughly constant after the jet break777The sideways expansion does not halt the radial progress of the jet (Granot, & Piran, 2012; Lamb et al., 2018); by assuming that it does, we can place a lower-limit on the Lorentz-factor of the second shell. then the Lorentz factor of the second shell is

[TABLE]

where $\Gamma_{2}\gtrsim 16$ for an observed collision time $t_{c}\sim 1$ day. The Lorentz factor of the forward shock at the collision is then $\Gamma_{1}(t)\gtrsim 8$.

We find that if the forward shock is refreshed when $\Gamma_{1}(t)=12$ and the resulting blast-wave has $12.5\times E_{{\rm k,iso}}$ of the initial outflow energy then the afterglow can account for the X-ray excess at $\sim 4$ days. The radio afterglow at $\sim 10$ days constrains the micro-physical parameter $\varepsilon_{e}\sim 0.3$, so as not to overproduce the radio flux. We assume throughout that the initial and final blast-wave have identical micro-physical parameters $\varepsilon_{B}$ and $\varepsilon_{e}$, electron index $p$, and $\theta_{j}$.

The early radio point at $\sim 0.1$ days requires a significant reverse shock. For this point to be forward shock dominated the X-ray and optical data constrain the characteristic synchrotron frequency to $\nu_{m}\sim 10^{12}$ Hz, much lower than the model estimate of $\nu_{m}\sim 3.5\times 10^{14}$ Hz. As $\nu_{m}\propto\Gamma^{4}\varepsilon_{B}^{1/2}n^{1/2}\varepsilon_{e}^{2}$, then the parameters that can successfully explain the X-ray and optical afterglow would need significantly lower values. Such lowered parameter values result in an afterglow that is inconsistent with the other observations and unphysical parameters in many cases. Following Harrison & Kobayashi (2013), the characteristic synchrotron frequency $\nu_{m}$ and the maximum flux $F_{\nu,{\rm max}}$ for the reverse shock can be found from the forward shock parameters. The reverse shock flux before and after the peak will scale following Kobayashi (2000); for the thin shell case and our parameters, the flux pre-peak will scale as $F_{\nu}\propto t^{5.7}$ and post peak $F_{\nu}\propto t^{-2.05}$. To accommodate the early radio detection, we need to use a magnetization parameter of $R_{B}\sim 8$. The model light-curve is shown in Figure 5, where we have taken an initial kinetic energy of $E_{{\rm k,iso}}=1.3\times 10^{51}$ erg and $\theta_{j}=0.033$, with all other parameters as discussed.

3.6 The Kilonova Modelling

The kilonova appears as an excess in the optical above the afterglow. From Figure 5, where the optical afterglow is shown as dotted lines, it is clear that all bands are in excess at $\sim 1$ day post-burst. The bluer bands ($g$, $r$, and $i$) follow the afterglow from $\sim 5$ days whilst the redder bands ($J$, $H$, and $K$) remain in excess until $\sim 10$ days post GRB.

Using two-component kilonova models from Kawaguchi et al. (2018), $K$-corrected to $z=0.16$, we find the model parameters via a $\chi^{2}$ minimisation fit to the data for the kilonova plus model afterglow. The kilonova is best described888The models have masses drawn from the parameter-grid $M_{\rm{dyn}}=[0.001,~{}0.002,~{}0.003,~{}0.005,~{}0.01]$ M*⊙, and $M_{\rm{pm}}=[0.01,~{}0.02,~{}0.03,~{}0.05,~{}0.1]$ M⊙*. by a secular ejecta (or post-merger wind driven by viscous and neutrino heating) with a mass $M_{\rm pm}=0.01$ M⊙, and a dynamic ejecta mass $M_{\rm dyn}=0.001$ M⊙. The density profile for each ejecta component is given by

[TABLE]

Here the top condition is for the secular ejecta, and the bottom condition for the dynamic ejecta. We find good fits for an upper-limit for the secular ejecta velocity, and lower-limits for the dynamic ejecta velocity, of $0.1-0.15c$. The function $\zeta(\theta)$ describes the angular distribution of the dynamic ejecta, and is given by

[TABLE]

where $\theta$ is the angle from the central axis.

The element abundances for the ejecta are determined following the results of r-process nucleosynthesis calculations by Wanajo et al. (2014) and assuming that the secular and dynamic ejecta have initially flat electron fraction $Y_{e}$ distributions ranging from 0.3 to 0.4 and from 0.1 to 0.4, respectively. Radiative transfer simulations were performed from 0.1 to 30 days resulting in a lightcurve with a statistical error in each band $\sim 0.1-0.2$ magnitudes.

The kilonova fit to the data depends on the afterglow subtraction, however, the precise details of the afterglow parameters are not crucial. As the optical afterglow is typically in the same spectral regime as the observed X-ray data for sGRBs, and supported by the similar spectral index between optical and X-rays at 0.1 and 10 days, then the optical afterglow will follow that at X-ray frequencies during the kilonova peak. The X-ray data extrapolated to the optical at $\sim 1-4$ days post-burst indicates that the afterglow contributes $\sim 10\%$. The typical photometric uncertainty is $\sim 10\%$, and the kilonova model uncertainty is $\sim 10\%$. Combining these uncertainties, and using the analytic scaling for luminosity with mass $L\propto M^{0.35}$ (e.g. Grossman et al., 2014), we can give limits on the mass estimates from the kilonova model fit of $\sim\pm 60\%$, however, we emphasize that both the masses and the uncertainties are model specific.

4 Discussion

We have shown that the afterglow of sGRB 160821B with extended X-ray emission until $\sim 300$ s post-burst exhibits a reverse shock at early times and a refreshed shock at late times. Early time observations at radio wavelengths require a reverse shock, while the complex light-curve at X-ray frequencies observed by Swift/XRT and XMM-Newton, combined with late time radio observations reveal a break at $\sim 0.35$ days and a re-brightening at $>1$ day. The jet is very narrow, at $\theta_{j}\sim 1.9$ degrees, and the slower second outflow episode that refreshes the forward shock carries significantly more energy than the initial outflow. However, the total combined energy of the jets, $E_{j}\sim 0.9\times 10^{49}$ erg, is consistent with the short GRB population (Fong et al., 2015).

Extended emission can be the result of a magnetar (e.g. Fan & Xu, 2006; Metzger et al., 2008; Bucciantini et al., 2012; Gompertz et al., 2013; Gibson et al., 2017), or energy dissipated within a jet launched due to mass fallback onto the central compact object (Fan et al., 2005; Rosswog, 2007; Kisaka, & Ioka, 2015; Kisaka et al., 2017); see also Barkov & Pozanenko (2011) for a two-component jet model. The refreshed shock at late times requires a second episode of jet activity and fallback accretion onto the central compact object supports both this late-time re-brightening and the extended emission. From the afterglow modelling, the second jet episode has a Lorentz factor of $\Gamma_{2}\sim 24$. Internal energy dissipation within such a low-$\Gamma$ jet is expected to be suppressed due to a large optical depth, see Lamb & Kobayashi (2016), however, any resulting emission will peak at X-ray frequencies and have a longer timescale than the initial dissipation timescale. Considering the energy required to refresh the forward shock, the efficiency of energy dissipation within the fallback launched jet is $\eta\sim 10^{-3}$, consistent with the expectation from a low-$\Gamma$ outflow (Lamb & Kobayashi, 2016). The fallback mass required to launch such an energetic second outflow can be estimated following Kisaka et al. (2017) giving a mass $\sim 2\times 10^{-3}$ M⊙.

As well as the EE and the refreshed shock, the afterglow reveals a reverse shock (the first confirmed reverse shock in an sGRB, see Lloyd-Ronning, 2018, who highlight the lack of observed reverse shocks in sGRBs); such a shock propagates into the colder and denser inner shell. To recreate the reverse shock emission we follow Lamb & Kobayashi (2019) and require a magnetization parameter of $R_{B}\sim 8$. Thus the magnetic field within the shell is much larger than the magnetic field induced by the forward shock. A high magnetic field indicates that the shell is endowed with primordial magnetic fields from the central engine.

In addition to these afterglow features, a kilonova is present at optical and near-IR frequencies. The best fitting model is one represented by a dynamic ejecta mass of $\sim 0.001$ M⊙ and a secular ejecta mass $\sim 0.01$ M⊙. The secular ejecta mass, required for the early blue excess, is consistent with the expectation of the mass-loss from a torus surrounding a massive neutron star (Fujibayashi et al., 2018; Fernández et al., 2019). However, the best-fit model from our parameter sample under-predicts the observed $g$-band emission at $\sim 2$ and $\sim 4$ days post-burst, this is likely due to the finite parameter spacing of the kilonova model samples. A small secular ejecta mass $\sim 0.01$ M⊙ and the low dynamic ejecta mass $\sim 0.001$ M⊙ may indicate that the remnant collapses to a black hole promptly after the merger (Kiuchi et al., 2009; Sekiguchi et al., 2016; Coughlin et al., 2018; Radice et al., 2018). In such a scenario the electron fraction, $Y_{e}$, will be lower. To test this, we compared the kilonova light-curve of the best-fit model with a model using a lower electron fraction distribution for the post-merger wind $Y_{e}=0.1-0.3$ as expected from a prompt collapse scenario. A comparison of the light-curves for these two scenarios was performed, the results indicate that the prompt collapse to a black hole, with a low-$Y_{e}$ and a higher velocity, will overproduce the red excess at late times and underproduce the early blue excess; see Figure 6. Thus, the observed blue emission in the early phase suggests the existence of a low opacity component, when interpreted as kilonova emission, and we can conclude that a very prompt collapse to a black hole is unlikely to explain the observed transient when considering the observed features. Note that the afterglow subtracted data at $\gtrsim 4$ days is typically brighter than the kilonova model we use, especially at K-, J-, $r$- and $g$-bands. This excess at bluer wavelengths is due to the afterglow subtraction, where the emission is afterglow dominated and the model afterglow slightly under-predicting the observed flux. The observed K- and J-band excesses ($\sim 4$ and $\sim 10$ days post-burst) have large associated errors, and the best-fit model is within 2$\sigma$ of each detection without considering the model uncertainty (see Fig. 5).

Of the five widely discussed GRBs with candidate kilonova contributions to their light-curves – GRBs 050709, 060614, 070809, 130603B and 150101B (Jin et al., 2016; Yang et al., 2015; Jin et al., 2019; Gompertz et al., 2018; Troja et al., 2018) – the kilonova in sGRB 160821B is the best sampled. At $\sim 0.011$ M⊙, the kilonova in sGRB 160821B has an ejecta mass towards the lower end of the range proposed for any of these other cases, and is consistent with the $<0.03$ M⊙ found by Kasliwal et al. (2017b). The kilonova following GW170817 had an ejecta mass $\sim 0.03-0.05$ M⊙ (e.g. Smartt et al., 2017; Pian et al., 2017), similar to the mass estimates for sGRB 130603B, $\sim 0.03$ M⊙ (e.g. Jin et al., 2016), whereas, GRB 050709, 060614, 070809 and 150101B have masses $\sim 0.05$, $0.13$, $0.015$ and $<0.004$ M⊙ respectively (Jin et al., 2016; Yang et al., 2015). However, we note that upper-limits implied by kilonova non-detection in some other sGRBs could indicate the existence of fainter kilonovae indicating still lower ejecta masses999The heating rates and therefore the estimated masses depend on the chosen nuclear mass formula (e.g. Barnes et al., 2016; Rosswog et al., 2017). For the very low $Y_{e}$ ejecta the r-process path passes close to the neutron-dripline in the nuclear chart, this is experimentally uncharted territory, and we rely on purely theoretical mass formulae. The amounts of trans-lead nuclei, important since they are efficient in releasing energy and their decay products are efficiently thermalizing with the ambient medium, depend quite sensitively on the chosen mass formula (e.g. Gompertz et al., 2018).

The best-fit kilonova model is consistent with the scenario where, following the merger, a massive neutron star survives for a short period (Fujibayashi et al., 2018). This scenario is similar to the case of GRB 170817A, for which various arguments point to a short-lived massive neutron star (e.g. Margalit & Metzger, 2017; Ai et al., 2018; Pooley et al., 2018, see Piro et al. 2019 for an alternative interpretation); however, the lower ejecta mass in sGRB 160821B could point to a more rapid collapse of the remnant massive neutron star. Extended emission was present in sGRB 160821B and is used to argue for significant mass fallback in this case, however, for GRB 170817A Swift/XRT did not begin observations until $\sim 15$ hours after the initial burst (Evans et al., 2017) and any EE would have long faded. The total energy in the jets in sGRB 160821B is lower than the energy required to drive the afterglow to GRB 170817A and, additionally, the required outflow structure is very different (e.g. Lamb et al., 2019). These differences, combined with the lower mass of the ejecta in sGRB 160821B when compared to GRB 170817A, could offer some clue as to the dynamical differences between mergers and sGRB phenomena. Understanding these differences may help explain the diversity in sGRB properties; especially amongst systems with a similar progenitor i.e. BNS mergers.

5 Conclusions

We have reported ground- and space-based optical and near-infrared monitoring of sGRB 160821B. We see clear evidence for red to blue evolution in the colour of the transient, indicative of a kilonova. The data-set presented here makes the kilonova in sGRB 160821B the best-sampled kilonova without a coincident gravitational wave signal. We find that a kilonova model with a dynamic ejecta mass $M_{\rm dyn}\sim 0.001$ M⊙, a velocity distribution $(0.15-0.9)c$, and a flat electron fraction distribution $Y_{e}=0.1-0.4$; and a secular ejecta with $M_{\rm pm}\sim 0.01$ M⊙, a velocity distribution $(0.025-0.15)c$, and $Y_{e}=0.3-0.4$ can best explain the observed emission, while the mass estimates have $\sim 60\%$ uncertainty. The blue excess, the mass of the dynamic and secular ejecta, and the electron fraction supports the existence of a short-lived massive neutron star that does not immediately collapse to a black hole.

We have also presented Swift and XMM-Newton observations of the event and combining with constraints from VLA radio observations find a complex afterglow with a radio-emitting reverse shock into a magnetised shell and a late-time, broadband, refreshed shock. The jet is very narrow with $\theta_{j}\sim 1.9$ degrees, and the second episode is significantly more energetic than the first. We find the prompt and extended emission, plus the early- and late-time re-brightening afterglow to be consistent with multiple accretion episodes onto the central compact object with the second episode consistent with a fallback mass of $\sim 0.002$ M⊙.

The authors thank the anonymous referee for helpful and constructive comments. GPL additionally thanks Alice Breeveld, Kunihito Ioka, Geoff Ryan, Graham Wynn, and Tomos Meredith for useful discussions and Yizhong Fan for helpful comments. Partly based on observations made with the Gran Telescopio Canarias (GTC), installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma; and with the Nordic Optical Telescope, operated by the Nordic Optical Telescope Scientific Association at the Observatorio del Roque de los Muchachos (program 51-504); and with the Italian Telescopio Nazionale Galileo (TNG) operated by the Fundación Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos (program A32TAC_5). 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. Based on data from the GTC Public Archive at CAB (INTA-CSIC). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. NRT, AJL, KW & BG have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no 725246, TEDE, PI Levan). AJL, JDL acknowledge support from STFC via grant ST/P000495/1. NRT, GPL acknowledge support from STFC via grant ST/N000757/1. EP acknowledges support from grant ASI/INAF I/088/06/0 JH was supported by a VILLUM FONDEN Investigator grant (project number 16599). DBM acknowledges support from the Instrument center for Danish astrophysics (IDA). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. A.C. acknowledges support from the National Science Foundation CAREER award #1455090. AdUP, CCT, ZC, LI and DAK acknowledge support from the Spanish research projects AYA2014-58381-P and AYA2017-89384-P, 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). AdUP and CCT acknowledge support from Ramón y Cajal fellowships (RyC-2012-09975, and RyC-2012-09984). ZC, LI and DAK acknowledge support from Juan de la Cierva Incorporación fellowships (JdCI-2014-21669, IJCI-2016-30940, and IJCI-2015-26153). PAE and KLP acknowledge support from the UK Space Agency. KEH and PJ acknowledge support by a Project Grant (162948–051) from The Icelandic Research Fund. SR has been supported by the Swedish Research Council (VR) under grant number 2016-03657_3, by the Swedish National Space Board under grant number Dnr 107/16 and by the research environment grant “Gravitational Radiation and Electromagnetic Astrophysical Transients (GREAT)” funded by the Swedish Research council (VR) under Dnr 2016- 06012.