Spitzer mid-infrared detections of neutron star merger GW170817 suggests synthesis of the heaviest elements

Abstract

1 INTRODUCTION

The discovery of gravitational waves from merging neutron stars, GW170817 (Abbott et al. 2017), offered the first opportunity to directly test the long-standing hypothesis of whether these are the sites of heavy element production (Lattimer & Schramm 1974). The discovery of long-lived infrared emission from GW170817 has provided unequivocal evidence that these are indeed prolific sites of r-process nucleosynthesis (Arcavi et al. 2017; Coulter et al. 2017; Cowperthwaite et al. 2017; Drout et al. 2017; Evans et al. 2017; Kasliwal et al. 2017; Smartt et al. 2017; Soares-Santos et al. 2017). The rapid photometric evolution to the redder wavebands and the sustained luminous infrared emission for a few weeks was consistent with predictions from a suite of kilonova/macronova models (Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010; Barnes & Kasen 2013; Kasen, Badnell & Barnes 2013; Tanaka & Hotokezaka 2013; Rosswog et al. 2017; Wollaeger et al. 2018). The photospheric infrared spectra showed possible evidence of lanthanides such as Neodymium (Kasen et al. 2017).

While there is an emerging consensus in the literature that ≈0.04–0.05 M_⊙ of heavy elements was synthesized and moving at ≈0.1–0.3c, there is much debate on which of the heavy elements were synthesized and whether the abundance distribution matches solar. The solar heavy element distribution has three distinct abundance peaks between atomic mass numbers 70 and 200: the first abundance peak spans 70–88, the second peak spans 120–140, lanthanides span 139–180, and the third peak spans 180–200. The photometric data has been explained both with and without the third r-process peak including the heaviest elements (Smartt et al. 2017; Waxman et al. 2017; Rosswog et al. 2018). The very red spectral energy distribution suggests a non-zero lanthanide fraction, which in-turn would suggest that elements at all three r-process peaks are synthesized (Kasen et al. 2017; Pian et al. 2017).

Despite an intensive campaign by a suite of telescopes worldwide, the infrared monitoring of GW 170817 came to a grinding halt three weeks post-merger due to the target becoming too close to the Sun. Here, we report late-time infrared photometry from the warm Spitzer Space Telescope (Werner et al. 2004), the only telescope that was able to collect infrared data at +43  d and  +74  d after merger despite the proximity to the sun. We use Spitzer observations at +264  d as a reference for image subtraction analysis. Our photometry is inconsistent with that reported in Villar et al. 2018.

2 OBSERVATIONS

We observed GW170817 thrice with the InfraRed Array Camera (IRAC; Fazio et al. 2004) aboard the warm Spitzer Space Telescope at the beginning and end of the first visibility window after explosion and again in the second visibility window (PID 13202, PI Kasliwal). Each epoch constituted a 10  h integration split into 30  s to minimize galaxy core saturation. Due to the larger data volume, each epoch was split into two back-to-back observing blocks. Each epoch had observations in both the 3.6 μm filter and 4.5 μm filter. The first epoch on 2017 Sep 29 was +43  d after merger, the second epoch on 2017 Oct 30 was +74  d after merger, and the third epoch on 2018 May 8 was +264  d after merger. Archival imaging of NGC 4993 also exists from 2014 (PID 10043, PI K. Sheth, we stacked 2014-10-12 and 2014-09-12 data). Data was reduced and mosaiced by the IRAC pipeline.

The IRAC point-spread function is complex: although its core is compact, significant light is scattered far from the position of an object in a complex, asymmetric pattern. The fixed detector/optical diffraction patterns in the IRAC PSF profiles complicate identification of sources in the vicinity of the bright galaxy core since they cannot be easily matched and subtracted between observations taken at different position angles. Specifically, the position angles were 106 deg for Epoch 1, 114  deg for Epoch 2, −65  deg for Epoch 3, and 108  deg for the archival observation. Given the complexity of the underlying galaxy background and the mismatched position angles, we undertook several independent data analysis methods that facilitated multiple consistency checks.

