Full Transport Model of GW170817-Like Disk Produces a Blue Kilonova

TL;DR

This paper presents the first 3D general relativistic neutrino transport simulations of a neutron star merger's accretion disk, demonstrating how neutrino interactions increase electron fraction and produce a blue kilonova.

Contribution

It introduces a comprehensive 3D GRRMHD simulation of GW170817-like mergers, revealing the impact of neutrino transport on ejecta composition and kilonova color.

Findings

Neutrino transport couples the disk's electron fraction globally.

Electron fraction reaches up to Y_e ~ 0.4 in polar outflows.

Nucleosynthesis predicts a blue kilonova from the disk wind.

Abstract

The 2017 detection of the inspiral and merger of two neutron stars in gravitational waves and gamma rays was accompanied by a quickly-reddening transient. Such a transient was predicted to occur following a rapid neutron capture (r-process) nucleosynthesis event, which synthesizes neutron-rich, radioactive nuclei and can take place in both dynamical ejecta and in the wind driven off the accretion torus formed after a neutron star merger. We present the first three-dimensional general relativistic, full transport neutrino radiation magnetohydrodynamics (GRRMHD) simulations of the black hole-accretion disk-wind system produced by the GW170817 merger. We show that the small but non-negligible optical depths lead to neutrino transport globally coupling the disk electron fraction, which we capture by solving the transport equation with a Monte Carlo method. The resulting absorption drives up the electron fraction in a structured, continuous outflow, with electron fraction as high as $Y_e\sim 0.4$ in the extreme polar region. We show via nuclear reaction network and radiative transfer calculations that nucleosynthesis in the disk wind will produce a blue kilonova.

Equations

$$Q_{\text{mri}}^{(\theta )}=\frac{2\pi b^{(\theta )}}{\Delta x^{(\theta )}\sqrt{w+b^2}\Omega },$$

$$Q_{\text{mri}}^{(\text{c})}=\frac{b}{b^{(\theta )}}{Q}_{\text{mri}}^{(\theta )}$$

$$Q_{\text{rad}}=\underset{\Omega }{\min }\left(\frac{\partial N}{\partial t}\frac{u}{J}\right),$$

$$\rho \ge \text{max}\left(\frac{{\rho }_0}{\lfloor \rho \rfloor }\frac{c^4}{G^2{M}_{BH}^2}\frac{1}{r^2},{\rho }_{\text{min}}\right),$$

Full text

I Introduction

In August, 2017, the inspiral and merger of a pair of neutron stars (GW170817) was jointly detected by gravitational wave detectors and electromagnetic telescopes around the world Abbott et al. (2017a). This detection confirms that such mergers are central engines of short gamma ray bursts Soares-Santos et al. (2017); Eichler et al. (1989); Narayan et al. (1992) and a site of r-process nucleosynthesis Abbott et al. (2017b); Rosswog et al. (2018), where the heaviest elements in our universe are formed Blinnikov et al. (1984); Lattimer and Schramm (1976); Lattimer et al. (1977); Côté et al. (2018).

The radioactive decay of r-process elements produces an optical and infra-red afterglow—the kilonova Metzger et al. (2010); Côté et al. (2018), which was observed clearly in the aftermath of GW170817 Abbott et al. (2017b). This afterglow is likely driven by at least two components Cowperthwaite et al. (2017); Tanvir et al. (2017); Tanaka et al. (2017): a “blue” kilonova driven by polar outflow Evans et al. (2017); Nicholl et al. (2017) and a “red” kilonova driven by equatorial outflow Chornock et al. (2017); Tanvir et al. (2017); Troja et al. (2017). These distinct components are believed to arise due to the different compositions of these outflows Cowperthwaite et al. (2017); Tanvir et al. (2017); Tanaka et al. (2017). Relatively neutron rich outflows with an electron fraction $Y_{e}\lesssim 0.25$ can produce lanthanides Mumpower et al. (2012); Lippuner and Roberts (2015), which are opaque to blue light Kasen et al. (2013); Fontes et al. (2015, 2019). Less neutron-rich outflows ($Y_{e}\gtrsim 0.25$) will produce nucleosynthetic yields which allow blue light to escape the photosphere Lippuner and Roberts (2015); Martin et al. (2015).

