The afterglow and kilonova of the short GRB 160821B

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 modeling 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 while expanding into a low-density medium. Analysis of the kilonova properties suggests a rapid evolution toward red colors, 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 Msun) and velocities (v>0.05 c) of lanthanide-rich ejecta are consistent with a binary neutron star merger progenitor.


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 a 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 et al. 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;. 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, 1976Freiburghaus et al. 1999;Roberts et al. 2011;Goriely et al. 2011;Grossman et al. 2014;Rosswog et al. 2014). Evidence of kilonovae associated to short GRBs was only recently found Yang et al. 2015;Jin et al. 2016). The timescales and red color of these first candidate kilonovae suggested that the merger ejecta were highly opaque, as predicted for the heaviest r-process elements (Barnes & Kasen 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;Abbott et al. 2017b). Thanks to the proximity of the event, the associated kilonova AT2017gfo was characterized in great detail (Chornock et al. 2017;Kasliwal et al. 2017b;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;Nicholl et al. 2017;Shappee 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 Perego et al. 2017;Metzger et al. 2018;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 (e.g. GRB150101B, GRB060505, GRB160821B, and GRB050709) have optical luminosities comparable to AT2017gfo. In particular, they showed that GRB150101B was a likely analogue to GW170817, characterized by a latepeaking 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 rprocess 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 180. NS merger ejecta with an electron fraction Ye 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;. Only the detection of such red component, peaking at IR wavelengths, is a clear signature of the production of the third-peak r-process elements. 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. The short duration GRB160821B, thanks to its low redshift and a rich multi-wavelength dataset, is an excellent testbed for kilonova searches and afterglow studies. The presence of a kilonova was discussed by Kasliwal et al. (2017a) and Jin et al. (2018), who did not find conclusive evidence for it due to the sparse dataset 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 dataset, 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 color, and interpret it as evidence for a lanthanide-rich kilonova emission.
Throughout the paper, we adopt a standard ΛCDM cosmology (Planck Collaboration et al. 2018). Unless otherwise stated, the quoted errors are at the 68% confidence level, and upper limits are at the 3 σ confidence level.
Its short duration and possible low redshift made it a prime target for kilonova searches, triggering an intense multi-wavelength campaign, as we describe in detail below. A log of the observations is reported in Table 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 days 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 2 (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 time-averaged spectrum, from 4 ks to 74 ks, is well described by a power-law model with photon index Γ=1.8±0.2 and negligible intrinsic absorption in addition to the Galactic value. Based on this model, we derive a counts-to-unabsorbed flux conversion factor of 4.4×10 −11 erg cm −2 ct −1 .

XMM-Newton
The X-ray afterglow was also observed with XMM-Newton at two 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. 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, we analyzed the two observations selecting only events with FLAG=0 and PATTERN≤4 and PATTERN≤12 for PN and MOS, respectively. The PN, MOS1 and MOS2 exposure times were, respectively: 24.6 ks, 27.1 ks, 27.1 ks for the first observation and 34.1 ks, 36.0 ks, 35.9 ks for the second observation.
The afterglow is detected with high significance during the first XMM epoch. Source spectra were extracted from a circular region with radius of 20 ′′ , and the background was estimated from two independent source-free boxes around the afterglow position. The spectrum was modeled with an absorbed power-law function with gas column density NH = 5.7 × 10 20 cm −2 fixed to the Galactic value (Willingale et al. 2013), and a photon index Γ = 1.88 ± 0.14. 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). As the source is too faint for spectroscopy, we derived the corresponding X-ray flux by assuming the same spectral parameters of the first epoch. Our results are reported in Table 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).

