The Combined Ultraviolet, Optical, and Near-infrared Light Curves of the Kilonova Associated with the Binary Neutron Star Merger GW170817: Unified Data Set, Analytic Models, and Physical Implications

The joint detection of gravitational waves and electromagnetic radiation from the binary neutron star merger GW170817 marks the beginning of a new era in observational astrophysics. The merger was detected and localized by the Advanced LIGO and Virgo detectors to a sky region of about 30 deg² at a distance of ≈24–48 Mpc, with inferred component masses of ≈1.36–1.60 and $\approx 1.17\mbox{--}1.36\,{M}_{\odot }$ (90% confidence ranges for the prior of low neutron star spins; Abbott et al. 2017a). A spatially coincident short-duration gamma-ray burst was detected with a delay of 1.7 s relative to the merger time (Abbott et al. 2017c; Goldstein et al. 2017; Savchenko et al. 2017). About 11 hr post-merger several groups (Abbott et al. 2017b; Coulter et al. 2017; Soares-Santos et al. 2017; Valenti et al. 2017a) independently detected an optical counterpart coincident with the quiescent galaxy NGC 4993 at a distance of 39.5 Mpc (Freedman et al. 2001).

Subsequently, multiple ground- and space-based observatories followed up the optical counterpart in the UV, optical, and NIR (hereafter UVOIR), extending to about 30 days post-merger when the location of the source near the Sun prevented further observations. These observations were published in multiple papers that appeared when the detection was publicly announced on 2017 October 16 (Andreoni et al. 2017; Arcavi et al. 2017; Coulter et al. 2017; Cowperthwaite et al. 2017; Díaz et al. 2017; Drout et al. 2017; Evans et al. 2017; Hu et al. 2017; Kasliwal et al. 2017; Lipunov et al. 2017; Pian et al. 2017; Pozanenko et al. 2017; Smartt et al. 2017; Tanvir et al. 2017; Troja et al. 2017; Utsumi et al. 2017; Valenti et al. 2017b). The various papers generally conclude that the UVOIR emission is due at least in part to a kilonova, a quasi-thermal transient powered by the radioactive decay of newly synthesized r-process nuclei and isotopes (Li & Paczyński 1998; Metzger et al. 2010; Roberts et al. 2011; Metzger & Berger 2012; Barnes & Kasen 2013; Tanaka & Hotokezaka 2013). In particular, there is general agreement that the observed light curves require at least two distinct components: a "blue" component that dominates the emission in the first few days, followed by a transition to a "red" component. This multi-component behavior is also seen in optical and NIR spectroscopic observations of the transient (Chornock et al. 2017; Nicholl et al. 2017; Pian et al. 2017; Shappee et al. 2017; Smartt et al. 2017). The blue emission is interpreted to be due to ejecta dominated by Fe-group and light r-process nuclei (atomic mass number $A\lesssim 140$), while the red emission is likely due to ejecta rich in lanthanides and heavy r-process material ($A\gtrsim 140$).

In Cowperthwaite et al. (2017), we modeled photometric data from the Dark Energy Camera (DECam), Swift/UVOT, Gemini, and the Hubble Space Telescope (HST) using the flexible light curve modeling code MOSFiT (Guillochon et al. 2017a). The analysis demonstrated that the UVOIR data cannot be explained by the radioactive decay of ⁵⁶Ni, nor with the associated opacity from Fe-peak elements alone. The data could be well matched by a kilonova model using r-process heating but required at least two distinct components (red and blue) with different opacities, masses, and velocities. A model with a third component (with a higher lanthanide fraction) fit the data equally well (Cowperthwaite et al. 2017). A similar conclusion was reached by several other groups modeling independent sets of observations (e.g., Kilpatrick et al. 2017; Tanaka et al. 2017b). However, given our limited data set, we were unable to break degeneracies between the two- and three-component models.

Following the publication of multiple data sets, we undertake here the first effort to aggregate, homogenize, and model all of the available UVOIR measurements. In total, the UVOIR data set includes 714 individual measurements from 46 different instruments. After collecting the data, we identify measurements that are clearly discrepant from the majority of similar observations, and where possible correct for systematic deviations in order to include as many photometric points as possible. The final unified data set includes 647 measurements. With this extensive data set we revisit the models first explored in Cowperthwaite et al. (2017) with a number of refinements to the physical setup; the model setup is available via the Open Kilonova Catalog⁶ (OKC).

The layout of this Letter is as follows. In Section 2, we discuss the various data sets and describe our approach to standardize the data. In Section 3, we present our model, including additional parameters designed to capture possible asymmetries in the ejecta geometry. We present the results of the model fits in Section 4 and explore their implications in Section 5.

Following the public announcement of the discovery and observations of GW170817, we aggregated the UVOIR photometry available in the literature, which we provide in this Letter and in the OKC. The data span from 0.45 days to 29.4 days post-merger and were collected with 46 instruments in 37 unique filters. This extensive data set represents a departure from most transient light curves, with over 20 observations taken each night on average with fairly complete color coverage during the duration of the event. For each published set of observations, we summarize the instruments and filters used, the details of the photometry methods, and any relevant notes in Table 1. All photometry is reported as AB magnitudes with no correction for Milky Way extinction.