Several mechanisms can produce these outflows Rosswog (2015); Fernández and Metzger (2016). Tidal ejecta typically produce a red component, while shock-driven, near-polar dynamical ejecta can potentially be blue Sekiguchi et al. (2015); Foucart et al. (2016); Radice et al. (2016); Bovard et al. (2017); Martin et al. (2018). Wind off of a remnant hypermassive, supramassive, or stable neutron star can also be blue Dessart et al. (2008); Perego et al. (2014); Martin et al. (2015). Finally, a remnant-disk system can drive a wind Sekiguchi et al. (2015); Perego et al. (2014); Ruffert et al. (1997); Popham et al. (1999); Shibata and Taniguchi (2006); Shibata et al. (2007); Surman et al. (2008); Metzger et al. (2008); Beloborodov (2008); Metzger et al. (2009); Lee et al. (2009); Fernández and Metzger (2013); Fernández et al. (2014); Janiuk (2014); Just et al. (2015); Richers et al. (2015); Foucart et al. (2015); Wu et al. (2016); Lippuner et al. (2017); Hossein Nouri et al. (2018); Siegel and Metzger (2018); Fernández et al. (2018); Foucart et al. (2018). For this last source, the composition is as-yet uncertain. Some studies show the disk wind to have an electron fraction ranging from $Y_{e}{\mathord{\sim}}0.2-0.4$ and thus produce a blue component Fernández and Metzger (2013); Fernández et al. (2014); Janiuk (2014); Just et al. (2015); Fernández et al. (2018); Foucart et al. (2018). Other work shows the disk wind to be uniformly composed of $Y_{e}{\mathord{\sim}}0.2$ material that produces only a red component Wu et al. (2016); Siegel and Metzger (2018).

We focus on the evolution of the post-merger disk. Until now, studies of the remnant disk wind have employed various approximations to the neutrino transport, neutrino-matter coupling, or magnetohydrodynamics (MHD). In this work, we present, for the first time, fully three-dimensional general-relativistic radiation magnetohydrodynamics (GRRMHD) simulations of a post-merger disk system with full neutrino transport using a Monte Carlo method.

We model a black hole accretion disk system which may have formed from the GW170817 merger Shibata et al. (2017). Magnetohydrodynamic turbulence Balbus and Hawley (1991) drives a wind Blandford and Payne (1982) off the disk. We find the electron fraction of this outflow ranges from $Y_{e}{\mathord{\sim}}0.2$ to $Y_{e}{\mathord{\sim}}0.4$. Moreover, we find that the composition of the outflow varies significantly with angle off of the midplane, suggesting that the observed character of the outflow depends heavily on viewing angle. Thus, a blue, wind-produced kilonova will be visible if the remnant is viewed close to the polar axis.

II Methods

We perform a GRRMHD simulation in full three dimensions with our code, $\nu\texttt{bhlight}$Miller et al. (2019). We assume a Kerr background metric, consistent with the relatively small disk mass compared to black hole mass. The radiation transport is treated via explicit Monte Carlo and the MHD is treated via finite volumes with constrained transport. The two methods are coupled via operator splitting.

We use the SFHo equation of state Steiner et al. (2013) as tabulated in O’Connor and Ott (2010, 2010–) and the neutrino-matter interactions described in Miller et al. (2019) and tabulated in Burrows et al. (2006). For initial data, we use parameters consistent with a remnant from GW170817 Shibata et al. (2017); Abbott et al. (2017a, 2019): an equilibrium torus Fishbone and Moncrief (1976) of mass $M_{d}=0.12$ $M_{\odot}$ and constant electron fraction $Y_{e}=0.1$ around a black hole of mass $M_{BH}=2.58$ $M_{\odot}$ and dimensionless spin $a=0.69$. We thread our torus with a single poloidal magnetic field loop such that the minimum ratio of gas to magnetic pressure is 100.

