The afterglow and kilonova of the short GRB 160821B

TL;DR

This paper presents multi-wavelength observations and modeling of GRB 160821B, revealing a structured jet and a kilonova with properties indicative of a binary neutron star merger, including ejecta mass and velocity estimates.

Contribution

It provides detailed broadband modeling of the afterglow and kilonova of a short GRB, highlighting the jet structure and ejecta characteristics associated with neutron star mergers.

Findings

Structured jet observed slightly off-axis.

Kilonova with rapid red evolution and low luminosity.

Ejecta properties suggest a neutron star merger progenitor.

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 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.

Tables

Table 1. Table 1: Observations of GRB160821B.

Date	T-T0	Telescope	Instrument	Exposure	Band	AB mag	Flux density44footnotemark: 4
MJD	(d)			(s)			($\mu$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$\pm$0.10	3.4$\pm$0.3
57622.017	0.080	GTC	OSIRIS	270	$i$	22.39$\pm$0.07	4.3$\pm$0.3
57622.020	0.083	GTC	OSIRIS	180	$z$	22.28$\pm$0.06	4.7$\pm$0.3
57623.006	1.07	WHT	ACAM	1440	$r$	23.83$\pm$0.25	1.2$\pm$0.3
57623.027	1.09	WHT	ACAM	1680	$z$	23.6$\pm$0.3	1.4$\pm$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	$K_s$	$>$23.3	1.7
57623.958	2.021	GTC	OSIRIS	800	$g$	25.67$\pm$0.15	0.22$\pm$0.03
57623.965	2.028	GTC	OSIRIS	720	$r$	25.12$\pm$0.12	0.36$\pm$0.04
57623.973	2.036	GTC	OSIRIS	450	$i$	24.56$\pm$0.12	0.58$\pm$0.06
57623.980	2.043	GTC	OSIRIS	420	$z$	24.31$\pm$0.17	0.72$\pm$0.11
57625.564	3.627	HST	WFC3	2484	$F606W$	26.02$\pm$0.06	0.157$\pm$0.009
57625.631	3.694	HST	WFC3	2397	$F160W$	24.53$\pm$0.08	0.57$\pm$0.04
57625.697	3.760	HST	WFC3	2397	$F110W$	24.82$\pm$0.05	0.44$\pm$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	$K_s$	24.0$\pm$0.4	0.9$\pm$0.3
57626.922	4.985	GTC	OSIRIS	800	$r$	26.49$\pm$0.20	0.101$\pm$0.019
57629.402	7.465	Keck I	MOSFIRE	145	$K_s$	$>$23.9	0.9
57630.321	8.383	Keck I	MOSFIRE	110	$K_s$	$>$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$\pm$0.3	0.028$\pm$0.008
57632.383	10.446	HST	WFC3	2397	$F110W$	26.9$\pm$0.4	0.07$\pm$0.02
57632.449	10.512	HST	WFC3	2397	$F160W$	26.6$\pm$0.3	0.08$\pm$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 $\pm$ 5
57623.06	1.13	VLA	–	3600	C	–
57632.01	10.07	VLA	–	6480	X	–
57639.04	17.10	VLA	–	6460	X	–
X-ray
57721.995	0.057	Swift	XRT	185	0.3 – 10 keV	–	0.15${}_{-0.06}{}^{}{}_{}{}^{+0.08}$
57722.001	0.063	Swift	XRT	566	0.3 – 10 keV	–	0.05${}_{-0.02}{}^{}{}_{}{}^{+0.03}$
57722.008	0.070	Swift	XRT	784	0.3 – 10 keV	–	0.05$\pm$0.02
57722.064	0.126	Swift	XRT	363	0.3 – 10 keV	–	0.08${}_{-0.03}{}^{}{}_{}{}^{+0.04}$
57722.068	0.130	Swift	XRT	396	0.3 – 10 keV	–	0.07${}_{-0.03}{}^{}{}_{}{}^{+0.04}$
57722.074	0.136	Swift	XRT	865	0.3 – 10 keV	–	0.037${}_{-0.014}{}^{}{}_{}{}^{+0.018}$
57722.132	0.195	Swift	XRT	862	0.3 – 10 keV	–	0.035$\pm$0.019
57722.222	0.285	Swift	XRT	1273	0.3 – 10 keV	–	0.025${}_{-0.010}{}^{}{}_{}{}^{+0.014}$
57722.264	0.327	Swift	XRT	584	0.3 – 10 keV	–	0.05$\pm$0.02
57722.278	0.340	Swift	XRT	1584	0.3 – 10 keV	–	0.029$\pm$0.011
57722.356	0.419	Swift	XRT	3908	0.3 – 10 keV	–	0.014$\pm$0.05
57722.962	1.024	Swift	XRT	9008	0.3 – 10 keV	–	(3.6$\pm$0.2)$\times$10-3
57724.725	2.327	Swift	XRT	8777	0.3 – 10 keV	–	4.6$\times$10-3
57625.879	3.942	XMM-Newton	EPIC/PN	10880	0.3 – 10 keV	–	(2.3$\pm$0.3)$\times$10-3
57631.913	9.976	XMM-Newton	EPIC/PN	22665	0.3 – 10 keV	–	(2.9$\pm$1.5)$\times$10-4
57737.570	15.172	Swift	XRT	26000	0.3 – 10 keV	–	10-3

Equations

$$\text{centering}\log E_{0,\mathrm{iso}}(\mathrm{erg})\text{@add@centering}$$

$$\log n({\mathrm{cm}}^{-3})$$

$$\log {ϵ}_e$$

$$\log {ϵ}_B$$

$$p$$

$${\theta }_{core}(\mathrm{rad})$$

$${\theta }_{view}(\mathrm{rad})$$

Full text

††pubyear: 2019††pagerange: The afterglow and kilonova of the short GRB 160821B– The afterglow and kilonova of the short GRB 160821B

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 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; 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; Roberts et al., 2011; Goriely 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 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; 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 et al., 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$ $\gtrsim$140. NS merger ejecta with an electron fraction $Y_{e}$ $\gtrsim$0.25 do not have enough neutrons to push the nuclear chain past the second r-process peak at $A$ $\approx$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$ $\gtrsim$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 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 $\Lambda$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 $\sigma$ 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_{0}$. Based on its duration of $T_{90}$ =0.48$\pm$0.07 s (Lien et al., 2016), it is classified as a short burst. According to the BAT GRB catalogue111https://swift.gsfc.nasa.gov/results/batgrbcat/, its spectrum is relatively soft, and best described ($\chi^{2}$=59 for 59 dof) by a power-law with a high-energy exponential cut-off at $\approx$50 keV. The time-averaged fluence is $(1.10\pm 0.10)\times 10^{-7}$ erg cm*-2* s*-1* in the 15-150 keV energy band. The Fermi GBM Burst catalogue222https://heasarc.gsfc.nasa.gov/W3Browse/fermi/fermigbrst.html reports a fluence of $(2.0\pm 0.2)\times 10^{-7}$ erg cm*-2* s*-1* over the broader the 10-1,000 keV energy band. Using this value and the redshift $z$ $\sim$ 0.1613 from the candidate host galaxy (see Sect. 2.2.6), we derive an isotropic-equivalent energy $E_{\gamma,iso}$ = (1.3 $\pm$ 0.3) $\times$ $10^{49}$ erg in the 10-1,000 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 multi-wavelength campaign, as we describe in detail below. A log of the observations is reported in Table 1.

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 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 repository333http://www.swift.ac.uk/xrt_products/ (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_{\rm H}=5.7\times 10^{20}$ cm*-2* (Willingale et al., 2013). The time-averaged spectrum, from 4 ks to 74 ks, is well described by an absorbed power-law function with photon index ${\Gamma=1.8\pm 0.2}$. Based on this model, we derive a counts-to-unabsorbed flux conversion factor of 4.4$\times$10*-11* erg cm*-2* ct*-1*. The derived X-ray fluxes are reported in Table 1. No evidence for intrinsic absorption was found at the galaxy’s redshift of z$$\sim0.1613, and we place a 3 $\sigma$ upper limit ${N_{{\rm H},z}<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 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. 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.

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 $\approx$55% for the first observation and $\approx$30% for the second observation. For our analysis we selected only events with FLAG=0 and PATTERN$\leq$4 and PATTERN$\leq$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 ″, 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 modeled the spectrum as described in Sect. 2.1.1 with an absorbed power-law function of photon index ${\Gamma=1.88\pm 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 $\sigma$ 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 hr 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 $\approx$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$\sigma$ 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 ($\approx$ $T_{0}$+2 d), a GTC (+OSIRIS) spectrum (3 $\times$ 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.

2.2.3 William Herschel Telescope (WHT)

Optical imaging with the Auxiliary port CAMera (ACAM) started on 2016-08-22 at 23:55 ($\approx T_{0}$ + 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*′′*. The data were reduced within IRAF555IRAF 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. following standard procedures (e.g., bias subtraction, flat-fielding, etc.). Aperture photometry was performed using the photutils666https://photutils.readthedocs.io/en/stable/ 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 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 AstroConda777AstroConda 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/ 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, determined from the optical frames, is R.A. (J2000) = 18:39:54.56, Dec (J2000) = +62:23:30.35 with a 1 $\sigma$ uncertainty of 0.04″, and lies 5.7″ from the center of a bright face-on spiral (Figure 1, which is likely the GRB host galaxy. At a redshift z$\sim$0.1613, this corresponds to a projected physical offset of 16.40$\pm$0.12 kpc.

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. (2017a). A possible detection was found during the first epoch at $T_{0}$+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. 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 ($\approx$3.5 $\sigma$) at the GRB position. A magnitude of $K_{s}$=22.12$\pm$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.

2.2.6 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 $\sim$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 $\lambda_{obs}$ $\approx$ 4329, 5814, 5759, and 5647 Å, associated with [O II], $H{\beta}$ and [O III] transitions, respectively, and identify the galaxy to be at $z$ = 0.1613$\pm$0.0004, consistent with the preliminary estimate of Levan et al. (2016). The bright line at 7622 Å, associated with the Balmer $H{\alpha}$ line, falls within the telluric A band. Lines properties were derived by modeling 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 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 flux of $\approx$26 $\mu$Jy, consistent within the uncertainties with the value reported in Fong et al. (2016). By performing a 2D Gaussian fit to the source region we obtain a position of R.A. (J2000) = 18:39:54.56, Dec (J2000) = 62:23:30.32 with an error of 0.07 arcsec and 0.08 arcsec, respectively. This is consistent with the optical position (Sect. 2.2.4). For the other epochs a $3~{}\sigma$ 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 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, 2015). The color evolution of GRB 160821B (Figure 2) is consistent with a similar behavior. By using a simple power-law model, $F_{\nu}\propto^{-\beta}$ to describe the data, we derive $\beta$=0.70$\pm$0.20 at $T_{0}+$2 hr, rather standard for an afterglow spectrum with $\nu_{m}<\nu<\nu_{c}$ (Granot & Sari, 2002). Data from the first night of observations show a flatter spectrum, $\beta$ $\approx$0.30, although with larger uncertainties. The optical/nIR counterpart displays a much redder color, $\beta\approx$1.4–1.8, between 2 and 4 days, and then return to a more typical value of $\beta$=1.1$\pm$0.30. Evidence of a possible spectral break at $\lambda$ $\approx$10,000 Å is seen in the GTC OSIRIS and CIRCE data at $T_{0}$+2 d, but no longer visible at later times. Whereas spectral breaks are observed in many GRB afterglows, the change in spectral slope $\Delta\beta$ $\gtrsim$1.0 implied by our observation is hard to reconcile with standard afterglow theory. Substantial color variation is often observed during the early ($\lesssim$103 s) afterglow phases (e.g. Melandri, et al., 2017), but is uncommon on longer timescales (Li, et al., 2018). The observed color evolution in GRB 160821B appears atypical for a GRB afterglow, and instead consistent with the onset of a kilonova peaking in the optical at $\lesssim$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).

No significant emission from the kilonova AT2017gfo was detected at X-ray or radio wavelengths (Troja et al., 2017; Hallinan 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 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)\propto E^{-p}$, and emit their energy via synchrotron radiation. The resulting broadband spectrum is characterized by four quantities: the self-absorption frequency $\nu_{a}$, the synchrotron frequency $\nu_{m}$, the cooling frequency $\nu_{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); 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 (Figure 3). The measured X-ray flux decreases from 2$\times$10*-10* erg cm*-2* s*-1* at 250 s post-burst to 2$\times$10*-12* erg cm*-2* s*-1* at 400 s, implying a temporal slope $\alpha$ $\approx$9, where $F_{X}\propto t^{-\alpha}$. 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 $\alpha_{1}=0.84\pm$0.08. The condition $\beta_{X}\gtrsim\beta_{\rm opt}\approx 0.7$ is satisfied for an electron spectral index $p\approx 2.3$ when the cooling frequency $\nu_{c}^{FS}$ lies within or above the X-ray band. Given its slow temporal evolution, $\nu_{c}\propto 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 $\alpha_{1}$ is shallower than the predicted value $3\beta/2$ $\approx$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 ($\lesssim$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 $\approx$4 d (3.5 d rest-frame) and no spectral evolution. The measured spectral index is $\beta=\Gamma-1$ =0.88$\pm$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 $\alpha_{2}\sim$2.3, and may indicate that a temporal break, likely a “jet-break” (Rhoads, 1999; Sari et al., 1999), indeed occurred after 3.5 d.

Radio observations show a fading afterglow since early times. For standard FS emission, this may indicate that $\nu_{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 $\beta_{OR}\approx$0.3 at 4 hr rules out such regime, and shows that $\nu_{m}^{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 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 $T_{90}\approx$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_{pk}^{RS}\gtrsim$26 $\mu$Jy and frequency $\nu_{pk}^{RS}\lesssim$6 GHz at 3.5 hrs. Given the predicted evolution of the RS characteristic frequency $\nu_{m}^{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\gtrsim$21 mag at $\approx$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_{pk}^{FS}\lesssim$50 $\mu$Jy, and synchrotron frequency, $\nu_{pk}^{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 d and 17 d therefore supports the presence of an earlier jet-break.

When viewed within the effective core opening angle $\theta_{core}$, that is the viewing angle $\theta_{view}\lesssim\theta_{core}$, the afterglow of a 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_{core}^{8/3}$, with deviations as $\theta_{view}$ approaches $\theta_{core}$ (Rossi et al., 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_{core}+\theta_{view})/0.2)^{8/3}$, where $E_{50}$ is the on-axis isotropic equivalent energy in units of $10^{50}$ erg and $n_{-3}$ is the circumburst medium number density in units of $10^{-3}$ cm*-3* (van Eerten et al., 2010). Then, using fiducial values $E_{50}=1$ and $n_{-3}=1$, the observed break in the light curve is consistent with a jet break if $\theta_{core}+\theta_{view}\approx 0.2$ rad.

3.2 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 $\theta_{core}$. We included in the fit the X-ray data, the late radio upper limits, and the first epoch of GTC optical/nIR observations ($\approx$ $T_{0}$+0.08 d). The rest of the optical/nIR data (from $T_{0}$+1 d to $T_{0}$+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}\approx$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. (2018a) and Troja et al. (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:

[TABLE]

[TABLE]

where $E_{\rm 0,iso}$ is the isotropic equivalent blastwave energy, $n$ the ambient density, $\epsilon_{e}$ and $\epsilon_{B}$ the shock microphysical parameters, and $\theta_{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 modeling 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 $\approx$3, that is $\eta_{\gamma}\approx$0.15.

Our best fit model and its uncertainty are shown in Figure 4. In this model early emission at $T_{0}+0.1$ d is dominated by a small patch around the line of sight with Lorentz factor $\sim 8$. The patch widens as the jet decelerates, through the jet break, until at $T_{0}+10$d whole jet is in sight and has a Lorentz factor of roughly $1.5$. Whereas at early times ($<T_{0}$+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 $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 5 reports the spectral energy distribution at three different epochs. Whereas at early times ($T_{0}$+2 hr) 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 modeled with a black-body spectrum of decreasing temperature, from $\approx$4,500 K at 2 d to $\approx$3,300 K at 4 d, and radius $R\approx$1015 cm.

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$ $\gtrsim$140). In Figure 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 color and temporal evolution of AT2017gfo, although less luminous by $\approx$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, $L_{\rm nIR}$ $\approx$2$\times$1039 erg s*-1*, and timescales, $t_{\rm pk}$ $\lesssim$ 3 d, imply a low ejecta mass $M_{ej}$ $\lesssim$0.006 $M_{\odot}$ and high velocity $v_{ej}$ $\gtrsim$0.05$c$ for an opacity $\kappa$ $\approx 10$ g cm*-3*. This result agrees well with the constraints $v\approx R_{bb}/t\approx$0.1-0.2 $c$, derived from simple black body fits (Figure 5). 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.

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 ($M_{ej}$ $\approx$0.01 $M_{\odot}$) 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 et al., 2017). 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 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 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 engine activity. 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 proto-magnetar (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\approx 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}$ $\approx$5$\times$1047 erg s*-1*, and its lifetime, $T$ $\approx$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_{0}\approx$9 ms and a dipolar magnetic field $B_{0}\approx$3$\times$1016 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 color evolution of GRB160821B supports the presence of a red kilonova evolving on timescales similar to AT2017gfo, although less luminous by a factor of $\approx$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\gtrsim$1. Similar distance scales are disfavored by the presence of a kilonova, for which luminosities $>>$1043 erg would lead to implausible properties of the merger ejecta. We consider more natural to link GRB160821B with the nearby bright galaxy (Figure 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″, the probability of a chance alignment between the GRB and the galaxy is $\approx$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 $M_{B}\approx$ -19.9 AB, approximately 0.7$L^{*}_{B}$ for a late-type galaxy (Type 3+4; 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 $W1$(3.4 $\mu$m) = 20.52$\pm$0.10 and $W2$(4.6 $\mu$m) = 21.0 $\pm$ 0.3 AB mag, from which we estimate a stellar mass log ($M/M_{\odot}$) $\approx$ 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]) $\gtrsim$ 1.5 $M_{\odot}$ yr*-1* (Kennicutt, 1998), consistent with the estimate SFR(FUV) $\approx$ 1.2 $M_{\odot}$ yr*-1* derived from the UV luminosity $M_{w2}$ $\approx$-18.5 AB mag (Hao et al., 2011). Line ratios, log ([OIII]/$H\beta$)$\approx$0.4 and log ([OIII]/[OII])$\approx$-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_{23}$ (Pagel et al., 1979) and $O_{32}$ indicators is degenerate, and we can only constrain the upper metallicity branch (Kobulnicky & Kewley, 2004), for which we find 12 + log(O/H) $\approx$ 8.7, consistent with a solar metallicity (Asplund et al., 2009).

The projected physical offset is $\approx$16.4 kpc from the galaxy’s center, 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 $\approx$300 Myr and a velocity dispersion $v_{disp}\approx$120 km s*-1*, a minimum kick velocity $v_{kick}\gtrsim$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).

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 modeling (Sect. 3.2) is also compatible with an origin in a globular cluster (Freire et al., 2001).

6 Conclusions

GRB160821B is a nearby ($z$ $\sim$ 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 $\approx$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.

7 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 modeling. 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 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.