Table 1.  Data Summary

Reference	Bands	Instruments	Telescopes	Photometry	Comments
Andreoni et al. (2017)	g,r,i,C	SkyMapper, 2k2k CCD, 1k2k CCD, NAOS-CONICA, VISIR	SkyMapper, Zadko, VIRT, VLT	image subtraction	Additional data to be published by authors
Arcavi et al. (2017)	V, g, r, i, z, w	Sinistro	LCO 1m/CTIO, SAAO, Siding Spring	image subtraction	Possible template contamination in V-, g-, r-, and i-band; w-band calibrated using r-band SDSS reference stars
Coulter et al. (2017)	B, V, g, r, i	E2V 4k4k CCD	Swope	PSF-fitting
Cowperthwaite et al. (2017)	u, g, r, i, z, Y	DECam	Blanco/CTIO,	image subtraction
Cowperthwaite et al. (2017)	$F336W$, $F475W$, $F625W$, $F775W$, $F850{LP}$, $F110W$, $F160W$, H, Ks	WFC3/UVIS, ACS/WFC, WFC3/IR, Flamingos-2	HST, Gemini-South	PSF-fitting
Díaz et al. (2017)	g, r, i	T80Cam	T80S/CTIO	PSF-fitting
Drout et al. (2017)	B, g, r, i, z, J1, J, H, Ks	IMACS, LDSS-3, FourStar, RetroCam	Magellan, du Pont	PSF-fitting	Used rotated image of galaxy as template
Drout et al. (2017)	U, V, g, I, J, H, Ks	EFOSC2, SOFI, LRIS	NTT, Keck-I	PSF-fitting
Evans et al. (2017)	UVW2, UVM2, U, B, V	UVOT	Swift	host count rate subtraction
Hu et al. (2017)	i	10k10k CCD	AST3-2	image subtraction	Possible template contamination in i-band
Valenti et al. (2017b)	r	Alta U47+	PROMPT5	image subtraction	Pre-existing template
Kasliwal et al. (2017)	$F225W$, $F336W$, B, g, V, r, R, i, I, z, u, J, H, Ks	Flamingos-2, GMOS, WIRC, SIRIUS, ANDICAM, NICFPS, VISIR, WFC3/UVIS	Gemini, Palomar, IRSF, CTIO 1.3m, APO 3.5m, VLT, HST	PSF-fitting, aperture photometry	Subtraction of median-filtered image to remove galaxy
Lipunov et al. (2017)	$B,V,R,W$	MASTER	OAFA, SAAO	image subtraction	Pre-existing template
Pian et al. (2017)	$B,V,g,r,R,i,I,z$	FORS2, ROS2, X-shooter, OmegaCam	VLT, VST, REM	PSF-fitting
Pozanenko et al. (2017)	LUM	4k4k CCD	RC-1000	image subtraction	LUM-band calibrated using r-band reference stars
Shappee et al. (2017)	$B,V,R,I$ , $g,r,i,z$	IMACS, LDSS-3	Magellan	synthetic photometry	Generated synthetic photometry from spectra
Smartt et al. (2017)	$g,r,i,z,y,J,H,K$	GFC, EFOSC2	Pan-STARRS, NTT, 1.5B	image subtraction	Pre-existing template
Smartt et al. (2017)	$U,g,r,i,z,J,H,K$	GROND	MPI/ESO 2.2m	image subtraction	Possible template contamination in GROND K-band
Tanvir et al. (2017)	F475W, r, F606W, i, F814W, $z,Y,J$, F110W, F160W, Ks	VIMOS, WFC-UVIS, FORS, DK1.5, VISTA, NOTCam, WFC-IR, HAWK-I	HST, VLT, HST, DK1, VISTA, NOT	aperture photometry	Local background subtraction; F110W calibrated to J-band.
Troja et al. (2017)	F275W, $B,V$, F475W, F606W, R, I, z, J, H, Ks, F110W, F160W	WFC-IR, WFC-UVIS, GMOS	HST, KMTNet, Gemini	image subtraction
Utsumi et al. (2017)	$V,R,g,r,i,z,J,H,K$	HSC, SIRIUS, MOA-II, MOACam, MOIRCS	B&C, IRSF, Tripol5, Subaru	PSF-fitting	MOACam R-band converted to standard R-band using empirical relationship

Download table as:  ASCIITypeset image

Thanks to the extensive observations from multiple telescopes there is significant redundancy of photometric measurements. This allows us to compare individual data sets to the bulk of the other observations and hence to homogenize and prune the data set. With this approach we find that some corrections are required for three data sets: gri-band data from Arcavi et al. (2017), some K_s-band data from Smartt et al. (2017), and i-band data from Hu et al. (2017). All of these data sets utilized image subtraction to isolate the flux of the transient. However, we find that for the specific filters listed above the resulting light curves were typically dimmer, and faded more rapidly, than the rest of the data. We interpret this as being due to residual emission from the transient in the reference templates, since in each case the template was obtained after the discovery of the source. Using the dates of the template images (I. Arcavi 2017, private communication; Hu et al. 2017; Smartt et al. 2017), we estimate the kilonova brightness for each filter and add this residual flux to the reported photometry. Specifically, we use estimated template magnitudes of: 20.8 (g), 20.9 (r), 20.3 (i), and 20.0 (z) mag to the Arcavi et al. (2017) data set; 19.4 (K_s, GROND data only) mag to the Smartt et al. (2017) data set; and 19.9 (i) mag to the Hu et al. (2017) data set. With these corrections the data are in good agreement with the photometry from other sources (to $\lesssim 0.2$ mag).

