The Emergence of a Lanthanide-Rich Kilonova Following the Merger of Two Neutron Stars

We report the discovery and monitoring of the near-infrared counterpart (AT2017gfo) of a binary neutron-star merger event detected as a gravitational wave source by Advanced LIGO/Virgo (GW170817) and as a short gamma-ray burst by Fermi/GBM and Integral/SPI-ACS (GRB170817A). The evolution of the transient light is consistent with predictions for the behaviour of a"kilonova/macronova", powered by the radioactive decay of massive neutron-rich nuclides created via r-process nucleosynthesis in the neutron-star ejecta. In particular, evidence for this scenario is found from broad features seen in Hubble Space Telescope infrared spectroscopy, similar to those predicted for lanthanide dominated ejecta, and the much slower evolution in the near-infrared Ks-band compared to the optical. This indicates that the late-time light is dominated by high-opacity lanthanide-rich ejecta, suggesting nucleosynthesis to the 3rd r-process peak (atomic masses A~195). This discovery confirms that neutron-star mergers produce kilo-/macronovae and that they are at least a major - if not the dominant - site of rapid neutron capture nucleosynthesis in the universe.


Introduction
When compact binary star systems merge, they release copious amounts of energy in the form of gravitational waves (GWs; Abbott et al. 2016Abbott et al. , 2017). If the system is either a binary neutron star (BNS) or a neutron star and stellar mass black hole (NSBH), the merger is expected to be accompanied by various electromagnetic phenomena. In particular, systems of this sort have long been thought to be the progenitors of short-duration gamma-ray bursts (short-GRBs; e.g., Eichler et al. 1989;Nakar 2007), while their neutron-rich ejecta should give rise to a socalled "kilonova" or "macronova" (KN/MN) explosion (Li & Paczyński 1998;Kulkarni 2005;Rosswog 2005;Metzger et al. 2010). Short-GRBs are bright and conspicuous high-energy events. However, since they are thought to be jetted systems, they are expected to be observed for only a subset of such mergers, as the most intense emission from a given merger will usually not intersect our line of sight. Although considerably fainter, KN/MN, which are powered by radioactive decay, emit more isotropically (e.g., Grossman et al. 2014) and peak later than short-GRB afterglows. Thus, they are generally considered to provide the best prospects for electromagnetic (EM) counterparts to GW detections (e.g., Metzger & Berger 2012;Kelley et al. 2013;Fernandez & Metzger 2016;Rosswog et al. 2017).
However, it has been argued that the high opacity of newly synthesized heavy elements in the KN/MN ejecta, particularly lanthanides and actinides, will render them faint in the optical, with emission instead appearing primarily in the near-infrared on timescales of several days Kasen et al. 2013;Tanaka & Hotokezaka 2013). This connects them closely to cosmic nucleosynthesis. The "rapid neutron capture," or "r-process," is responsible for about half of the elements heavier than iron and had traditionally been attributed to core collapse supernovae (Burbidge et al. 1957). A number of recent studies, however, have disfavored supernovae as their conditions were found unsuitable for producing at least the heaviest elements of the "platinum peak" near atomic mass A=195. At the same time, neutron-star mergers have gained increasing attention as a major r-process production site. Lattimer & Schramm (1974) first discussed such compact binary mergers as an r-process site, and since the first nucleosynthesis calculations (Rosswog et al. 1998;Freiburghaus et al. 1999) a slew of other studies (e.g., Goriely et al. 2011;Korobkin et al. 2012;Just et al. 2015;Mendoza-Temis et al. 2015) have confirmed their suitability for the production of the heaviest elements in the Universe.
To date, the most compelling evidence in support of this scenario was provided by the observation of excess infrared light (rest frame 1.2 l~μm) at the location of a short-GRB about a week (in the rest frame) after the burst occurred (GRB 130603B; Tanvir et al. 2013;Berger et al. 2013). Subsequent work has uncovered possible "kilonova" components in several other short-GRBs (Jin et al. , 2016Yang et al. 2015), although other late-time emission processes cannot be ruled out, and the fact that in these instances the excess was in the rest-frame optical bands suggested that it did not originate in lanthanide-rich ejecta.
While the initial focus has been on the extremely neutronrich, low electron fraction (Y e ), "tidal" ejecta component, recent studies (Perego et al. 2014;Wanajo et al. 2014;Just et al. 2015;Radice et al. 2016) have highlighted that this material is likely complemented by higher Y e material that still undergoes r-process nucleosynthesis, but does not produce the heaviest elements (such as gold or platinum) in the third r-process peak. This higher Y e material results from either shocks, neutrinodriven winds, and/or the unbinding of the accretion torus that is formed in the merger. Being free of lanthanides, this material possesses lower opacities and produces earlier and bluer optical transients (e.g., Metzger & Fernández 2014;Kasen et al. 2015). Geometrically, the low-Y e , high-opacity matter is ejected preferentially in the binary orbital plane, while the higher-Y e , low-opacity ejecta are concentrated toward the binary rotation axis. Existing numerical studies suggest that dynamical ejecta have higher velocities ( c 0.1 ; > e.g., Kasen et al. 2015;Rosswog et al. 2017) and could-if viewed edge-on-obscure the windtype ejecta. Therefore, significant viewing-angle effects are expected for the EM signatures of neutron-star mergers.
Here we present the optical and infrared light curve of an explosive transient seen in the hours and days following the detection of a BNS merger by Advanced LIGO/Virgo. We also present optical and near-infrared spectra of the transient. The data show a marked color change from blue to red on a timescale of days as well as conspicuous spectral features, strongly indicative of a KN/MN showing both rapidly evolving blue and more slowly evolving red components.
We use AB magnitudes throughout and, except where otherwise stated, correct for Milky Way foreground extinction according to A 0.338 V = mag from Schlafly & Finkbeiner (2011).

