Spiral-wave wind for the blue kilonova

TL;DR

This paper proposes a new mechanism involving spiral density waves in neutron star merger remnants that can produce the blue kilonova observed in GW170817, explaining its early brightness and nucleosynthesis.

Contribution

The study introduces a novel spiral-wave wind mechanism from neutron star merger remnants, supported by detailed numerical simulations, to explain the blue kilonova and r-process element production.

Findings

Spiral density waves generate a wind of about 0.01 solar masses at 0.2c.

Ejected material has electron fraction mostly above 0.25.

The combined ejecta can reproduce observed light curves and r-process abundances.

Abstract

The AT2017gfo kilonova counterpart of the binary neutron star merger event GW170817 was characterized by an early-time bright peak in optical and UV bands. Such blue kilonova is commonly interpreted as a signature of weak $r$-process nucleosynthesis in a fast expanding wind whose origin is currently debated. Numerical-relativity simulations with microphysical equations of state, approximate neutrino transport, and turbulent viscosity reveal a new mechanism that can power the blue kilonova. Spiral density waves in the remnant generate a characteristic wind of mass ${\sim}10^{-2}~M_{\odot}$ and velocity ${\sim}0.2$c. The ejected material has electron fraction mostly distributed above $0.25$ being partially reprocessed by hydrodynamic shocks in the expanding arms. The combination of dynamical ejecta and spiral-wave wind can account for solar system abundances of $r$-process elements and early-time observed light curves.

Full text

1. Introduction

The observation of the kilonova (kN) AT2017gfo (Coulter et al., 2017; Chornock et al., 2017; Nicholl et al., 2017; Cowperthwaite et al., 2017; Tanvir et al., 2017; Tanaka et al., 2017) associated to the binary neutron star (BNS) merger GW170817 (Abbott et al., 2017) provided evidence that the ejection of neutron-rich matter from compact binary mergers is a primary site for $r$-process nucleosynthesis (Lattimer and Schramm, 1974; Li and Paczynski, 1998; Kulkarni, 2005; Rosswog, 2005; Metzger et al., 2010; Roberts et al., 2011; Kasen et al., 2013). In this scenario, the electromagnetic UV/optical/NIR transient is powered by the radioactive decay of the freshly synthesized elements. The NIR luminosity of AT2017gfo peaked at several days after the merger (Chornock et al., 2017), and it is consistent with expectations that the opacities of expanding $r$-process material are dominated by the opacities of lanthanides and possibly actinides (Kasen et al., 2013). The UV/optical luminosity peaked instead in less than one day after the merger (Nicholl et al., 2017), and it originates from ejected material that experienced only a partial $r$-process nucleosynthesis (Martin et al., 2015). A fit of AT2017gfo light curves to a semianalytical two-components spherical model indicates a lanthanide poor (rich) blue (red) component of mass $2.5\times 10^{-2}M_{\odot}$ ($5.0\times 10^{-2}M_{\odot}$) and velocity $0.27$c ($0.15$c) (Cowperthwaite et al., 2017; Villar et al., 2017) (See however (Waxman et al., 2018) for an alternative interpretation.) Similar results are obtained using more sophisticated 1D simulations of radiation transport along spherical shells of mass ejecta (Tanvir et al., 2017; Tanaka et al., 2017).

Numerical relativity (NR) simulations produce dynamical ejecta of a few times $10^{-3}{\rm M_{\odot}}$ with velocities distributed around ${\sim}0.1{-}0.3$c  (Hotokezaka et al., 2013; Bauswein et al., 2013; Radice et al., 2018a). Dynamical ejecta are characterized by a range of electron fractions $0.05\lesssim Y_{e}\lesssim 0.4$; with larger values distributed towards polar regions above the remnant (as part of the shocked component) and lower values across the equatorial plane. These properties are largely independent of the NS equation of state (EOS) (Sekiguchi et al., 2015; Radice et al., 2018a). Additional ejecta from the disk are expected on longer timescales (Perego et al., 2014; Just et al., 2015; Kasen et al., 2015; Metzger and Fernández, 2014; Wu et al., 2016; Siegel and Metzger, 2017; Fujibayashi et al., 2018; Miller et al., 2019); disk mass and composition depend on the binary mass and EOS (Radice et al., 2018b; Perego et al., 2019). Neutrino irradiation can unbind ${\sim}5$% of the disk mass with $Y_{e}>0.25$ and velocities ${\lesssim}0.08$c from the polar region (Perego et al., 2014; Martin et al., 2015). A significant fraction of the disk mass, up to 40%, can be ejected on time scales ${\gtrsim}100$ ms due to magnetic-field induced viscosity and/or nuclear recombination, (Dessart et al., 2009; Fernández et al., 2015; Wu et al., 2016; Lippuner et al., 2017; Siegel and Metzger, 2017; Fujibayashi et al., 2018; Radice et al., 2018c; Fernández et al., 2019; Miller et al., 2019). These secular ejecta are expected to have velocities ${\lesssim}0.05{-}0.1$c and electron fraction in the broad range $0.1\lesssim Y_{e}\lesssim 0.5$, where lower (higher) values are found for black-hole (long-lived NS) remnant. If present, the secular ejecta might give the dominant contribution to the kN on timescales of days to months (Fahlman and Fernández, 2018).