We additionally exclude two data sets from our model fitting: the r-band data set from Pozanenko et al. (2017), which was obtained in the LUM filter but calibrated to r-band reference stars; and the w-band from Arcavi et al. (2017), which was similarly calibrated using r-band reference stars. Because the kilonova colors differ so drastically from the comparison stars (see, e.g., Cowperthwaite et al. 2017), these calibrations are unreliable.

Due to the fact that the observations conducted by the Swift UV/Optical Telescope (UVOT) were publicly available, three papers presented independent analyses and photometry of these data (Cowperthwaite et al. 2017; Drout et al. 2017; Evans et al. 2017). However, in our homogenized data set we only use the photometry presented by the Swift team (Evans et al. 2017) without alteration. Early photometry is largely consistent among the three papers to within ≈0.2 mag, although the reported observation times differ by several hours due to different choices of time binning.

Similarly, several teams independently analyzed some Gemini-South FLAMINGOS-2 data (Cowperthwaite et al. 2017; Kasliwal et al. 2017; Troja et al. 2017), some NTT EFOSC2 data (Drout et al. 2017; Smartt et al. 2017), and some HST/WFC3 data (Tanvir et al. 2017; Troja et al. 2017). All of the measurements are listed in Table 3 but marked as repeated observations. The HST/WFC3/F110W data from Tanvir et al. (2017) are re-calibrated to ground-based J-band photometry, so we use the data for these epochs from Troja et al. (2017). For all other epochs with multiple analyses of the same data we take a weighted average of the rep-orted photometry for use in the model fitting, excluding outliers (see below); we report the averaged values in Table 3.

Finally, we identify individual outlying data points through visual inspection and comparison. In total, we find 15 such data points. Three of these are photometry of common data analyzed by multiple teams, so we simply exclude these points from our averaged photometry. We include the 12 other outliers in our modeling, but specifically identify these outliers in Table 3.

The combined data set is listed in Table 3. This table includes the MJD date and phase of each observation; the instrument, telescope, and filter combination; our corrected magnitudes and uncertainties; the correction applied to the original magnitudes (where applicable); a reference to the original paper; and a note indicating if the data were excluded from modeling ("X"), were included in modeling ("*"), represent a repeated reduction of the same observations ("R"), are averaged values from repeated observations ("A"), or are marked as outliers ("O"). We request that any use of the data in this table includes appropriate citation to the original papers, as well as to our compilation.

To properly model this extensive and heterogeneous data set we use the appropriate transmission curve (or close equivalent) for each filter, instrument, and telescope combination.⁷

Photometric modeling of the host galaxy, NGC 4993, suggests that the host environment contributes minimal extinction (Blanchard et al. 2017).⁸ We therefore only include a correction for Milky Way extinction, with $E(B-V)\,=0.105$ mag (Schlafly & Finkbeiner 2011).

In this section, we outline the analytical kilonova model first introduced in Metzger (2017) and implemented in MOSFiT by Villar et al. (2017). This model was also used in Cowperthwaite et al. (2017) to model our own set of observations.

Following decompression from high densities, seed nuclei within the neutron-rich ejecta from a BNS merger undergo rapid neutron capture (r-process) nucleosynthesis (Li & Paczyński 1998; Metzger et al. 2010), and it is the radioactive decay of these freshly synthesized nuclei that powers the kilonova (Metzger 2017). Unlike SNe, which are powered primarily by the radioactive decay of one species (⁵⁶Ni) and therefore undergo exponential decline in their bolometric light curves, kilonovae are powered by the decay of a wide range of r-process nuclei with different half-lives, leading to a power-law decay. At very early times (first few seconds), the energy generation rate is roughly constant as neutrons are consumed during the r-process, but subsequently the r-process freezes out and the energy generation rate approaches a power-law decay, $\propto {t}^{-\alpha }$ with $\alpha \approx 1.3$ (Metzger et al. 2010). The temporal evolution of the radioactive heating rate can be approximated by the parameterized form (Korobkin et al. 2012):

$$\begin{eqnarray}{L}_{\mathrm{in}}(t) & = & 4\,\times \,{10}^{18}{M}_{\mathrm{rp}}\\ & & \times \,{\left[0.5-{\pi }^{-1}\arctan \left(\displaystyle \frac{t-{t}_{0}}{\sigma }\right)\right]}^{1.3}\mathrm{erg}\ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ 1 }$$