III Outflow Properties

Our disk drives a wind consistent with other GRMHD simulations of post-merger disks Fernández and Metzger (2013); Fernández et al. (2014); Janiuk (2014); Just et al. (2015); Fernández et al. (2018); Foucart et al. (2018); Wu et al. (2016); Siegel and Metzger (2018), which expands outward from the disk in polar lobes as shown in figure 1. We record material crossing a sphere of radius $r\sim 10^{3}$ km. Figure 2 bins outflow material in both electron fraction $Y_{e}$ and in angle off the equator, $|90^{\circ}-\theta_{\text{bl}}|$ for Boyer-Lindquist angle $\theta_{\text{bl}}$, integrated in time. The 90% confidence interval for the viewing angle for GW170817 Finstad et al. (2018) is bounded by the dashed lines.

We choose two regions, one close to the midplane, and one far from it, highlighted in the red and blue rectangles. We bin the electron fraction in these regions in the red and blue histograms. Regardless of electron fraction, ejected material has an average entropy, $s$, of about 20 $k_{b}/$baryon and an average radial velocity (as measured at a radius of 1000 km) of about 0.1$c$.

The electron fraction depends on angle off of the midplane and this dependence persists through time. The right panel of figure 3 shows the average electron fraction of gravitationally unbound material passing through a surface at $t\sim 10^{3}$ km as a function of angle off the equator and time. For any given time, larger $|90^{\circ}-\theta_{\text{bl}}|$ correlates with larger $Y_{e}$.

We use the nuclear reaction network SkyNet Lippuner and Roberts (2017) to compute nucleosynthetic yields on tracer particles advected with gravitationally unbound material. We start the network calculation when the tracer reaches $T{\mathord{\sim}}10$ GK and we assume a nuclear statistical equilibrium (NSE) composition at that time. The network is run up to $t=10^{9}$ s assuming homologous expansion ($\rho\propto t^{-3}$) and uses the same nuclear physics inputs as in Lippuner et al. (2017); Roberts et al. (2017), namely: 8000 nuclides and 140,000 nuclear reactions, including fission, with rates from Cyburt et al. (2010); Frankel and Metropolis (1947); Panov et al. (2010); Mamdouh et al. (2001); Wahl (2002); Fuller et al. (1982); Oda et al. (1994); Langanke and Martínez-Pinedo (2000).

Figure 4 plots nucleosynthetic yields. We plot three angular cuts: in red for material near the midplane, in blue for material near the poles, and in black for material within the viewing angle for GW170817 Finstad et al. (2018). We sketch out the range of possible yields in gray. The second, rare-earth and third peaks are suppressed by up to a factor of 100 with respect to the first peak in the polar regions.

IV Outflow Mass

The left panel of figure 3 shows the total mass in the outflow as a function of time. Due to computational cost, we did not run our simulation for long enough to observe the total amount of mass that becomes gravitationally unbound. As a lower bound, we report the amount of material with Bernoulli parameter $B_{e}>0$ Novikov and Thorne (1973) at a radius greater than 125 gravitational radii (${\mathord{\sim}}500$ km) at the end of the simulation (${\mathord{\sim}}127$ ms). (This includes material that has already left the domain.) We find this to be about $4.33\times 10^{-3}M_{\odot}$ and the ratio of mass in the outflow to accreted mass is about 9% by this time in the simulation. About 18% of this outflow has an electron fraction of $Y_{e}\geq 0.275$ and about 14.5% is within the expected range of viewing angles for GW170817.

V Neutrino Transport

A characterization of the importance of neutrino absorption is the neutrino absorption optical depth $\tau$ of the disk. $\tau\ll 1$ implies free-streaming and $\tau\gg 1$ implies no neutrino can escape. At relatively early times ($t\lesssim 30$ ms), we find $\tau\sim 10$. In this phase, $Y_{e}$ evolution is dominated by emission of electron neutrinos in the core of the disk and their absorption in the corona. At later times ($t\gtrsim 30$ ms), the disk achieves a quasistationary state with $\tau\sim 0.1$. Although this later stage is emission dominated, reaching it requires properly treating absorption.