KN light curve models need to account for multiple ejecta (dynamical, wind, viscous, etc.), for the anisotropy of the ejecta composition, and for the irradiation among the ejecta components to fully explain AT2017gfo. Indeed, outflow properties inferred for AT2017gfo using multi-components and 2D kN models including these effects are broadly compatible with the results from simulations, e.g. (Perego et al., 2017; Kawaguchi et al., 2018). The early blue kN however, remains a challenging aspect to model. Both semi-analytical and radiation transport models require ejecta properties different from those found in simulations. In particular, simulations cannot produce ejecta with the large velocities and electron fraction inferred from the electromagnetic data (Fahlman and Fernández, 2018).

There exist alternative explanations of the blue kN based on the interaction between a relativistic jet and the ejecta (Lazzati et al., 2017; Bromberg et al., 2018; Piro and Kollmeier, 2017) but simulations show that successful jets do not deposit a sufficient amount of thermal energy in the ejecta for this mechanism to work (Duffell et al., 2018). Other possibilities include the presence of highly magnetized winds (Metzger et al., 2018; Fernández et al., 2019), or the presence of the so-called viscous-dynamical ejecta (Radice et al., 2018d). However, both models rely on the development of large-scale strong magnetic fields. Here, we identify a new generic hydrodynamics-driven mechanism that works in self-consistent ab-initio simulations and does not require the presence of a strong ordered magnetic field.

2. Method

We perform 3+1 NR simulations of two binaries with mass $M=(1.364+1.364){\rm M_{\odot}}$ and NS described by the microphysical EOS HS(DD2) (Typel et al., 2010; Hempel and Schaffner-Bielich, 2010) and LS220 (Lattimer and Swesty, 1991). The simulations include the merger and the remnant evolution for a timescale of at least 30 ms and up to 100 ms depending on the binary. The results presented here are representative cases producing a long-lived NS remnant (DD2) and a short-lived NS (LS220) from a larger set of simulations that will be presented elsewhere.

We use the WhiskyTHC code (Radice and Rezzolla, 2012; Radice et al., 2014a, b, 2018c) with the approximate neutrino transport scheme developed in (Radice et al., 2016a, 2018a). The simulations treat turbulent viscosity using the general-relativistic large eddy simulations method (GRLES) (Radice, 2017). The interactions between the fluid and neutrinos are treated with a leakage scheme in the optically thick regions (Ruffert et al., 1996; Neilsen et al., 2014) while free-streaming neutrinos are evolved according to the M0 scheme discussed in Ref. (Radice et al., 2018a). The turbulent viscosity in the GRLES is parametrized as $\sigma_{T}=\ell_{\rm mix}c_{s}$, where $c_{s}$ is the sound speed and $\ell_{\rm mix}$ is a free parameter that depends on the intensity of the turbulence. We perform two groups of simulations in this work with $\sigma_{T}$ either set to zero, or prescribed as a function of the rest-mass density as in (Perego et al., 2019) using the results of (Kiuchi et al., 2018). We perform simulations with the same grid setup as in Ref. (Radice et al., 2018a). In particular, the adaptive mesh refinement grids have seven 2:1 refinement levels with finest linear resolutions of $h=246,185,123$m, which are labelled LR, SR and HR. Each model was evolved at least at two different resolutions (LR and SR).

The ejecta are calculated on coordinate spheres at $r=294$ km employing the geodesic criterion for the dynamical ejecta (Radice et al., 2018a). For the wind we use the Bernoulli criterion, which is appropriate for steady-state flow, assuming $(\partial_{t})^{a}$ is an approximate Killing vector (see e.g. (Kastaun and Galeazzi, 2015)). The Bernoulli calculation is started after the ejecta mass computed with the geodesic criterion has saturated to its final value. From the fluid’s stress energy tensor, we compute the angular momentum density flux $J_{r}=T_{ra}(\partial_{\phi})^{a}$, where $\phi$ is the cylindrical angular coordinate; angular momentum is conserved if $(\partial_{\phi})^{a}$ is a Killing vector. $r$-process nucleosynthesis yields are computed using the method detailed in (Radice et al., 2018a).

