Spectroscopic identification of r-process nucleosynthesis in a double neutron star merger

References

1 Optical/NIR imaging

Our first observations of the field of SSS17a were carried out with the 60-cm robotic telescope REM located at the ESO La Silla Observatory (Chile) in the g, r, i, z and H bands starting on 2017 Aug 18 at 01:29:28 UT (i.e. 12.8 hours after the GW event). The field was included in the selection we made to carry out targeted observations of catalogued galaxies in the LVC skymap aimed at searching for an optical/NIR counterpart of the GW event starting on 2017 Aug 17 at 23:11:29 UT (i.e. 10.5 hours after the GW event). Following this first detection, we started an extensive follow-up campaign of optical/NIR imaging carried out with an almost daily cadence from about 1.5 to 15.5 days after the time of the GW trigger. These observations were performed using the ESO VLT telescopes equipped with the X-shooter acquisition camera, the FORS2 instrument, and the ESO VST equipped with OmegaCam instrument. The complete log of our photometric observations is reported in Extended Data Table 1. The optical/NIR light curves are shown in Figure 1. Concerning REM and FORS2 imaging, data reduction was carried out following the standard procedures: subtraction of an averaged bias frame and division by a normalized flat frame. The astrometric solution was computed against the USNO-B1.0 catalogue (http://www.nofs.navy.mil/data/fchpix/). Aperture photometry was performed using SExtractor and the PHOTOM package part of the Starlink software distribution (http://starlink.eao.hawaii.edu/starlink). The photometric calibration was achieved by observing Landolt standard fields and the Pan-STARRS catalogue (https://panstarrs.stsci.edu). In order to minimize any systematic effect, we performed differential photometry with respect to a selection of local isolated and non-saturated reference stars. As shown in Extended Data Figure 1, the transient is embedded in the host galaxy light, so that the background around the transient position is highly inhomogeneous, making accurate photometry measurements arduous. In order to minimize the effect of flux contamination from the host light, we fitted it with an analytical profile. The result obtained from the fit was then subtracted from the image in a neighborhood of the transient. This procedure was repeated for each frame. After this subtraction, the background around the transient position is much more uniform, enabling accurate photometric measurements. A dedicated procedure was applied for the reduction and analysis of the wide-field images obtained with the VLT Survey Telescope (VST). The telescope is equipped with OmegaCam , a camera with one square degree field of view (FOV) matched by 0.21 arcsec pixels scale. Data have been processed with a dedicated pipeline for the VST-OmegaCAM observations (dubbed VST-tube). The pipeline searches for new data in the ESO Data archive and, if available, automatically downloads and processes them performing the following main steps: pre-reduction; astrometric and photometric calibration; mosaic production. The OT magnitude, in the AB system, is the PSF fitting magnitude measured on the image after subtracting a model of the galaxy obtained fitting the isophotes with the IRAF/STSDAS task ELLIPSE . The reference catalog used for the absolute photometric calibration is the APASS DR9.

2 FORS2 spectroscopic observations

FORS2 spectra were acquired with the 600B and 600RI grisms, covering the 3500–8600 Å wavelength range. We used in all cases a $1^{\prime\prime}$ slit, for an effective resolution of $R\sim 800-1000$. Spectral extraction was performed with the IRAF software package (IRAF is the Image Reduction and Analysis Facility made available to the astronomical community by the National Optical Astronomy Observatories, which are operated by AURA, Inc., under contract with the US National Science Foundation. It is available at http://iraf.noao.edu.). Wavelength and flux calibration of the spectra were accomplished using helium-argon lamps and spectrophotometric stars. A check for slit losses was carried out by matching the flux-calibrated spectra to our simultaneous photometry (see Extended Data Table 1 and Extended Data Table 2). This shows that the derived spectral shape is robust.

3 X-shooter spectroscopic observations

The cross-dispersed echelle spectrograph, X-shooter, mounted on the VLT, was used to observe the optical/near-infrared counterpart of GW170817. The observing campaign started on the night following the discovery and continued until the source had faded below the detection limit (see Extended Data Table 2) of X-shooter. The observations were carried out using a standard ABBA nodding pattern. Similar position angles of the slit were used for all observations. The position of the slit on the source is shown in Extended Data Figure 1.

The spectroscopic data obtained with X-shooter were managed with the Reflex interface and reduced using version 2.9.3 of the X-shooter pipeline. The reduction cascade consists of bias subtraction, order tracing, flat fielding, wavelength calibration, flux calibration using the spectrophotometric standard EG274 , background subtraction and order rectification – all carried out using the nightly obtained calibration files. A refinement to the wavelength solution was obtained by cross correlating the observed sky spectra with a synthetic sky spectrum, leading to a wavelength solution more accurate than 1 km s-1. Because X-shooter is a cross-dispersed echelle spectrograph, the individual echelle orders are curved across each detector and a rectification algorithm, which correlates neighboring pixels, must be employed. A sampling of 0.2/0.2/0.6 Å per pixel (in the UVB, VIS, and NIR arms, respectively) in the rectified image was chosen to minimize this correlation while conserving the maximal resolving power. The effective resolving power, $R$, of each observation was obtained from fits to unsaturated telluric absorption lines and yielded mean values of 4290/8150/5750 in the UVB/VIS/NIR arms, respectively. This is better than nominal values, owing to a seeing PSF being narrower than the slit width. Immediately following the observations each night, telluric standard stars were observed at an airmass comparable to the target from which the atmospheric transmission spectrum was obtained using Molecfit. Host continuum contamination is visible as a faint background gradient along the slit. An effort has been made to minimize this contamination by using the background regions closest to the target. The images are combined in nightly sets using a weighting scheme based on a moving background variance measure wide enough to avoid it being pixel based and therefore unsuitable for Poisson-noise dominated images. For a subset of the observations, the signal-to-noise (S/N) in the spectral trace is large enough to build a model of the spectral line-spread function to employ an optimal extraction algorithm , but for the majority of the data, an aperture covering the entire trace is used. To establish an accurate flux calibration, slit loss corrections were calculated using the average seeing FWHM of the nightly observations along with the theoretical wavelength dependence of seeing . The slit losses are obtained by integrating a synthetic 2D PSF over the width of the slits and corrections are made accordingly.

4 Foreground dust extinction

We have estimated the intervening dust extinction toward the source using the Na I D line doublet at 5896 Å. Based on the strength of the line in our Galaxy we derive E($B-V$) = 0.09 mag using component D1, E($B-V$) = 0.05 mag using component D2, and E($B-V$) = 0.06 mag using the sum. The Galactic extinction is thus limited to E($B-V$) $<$ 0.1 mag. Similar upper limits on E($B-V$) are obtained from the upper limits on the equivalent widths of the undetected K I 7699 Å absorption line (EW $<0.025$ Å) and undetected 8620 Å diffuse interstellar band (EW $<0.04$ Å). These estimates and limits are marginally consistent with the value of E($B-V$) = 0.11 mag obtained from COBE/DIRBE maps covering that sky region.

5 Spectrum analysis and interpretation

The first epoch X-shooter spectrum was fit with a black-body with temperature of $5000\pm 200$ K. The main deviations from this fit are two absorption-like lines at 8100 and 12300 Å, that evolve with time and become more pronounced in the second spectrum. Altogether, all deviations from a black-body in the first spectrum are below $\sim$10% from 3500 Å to 20000 Å, indicating that the fit is very satisfactory. Moreover, the expansion speed of $0.2c$ we derive from the black-body radius at the epoch of the first spectrum (1.5 days) is compatible with the width of the absorption lines we observe in the second spectrum ($\Delta\lambda/\lambda\sim 0.1-0.2$), confirming that the black-body emission in the first spectrum is highly efficient.

The first 4 X-shooter spectra were compared with kilonova models from Tanaka et al. (2017). The model uses atomic structure calculations for Se (Z = 34), Ru (Z = 44), Te (Z = 52), Ba (Z = 56), Nd (Z = 60), and Er (Z = 68) to construct the atomic data for a wide range of r-process elements. By using two different atomic codes, they confirmed that the atomic structure calculations returned uncertainties in the opacities by a factor of up to $\sim$2. Thereafter, they apply multiwavelength radiative transfer simulations to predict a possible variety of kilonova emission. For each model, the abundance is assumed to be homogeneous in the ejecta, However, a high-$Y_{e}$ component should preferentially dominate near the polar region and low-$Y_{e}$/dynamical component develops in the equatorial region. For each model, the energy release is similar to a power-law ($t^{-1.3}$) owing to the sum of the radioactive decays of various nuclei with different lifetimes. The efficiency of the energy deposition is also taken into account, and the energy deposition rate is somewhat steeper than $t^{-1.3}$ because the gamma-rays can escape without depositing energy.

We emphasize that we have not attempted a real fit of this model to our X-shooter spectra, but have rather looked into an interpretation that was in reasonable agreement. The match is satisfactory only for the first X-shooter spectrum, and not completely satisfactory for the following three. For this reason, we refrained from deriving a light curve model. Infact, in principle, one may fold the synthetic spectral model with the sensitivity curve of any given broad-band filter and integrate the flux in the corresponding band to compare with the observed one. However, the result may be misleading independent of how persuasive it is at face value. The spectral comparison allows one to appreciate in which wavelength ranges the model is effective and in which ones it fails. Integration of the model over a broad wavelength interval cancels the spectral ”memory” and prevents a critical judgment. In other words, since the spectral model is not completely satisfactory, the comparison of synthetic and observed photometry is not significant, although it may appear good.

6 Description of the spectral evolution

The first X-shooter spectrum obtained at $t=1.5$ d after the GW trigger shows an almost featureless, moderately blue continuum. The overall spectral energy distribution is similar to that of early, broad line core collapse SNe. While in general at this relatively low temperature ($\sim$5000 K) SNe typically show strong broad features using the supernova spectral classification tool GELATO a good match is obtained with the early spectra of the type Ib SN2008D/XRF080109. As shown in Extended Data Figure 2, the X-shooter extended spectral range displays, by comparison with the black-body fit (dotted line) the presence of some large scale modulations that are suggestive of multi-component contributions already suggestive of a kilonova event.

In the next two days the spectrum shows a very rapid evolution. The continuum temperature rapidly drops to about 3300K and broad features emerges, with peaks at 10700 Å and 16000 Å. The broad features point to very high expansion velocity and the rapid evolution to a low ejected mass. The combined spectral properties and evolution are unlike those of any known SN types and instead they are very similar to the predicted outcomes of kilonova models.

In the following week the temperature derived from the optical continuum seems to remain roughly constant while the peak at 10700 Å drifts to longer wavelengths (11200 Å at day 6) and decreases in intensity until, at ten days from discovery, the dominant feature in the spectrum is a broad emission centered at about 21000 Å.

7 Host emission analysis

Extending 3–10′′ (0.6 – 2.0 kpc in projection) from the position of the GW counterpart are emission lines formed in the host. The lines are identified as [O II]$\lambda 3726,3729$, H$\beta$, [O III]$\lambda 4959,5007$, H$\alpha$, [N II]$\lambda 6549,6583$ and [S II]$\lambda 6717,6731$, and they exhibit both spatial and velocity structure along the extent of the slit, as shown in Extended Data Figure 3.

From the brightest blob of emission, centered at 6′′ (1.2 kpc in projection) from the source, we measure a receding velocity of $247\pm 15\,{\rm km}\,{\rm s}^{-1}$ relative to the host nucleus (adopting a systemic velocity of NGC 4993 of $2916\pm 15\,{\rm km}\,{\rm s}^{-1}$). Along the spatial direction of the slit, closer to the source, the emission line centroids become more blue-shifted, approaching a recession velocity of $100\,{\rm km}\,{\rm s}^{-1}$ relative to the NGC 4993 systemic velocity. The velocity range ($150\,{\rm km}\,{\rm s}^{-1}$) of the line emission along the slit indicates coherent motion of the gas along the slit. This is further supported by the dust lanes superposed on the host nucleus. The presence of spiral arms was also noted by. A strong [N II]$\lambda$6583 relative to H$\alpha$ combined with a weak H$\beta$ relative to [O III]$\lambda$5007 indicates a radiation field dominated by AGN activity, as also reported previously and supported by the presence of a central radio source. Using the Balmer decrement, the inferred extinction at the position of the line emission is E($B-V$) $=0.21\pm 0.21$.

8 Off–beam jet scenario

GRB170817A had a fluence of $2.2\times 10^{-7}$ erg cm-2 in the 10-1000 keV energy range as observed by the GBM which, at a distance of 40 Mpc, corresponds to a $\gamma$–ray isotropic equivalent energy $E_{\rm iso}\sim 4.3\times 10^{46}$ erg. The peak energy is $E_{\rm peak}=128\pm 48$ keV. The observed $E_{\rm iso}$ is three to four orders of magnitude smaller than the average energy of short GRBs with known redshift.

For illustration let us consider a very simple model: a uniform conical jet of semi-aperture angle $\theta_{\rm jet}$ observed off–beam, i.e at a viewing angle $\theta_{\rm view}>\theta_{\rm jet}$. In this case larger bulk Lorentz factors $\Gamma$ correspond to larger de–beaming factors $b=E_{\rm iso}(0^{\circ})/E_{\rm iso}(\theta_{\rm view})$ for a fixed $\theta_{\rm view}$. Given the small distance of 40 Mpc, and a likely luminosity function decreasing with increasing luminosity (e.g. ), we can assume that the on–axis luminosity of this burst belongs to the low–luminosity tail. For this reason we assume $E_{\rm iso}(0^{\circ})=10^{50}$ erg. Therefore $b=2500$. The probability of a jet oriented at an angle $<\theta_{\rm view}$ is $P(<\theta_{\rm view})=1-\cos\theta_{\rm view}$. A probability of at least $P>10$% implies $\theta_{\rm view}>26^{\circ}$. An off-axis viewing angle larger than $\sim 30^{\circ}$ is also suggested by the expected rate of joint GW and Fermi-GBM detection rescaled to the actual observations. Combining Eq. 2 and 3 from it is possible to estimate the observed energy $E_{\rm iso}$ and peak energy $E_{\rm peak}$ as a function of $\theta_{\rm view}$ and $\Gamma$ for a given $\theta_{\rm jet}$. With $\theta_{\rm view}=30^{\circ}$, $b=2500$ ($E_{\rm iso}(0^{\circ})=10^{50}$ erg) requires $\Gamma=10$ for $\theta_{\rm jet}=10^{\circ}$. The latter is within the currently few estimates of short GRB opening angles and $\Gamma\sim 10$ is within the dispersion of the $\Gamma-E_{\rm iso}$ relation for $E_{\rm iso}(0^{\circ})\sim 10^{50}$ erg. With these values $E_{\rm peak}(0^{\circ})$ turns out to be $\sim$2 MeV. The corresponding comoving frame peak energy would be $\sim$100 keV. If photons with much larger energies are absorbed by pair production we should expect (as observed at $30^{\circ}$) a spectral cutoff at $\sim$650 keV which is larger than the observed peak energy reported by the GBM. Though these values of $E_{\rm peak}(0^{\circ})$ and $E_{\rm iso}(0^{\circ})$ are consistent with those observed in short GRBs, they locate this burst relatively far from the possible spectral-energy correlations of short GRBs.

Extended Data Figure 4 shows the predicted afterglow light curves at 6 GHz, $R$ band and 1 keV. The filled circle shows the X–ray flux at 15 days. The arrows show two representative radio upper limits: at 8.65 days (obtained by co-adding six e-MERLIN observations at 5 GHz) and at 20 days (obtained with MeerKAT at 1.5 GHz). For the model curves the assumed parameters are: $\theta_{\rm jet}=10^{\circ}$, $\theta_{\rm view}=30^{\circ}$, isotropic equivalent kinetic energy $E_{\rm k,iso}=10^{50}$ erg, $\Gamma=10$, a uniform density ISM with $n=2\times 10^{-3}$ cm-3 and standard micro-physical parameters at the shock i.e. $\epsilon_{\rm e}=0.1$, $\epsilon_{\rm B}=0.01$ and electrons’ energy injection power law index $p=2.1$. Standard afterglow dynamics and radiation codes are used. As can be seen the R flux is always below $2\times 10^{-5}$ mJy, corresponding to R$>$28, and therefore orders of magnitude lower than the kilonova emission.

References