Figure 5 shows this transition. We plot for both phases the $\sinh_{10}^{-1}$ of the rate $\lambda=\partial N/\partial t$ of emitted neutrinos minus the rate of absorbed neutrinos for electron neutrinos ($x<0$) and electron antineutrinos ($x>0$), where we define $\sinh^{-1}_{10}$ as the inverse of $\sinh_{10}(x):=(10^{x}-10^{-x})/2$ such that for $|x|\gtrsim 10$, $\sinh^{-1}_{10}(x)\to\text{sign}(x)\log_{10}(|2x|)$. Red and orange imply $Y_{e}$ is decreasing due to neutrino interactions. Blue and purple imply it is increasing. Figure 6 shows the resulting change in the electron fraction in the Lagrangian frame for each phase.

VI Electromagnetic Counterpart

We compute spectra from the kilonova assuming spherically symmetric outflow composed of nucleosynthetic yields produced in material with $|90^{\circ}-\theta_{\text{bl}}|\leq 15^{\circ}$ and $|90^{\circ}-\theta_{\text{bl}}|\geq 50^{\circ}$. For comparison, we compute spectra for an outflow with solar-like abundances such as those reported in Siegel and Metzger (2018). For the former, we assume an outflow mass of $M_{e}=10^{-2}$ $M_{\odot}$, consistent with our results. For the latter, we assume an outflow mass of $M_{e}=2\times 10^{-2}$ $M_{\odot}$, consistent with Siegel and Metzger (2018). We use a mean radial velocity of $0.1c$.

To compute spectra for each model, we simulate radiative transfer with the Monte Carlo code SuperNu Wollaeger and van Rossum (2014); Wollaeger et al. (2019), using opacity from the LANL suite of atomic physics codes (Fontes et al., 2015). We use a complete suite of lanthanide opacities (Fontes et al., 2019), and some representative wind opacities (Wollaeger et al., 2018). These calculations do not explore the effect of aspherical morphology or uncertainties in r-process heating or composition.

Figure 7 shows computed spectra for several epochs after merger. At early times, the polar outflow produces more luminous spectra peaked at a blue wavelength, consistent with a blue kilonova. Differences in these early-time spectra amount to about a 2 magnitude difference in brightness between polar and equatorial outflows. At late times, the more neutron-rich outflows are more luminous and peaked at long wavelengths, consistent with a red kilonova. The luminosity peaks at ${\mathord{\sim}}3\times 10^{41}$ erg$/$s after ${\mathord{\sim}}0.3$ days for the $|90^{\circ}-\theta_{\text{bl}}|\leq 15^{\circ}$ outflow, at $3\times 10^{40}$ erg$/$s after ${\mathord{\sim}}2$ days for the $|90^{\circ}-\theta_{\text{bl}}|\geq 50^{\circ}$ outflow, and $4\times 10^{40}$ erg$/$s after ${\mathord{\sim}}4$ days for the solar-like outflow.

VII Outlook

We have explored a possible disk-driven outflow from the remnant of the GW170817 merger using the first full transport GRRMHD simulations of a post-merger accretion disk system. We calculate nucleosynthetic yields and spectra of the electromagnetic counterpart that would be observed given these yields. These spectra indicate a blue kilonova, as viewed off the midplane, and a red one, as viewed from the midplane We find about 9% of accreted mass ends up in this outflow.

The range of electron fractions in our outflow is consistent with Fernández and Metzger (2013); Fernández et al. (2014); Just et al. (2015) and disagrees with Wu et al. (2016); Siegel and Metzger (2018). The former include neutrino absorption but approximate magnetic fields with a viscous prescription, while Siegel and Metzger (2018) uses GRMHD but includes neutrino absorption only approximately in post-processing. Our full treatment allows us to conclude that neutrino absorption is critical to attaining this range of $Y_{e}$’s.