3. Results

The key dynamical feature of relevance here is the development of spiral arms in the remnant (Shibata and Uryu, 2000; Shibata and Taniguchi, 2006; Bernuzzi et al., 2014; Kastaun and Galeazzi, 2015; Bernuzzi et al., 2016; East et al., 2016; Paschalidis et al., 2015; Radice et al., 2016b; Lehner et al., 2016). The hydrodynamic instability is monitored by a decomposition in Fourier modes $e^{-{\rm i}m\phi}$ of the Eulerian rest-mass density on the equatorial plane [see Eq. (1) of (Radice et al., 2016b)] and characterized by the development of a $m=2$ followed by a $m=1$ mode (East et al., 2016; Paschalidis et al., 2015; Radice et al., 2016b; Lehner et al., 2016; Bernuzzi et al., 2014; Kastaun and Galeazzi, 2015). In the short-lived remnant (LS220) the $m=1$ mode is subdominant with respect to the $m=2$, and it reaches a maximum close to the collapse (Bernuzzi et al., 2014). Instead, in the long-lived remnant (DD2) the $m=1$ becomes the dominant mode at $\sim$20 ms and persists throughout the remnant’s lifetime, while the $m=2$ efficiently dissipates via gravitational-wave emission (Bernuzzi et al., 2016; Radice et al., 2016b). Considering the turbulent viscosity effect, we find that the $m=2$ mode is suppressed more rapidly in presence of viscosity than without viscosity. By contrast, the $m=1$ modes are not significantly affected by viscosity. The spiral arms propagate from the remnant NS into the disk and transport angular momentum outwards as shown in Fig. 1. Such global density waves are a generic and efficient mechanism to redistribute energy and eventually deplete accretion disks (Goodman and Rafikov, 2001; Rafikov, 2016; Arzamasskiy and Rafikov, 2018). Crucially, we find that both the $m=1$ and $m=2$ modes generate a spiral-wave wind from the disk’s outer layers that is distinct from the dynamical ejecta, see Fig. 2.

The long lived NS remnant (DD2) develops a spiral-wave wind more massive than the dynamical ejecta, as shown also in Fig. 2. The spiral-wave wind mass is larger the longer the remnant survives and the more massive the disks are. It continues as long as as the remnant does not collapse and the spiral modes persist. Thus, binary mass asymetry can enhance the spiral-wave wind as we find in simulations discussed elsewhere [In Prep.]. The inclusion of turbulent viscosity alters all the ejecta masses with an additional component and, for the viscosity parametrization we have considered, it enhances the DD2 spiral-wave wind mass by ${\sim}25$%. The viscosity effect is larger than resolution effects. Comparing data at different grid resolutions we find that the largest variation is in the wind mass. The relative variation of mass from data pairs at increasing resolutions is ${\sim}+15\%$ (LR-SR) and ${\sim}+8\%$ (SR-HR). Hence, finite grid effects tend to increase mass. A similar analysis on the average electronfraction and velocity indicate variations below $4\%$.

The spiral-wave wind has an angular distribution of mass similar to the dynamical ejecta with material mostly confined to the orbital plane, as shown by the histograms in Fig. 2. On the contrary, the velocity profiles show a drastic difference between the two ejecta components. While the dynamical ejecta has a broad velocity distribution (Hotokezaka et al., 2013; Bauswein et al., 2013; Radice et al., 2018a), the spiral-wave wind velocity is narrowly distributed around $0.2$c in the case of a long-lived remnant (DD2). The spiral-wave wind from the short-lived remnant (LS220) has a broader velocity distribution extending down to $0.1$c. This is due to the spiral-wave shutting down and the disk transition to a more steady accretion. As a consequence, the spiral-wave wind ceases but ejecta continue as a slower disc wind driven by nuclear recombination solely. The electron fraction of the spiral-wave wind has a narrower distribution than the dynamical ejecta in both cases. But because disks around NS remnants are less compact, colder, and optically thicker than those around black holes (Perego et al., 2019), the outer layers of the DD2 disk have a lower $Y_{e}$ than the LS220 disk and so does the spiral-wave wind coming from those layers. While the spiral-wave wind is generic in its hydrodynamics origin, the quantification of its properties relies on the accurate microphysics and neutrino treatment in our simulations.