where ${M}_{\mathrm{rp}}$ is the mass of the r-process ejecta, and ${t}_{0}=1.3\,{\rm{s}}$ and $\sigma =0.11\,{\rm{s}}$ are constants. Our chosen input luminosity described above neglects any contribution from fall-back accretion on the newly formed remnant. Hydrodynamical simulations suggest that disk winds prevent the fall-back material from reaching the remnant on timescales $\gtrsim 100$ ms (Fernández & Metzger 2013; Metzger 2017); however, some contribution to the bolometric light curve from fall-back accretion is possible on longer (days to weeks) timescales.

Although ${L}_{\mathrm{in}}$ provides the total power of radioactive decay (shared between energetic leptons, γ-rays, and neutrinos), only a fraction ${\epsilon }_{\mathrm{th}}\lt 1$ of this energy thermalizes within the plasma and is available to power the kilonova (Metzger et al. 2010). The thermalization efficiency decreases as the ejecta become more dilute with time, in a manner that can be approximated analytically as (Barnes et al. 2016)

$$\begin{eqnarray}&&{\epsilon }_{\mathrm{th}}(t)=0.36\left[{e}^{-{at}}+\displaystyle \frac{\mathrm{ln}(1+2{{bt}}^{d})}{2{{bt}}^{d}}\right],\end{eqnarray} \tag{ 2 }$$

where a, b, and d are constants of order unity that depend on the ejecta velocity and mass. We use an interpolation of Table 1 of Barnes et al. (2016) for these values.

Assuming that the energy deposition is centrally located and the expansion is homologous, we can use the formalism originally outlined in Arnett (1982) to compute the observed bolometric luminosity (Chatzopoulos et al. 2012):

$$\begin{eqnarray}&&{L}_{\mathrm{bol}}(t)=\exp \left(\displaystyle \frac{-{t}^{2}}{{t}_{d}^{2}}\right)\times {\int }_{0}^{t}{L}_{\mathrm{in}}(t){\epsilon }_{\mathrm{th}}(t)\exp ({t}^{2}/{t}_{d}^{2})\displaystyle \frac{t}{{t}_{d}}{dt},\end{eqnarray} \tag{ 3 }$$

where ${t}_{{\rm{d}}}\equiv \sqrt{2\kappa {M}_{\mathrm{rp}}/\beta {vc}}$, κ is the gray opacity, and $\beta =13.4$ is a dimensionless constant related to the ejecta mass geometric profile. We note that the assumption of a centrally concentrated power source is not necessarily true for kilonovae, as here we assume that the ejecta consists entirely of radioactive r-process material. Relaxation of this assumption should be explored in future work.

We explore multi-component models in which each component has a different opacity corresponding to theoretical expectations for different ejecta compositions. The opacity is largely determined by the fraction of lanthanides in the ejecta, with lanthanide-poor ejecta having a typical opacity of $\kappa \approx 0.5$ cm² g⁻¹ and lanthanide-rich ejecta having a typical opacity of $\kappa \approx 10$ cm² g⁻¹ (Tanaka et al. 2017a). A larger opacity results in a slower light curve evolution and a shift of the spectral energy distribution peak to redder wavelengths. We specifically explore a model with two components ("blue," $\kappa =0.5$ cm² g⁻¹ and "red," κ left as a free parameter), and with three components ("blue," $\kappa =0.5$ cm² g⁻¹; "purple," $\kappa \,=3$ cm² g⁻¹; and "red," $\kappa =10$ cm² g⁻¹; Tanaka et al. 2017a). The purple component corresponds to ejecta with a low, but non-negligible, lanthanide fraction. Each component of the multi-component model is evolved independently, accounting for the unique opacities and therefore diffusion timescales.

To model the multi-band light curves, we assume that each component has a blackbody photosphere with a radius that expands at a constant velocity (${v}_{\mathrm{phot}}\equiv v$, where v is the ejecta velocity). At every point in time, the temperature of each component is defined by its bolometric luminosity and radius, using the Stefan–Boltzmann law. However, when the ejecta cool to a critical temperature (T_c) the photosphere recedes into the ejecta and the temperature remains fixed. The full spectral energy distribution (SED) of the transient is given by the sum of the blackbodies representing each component. The blackbody approximation and temperature floor behavior have both been seen in more sophisticated simulations (Barnes & Kasen 2013); the temperature floor may relate to the first ionization temperature in lanthanide species. The analytic form of the blackbody behavior is