The electromagnetic counterpart we compute is incomplete, as it is sourced only by the disk outflow. Depending on the equation of state, it is possible that a transient remnant neutron star supported by differential or rigid rotation existed for some time before collapse to a black hole. For some systems, a long-lived neutron star may be the remnant Rosswog (2015); Fernández and Metzger (2016). This remnant can produce its own separate, potentially massive outflow Fryer et al. (1996, 2006); Fryer (2009). This outflow and shock-driven dynamical ejecta Sekiguchi et al. (2015); Foucart et al. (2016); Radice et al. (2016); Bovard et al. (2017); Martin et al. (2018) can also contribute to a blue kilonova Perego et al. (2014); Lippuner et al. (2017); Martin et al. (2015).

For a generic binary neutron star merger, our results imply that a blue kilonova does not necessarily require a short- or long-lived neutron star remnant. However, these additional sources may be required to explain both the total mass and the velocity of the source of the blue component of the afterglow of GW170817 as determined by early light curve models Evans et al. (2017); Cowperthwaite et al. (2017); Nicholl et al. (2017); Tanvir et al. (2017); Tanaka et al. (2017); Troja et al. (2017). Combining our results with other potential sources for a blue component and the red kilonova from tidal ejecta suggests a three (or more) component kilonova model, such as the ones described in Kasliwal et al. (2017); Perego et al. (2017).

In a black hole-neutron star merger, only the tidal ejecta and accretion disk are present. An important observational implication of our model is that this disk-wind system is sufficient to produce a blue kilonova. This is in contrast to Wu et al. (2016); Siegel and Metzger (2018), which would imply that black hole-neutron star mergers only produce a red kilonova.

Another important implication of our model is that accurately capturing the early transient phase of the disk, when optical depths are relatively large, is critical to correctly predicting the long-term outflow. Unfortunately, initial conditions are a source of uncertainty in kilonova disk modeling. A hot hypermassive or supramassive neutron star can emit its own neutrino flux, which can reset the electron fraction of the disk. Even in the absence of a hot remnant, the seed magnetic field is uncertain in both strength and topology. As the community moves forward more attention should be paid to both the initial transient phase of the disk and the initial conditions that drive this early phase.

VIII Acknowledgements

We thank Francois Foucart, Daniel Siegel, Ingo Tews, Patrick Mullen, Roseanne Cheng, and especially Ben Prather for their insight. We also thank our anonymous reviewer for their thoughtful questions and suggestions.

We gratefully acknowledge support from the U.S. Department of Energy (DOE) Office of Science and the Office of Advanced Scientific Computing Research via the SciDAC4 program and Grant DE-SC0018297, from the U.S. NSF grant AST-174267, and from the U.S. DOE through Los Alamos National Laboratory (LANL). This work used resources provided by the LANL Institutional Computing Program. Additional funding was provided by the LDRD Program and the Center for Nonlinear Studies at LANL under project number 20170508DR. LANL is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. DOE (Contract No. 89233218CNA000001). Authorized for unlimited release under LA-UR-19-22623.

Appendix A Resolution and grid

We use a radially logarithmic, quasi-spherical grid in horizon penetrating coordinates, as first introduced in Gammie et al. (2003), with $N_{r}\times N_{\theta}\times N_{\phi}=192\times 128\times 64$ zones out to a radius of ${\mathord{\sim}}4\times 10^{3}$ km. Our grid focuses resolution at the midplane and we are roughly three times more resolved at the midplane than a grid with uniform resolution in $\theta$ for an effective resolution of $N_{\theta}^{\text{eff}}\sim 384$.