Observations
The discovery of GW170817 by LIGO/Virgo was announced to electromagnetic follow-up partners shortly after the trigger time of 12:41:04 UT on 2017 August 17 (LIGO & Virgo Collaboration 2017a). The potential importance of this event was immediately realized due to its temporal and (within the large error bounds) spatial coincidence with a short-GRB (170817A) detected by Fermi/Gamma-ray Burst Monitor (GBM) at 12:41:06.47 UT (von Kienlin et al. 2017;Goldstein et al. 2017) and also INTEGRAL/SPI-ACS (Savchenko et al. 2017a(Savchenko et al. , 2017b. The existence of a short gamma-ray signal could be interpreted as requiring a close to pole-on viewing angle, but the absence of a normal GRB afterglow in subsequent monitoring (e.g., in X-rays; Evans et al. 2017) instead suggests the possibility of some kind of off-axis emission mechanism, such as may be produced by a shocked cocoon around the primary jet (e.g., Gottlieb et al. 2017;Lazzati et al. 2017).

Imaging
We triggered observations with the European Southern Observatory (ESO) Visible and Infrared Survey Telescope for Astronomy (VISTA; Sutherland et al. 2015) covering two fields within the GW error region and containing high densities of galaxies in the plausible distance regime to have produced such a signal (LIGO & Virgo Collaboration 2017b). Observations began in Chilean twilight at 23:24 UT using the Y (1.02 μm), J (1.25 μm) and K s (2.15 μm) filters. In the second field we identified a bright new point source, visible in all three filters, which was not apparent in prior imaging of the field obtained as part of the VISTA Hemisphere Survey (McMahon et al. 2013). These images were processed using a tailored version of the VISTA Data Flow System that follows the standard reduction path described in González-Fernández et al.  Figure 1). Contemporaneous observations made independently with several optical telescopes also revealed a new source at this location (Allam et al. 2017;Coulter et al. 2017;Valenti et al. 2017), which was designated AT2017gfo (also referred to as SSS17a or DLT17ck).
Subsequently, we monitored AT2017gfo with VISTA at roughly nightly cadence until the field became too difficult to observe due to its proximity to the Sun, after ∼25 days. At later epochs, observations were restricted to the K s -band, which is least affected by twilight observing.
Additionally, we imaged the field with the ESO Very Large Telescope (VLT), the Hubble Space Telescope (HST), the Nordic Optical Telescope (NOT), and the Danish 1.5 m Telescope (DK1.5), including optical observations (a full list of observations and description of photometric measurements is given in Table 1). VLT observations were taken with VIMOS and HAWK-I in the optical (r, z), and infrared (K s ) bands, respectively. Observations were processed through esorex in a standard fashion. HST observations were obtained in the optical (F475W, F606W, and F814W) and infrared (IR; F110W and F160W) and reduced using astrodrizzle to combine, distortion correct, and cosmic-ray reject individual images. The images were ultimately drizzled to plate scales of 0 025 pixel −1 (for UVIS) and 0 07 pixel −1 (for the IR).
For each image, the light from the host galaxy was modeled and subtracted using custom routines in order to aid photometry of the transient, which was performed using the GAIA software. 27 The ground-based J-and K s -bands were calibrated to the 2MASS 28 stars in the field, while the Y-band was calibrated via the relations given in González-Fernández et al. (2017). The optical filters were calibrated to the Pan-STARRS 29 scale. The HST photometry used the standard WFC3 calibrations, 30 apart from the F110W observations, which were also calibrated to the J-band to aid comparison with the other J-band photometry.
Over the first several days AT2017gfo exhibited marked color evolution from blue to red (Figures 2 and 3). Following a slow rise within the first day or so, the optical light declined rapidly from a peak in the first 36 hr, and proceeded to follow an approximately exponential decline (half-life in r-band 40 » hr). The Y-and J-band light curves track each other closely, and again decline following a peak in the first ∼36 hr. By contrast, the K s -band exhibits a much broader peak than the optical, varying by only 20 » % in flux from about 30 hr to 6 days post-merger.
Although there is some evidence for dust lanes in the galaxy, its early-type nature and the absence of host absorption lines (Section 2.2) suggests little dust extinction. Furthermore, the transient is located away from these obviously dusty regions (see Levan et al. 2017 for details of host morphology and transient location). This is supported by the linear polarimetry of the transient, which shows very low levels of polarization (Covino et al. 2017), implying a line of sight dust column in the host galaxy of )and linear polarization). Thus we only correct the photometry for dust extinction in the Milky Way. The measured peak apparent magnitudes are Y 17.22 0 = and K 17.54 0 = . The distance to NGC 4993 is not well established (Hjorth et al. 2017). The heliocentric velocity is 2930 km s −1 (z 0.0098; » Levan et al. 2017), and here we take the distance