We first describe our preferred method that yielded the highest Signal to Noise (S/N) detections presented in Table 1. To remove the flux from the bright galaxy, we employ the software package GALFIT (Peng et al. 2010), interfaced via a custom python wrapper for galaxy subtraction (Perley et al. 2016). We use the post-Basic Calibrated Data (PBCD) images from the Spitzer archive for each set of observations (split into two observing blocks for each Epoch). We used the first reference observation to fit a Sersic model for the host galaxy profile (simultaneously fitting the 20 brightest stars within the fitting box). The Point Response Function (PRF) files for both Spitzer/IRAC bands were downloaded from the Spitzer website and used as the PSFs, and our fit is restricted to a 58 x 58 arcsec box around the galaxy centroid. A reasonable fit is obtained for this model and most of the galaxy light and foreground starlight is effectively removed (Fig. 1). We see some asymmetric residuals which visually match the structure of the tidal ’shells’ visible in HST imaging of the field as well as some residuals from much fainter foreground stars and background galaxies. At 4.5 μm, we measure a Sersic index of 4.24, a half-light radius of 10.74 arcsec, an integrated magnitude of 12.29 (AB), and an axial ratio of 0.849. At 3.6 μm, we measure a Sersic index of 4.65, a half-light radius of 12.38 arcseconds, axial ratio of 0.838 and integrated magnitude of 11.65 AB. We then repeat this procedure for all other PBCD images, but freeze the fundamental galaxy fit parameters to their values above to ensure a consistent subtraction across every image. We allow the PA and centroid location to vary to allow for astrometric inaccuracies. This effectively cleans the galaxy light and its complex scattering pattern consistently for all images, with the exception of the faint tidal rings.

Figure 1.

Open in new tabDownload slide

Panel 1: Combined Spitzer 3.6 μm and 4.5 μm image, depicting that the faint transient GW 170817 is buried in the bright host galaxy NGC 4993. Panel 2: Subtracting the galaxy light by fitting a GALFIT model clearly shows the red transient in the first epoch image, +43 d after merger. Panel 3a: Proper image subtraction of Epoch 3 reference data from Epoch 1 using the ZOGY algorithm boosts the S/N of our detection of GW170817. Panel 3b: ZOGY subtraction of Epoch 3 reference data from Epoch 2 yields a marginal detection of GW170817. The orientation of all four panels is such that North is up and East is left. The dimensions of the panels are 2.75′ × 2.75′, 1.38′ × 1.38′, 0.69′ × 1.38′, and 0.69′ × 1.38′.

To remove this residual light, and other variations not well-fit by the Sersic model, we used the proper image subtraction routine for optimal transient detection and photometry described in Zackay, Ofek & Gal-Yam 2016 (ZOGY). We used a bright nearby star, 2MASS 13094158–2323149 (Skrutskie et al. 2006), as our reference for the ZOGY PSF and measure the flux-correction factors from the ZOGY output of the difference-image PSF. Following the guidelines in the IRAC Handbook, our correction factors are normalized to a 20 pix radius which includes 100 per cent of the flux (see Table 2). As a consistency check, we compute three apertures (radii of 4, 6, and 10 pix) that match those in the IRAC Handbook and get similar correction factors.

We applied ZOGY to subtract Epoch 3 images from each of the Epoch 1 and Epoch 2 images (Fig. 1). The advantage of using Epoch 3 as a reference for image subtraction over archival data is that it is both comparable in depth to Epoch 1/Epoch 2 and devoid of transient light. As an additional consistency check, we also applied ZOGY to subtract archival Spitzer data of NGC 4993 taken in 2014 from the new data as the archival images are better matched in position angle than the Epoch 3 reference. Additionally, we applied ZOGY to subtract Epoch 2 from Epoch 1 as a consistency check on the difference in flux.

All methods reveal a source at 4.5 μm in Epoch 1 and Epoch 2. All methods reveal no source at 3.6 μm in either Epoch 1 or Epoch 2. Each of these difference images were subject to both PSF-fit and aperture photometry tasks. We summarize the photometry and 3σ upper limits measured from our ZOGY subtractions via PSF-fit photometry in Table 1. (Paranthetical errors added to the measured magnitudes are on account of our analysis of the noisy residuals from stars in the subtracted image suggesting that there may be an additional systematic error of 0.05 mag). We note that our Epoch 1 photometry is inconsistent and brighter by 1 magnitude compared to that reported in Villar et al. 2018. We undertake the following consistency checks.

• We get consistent magnitudes for aperture photometry and PSF photometry. For a 2.5 pix aperture, the aperture magnitudes are 21.99 ± 0.04 for Epoch 1 and 24.14 ± 0.30 AB mag for Epoch 2, consistent with the results from PSF-fit photometry in Table 1.

• The sum of PSF-fit fluxes of the (Epoch 1 – Epoch 2) and (Epoch 2 – Epoch 3) difference images equal the flux in the (Epoch 1 – Epoch 3) difference image. Specifically, the sum of the measured flux in the first two difference images is (5.47 ± 0.14 μJy) +(1.04 ± 0.21 μJy), which is consistent with the measured flux in the (Epoch 1 – Epoch 3) difference image of (6.39 ± 0.21 μJy).