Gran Telescopio Canarias (GTC)
We observed the GRB afterglow (PI: 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 hr after the trigger, in order to characterize any spectral evolution of the GRB counterpart. A single epoch of nIR imaging in the J, H, and Ks 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. The resulting values are reported in Table 1.
On August 23 (≈ T0+2 d), a GTC (+OSIRIS) spectrum (3 × 1500s) with the R1000B grism and a 1. ′′ slit covering the 3,700,Å-7,500,Å 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. dard procedures (e.g., bias subtraction, flat-fielding, etc.). Aperture photometry was performed using the photutils 5 package, and calibrated to nearby Pan-STARRS sources (Chambers, et al. 2016). A 20% systematic uncertainty was added to our measurements to account for contaminating light from a nearby galaxy. Our measurements are reported in Table 1.

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 6 were used to align, drizzle and combine the exposures into the final images. The resulting plate scales were 0.067 ′′ /pixel for the WFC3/IR images and 0.033 ′′ /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 two epochs, whereas is no longer visible in subsequent visits. Its position lies 5.7 ′′ from the center of a bright faceon spiral (Figure 1), which is likely the GRB host galaxy.

Keck
Observations with the MOSFIRE instrument on the Keck I telescope were taken at three different epochs, as previously reported by Kasliwal et al. (2017a). A possible detection was found during the first epoch at T0+4.3 d, whereas the source was undetected at later times. We retrieved the archival data, and independently analyzed them using standard procedures for CCD data reduction. 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 Ks=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. (2017a). We used the offsets from Blanton & Roweis (2007) to convert the 2MASS Vega magnitudes to the AB system, as quoted in Table 1.

Discovery Channel Telescope (DCT)
The DeVeny spectrograph on the 4.3-meter DCT was used on 2017 March 19 with the 300 g mm −1 grating in the first order and a 1. ′′ 5 slit to obtain a spectrum for the GRB host galaxy, covering the 3,600,Å-8,000Å range at a dispersion of ∼2.2Å per pixel. Standard IRAF procedures were used for reduction and calibration.
The final spectrum is shown in Figure 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. The bright line at 7622Å, associated with the Balmer Hα line, falls within the telluric A band. Lines properties were derived by modeling them with Gaussian functions using the splot task in IRAF.

Radio
Radio observations were carried out with the Karl J. Jansky Very Large Array (JVLA) at four different epochs, on Aug 22 and 23, 2016 (project code: 15A-235; PI: Fong), and on Sep 1 and 8, 2016 (project code: 16B-386; PI: Gompertz) in B configuration at the center frequencies of 6 GHz with a bandwidth of 2 GHz (the former two epochs) and 10 GHz with a bandwidth of 4 GHz (the latter two epochs). The primary calibrator was 3C286 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 position and flux consistent with those quoted in Fong et al. (2016). For the other epochs a 3 σ flux density upper limit is provided.

DISENTANGLING THE AFTERGLOW AND KILONOVA EMISSION
The kilonova AT2017gfo was characterized by a quasithermal spectrum peaking in the optical/UV and rapidly evolving toward redder wavelengths (e.g. Pian et al. 2017;Drout et al. 2017), in overall agreement with kilonova models (e.g. Kasen et al. 2017;Kasen et al. 2015). The color evolution of GRB 160821B (Figure 2) is consistent with a similar behavior. By using a simple power-law model, Fν ∝ −β to describe the data, we derive β=0.70±0.20 at T0+2 hr, 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 afterglow displays a much redder color, β ≈1.4-1.8, between 2 and 4 days, 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 T0+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. Such color evolution is uncommon for a GRB afterglow, and appears 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 . No significant emission from the kilonova AT2017gfo was detected at X-ray or radio wavelengths (Troja et al. 2017;Hallinan et al. 2017;Margutti et al. 2018;Mooley et al. 2018b). 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 dataset 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) traveling 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 (

Basic Constraints to the Afterglow
The early X-ray afterglow displays a bright and rapidly fading light curve (Figure 3). The measured X-ray flux decreases from 2×10 −10 erg cm −2 s −1 at 250 s post-burst to 2×10 −12 erg cm −2 s −1 at 400 s, implying a temporal slope α ≈9, where FX ∝ 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 modeling.
After the first orbit, the X-ray afterglow follows a simple power-law decay with slope α1 = 0.84±0.08. The condition βX βopt ≈ 0.7 is satisfied for an electron spectral index p ≈ 2.3 when the cooling frequency ν F S c 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 color 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 α1 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) and Troja et al. (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 dataset. They are instead consistent with an uninterrupted power-law decay up to ≈3.5 d and no spectral evolution. The measured spectral index is β = Γ − 1 =0.88±0.14, 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 ∼2.3, and may indicate that a temporal break, likely a "jet-break" (Rhoads 1999;Sari et al. 1999), indeed occurred after T0+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 (Troja et al. 2019;D'Avanzo et al. 2018;Lyman et al. 2018;Margutti et al. 2018). However, for GRB160821B the flat radio-to-optical spectral index βOR ≈0.3 at 4 hr rules out such regime, and shows that ν F S m ≫ 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 a 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, favoring the RS scenario. For a short GRB with T90 ≈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 RS pk 26 µJy and frequency ν RS pk 6 GHz at 3.5 hrs. Given the predicted evolution of the RS characteristic frequency ν RS m ∝ 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 (Figure 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 F S pk 50 µJy, and synchrotron frequency, ν RS pk 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 d and 17 d therefore supports the presence of an earlier jet-break at ≈3.5 d.

Broadband Afterglow Modeling
Based on the preliminary constraints discussed in Sect. 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 observations (≈ T0+0.08 d). The rest of the optical/nIR data (from T0+1 d to T0+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 host V

≈0.
A small amount of extinction may affect our estimates of the spectral slope, but not the evident color evolution (Figure 2).
We followed the same procedure used in Troja et al. where E0,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. Our best fit model and its uncertainty are shown in Figure 4. Whereas at early times (< T0+0.1 d) it provides a good representation of the broadband dataset, at later times it does not naturally account for the drastic color change of the optical/nIR data. It therefore underpredicts the observed emission at nIR and, to a less extent, optical wavelengths.
In Figure 4 we also overplot the template light curves for AT2017gfo in the F 606W and F 110W 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.

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 180). In Figure 5, 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 ). The observed emission resembles more closely the color 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.
The complexity of these systems, in particular the poorly known nuclear physics involved, lead to large uncertainties in the modeling 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, LnIR ≈2×10 39 erg s −1 , and timescales, t pk 3 d, imply a low ejecta mass Mej 0.006 M⊙ and high velocity vej 0.05c for an opacity κ ≈ 10 g cm −3 . The ejecta mass is substantially lower than the values inferred for other GRB kilonovae (Jin et al. 2016;Yang et al. 2015;Tanvir et al. 2013), and comfortably within the range of dynamical ejecta from double NS mergers.  Galama et al. 1998) is also shown for comparison. Errors are 1 σ, downward triangles are 3 σ upper limits. For plotting purposes, optical r and z data were rescaled using the observed colors (Fig. 2) in order to match the F606W and F110W filters, respectively.
Our analysis also finds evidence for an early blue excess, although with larger uncertainties. It is suggestive that the luminosity and timescale of this blue component are consistent with the early optical emission in AT2017gfo (Figure 4, top right panel). The blue color and early onset require a larger mass (Mej ≈0.01 M⊙) of lanthanide-poor material, produced, for example, by the merger remnant.