Matter in the spiral-wave wind undergoes $r$-process nucleosynthesis, and produces predominantly elements up to the second peak (mass number $A<130$), see Fig. 3. The combined nucleosynthesis in the dynamical ejecta and the spiral-wave wind reproduces the solar abundances to within the uncertainties due to nuclear physics. The radioactive decay in the spiral-wave wind contributes to a blue day-long kN emission similar to the neutrino wind and viscous ejecta (Perego et al., 2014; Martin et al., 2015; Metzger and Fernández, 2014; Miller et al., 2019). But in comparison to the latter, the spiral-wave wind is distributed closer to the equatorial plane, it is faster and more massive.

We calculate light curves in different photometric bands by postprocessing the simulation data with the anisotropic multi-component model of (Perego et al., 2017). In order to emulate the spiral-wave wind from different BNS, the DD2 spiral-wave wind data are extracted every $10$ ms until the end of the simulation (${\sim}90$ ms) and then linearly extrapolated to $250$ ms. The LS220 simulation has instead a complete ejecta, since both the dynamical and the spiral-wave wind have terminated at the end of our simulation. We stress that we do not include additional ejecta components to the ones extracted from the simulations, although we expect additional material to be unbound due to viscous processes and nuclear recombination on even longer timescales (Radice et al., 2018c).

When comparing our results to the early emission of AT2017gfo in Fig. 4, we find good agreement between the observed luminosities in the high frequency bands and our kN model informed by the DD2 simulation with spiral-wave wind masses ${\sim}0.75{-}1.25\times 10^{-2}M_{\odot}$. By contrast, the LS220 simulation does not produce enough ejecta to explain the observations with this light-curve model. Explaining the low frequency bands with the DD2 data would require a more massive spiral-wave wind with mass ${\gtrsim}2\times 10^{-2}M_{\odot}$, implying a remnant lifetime of ${\gtrsim}200$ ms. However, a more massive spiral-wave wind is incompatible with the early emission for the considered simulations. Late-time luminosities (peaking at $t\approx 3{-}10$ days), could be explained by a combination of spiral-wave wind and viscous ejecta from the disintegration of the disk. These results have uncertainties related to our simplified calculation of the kilonova light curves which is expected to be less accurate at late times when absorption features and deviations from local thermodynamics equilibrium become more relevant, e.g. (Smartt et al., 2017). Indeed, time- and energy-dependent modeling of the photon radiation transport will be needed to model more robustly the kN emission, and quantitatively reproduce the observed spectra (Kasen et al., 2017; Tanaka et al., 2017; Miller et al., 2019; Bulla, 2019). Furthermore, all current kilonova models suffer systematic uncertainties in nuclear (e.g. mass models, fission fragments and $\beta$-decay rates) and atomic (e.g. detailed wavelength dependent opacities for $r$-process element) physics (Eichler et al., 2015; Rosswog et al., 2017; Gaigalas et al., 2019).

4. Conclusion

Standard kN models applied to the early AT2017gfo light curve are in tension with ab-initio simulations conducted so far. While alternative interpretations have been proposed, they are either disfavored by current simulations and observations (e.g. jets) (Bromberg et al., 2018; Duffell et al., 2018), or require the presence of large-scale strong magnetic fields which might not be formed in the postmerger (Metzger et al., 2018; Fernández et al., 2019; Radice et al., 2018d; Ciolfi et al., 2019). We identified a robust dynamical mechanism for mass ejection that explains early-time observations without requiring any fine-tuning. The resulting nucleosynthesis is complete and produces all $r$-process elements in proportions similar to solar system abundances. Methodologically, our work underlines the importance of employing NR-informed ejecta for the fitting of light-curves. Further work in this direction should include better neutrino-radiation transport and magnetohydrodynamic effects (Siegel and Metzger, 2017; Fujibayashi et al., 2018; Radice et al., 2018c, a; Miller et al., 2019).

SB and NO acknowledge support by the EU H2020 under ERC Starting Grant, no. BinGraSp-714626. DR acknowledges support from a Frank and Peggy Taplin Membership at the Institute for Advanced Study and the Max-Planck/Princeton Center (MPPC) for Plasma Physics (NSF PHY-1804048). Computations were performed on the supercomputer SuperMUC at the LRZ Munich (Gauss project pn56zo); on the supercomputers Bridges, Comet, and Stampede (NSF XSEDE allocation TG-PHY160025); on NSF/NCSA Blue Waters (NSF AWD-1811236); on supercomputer Marconi at CINECA (ISCRA-B project number HP10BMHFQQ).