• If we increase the Gaussian FWHM of the PSF-fit to 3.5 pix and apply the appropriate correction factor, we measure a magnitude of 21.93 ± 0.06, consistent with the 2.8 pix FWHM measurement at Epoch 1.

• We get consistent fluxes for Epoch 1 using either the shallow archival reference or the deeper Epoch 3 reference. The PSF-magnitude of Epoch 1 in the archival difference is 21.79 ± 0.09 AB mag, consistent with the late-time difference albeit with larger error bars.

• We get consistent fluxes for Epoch 1 if we directly apply ZOGY to subtract Epoch 3 without first applying the GALFIT-model. We derive 21.94 ± 0.25 mag. The subtraction is noisier by direct subtraction, hence, we prefer the two-step method described above.

• We re-do aperture corrections with a different sky annulus (5–7 pix) and scaling the ZOGY PSF to the standard PRF after re-normalizing the sky. We also take into account colour corrections for this red source by multiplying the measured 4.5 μ m flux by 1.024 (and 3.6 μ m flux by 1.0614) . This gives 21.92 ± 0.09 mag at Epoch 1 and 23.94 ± 0.4 mag at Epoch 2, consistent with Table 1.

All aperture corrections are summarized in an Appendix Table available online.

Converting to flux density, we get F_ν  = 6.43 × 10⁻²⁹ erg s⁻¹ cm⁻² Hz⁻¹ at Epoch 1 and F_ν  = 1.04 × 10⁻²⁹ erg s⁻¹ cm⁻² Hz⁻¹ at Epoch 2. Now, Δν–-L_ν would be a strict lower limit on the total bolometric luminosity. If we assume a power-law ν–L_ν approximation to bolometric, the assumed correction factor is the ratio between the central frequency and bandwidth i.e. a multiplicative factor of 4.3 (since Channel 2 of Spitzer/IRAC spans 3.955 μm to 5.015 μm).

At this late phase, we expect optically thin, nebular conditions and a blackbody approximation with a photosphere is unlikely to be applicable. Nevertheless, we proceed with blackbody calculations as another way to estimate the bolometric correction. The observed Spitzer/IRAC colour ([4.5] – [3.6]) of 1.3  mag suggests a blackbody temperature of 420 K at Epoch 1 (the Epoch 2 colour is not constraining). This suggests a multiplicative bolometric correction factor of ≈16. In the rest of the paper, we assume a ν–L_ν approximation to the bolometric luminosity of 7.8 × 10³⁸ erg s⁻¹ at Epoch 1 and 1.3 × 10³⁸ erg s⁻¹ at Epoch 2.

We check whether synchrotron emission could contribute to the observed flux. Assuming the spectral index presented in Mooley et al. (2018), and a flux density of 44  μJy at 3 GHz measured at the same phase, we estimate that the synchrotron contribution at 4.5μm would be 1.1 × 10⁻³⁰ erg s⁻¹ cm⁻² Hz⁻¹ at Epoch 1. This is ≈60 times smaller than the observed flux density and hence, we conclude that the synchrotron contribution is negligible.

3 IMPLICATIONS ON ABUNDANCES OF R-PROCESS ELEMENTS

At the epochs of the Spitzer observations (t ≳ 40 d) the ejecta of kilonovae are expected to be optically thin to optical/infrared photons. The bolometric luminosity should then be independent of viewing angle and follow the instantaneous radioactive heating rate, |$L(t) \approx M_{\rm ej} \dot{\epsilon }(t) f(t)$|⁠, where M_ej is the ejecta mass, |$\dot{\epsilon }(t)$| the radioactive power per gram, and f(t) the efficiency with which radioactive energy is thermalized. The late-time Spitzer data can thus be used to derive constraints on the ejecta mass and composition that are independent of the complex ejecta opacity and geometry. The main limitation is the uncertain bolometric corrections.