Spectroscopy
We observed AT2017gfo with the Multi Unit Spectroscopic Explorer (MUSE) integral field spectrograph on the VLT, which provides optical spectroscopy of both the transient and also the surrounding galaxy (a more detailed description of these data and the analysis of the environment is presented in Levan et al. 2017).
Later spectroscopy was obtained with HST using the Wide-Field Camera 3 Infrared channel (WFC3-IR), with both available grisms, G102 and G141. These observations were pre-reduced by the WFC3 pipeline. The pipeline products were astrometrically calibrated and flat-field corrected, and the diffuse sky background subtracted, using the python-based package grizli. 31 The significant background contamination, caused by the bright host galaxy, was fitted with a two-dimensional polynomial model in a region around the target spectrum, then subtracted using astropy (Astropy Collaboration et al. 2013). The grizli package was then used to optimally extract and combine the spectra from individual exposures. We confirmed these features are robust by comparing the results to extractions from the standard aXe software.
The spectroscopic observations are summarized in Table 2, and the spectra are plotted in Figure 4. The first spectrum at roughly 1.5 days post-merger peaks around 0.6 μm in the optical. The continuum is smooth, with only weak troughs around 0.55 μm, 0.58 μm, 0.75 μm, and 0.8 μm, with a more pronounced break at 0.7 μm. Subsequently, the HST spectra monitor the behavior in the near-infrared, and show that by five days the spectrum is dominated by a prominent peak at ∼1.1 μm. Lesser peaks are apparent at ∼1.4 μm and ∼1.6 μm, and a weak peak at ∼1.22 μm. The breadth of the features is reminiscent of broadline supernova spectra (e.g., Hjorth et al. 2003), and their positions, particularly that of the ∼1.1 μm peak, matches qualitatively the model spectra of Kasen et al. (2013), which adopted opacity based on the lanthanide neodymium. These features appear to be present through the sequence, although they diminish in significance and move toward slightly longer wavelengths. This is consistent with the photosphere moving  deeper with time to slower-moving ejecta as the faster-moving outer layers cool and recombine. Overall, the spectra match well those seen in the extensive ground-based spectroscopic sequence of (Pian et al. 2017), although the absence of atmospheric absorption, compared to ground-based spectra, is particularly beneficial in revealing clearly the 1.4 μm feature.