$$\begin{eqnarray}&&{T}_{\mathrm{phot}}(t)=\max \left[{\left(\displaystyle \frac{L(t)}{4\pi {\sigma }_{\mathrm{SB}}^{2}{v}_{\mathrm{ej}}^{2}{t}^{2}}\right)}^{1/4},{T}_{c}\right],\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}{R}_{\mathrm{phot}}(t)=\left\{\begin{array}{ll}{v}_{\mathrm{ej}}t & {\left(\displaystyle \frac{L(t)}{4\pi {\sigma }_{\mathrm{SB}}^{2}{v}_{\mathrm{ej}}^{2}{t}^{2}}\right)}^{1/4}\gt {T}_{c}\\ {\left(\displaystyle \frac{L(t)}{4\pi {\sigma }_{\mathrm{SB}}{T}_{{\rm{c}}}^{4}}\right)}^{1/2} & {\left(\displaystyle \frac{L(t)}{4\pi {\sigma }_{\mathrm{SB}}^{2}{v}_{\mathrm{ej}}^{2}{t}^{2}}\right)}^{1/4}\leqslant {T}_{c}\end{array}\right..\end{eqnarray} \tag{ 5 }$$

3.1. Asymmetric Model

3.2. Fitting Procedure

We model the combined data set using the light curve fitting package MOSFiT (Guillochon et al. 2017a; Nicholl et al. 2017; Villar et al. 2017), which uses an ensemble-based Markov Chain Monte Carlo method to produce posterior predictions for the model parameters. The functional form of the log-likelihood is

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{n}\left[\displaystyle \frac{{({O}_{i}-{M}_{i})}^{2}}{{\sigma }_{i}^{2}+{\sigma }^{2}}-\mathrm{ln}(2\pi {\sigma }_{i}^{2})\right]-\displaystyle \frac{n}{2}\mathrm{ln}(2\pi {\sigma }^{2}),\end{eqnarray} \tag{ 6 }$$

where O_i, M_i, and ${\sigma }_{i}$, are the $i\mathrm{th}$ of n observed magnitudes, model magnitudes, and observed uncertainties, respectively. The variance parameter σ is an additional scatter term, which we fit, that encompasses additional uncertainty in the models and/or data. For upper limits, we use a one-sided Gaussian penalty term.

For each component of our model there are four free parameters: ejecta mass (${M}_{\mathrm{ej}}$), ejecta velocity (${v}_{\mathrm{ej}}$), opacity (κ), and the temperature floor (T_c). We use flat priors for the first three parameters, and a log-uniform prior for T_c (which is the only parameter for which we consider several orders of magnitude). In the case of the asymmetric model, we assume a flat prior for the half-opening angle (θ).

For each model, we ran MOSFiT for approximately 24 hr using 10 nodes on Harvard University's Odyssey computer cluster. We utilized 100 chains until they reached convergence (i.e., had a Gelman–Rubin statistic $\lt 1.1;$ Gelman & Rubin 1992). We use the first $\simeq 80$% of the chain as burn-in. We compare the resulting fits utilizing the Watanabe–Akaike Information Criteria (WAIC; Watanabe 2010; Gelman et al. 2014), which accounts for both the likelihood score and number of fitted parameters for each model.

We fit three different models to the data: a spherical two-component model, a spherical three-component model, and an asymmetric three-component model. The results are shown in Figures 1–5 and summarized in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Kilonova Model Fits

Model	${M}_{\mathrm{ej}}^{\mathrm{blue}}$	${v}_{\mathrm{ej}}^{\mathrm{blue}}$	${\kappa }_{\mathrm{ej}}^{\mathrm{blue}}$	T blue	${M}_{\mathrm{ej}}^{\mathrm{purple}}$	${v}_{\mathrm{ej}}^{\mathrm{purple}}$	${\kappa }_{\mathrm{ej}}^{\mathrm{purple}}$	T purple	${M}_{\mathrm{ej}}^{\mathrm{red}}$	${v}_{\mathrm{ej}}^{\mathrm{red}}$	${\kappa }_{\mathrm{ej}}^{\mathrm{red}}$	T red	σ	θ	WAIC
2-Comp	${0.023}_{0.001}^{0.005}$	${0.256}_{0.002}^{0.005}$	(0.5)	${3983}_{70}^{66}$	⋯	⋯	⋯	⋯	${0.050}_{0.001}^{0.001}$	${0.149}_{0.002}^{0.001}$	${3.65}_{0.28}^{0.09}$	${1151}_{72}^{45}$	${0.256}_{0.004}^{0.006}$	⋯	−1030
3-Comp	${0.020}_{0.001}^{0.001}$	${0.266}_{0.008}^{0.008}$	(0.5)	${674}_{417}^{486}$	${0.047}_{0.002}^{0.001}$	${0.152}_{0.005}^{0.005}$	(3)	${1308}_{34}^{42}$	${0.011}_{0.001}^{0.002}$	${0.137}_{0.021}^{0.025}$	(10)	${3745}_{75}^{75}$	${0.242}_{0.008}^{0.008}$	⋯	−1064
Asym. 3-Comp	${0.009}_{0.001}^{0.001}$	${0.256}_{0.004}^{0.009}$	(0.5)	${3259}_{306}^{302}$	${0.007}_{0.001}^{0.001}$	${0.103}_{0.004}^{0.007}$	(3)	${3728}_{178}^{94}$	${0.026}_{0.002}^{0.004}$	${0.175}_{0.008}^{0.011}$	(10)	${1091}_{45}^{29}$	${0.226}_{0.006}^{0.006}$	${66}_{3}^{1}$	−1116

Download table as:  ASCIITypeset image

For the spherical two-component model we allow the opacity of the red component to vary freely. This model has a total of eight free parameters: two ejecta masses, velocities and temperatures, one free opacity, and one scatter term. We find best-fit values of ${M}_{\mathrm{ej}}^{\mathrm{blue}}={0.023}_{-0.001}^{+0.005}$ ${M}_{\odot }$, ${v}_{\mathrm{ej}}^{\mathrm{blue}}={0.256}_{-0.002}^{+0.005}c$, ${M}_{\mathrm{ej}}^{\mathrm{red}}\,={0.050}_{-0.001}^{+0.001}$ ${M}_{\odot }$, ${v}_{\mathrm{ej}}^{\mathrm{red}}={0.149}_{-0.002}^{+0.001}c$, and ${\kappa }^{\mathrm{red}}\,={3.65}_{-0.28}^{+0.09}$ cm² g⁻¹. Although the model provides an adequate fit, it predicts a double-peaked structure in the NIR light curves at ≈2–5 days that is not seen in the data (Figure 5).

Our best-fitting model, the spherical three-component model, has a total of 10 free parameters: three ejecta masses, velocities and temperatures, and one scatter term. The best-fit values are ${M}_{\mathrm{ej}}^{\mathrm{blue}}\,={0.020}_{-0.001}^{+0.001}$ ${M}_{\odot }$, ${v}_{\mathrm{ej}}^{\mathrm{blue}}={0.266}_{-0.008}^{+0.008}c$, ${M}_{\mathrm{ej}}^{\mathrm{purple}}={0.047}_{-0.002}^{+0.001}$ ${M}_{\odot }$, ${v}_{\mathrm{ej}}^{\mathrm{purple}}={0.152}_{-0.005}^{+0.005}c$, ${M}_{\mathrm{ej}}^{\mathrm{red}}={0.011}_{-0.001}^{+0.002}$ ${M}_{\odot }$, and ${v}_{\mathrm{ej}}^{\mathrm{red}}\,={0.137}_{-0.021}^{+0.025}c$. The parameters in this model are overall comparable to the two-component model in terms of the ejecta masses and velocities of the bluer and redder components, but here the ejecta in the redder component are distributed among the purple and red components. This model underpredicts some of the optical data at $\lesssim 1$ day and overpredicts the late time ($\gtrsim 15$ days) $K,{K}_{{\rm{s}}}$-band data; however, these deviations are less significant than for the two-component model. We additionally explored a version of this model in which the three opacities were allowed to vary freely, but found that these values fell close to our fixed values and did not significantly improve the fit.

Finally, the three-component model with an asymmetric ejecta distribution has a total of 11 free parameters: three ejecta masses, velocities and temperatures, one scatter term, and the opening angle. We find best-fit values of ${M}_{\mathrm{ej}}^{\mathrm{blue}}={0.009}_{-0.001}^{+0.001}$ ${M}_{\odot }$, ${v}_{\mathrm{ej}}^{\mathrm{blue}}\,={0.256}_{-0.004}^{+0.009}c$, ${M}_{\mathrm{ej}}^{\mathrm{purple}}={0.007}_{-0.001}^{+0.001}$ ${M}_{\odot }$, ${v}_{\mathrm{ej}}^{\mathrm{purple}}={0.103}_{-0.004}^{+0.007}c$, ${M}_{\mathrm{ej}}^{\mathrm{red}}={0.026}_{-0.002}^{+0.004}$ ${M}_{\odot }$, ${v}_{\mathrm{ej}}^{\mathrm{red}}={0.175}_{-0.008}^{+0.011}c$, and $\theta ={66}_{-3}^{+1}$ degrees. This model overpredicts the intermediate time (≈5 days) optical photometry and underpredicts the early NIR photometry. Although this model has additional freedom due to the opening angle, the ejecta masses become linked through this additional parameter. Due to the simplicity of the asymmetric model, we do not take the derived parameters and uncertainties at face value, and instead use them as a guide for the effects of asymmetry. We find that an asymmetric ejecta distribution leads to masses that are ≈50% lower than in the spherical case.

We note that the inferred value of θ is consistent with the blue component being visible at an orbital inclination angle of ≈20°–50°, as inferred from a comparison of the GW waveform to the source distance and from an analysis of the radio and X-ray data in the context of an off-axis jet (Abbott et al. 2017d; Alexander et al. 2017; Guidorzi et al. 2017; Hallinan et al. 2017; Margutti et al. 2017; Murguia-Berthier et al. 2017). The relatively large angle is also consistent with the low polarization found by Covino et al. (2017).

Our spherical three-component model realization of highest likelihood (the "best fit") is shown with the complete data set in Figure 1, and its corresponding corner plot is shown in Figure 2. Overall, the model provides a good fit to the complete data set. We find that most parameters are constrained to within ≲10%. The true errors in our models are likely larger, suggesting that the uncertainty is likely dominated by systematic effects (e.g., uncertainty in thermalization efficiency, heating rate, etc.).

We show the individual filters with each of the three components (and their sum) in Figure 4. We find that the blue component dominates across all bands at ≲2–3 days, while the purple component dominates at later times. Because of its low ejecta mass, the reddest component is sub-dominant at all times but contributes necessary flux to the redder bands at late times.

We explore the color evolution of our model compared to that of the kilonova in Figure 3, and again find that the model largely recovers the rapid color evolution, although it slightly deviates from the observed NIR colors at $\gtrsim 12$ days. Finally, we show specific representative filters (r, H, K_s) with a comparison of all three models in Figure 5. Although the differences are subtle, the three-component model provides a statistically better fit to the overall light curves. We stress that the overall success of all three models is remarkable given the extensive scope of the data in time and wavelengths, and the simplifying assumptions in our analytic approach.

Our best-fit three-component model, dominated by an intermediate purple component, is consistent with previous findings (e.g., Cowperthwaite et al. 2017; Chornock et al. 2017; Nicholl et al. 2017). Compared to our previous modeling presented in Cowperthwaite et al. (2017), both the blue and purple ejecta masses and the purple velocity increased by ≈40%. The other parameters remained within $\approx 1\sigma$ of the previously reported values. The uncertainties on the fitted parameters have decreased by ≈10%–50% due to the dramatic increase in the number of data points. Our inferred total ejecta mass of $\approx 0.078$ ${M}_{\odot }$, somewhat higher than the values inferred by several groups based on their individual subsets of the data set we modeled here ($\approx 0.02\mbox{--}0.06\,{M}_{\odot };$ Kasliwal et al. 2017; Kilpatrick et al. 2017; Tanaka et al. 2017b). Additionally, modeling of the optical and NIR spectra indicates that the early blue emission is best described by material with a gradient of lanthanide fraction, with the fraction increasing with time (Chornock et al. 2017; Nicholl et al. 2017). This is consistent with our findings that the purple component begins to dominate the UVOIR light curves at ≈2–3 days post-merger.

The inferred high velocity of the blue ejecta is most naturally explained by relatively proton-rich (high electron fraction, Y_e) polar dynamical ejecta created by the shock from the collision between the merging neutron stars (e.g., Oechslin & Janka 2006; Bauswein et al. 2013; Radice et al. 2016; Sekiguchi et al. 2016). In this scenario, the inferred high ejecta mass ($\approx 0.02\,{M}_{\odot }$) is indicative of a small neutron star radius of ≲12 km when compared to the results of numerical simulations (Bauswein et al. 2013; Hotokezaka et al. 2013; see also Nicholl et al. 2017). Alternatively, the blue ejecta could arise from a neutrino-heated outflow from a hyper-massive neutron star (e.g., Rosswog & Ramirez-Ruiz 2002; Dessart et al. 2009), although the high mass and velocity of the blue ejecta greatly exceed the expectations from a standard neutrino wind and would likely require additional acceleration of the wind by strong magnetic fields (e.g., Metzger et al. 2008).

The red ejecta component could in principle originate from the dynamically ejected tidal tails in the equatorial plane of the binary (e.g., Rosswog et al. 1999; Hotokezaka et al. 2013), in which case the high ejecta mass would require a highly asymmetric merger with a binary mass ratio of $q\lesssim 0.8$ (Hotokezaka et al. 2013). However, the velocity of this component (≈0.1c) is much lower than those typically found in simulations of NS mergers with extreme mass ratios (≈0.2–0.3c; Kilpatrick et al. 2017) potentially disfavoring this explanation. Additionally, our large mass estimate is on the upper end of the dynamical ejecta mass estimated by Abbott et al. (2017e), suggesting that not all of this mass is dynamically ejected.

A more promising source for the red and purple ejecta components is a delayed outflow from the accretion disk formed in the merger (Metzger et al. 2009; Fernández & Metzger 2013; Perego et al. 2014; Just et al. 2015; Siegel & Metzger 2017), for which the outflow velocity is expected to be $\approx 0.03\mbox{--}0.1c$. The relatively high neutron abundance of this matter (${Y}_{e}\lesssim 0.25\mbox{--}0.3$ as needed to synthesize lanthanide nuclei) would be consistent with the moderate amount of neutrino irradiation of the outflow from a black hole accretion disk (Just et al. 2015) but would disfavor a particularly long-lived ($\gtrsim 100$ ms) hyper-massive or supra-massive neutron star remnant (Metzger & Fernández 2014; Murguia-Berthier et al. 2014; Kasen et al. 2015; Lippuner et al. 2017; see also Margalit & Metzger 2017). In this context, the properties of the red/purple ejecta provide evidence for a relatively prompt formation of a black hole remnant.

The asymmetric model indicates a half-opening angle for the blue component of $\theta \approx 66^\circ$. This is consistent with the blue component being visible given the inclination angle of the system inferred both from a comparison of the GW waveform and the distance of the event, and from off-axis jet models of the radio and X-ray light curves (≈20°–50°; Abbott et al. 2017d; Alexander et al. 2017; Margutti et al. 2017). Our simple asymmetric model suggests that the total ejecta mass may be ≈50% smaller than inferred in the spherical model. The effects of other simplifying assumptions, such as the blackbody SED and constant opacities as a function of time and wavelength, should be explored in future work.

Finally, we compare our inferred total ejecta mass to the amount necessary to reproduce the Milk Way r-process production rate using the updated BNS merger rate inferred from Advanced LIGO of ${R}_{0}={1500}_{-1220}^{+3200}$ Gpc⁻³ yr⁻¹ (Abbott et al. 2017d) following a similar methodology as Cowperthwaite et al. (2017) and Kasen et al. (2017). For light r-process nuclei, the primary source of ejecta in our three-component model, the inferred Milky Way production rate is ${\dot{M}}_{\mathrm{rp},{\rm{A}}\lesssim 140}\approx 7\,\times {10}^{-7}\,{M}_{\odot }$ yr⁻¹ (Qian 2000). Combining this with the BNS rate and density of Milky Way–like galaxies ($\approx 0.01$ Mpc⁻³), we estimate the Milky Way rate of BNS mergers as ${R}_{\mathrm{MW}}\,\approx 150$ Myr⁻¹. Thus, the average ejecta mass necessary for a blue/purple kilonova is ${\dot{M}}_{\mathrm{rp},{\rm{A}}\lesssim 140}/{R}_{\mathrm{MW}}\approx 5\times {10}^{-3}$ ${M}_{\odot }$, with an uncertainty of about a factor of $\approx 5$ due to the large range of R₀. For heavy r-process elements (our red component), the Milky Way inferred production rate is ${\dot{M}}_{\mathrm{rp};{\rm{A}}\gtrsim 140}\approx {10}^{-7}\,{M}_{\odot }$ yr⁻¹ (Bauswein et al. 2014). The average ejecta mass necessary for a red kilonova is therefore ${\dot{M}}_{\mathrm{rp},{\rm{A}}\lesssim 140}/{R}_{\mathrm{MW}}\approx 7\times {10}^{-4}\,{M}_{\odot }$, again with an uncertainty of about a factor of 5. In both cases, this order of magnitude estimate is about a factor of 10 times smaller than our estimated ejecta masses for this event, although the rate errors (and potentially lower ejecta masses in the asymmetric case) are large enough to account for the discrepancy.⁹ However, we note that the ratio of red to blue/purple ejecta masses in our model, $\approx 0.16$, is in good agreement with the relative production rates of $A\gtrsim 140$ and $A\lesssim 140$ nuclei in the Milky Way.

If the BNS merger rate from future events is shown to be at the high end of the current estimates, the results inferred here would indicate that a large fraction of synthesized r-process material may remain in the gas phase within the interstellar medium or escape the galaxy entirely via galactic winds (Shen et al. 2015). It may also suggest that the kilonova in GW170817 is an outlier in terms of total r-process material produced. Future events will clarify the population parameters of kilonovae.

We presented the first effort to aggregate, homogenize, and uniformly model the complete UV, optical, and NIR data set for the electromagnetic counterpart of the binary neutron star merger GW170817, allowing us to better determine the likely combinations of parameters responsible for the observed kilonova. We are able to remove systematic offsets from several data sets and to identify outlying data points, providing the community with cleaned and uniform photometry for future analyses. Our key findings are as follows:

• 1.  
  We present 647 photometric measurements from the kilonova accompanying the binary neutron star merger GW170817, spanning from 0.45 to 29.4 days post-merger and providing nearly complete color coverage at all times. We make the homogenized data set available to the public in Table 3, in the OKC, and through https://kilonova.org/.
• 2.  
  The kilonova UVOIR light curves are well fit by a spherically symmetric, three-component model with an overall ejecta mass of $\approx 0.078\,{M}_{\odot }$, dominated by light r-process material ($A\lt 140$) with moderate velocities of ≈0.15c.
• 3.  
  We find evidence for a lanthanide-free component with mass and velocity of $\approx 0.020\,{M}_{\odot }$ and ≈0.27c, respectively. This component is indicative of polar dynamical ejecta, and hence a BNS origin (instead of NS–BH). The large ejecta mass implies a small neutron star radius of $\lesssim 12$ km.
• 4.  
  The mass and velocities of the purple/red components are consistent with a delayed outflow from an accretion disk formed in the merger. This disfavors a long-lived ($\gtrsim 100$ ms) hyper-massive neutron star remnant and provides evidence for relatively prompt formation of a black hole remnant.
• 5.  
  The asymmetric model extension implies that the total ejecta mass may be up to a factor of 2 times lower than for the symmetric model.
• 6.  
  Given the large uncertainties in BNS merger rates, we find that the r-process production rates are comfortably above the Galactic production rate, consistent with the idea that BNS mergers are the dominant source of r-process nucleosynthesis in the universe.

The sheer size of the data set for this event, which was the subject of unprecedented follow-up efforts by the observational astronomy community, represents a departure from typical transient events, allowing for more detailed modeling than typically feasible. Although future observing runs of Advanced LIGO/Virgo will lead to many more kilonova detections, it is likely that this event will remain one of the best-observed objects for years to come due to its vicinity and hence ease of follow-up. Thus, the broad UVOIR data set collected by multiple teams, and aggregated and homogenized here, will be an invaluable resource to explore questions about kilonova phenomenology that may be otherwise intractable using more sparsely sampled data.

We thank the anonymous referee and the larger community for valuable feedback on this work. The Berger Time-Domain Group at Harvard is supported in part by the NSF through grant AST-1714498, and by NASA through grants NNX15AE50G and NNX16AC22G. V.A.V. acknowledges support by the National Science Foundation through a Graduate Research Fellowship. This research has made use of NASA's Astrophysics Data System.

Software: astrocats (Guillochon et al. 2017b), matplotlib (Hunter 2007), MOSFiT (Guillochon et al. 2017a), numpy (Van Der Walt et al. 2011), scipy (Jones et al. 2001).