An Appendix Figure available online shows calculations of the radioactive power |$\dot{\epsilon }(t)$| derived from detailed r-process nuclear reaction networks for outflows with a range of physical conditions (initial electron fractions Y_e = 0.05–0.5, expansion velocity of 0.2c, ejecta mass of 0.05 M_⊙ Rosswog et al. 2018). At +43 d, the radioactive power ranges from |$\dot{\epsilon } \approx 0.5 - 2.5 \times 10^{8}\, {\rm erg\, s^{-1}\, g^{-1}}$|⁠. Adopting the νL_ν luminosity at epoch 1 of |$L_{43} = 7.8 \times 10^{38}\, {\rm erg\, s^{-1}}$| and using an efficiency factor f = 0.1 (appropriate for t_e ≈ 30 days) implies an ejecta mass of |$M_{\rm ej} \approx 1.6 - 7.8 \times 10^{-2}\, M_\odot$|⁠. Within large uncertainties, the mass range is consistent with that inferred from analysis of early time observations of GW170817 (Coulter et al. 2017; Drout et al. 2017; Evans et al. 2017; Kasliwal et al. 2017; Smartt et al. 2017; Soares-Santos et al. 2017; Cowperthwaite et al. 2017; Arcavi et al. 2017), and provides additional evidence that the neutron star merger produced a large quantity of radioactive ejecta.

Between the two epochs of Spitzer observations, the luminosity dropped by a factor L₁/L₂ ≈ 6.2 corresponding to a power-law L ∝ t^{−3.4 ± 0.2}. This is steeper than the L ∝ t^−7/3 dependence of statistical distribution of isotopes with power |$\dot{\epsilon } \propto t^{-4/3}$| with inefficient thermalization f(t) ∝ t⁻¹. Alternately, the observed decline can be explained if the efficiency drops even more rapidly, f(t) ∝ t⁻², as suggested by Waxman et al. (2017) (although such a steep dependence of f(t) is not consistent with the numerical thermalization calculations of (Barnes et al. 2016)). Based on late-time optical data, Waxman et al. 2017 and Arcavi 2018 also suggested a similarly steep late-time power-law slope of t⁻³.

It is possible that the decline in luminosity between the two Spitzer epochs is a result of the spectral energy progressively moving out of 4.5μm band, such that the bolometric correction increases with time. If such a colour evolution occurred, the spectrum must have moved redward of |$5\, \mu$|m, as the upper limits in the 3.6μm band rule out a substantial increase of the flux at bluer wavelengths.

If the late time radioactivity is indeed dominated by a single isotope, this provides constraints on the ejecta composition. For merger outflows with electron fractions Y_e ≲ 0.25 the nucleosynthesis proceeds to the 3rd r-process peak (Appendix Figure, available online) and the radioactive power |$\dot{\epsilon }(t)$| steepens at times t ≳ 40 days to a decline rate consistent with the two Spitzer epochs (Fig. 2). For electron fractions Y_e ≳ 0.25, in contrast, the r-process stalls at the first or second r-process peak and the heating rate is flatter at late times due to the presence of long-lived radioisotopes. Thus, the Spitzer data provides conditional evidence that GW170817 produced the heaviest r-process elements.

Figure 2.

Open in new tabDownload slide

Comparing early-time bolometric data (circles, Kasliwal et al. 2017) and late-time Spitzer detections (stars, this paper) with the predicted radioactive luminosity as a function of time (lines). The dashed coloured lines show a luminosity L  = M_ej  |$\dot{\epsilon }(t)$| f(t), where the ejecta mass M_ej  = 0.05 M_⊙, the thermalization efficiency f(t) is from Kasen & Barnes 2018, and the radioactive power |$\dot{\epsilon }(t)$| is from the detailed nuclear reaction network calculations of Rosswog et al. 2018. |$\dot{\epsilon }(t)$| explores a range of electron fraction Ye and expansion velocity from 0.1c to 0.4c. Outflows with Ye<0.25 synthesize the heaviest r-process elements in the second-peak and third-peak and show a steeper late time decline, whereas those with Ye≳0.25 produce relatively lighter elements and have a shallower decline due to the presence of longer lived radioactive isotopes. Also shown is the power law inferred from early-time data (grey solid line) and an analytic estimate of beta decay rates assuming a statistical distribution (magenta solid line; Hotokezaka, Sari & Piran 2017).

Another simple check to this inference is to compare the bolometric light curve to the electron heating rates calculated based on the solar abundance pattern (Fig. 3). The Spitzer detections cannot be explained only by radioactive decay of elements in the first abundance peak as none of them have half-life between between 10–100 d. Abundant elements with relevant half-life include ⁸⁹Sr, ¹²⁵Sn, ¹³¹I, ¹⁴⁰La, ¹⁴¹Ce, ¹⁴³Pr, ¹⁴⁴Ce, ¹⁵⁶Eu, ¹⁸⁸Re, and ¹⁸⁸W. Thus, while the early ground-based data can be explained by many different subsets of r-process elements with different relative abundances, the late-time Spitzer data require the presence of the heaviest r-process elements.

Figure 3.

Open in new tabDownload slide