Interpretation
A natural question is whether any of the light could be due to a synchrotron afterglow, as is generally seen in GRBs. In particular, the absence of early X-ray emission (for 40 Mpc distance, L 5.24 10 X 40 <´erg s −1 at 0.62 days after the trigger; Evans et al. 2017) argues that any afterglow must be faint. A simple extrapolation of the early X-ray limit, assuming conservatively that F 1 n µ n -, gives J 19.9 > . This would at most be a minor contribution to the light observed at early times, so we neglect it here.   We currently lack KN/MN model predictions based on a complete set of likely elements present, and so conclusions are necessarily preliminary. From the large width of the bumps and troughs in the spectrum, which have roughly 0.1 l l D~, we may infer a characteristic ejecta velocity of up to v c 0.1 , assuming that the width is at least partly due to Doppler spreading (see Figure 4). Using this value of the velocity and the light-curve rise time (as well as the decay time of the optical light curves), the ejecta mass M is approximately (Arnett 1980;Metzger et al. 2010) where κ is the opacity. This would suggest that only 10 50 erg of kinetic energy are in the ejecta, despite an energy input of 10 53 erg during the merger. The observed peak isotropic bolometric luminosity of ∼few×10 41 erg s −1 (integrating between u and K, making use of the Swift/UVOT data in Evans et al. 2017) is much higher than predicted for diffusion through an expanding medium following this initial energy input. Continued powering from radioactive decay is required in order to explain the observations, and is consistent with the much slower-decaying infrared light curve. Parametrizing the total heating output of radioactive decay as fMc 2  º (e.g., Metzger et al. 2010), we can estimate f as 10 erg s 0.005 .
The fact that the counterpart was bright, even in the UV, in the first ∼24 hr after the merger (Evans et al. 2017), indicates a high-mass wind with a high Y e , and hence comparatively lowopacity ejecta. This component is likely also dominating the optical emission at early times.
On the other hand, the relatively rapid decline in the J-band compared to the K s -band light suggests that the latter must be dominated, at least from a few days post-merger, by emission from lanthanide-rich dynamical ejecta, in which nucleosynthesis has proceeded to the third r-process peak.