Effects of a long-lived NS
The merger of two NSs can lead either to a stellarmass BH or to a hypermassive highly magnetized NS (Giacomazzo & Perna 2013). The latter is thought to significantly affect both the kilonova colors and the afterglow evolution through its continuous energy injection and strong neutrino irradiation (Kasen et al. 2015;Gao et al. 2015;Lippuner et al. 2017;Radice et al. 2018). Indeed, the red colors of GW170817/AT2017gfo and its smooth afterglow light curve, mostly consistent with a standard FS emis- sion, 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 Metzger et al. 2018;Yu et al. 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 (Figure 3) often interpreted as a signature of a magnetar engine (Fan & Xu 2006;Troja et al. 2007;Liang 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, LX ≈5×10 47 erg s −1 , 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 P0 ≈9 ms and a dipolar magnetic field B0 ≈3×10 16 G, in overall agreement with the assumption of a magnetar. 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 color evolution of GRB160821B supports the presence of a red kilonova evolving on timescales 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.

ENVIRONMENT
GRB160821B is associated to a bright late-type galaxy (Figure 1). Its spectrum and disturbed morphology are consistent with a star-forming spiral galaxy, possibly interacting with a minor satellite system. At the measured offset of 5.7 ′′ , the probability of a chance alignment between the GRB and the galaxy is ≈2%. However, we find that the overall galaxy properties suggest an environment typical of short GRBs, thus strenghtening the case for a physical association.
The galaxy's absolute B-band magnitude is MB ≈ -19.9 AB, approximately 0.7L * B for a late-type galaxy (Zucca et al. 2009), and its colors are consistent with those of an irregular galaxy (Fukugita et al. 1995). It is detected in the WISE bands with W 1(3.4 µm) = 20.52±0.10 and W 2(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 [OII] line luminosity gives a star formation rate SFR([O II]) 1.5 M⊙ yr −1 (Kennicutt 1998), consistent with the estimate derived from the UV luminosity, Mw2 ≈-18.5 AB mag. Line ratios, log ([OIII]/Hβ)≈0.4 and log ([OIII]/[OII])≈-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 R23 (Pagel et al. 1979) and O32 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).
At a redshift of z=0.161, the projected offset of 5.7 ′′ corresponds to ≈16 kpc, 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 et al. 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 −1 , a minimum kick velocity v kick 80 km s −1 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 −1 are needed to explain the large offsets of some short GRBs (Behroozi et al. 2014).

CONCLUSIONS
GRB160821B is a nearby (z=0.161) 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 toward red colors, 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 dataset 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 R. Clavero and S. Dichiara for assistance, B. Metzger and J. Barnes for help in the preliminary modeling. 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. 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). RSR acknowledges support by Italian Space Agency (ASI) through the Contract n. 2015-046-R.0 and by AHEAD the European Union Horizon 2020 Programme under the AHEAD project (grant agreement n. 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 N. 2017-14-H.0). GR acknowledges the support of the University of Maryland through the Joint Space Science Institute Prize Postdoctoral Fellowship.