Fraction of electron heating contributed by various sets of elements as a function of time using the solar abundance pattern including the first peak. While the first peak dominates at early-time, our detection of late-time emission requires elements in the second and third peak.

One caveat here is that we cannot rule out an abundance distribution that cuts off at lanthanum (A = 140). However, nuclear reaction network calculations show that only a narrow range in Ye near 0.25 produce lanthanides but not third peak elements (e.g. Lippuner & Roberts 2015; Holmbeck et al. 2018; Rosswog et al. 2018). Simulations of NS mergers generally show that the Ye of the total ejecta has a broad distribution, and is not narrowly spiked around a single value like Ye  = 0.25 (e.g. Perego et al. 2014; Sekiguchi et al. 2016; Siegel & Metzger 2018; Fernández et al. 2019).

Finally, we note that alpha decay and spontaneous fission of the heaviest nuclei could further enhance the late-time infrared emission. Specifically, Zhu et al. 2018 predict that Californium-254 could boost the mid-IR emission to −13.2 mag at +50 days for a specific model. However, our measured absolute magnitude at 4.5 μm of −11.1 mag AB at +43 d is fainter than their predicted contribution from spontaneous fission of Californium-254 by 2 mag and decays more steeply than their prediction. Perhaps, a lower Californium-254 abundance and different bolometric corrections at the two Spitzer epochs could resolve this discrepancy.

4 CONCLUSIONS

In summary, the Spitzer 4.5 μm observations of GW170817 are the latest and reddest detections of the kilonova emission: 21.88 ± 0.04 mag at +43 d and 23.86 ± 0.22 mag at +74  d. The inferred luminosity suggests a broad composition of r-process elements including the heaviest elements in the second and third abundance peak.

We conclude that mid-IR observations are a critical diagnostic of r-process nucleosynthesis and directly constrain the relative composition. This information cannot be gleaned from the near-IR bands that are accessible from the ground. Currently, warm Spitzer is planned to be online during LIGO-Virgo’s third observing run in 2019. We hope that Spitzer remains online for the entire duration of the observing run. Photometry with Spitzer, in particular, the late-time bolometric luminosity and slope would uniquely constrain the abundance of the heaviest elements. The warm Spitzer analysis is currently limited only by the uncertain bolometric corrections as only two bands are available. This uncertainty would be alleviated by future measurements of the full spectral energy distribution in the mid-IR.

The James Webb Space Telescope (JWST) is planned to be launched soon after the LIGO–-Virgo interferometers plan to attain full sensitivity. Monitoring the SED evolution from 1–25 μm would rule out a chunk of parameter space in heating rates. Specifically, the JWST/NIRCAM F444W filter would be able to detect a GW170817-like event at +74  d as far out as 440  Mpc in less than 10 ks. Furthermore, JWST sensitivity is well-matched to obtaining spectra in the nebular phase that would be a direct diagnostic of the nuclear composition. The JWST/NIRSPEC G395M/F290LP instrument could get a R ≈ 1000 spectrum spanning 2.9–5.1 μm for a GW170817-like event at +43  d as far out as 92 Mpc. The JWST/MIRI F1000W could image at 10μm out to 183 Mpc.

Our Spitzer data decline for GW170817 indicates that elements with a half-life around 14 d could dominate at late-time. Only a handful of the heaviest elements synthesized by the r-process have half-life between 10 d–30 d, e.g. ¹⁴⁰Ba, ¹⁴³Pr, ¹⁴⁷Nd, ¹⁵⁶Eu, ¹⁹¹Os, ²²³Ra, ²²⁵Ra, ²³³Pa, _and²³⁴Th. Given the lower velocities at late-time and lower line-blending, we may even be able to directly read offline identifications and abundances from future JWST spectra.

SUPPORTING INFORMATION

appendix_table.txt

Figure2.pdf

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

This work was supported by the GROWTH (Global Relay of Observatories Watching Transients Happen) project funded by the National Science Foundation Partnership in International Research Program under NSF PIRE grant number 1545949. MMK and DK acknowledge stimulating discussions at KITP; this research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. DK is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract number DE-AC02–05CH11231 and DE-SC0017616, and by a SciDAC award DE-SC0018297. This research was supported in part by the Gordon and Betty Moore Foundation through Grant GBMF5076. MMK and EOO thank the United States-Israel Binational Science Foundation for BSF 2016227. We thank the anonymous referee for constructive feedback. We thank B. Metzger, J. L. Barnes, D. Siegel, G. Martinez-Pinedo, M. Wu, and R. Surman for valuable discussions.

REFERENCES