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

When compact binary star systems merge, they release copious amounts of energy in the form of gravitational waves (GWs; Abbott et al. 2016, 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 so-called "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 (Barnes & Kasen 2013; 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 $\lambda \sim 1.2$ μ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. 2015, 2016; Yang 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 neutron-rich, 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, neutrino-driven 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 ($\gt 0.1c;$ e.g., Kasen et al. 2015; Rosswog et al. 2017) and could—if viewed edge-on—obscure the wind-type 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}_{{\rm{V}}}=0.338$ mag from Schlafly & Finkbeiner (2011).

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

2.1. 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}_{{\rm{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. (2017) but allows for quick processing of data by using the most current set of calibration frames (mainly flat fields) available at the time. The sky location of the transient was R.A. (2000) = 13:09:48.09, decl.(2000) = −23:22:53.3, approximately 10'' from the center of the S0 galaxy NGC 4993 (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).

Zoom In Zoom Out Reset image size

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}_{{\rm{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 0025 pixel⁻¹ (for UVIS) and 007 pixel⁻¹ (for the IR).

Table 1.  Optical and Near-IR Photometry of AT2017gfo

${\rm{\Delta }}t$ (d)	${t}_{\exp }$ (s)	Telescope/Camera	Filter	$\text{Mag}\,(\text{AB})$
(1)	(2)	(3)	(4)	(5)
8.116	520	HST/WFC3-UVIS	F475W	23.14 ± 0.02
11.300	520	HST/WFC3-UVIS	F475W	24.08 ± 0.05
11.411	600	HST/WFC3-UVIS	F475W	23.96 ± 0.05
1.44	30	VLT/FORS	r	17.69 ± 0.02
2.44	10	VLT/FORS	r	18.77 ± 0.04
3.45	60	VLT/FORS	r	19.28 ± 0.01
4.46	240	VLT/VIMOS	r	19.86 ± 0.01
5.44	20	VLT/FORS	r	20.39 ± 0.03
8.46	600	VLT/VIMOS	r	21.75 ± 0.05
9.46	600	VLT/VIMOS	r	22.20 ± 0.04
10.46	1200	VLT/VIMOS	r	22.45 ± 0.07
11.44	360	HST/WFC3-UVIS	F606W	23.09 ± 0.03
12.44	1200	VLT/VIMOS	r	23.12 ± 0.31
2.459	150	DK1.5	i	18.37 ± 0.03
11.428	560	HST/WFC3-UVIS	F814W	22.32 ± 0.02
2.461	150	DK1.5	z	18.01 ± 0.13
4.451	240	VLT/VIMOS	z	18.73 ± 0.01
8.443	400	VLT/VIMOS	z	20.28 ± 0.03
9.445	400	VLT/VIMOS	z	20.85 ± 0.04
9.462	60	VLT/FORS	z	20.69 ± 0.11
13.440	480	VLT/VIMOS	z	22.30 ± 0.28
19.463	720	VLT/VIMOS	z	23.37 ± 0.48
0.49	120	VISTA/VIRCAM	Y	17.46 ± 0.01
1.47	120	VISTA/VIRCAM	Y	17.23 ± 0.01
2.47	120	VISTA/VIRCAM	Y	17.51 ± 0.02
3.46	120	VISTA/VIRCAM	Y	17.76 ± 0.01
4.46	120	VISTA/VIRCAM	Y	18.07 ± 0.02
6.47	120	VISTA/VIRCAM	Y	18.71 ± 0.04
7.47	120	VISTA/VIRCAM	Y	19.24 ± 0.07
8.46	120	VISTA/VIRCAM	Y	19.67 ± 0.09
9.46	120	VISTA/VIRCAM	Y	20.09 ± 0.14
0.48	120	VISTA/VIRCAM	J	17.88 ± 0.03
0.51	120	VISTA/VIRCAM	J	17.82 ± 0.03
1.46	120	VISTA/VIRCAM	J	17.45 ± 0.01
2.46	120	VISTA/VIRCAM	J	17.66 ± 0.02
3.46	120	VISTA/VIRCAM	J	17.86 ± 0.02
4.46	120	VISTA/VIRCAM	J	18.08 ± 0.03
4.79	298	HST/WFC3-IR	F110W	18.26 ± 0.01
6.47	120	VISTA/VIRCAM	J	18.74 ± 0.04
7.24	298	HST/WFC3-IR	F110W	19.06 ± 0.01
7.46	120	VISTA/VIRCAM	J	19.07 ± 0.08
8.45	120	VISTA/VIRCAM	J	19.69 ± 0.09
9.45	120	VISTA/VIRCAM	J	20.06 ± 0.14
10.46	120	VISTA/VIRCAM	J	20.94 ± 0.35
10.55	298	HST/WFC3-IR	F110W	20.82 ± 0.02
11.46	120	VISTA/VIRCAM	J	21.16 ± 0.40
4.923	298	HST/WFC3-IR	F160W	18.063 ± 0.03
9.427	298	HST/WFC3-IR	F160W	19.600 ± 0.06
10.619	298	HST/WFC3-IR	F160W	20.279  ± 0.09
0.47	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	18.62 ± 0.05
0.50	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	18.64 ± 0.06
1.32	360	NOT/NOTCam	${K}_{{\rm{s}}}$	17.86 ± 0.22
1.46	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	17.77 ± 0.02
2.45	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	17.67 ± 0.03
3.45	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	17.54 ± 0.02
4.45	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	17.60 ± 0.02
6.46	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	17.84 ± 0.03
7.45	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	17.95 ± 0.04
8.45	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	18.25 ± 0.03
9.45	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	18.49 ± 0.05
10.45	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	18.74 ± 0.06
12.46	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	19.34 ± 0.08
14.46	120	VISTA/VIRCAM	${K}_{{\rm{s}}}$	20.02 ± 0.13
17.45	780	VLT/HAWK-I	${K}_{{\rm{s}}}$	20.77 ± 0.13
20.44	1140	VLT/HAWK-I	${K}_{{\rm{s}}}$	21.58 ± 0.06
21.44	1320	VLT/HAWK-I	${K}_{{\rm{s}}}$	21.46 ± 0.08
25.44	600	VLT/HAWK-I	${K}_{{\rm{s}}}$	22.06 ± 0.22

Note. Column (1) is the start time of the observation with respect to the GW trigger time (LIGO & Virgo Collaboration 2017a). Column (5) gives AB magnitudes corrected for Galactic foreground extinction (see the text).

Download table as:  ASCIITypeset images: 1 2

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.²⁷ The ground-based J- and ${K}_{{\rm{s}}}$-bands were calibrated to the 2MASS²⁸ 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²⁹ scale. The HST photometry used the standard WFC3 calibrations,³⁰ 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 $\approx 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}_{{\rm{s}}}$-band exhibits a much broader peak than the optical, varying by only $\approx 20$% in flux from about 30 hr to 6 days post-merger.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 $E(B-V)\lesssim 0.2\,\text{mag}$ (assuming a Milky-Way-like relation between $E(B-V)$ and linear polarization). Thus we only correct the photometry for dust extinction in the Milky Way. The measured peak apparent magnitudes are ${Y}_{0}=17.22$ and ${K}_{0}=17.54$.

The distance to NGC 4993 is not well established (Hjorth et al. 2017). The heliocentric velocity is 2930 km s⁻¹ ($z\approx 0.0098;$ Levan et al. 2017), and here we take the distance to be d = 40 Mpc (distance modulus $\mu =33.01$). Thus, the peak absolute magnitudes from our measurements are ${M}_{Y,0}=-15.79$ and ${M}_{K,0}=-15.47$.

2.2. 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.³¹ 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.

Zoom In Zoom Out Reset image size

Table 2.  Optical and Near-IR Spectroscopy of AT2017gfo

${\rm{\Delta }}t$ (d)	${t}_{\exp }$ (s)	Telescope/Camera	Coverage (μm)
(1)	(2)	(3)	(4)
1.47	2600	VLT/MUSE	0.48–0.93
4.86	1812	HST/WFC3-IR	0.8–1.15 (G102)
4.93	1812	HST/WFC3-IR	1.08–1.7 (G141)
7.27	1812	HST/WFC3-IR	0.8–1.15 (G102)
9.43	1812	HST/WFC3-IR	1.08–1.7 (G141)
10.52	1812	HST/WFC3-IR	0.8–1.15 (G102)
10.65	1812	HST/WFC3-IR	1.08–1.7 (G141)

Note. Column (1) contains the start time of the observation with respect to the GW trigger time.

Download table as:  ASCIITypeset image

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}_{{\rm{X}}}\lt 5.24\times {10}^{40}$ erg s⁻¹ 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}_{\nu }\propto {\nu }^{-1}$, gives $J\gt 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 ${\rm{\Delta }}\lambda /\lambda \sim 0.1$, we may infer a characteristic ejecta velocity of up to $v\sim 0.1c$, 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)

$$\begin{eqnarray*}&&M\sim 5\times {10}^{-3}\,{M}_{\odot }\left(\displaystyle \frac{0.1\,{\rm{g}}\,{\mathrm{cm}}^{-2}}{\kappa }\displaystyle \frac{v}{0.1\,c}\right),\end{eqnarray*}$$

where κ is the opacity. This would suggest that only 10⁵⁰ erg of kinetic energy are in the ejecta, despite an energy input of $\sim {10}^{53}$ erg during the merger.

The observed peak isotropic bolometric luminosity of ∼few × 10⁴¹ erg s⁻¹ (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 $\epsilon \equiv {{fMc}}^{2}$ (e.g., Metzger et al. 2010), we can estimate f as

$$\begin{eqnarray*}&&f\sim {10}^{-6}\displaystyle \frac{{L}_{\mathrm{peak}}}{{10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}\displaystyle \frac{0.005{M}_{\odot }}{M}.\end{eqnarray*}$$

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 low-opacity 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}_{{\rm{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.

3.1. 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 ³² (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, 2015b, 2017). Two-component axisymmetric outflow consists of neutron-rich toroidal dynamical ejecta (Rosswog et al. 2014), 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}_{\mathrm{dyn}}=0.002\ {M}_{\odot }$, ${v}_{\mathrm{dyn}}=0.2\ c$, wind parameters ${m}_{\mathrm{wind}}=0.015\ {M}_{\odot }$, ${v}_{\mathrm{wind}}=0.08\ c$, and orientation angle $\theta =20^\circ$. 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.

Zoom In Zoom Out Reset image size

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 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}_{{\rm{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}_{{\rm{s}}}$ for longer. Higher values of the dynamical ejecta mass ${m}_{\mathrm{dyn}}=0.013\ {M}_{\odot }$ lead to a better fit in the ${K}_{{\rm{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 (b) shows the spherically symmetric wind-only model with varying mass ${m}_{\mathrm{wind}}=0.01\mbox{--}0.1\ {M}_{\odot }$. 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}_{{\rm{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}_{\mathrm{dyn}}=0.002\ {M}_{\odot }$). 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 late-time 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}_{{\rm{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.

3.2. Comparison to Other Claimed Kilonovae

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}_{e}\lt 0.25$) 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, lanthanide-rich ejecta. The second component avoids strong r-process via a raised electron fraction (${Y}_{e}\gt 0.25$) 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 heavy-element budget of the Universe.

Finally, we note that if this system was moderately close to being viewed pole-on (e.g., $\lesssim 30^\circ$), 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.

We thank the staff at ESO, both at Paranal and Garching, for their expert and enthusiastic support of the observations reported here.

We thank the staff at STScI, in particular Tricia Royle, Alison Vick, Russell Ryan, and Neill Reid, for their help in implementing such rapid HST observations.

The observations with VISTA were gathered by the ESO VINROUGE Survey (198.D-2010). Observations also used data from the VISTA Hemisphere Survey (VHS: 179.A-2010).

HST observations were obtained using programs GO 14771 (PI: Tanvir), GO 14804 (PI: Levan), and GO 14850 (PI: Troja).

VLT observations were obtained using programs 099.D-0688, 099.D-0116, and 099.D-0622.

N.R.T., K.W., P.T.O., J.L.O., and S.R. acknowledge support from STFC.

A.J.L., D.S., and J.D.L. acknowledge support from STFC via grant ST/P000495/1.

N.R.T. and A.J.L. have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (grant agreement No. 725246, TEDE, Levan).

I.M. acknowledges partial support from the STFC.

Ad.U.P., C.T., Z.C., and D.A.K. acknowledge support from the Spanish project AYA 2014-58381-P. Z.C. also acknowledges support from the Juan de la Cierva Incorporación fellowship IJCI-2014-21669, and D.A.K. from Juan de la Cierva Incorporación fellowship IJCI-2015-26153.

J.H. is supported by a VILLUM FONDEN Investigator grant (project number 16599).

P.D.A., S.C., and A.M. acknowledge support from the ASI grant I/004/11/3.

S.R. has been supported by the Swedish Research Council (VR) under grant No. 2016-03657_3, by the Swedish National Space Board under grant No. Dnr. 107/16, and by the research environment grant "Gravitational Radiation and Electromagnetic Astrophysical Transients (GREAT)" funded by the Swedish Research council (VR) under Dnr 2016-06012.

P.A.E. acknowledges UKSA support.

The VISTA observations were processed by C.G.F. at the Cambridge Astronomy Survey Unit (CASU), which is funded by the UK Science and Technology Research Council under grant ST/N005805/1.

This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under contract No. DE-AC52-06NA25396.

Based on observations made with the Nordic Optical Telescope (program 55-013, PI Pian), operated by the Nordic Optical Telescope Scientific Association at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofísica de Canarias.

Facilities: HST(WFC3) - Hubble Space Telescope satellite, VISTA(VIRCAM) - , VLT(MUSE - , HAWK-I - , VIMOS - , FORS), NOT(NOTCam) - .