Comparison to Theoretical Models
We compare our observations to the two-component models developed in Wollaeger et al. (2017). These models are computed using the multidimensional radiative Monte Carlo code SuperNu 32 (Wollaeger et al. 2013;Wollaeger & van Rossum 2014;van Rossum et al. 2016) with the set of multigroup opacities produced by the Los Alamos suite of atomic physics codes (Fontes et al. 2015a(Fontes et al. , 2015b. Twocomponent axisymmetric outflow consists of neutron-rich toroidal dynamical ejecta , and a slower spherically symmetric homologous outflow with higher electron fraction, broadly referred to as "wind." The r-process nucleosynthesis and radioactive heating are computed using the nuclear network code WinNet (Thielemann et al. 2011;Korobkin et al. 2012;Winteler et al. 2012) with a reaction rates compilation for the finite range droplet model (FRDM; Möller et al. 1995;Rauscher & Thielemann 2000). Coordinate-and time-dependent thermalization of nuclear energy is calculated using empirical fits developed in Barnes et al. (2016) and Rosswog et al. (2017).
The models are characterized by five parameters: mass and velocity of the dynamical ejecta; mass and velocity of the wind outflow; and inclination angle, which characterizes the remnant orientation. Below we explore a range of these parameters in comparison with the photometric and spectral observations. Figure 5 shows the photometry compared to a few models with varying individual parameters relative to the baseline model with dynamical ejecta parameters m M 0.002 , and orientation angle 20 q = . Figure 5(b) shows the observed spectrum compared to the synthetic spectrum of the baseline model. We conclude that it provides a reasonable fit given the uncertainties in our modeling.
Notice that the wind composition here is moderately neutron rich, with the initial electron fraction Y e =0.27 (denoted as "wind 2" in Wollaeger et al. 2017). Such neutron richness produces a composition of elements grouped around the first r-process peak, and, unlike models with a higher electron fraction (e.g., "wind 1" in Wollaeger et al. 2017), supplies sufficient nuclear heating to explain the observed early . VLT/MUSE and HST grism spectra at five epochs (days post-merger labeled). The later HST observations have been rebinned to reduce the noise. G141 grism spectra are plotted in a lighter line in order to distinguish them from the G102 spectra. The spectra are scaled to match our photometric observations, but have not been corrected for Galactic foreground extinction. Note that since the flux density axis here plots F l , the slopes of the spectra are not directly comparable to Figure 3 emission in the optical bands. Lanthanides in this composition are synthesized only in trace amounts and do not have any noticeable impact on the opacity.
Panel (a) in Figure 5 compares photometric observations in the optical rY-bands and near-IR J K s -bands to the light curves increasing the dynamical ejecta mass, with other parameters set to their default values. Although the fit is not perfect, particularly at later times, the evolution to ∼4 days is reasonably well reproduced, and in J and K s for longer. Higher values of the dynamical ejecta mass m M 0.013 dyn =  lead to a better fit in the K s -band near the peak, however the peak epoch shifts to much later time compared to the observed value. A higher dynamical ejecta mass also produces dimmer light curves in all bands at early times.
Panel ( . The wind-only model qualitatively captures the behavior in the rYJ-bands, but due to the absence of lanthanides it underproduces light in the K s -band. This demonstrates the need to include a secondary, neutron-rich outflow with lanthanides, which can redistribute the emission into the infrared bands. Panel (c) shows the impact of adding a small amount of neutron-rich dynamical ejecta (m M 0.002 dyn =  ). We can see that the infrared bands are reproduced fairly well; however, the addition of a highly opaque component leads to the rapid decay in the rY-bands at late times when compared to the observations. On the other hand, our models only explore a limited parameter space in terms of the composition; the latetime behavior in the optical bands can be cured by tuning the composition of the neutron-rich component. Since our intent here is only to demonstrate viability of the red kilonova hypothesis, adjusting the composition is beyond the scope of this paper.
Panel (d) shows the effect of remnant orientation. Notice that the K s -band becomes insensitive to the orientation after t=6days, indicating that at this epoch the remnant is transparent to the infrared emission and the photosphere disappears. Emission in the optical bands remains sensitive to the orientation, even at t=10days. Nevertheless, in these conditions the local thermal equilibrium (LTE) approximation may not be applicable anymore, so we stop our simulations beyond this epoch.
Panel (e) shows that higher values of nuclear heating (a possibility pointed out in Barnes et al. 2016;Rosswog et al. 2017) lead to only a marginal increase in brightness due to small mass of the ejecta and inefficient thermalization in the dilute dynamical ejecta at late times.

Comparison to Other Claimed Kilonovae
The KN/MN associated with GRB 130603B was observed at 6.94 days rest frame post-burst (corresponding to 7.0 days at 40 Mpc), with an inferred absolute magnitude M j = 15.35 0.2 - (J=17.66 at 40 Mpc). This is roughly a factor of three greater than the luminosity in the J-band at the equivalent epoch for the kilonova accompanying GW170817/ GRB 170817A, and could indicate a higher mass of dynamical ejecta, or additional energy injection from the central remnant (see Kisaka et al. 2016;Gao et al. 2017), in that case.
The candidate KN/MNe discussed by Yang et al. (2015) and Jin et al. (2015Jin et al. ( , 2016 are more difficult to disentangle from the afterglow contribution, but have absolute AB magnitudes (roughly rest frame r-band) around −14 to −15 in the range 3-10 days post-burst, which is again in excess of the emission from AT2017gfo.
These comparisons show that some diversity is to be expected, but it bodes well for the detection of dynamically driven emission components in BNS events at the distances accessible with the advanced GW arrays.

Discussion and Conclusions
Our densely sampled optical and near-IR light curves have revealed the emergence of a red kilonova following the merger of two neutron stars in a galaxy at ∼40 Mpc.
Our modeling of the multi-band light curves indicates the presence of at least two emission components: one with high opacity, and one with low opacity. The former is interpreted as being the "tidal part" of the dynamical ejecta that carry the original, very low electron fraction (Y 0.25 e < ) and result in "strong r-process" producing lanthanides/actinides. This conclusion is supported by near-IR spectroscopy that shows characteristic features expected for high-velocity, lanthaniderich ejecta. The second component avoids strong r-process via a raised electron fraction (Y 0.25 e > ) and may arise from different mechanisms such as neutrino-driven winds and/or the unbinding of accretion torus material. In either case, the ejecta are exposed to high-temperature/high neutrino irradiation conditions for much longer, which drives them to be more proton-rich. Taken together, this lends strong observational support to the idea that compact binary mergers not only produce the "strong r-process" elements, as previously suspected, but also elements across the entire r-process range.
Although the detection of this event in the Advanced LIGO/ Virgo O2 science run is encouraging for future detection rates, the fact that we have not previously seen a similar electromagnetic phenomenon in the low redshift universe is an indication of their rarity. For example, in over 12 years of operation Swift has only located one short-GRB which could potentially be associated with a host galaxy within 150 Mpc, and hence might have been comparable to the AT2017gfo event (Levan et al. 2008). In that case, no counterpart was found despite deep optical and near-IR follow-up that would have easily seen a transient as bright as AT2017gfo, unless it were heavily dust obscured.
The arguments for BNS and NS-BH mergers as heavy r-process nucleosynthesis factories (Vangioni et al. 2016;Rosswog et al. 2017), including from r-process-enriched dwarf galaxies (Beniamini et al. 2016) and the terrestrial abundance of plutonium-244 (Hotokezaka et al. 2015), are broadly in agreement with other observational constraints: from radio observations of Galactic double neutron-star binaries (e.g., O'Shaughnessy & Kim 2010), from the rate and beaming-angle estimates of short-GRBs (Fong et al. 2012), and from population synthesis models of binary evolution (Abadie et al. 2010, and references therein).
A single observed merger during the Advanced LIGO/Virgo O2 science run is consistent with this rate, and likely also consistent with the absence of previous serendipitous kilonova observations. On the other hand, the lack of Swift observations of other γ-ray bursts like this one places an upper limit on the rate of similar events. Future observations will pin down the rate of such events and their typical yields much more precisely, thus establishing their contribution to the heavyelement budget of the Universe.
Finally, we note that if this system was moderately close to being viewed pole-on (e.g., 30  ), as may be suggested by the detection of γ-rays, more highly inclined systems could appear fainter in the optical due to the wind component being obscured by more widely distributed lanthanide-rich ejecta. If this is the case, then near-IR observations could be critical for their discovery. The depth of our short VISTA observations is such that a similar transient would have been easily seen to ∼3 times the distance of NGC 4993, and a more favorable sky location (allowing longer exposures) would have allowed searches to the full BNS detection range (≈200 Mpc) expected for Advanced LIGO at design sensitivity.