In the region where opacities are nonvanishing (roughly 125 km), we use more than 38 million Monte Carlo radiation packets at every time step, resulting in an average packet density of more than 20 packets per finite volume cell. We track Lagrangian fluid elements via more than $10^{6}$ tracer particles, of which about 10% end up gravitationally unbound by the final time of ${\mathord{\sim}}127$ ms. We initialize our tracers so that, at the initial time, they roughly uniformly sample non-atmosphere regions by volume, as described in Miller et al. (2019) and suggested in Bovard and Rezzolla (2017).

The focusing effect of our grid reduces the number of grid points required to resolve the magneto-rotational instability (MRI) Balbus and Hawley (1991). Following the definition in Sano et al. (2004), we define a quality factor

[TABLE]

for the MRI to be the number of grid points per minimum unstable MRI wavelength inside the disk. Here $b^{(\theta)}$ is the $\theta$-component of the magnetic field four-vector, $\Delta x^{(\theta)}$ is grid spacing in the $\theta$ direction, $w$ is the enthalpy of the fluid, and $b^{2}=b^{\mu}b_{\mu}$ is total magnetic field strength. Following Siegel and Metzger (2018), we also define

[TABLE]

for $b=\sqrt{b^{\mu}b_{\mu}}$, which uses the strength of the magnetic field in the comoving frame, rather than the lab $\theta$-direction. Unfortunately, the nuclear equation of state we use Steiner et al. (2013) makes our enthalpy larger than for an equivalent disk with an ideal gas equation of state Porth et al. (2019), and this larger enthalpy makes resolving the MRI more challenging.

We plot both $Q_{\text{mri}}^{(\theta)}$ and $Q_{\text{mri}}^{(c)}$ in the mid-plane averaged over $\phi$ for several times in figure 8. On average, we find $Q^{(\text{c})}_{\text{mri}}\gg 10$ for all time. At early times, $Q^{(\theta)}_{\text{mri}}\gtrsim 10$. As the disk evolves, this quality factor drops on average to a minimum of about $Q^{\theta}_{\text{mri}}\gtrsim 2$ at $t=10^{4}GM_{BH}/c^{3}=127$ ms. We note that high-order spatial reconstructions such as the WENO-5 Tchekhovskoy et al. (2007) reconstruction we use may effectively improve this quality factor Porth et al. (2019). For various technical issues related to resolving the MRI in global simulations, see Porth et al. (2019).

Appendix B Resolution in the radiation sector

We define the Monte Carlo quality factor

[TABLE]

minimized over the simulation domain $\Omega$. $N$ is the number of emitted Monte Carlo packets, $u$ is gas internal energy density by volume, and $J$ is the total frequency and angle integrated neutrino emissivity. $Q_{\text{rad}}$ roughly encodes how well resolved the radiation field is, with $Q_{\text{rad}}=1$ a marginal value. In our simulation, we find $Q_{\text{rad}}\geq 100$.

Appendix C Artificial Atmosphere Treatment

Since our Eulerian code cannot handle vacuum, we demand densities $\rho$ obey

[TABLE]

where $\rho_{0}=10^{-5}$ is a unitless, simulation-dependent parameter, $\lfloor\rho\rfloor=1.1\times 10^{13}$ g$/$cm3 is the code unit for density and $\rho_{\text{min}}=1.6\times 10^{2}$ g$/$cm3 is the minimum density in our tabulated equation of state. We set our initial atmosphere regions to nearly virial temperatures to prevent the atmosphere from falling back onto the disk. We track our artificial atmosphere with a passive scalar and ensure it does not contribute to any reported quantities such as outflow mass and electron fraction.

As the disk-wind system evolves, outflow will displace the artificial atmosphere by pushing it through the outer boundary of the domain. There is no artificial atmosphere remaining after about 63 ms. The total amount of atmosphere displaced in this way before the wind completely displaces it is roughly $10^{-5}M_{\odot}$, or ${\mathord{\sim}}10^{-4}$ the mass of the disk and ${\mathord{\sim}}10^{-2}$ the mass of the outflow. We find that the radial momentum flux in atmosphere regions is always more than three orders of magnitude less than in wind regions, giving us confidence that our artificial atmosphere does not interfere with the dynamics of the outflow.