Three-dimensional GRMHD simulations of the remnant accretion disks from neutron star mergers: outflows and r-process nucleosynthesis

TL;DR

This paper presents the first 3D GRMHD simulations of neutron star merger accretion disks, revealing turbulence-driven outflows that produce heavy r-process elements, highlighting their significance in nucleosynthesis.

Contribution

It introduces comprehensive 3D GRMHD simulations with realistic physics, demonstrating the role of turbulence and magnetic fields in disk outflows and r-process element production.

Findings

20% of disk mass ejected in outflows

Outflows reach velocities up to 0.1c

Robust production of second and third r-process peaks

Abstract

The merger of binary neutron stars, or of a neutron star and a stellar-mass black hole, can result in the formation of a massive rotating torus around a spinning black hole. In addition to providing collimating media for gamma-ray burst jets, unbound outflows from these disks are an important source of mass ejection and rapid neutron capture (r-process) nucleosynthesis. We present the first three-dimensional general-relativistic magnetohydrodynamic (GRMHD) simulations of neutrino-cooled accretion disks in neutron star mergers, including a realistic equation of state valid at low densities and temperatures, self-consistent evolution of the electron fraction, and neutrino cooling through an approximate leakage scheme. After initial magnetic field amplification by magnetic winding, we witness the vigorous onset of turbulence driven by the magneto-rotational instability (MRI). The disk quickly reaches a balance between heating from MRI-driven turbulence and neutrino cooling, which regulates the midplane electron fraction to a low equilibrium value $Y_\text{e} \approx 0.1$. Over the 380 ms duration of the simulation, we find that a fraction $\approx 20\%$ of the initial torus mass is unbound in powerful outflows with velocities $v \approx 0.03-0.1\,c$ and electron fractions $Y_\text{e} \approx 0.1-0.25$. Post-processing the outflows through a nuclear reaction network shows the production of a robust second and third peak r-process. Though broadly consistent with the results of previous axisymmetric hydrodynamical simulations, extrapolation of our results to late times suggests that the total ejecta mass from GRMHD disks is significantly higher. Our results provide strong evidence that post-merger disk outflows are an important site for the r-process.

Figures

Tables

Table 1. Table 1: Initial configuration: BH mass and dimensionless spin, torus mass, inner and outer torus radius, radius at maximum density, specific entropy, electron fraction, and maximum magnetic-to-fluid pressure ratio.

$M_{\mathrm{BH}}$	${\chi }_{\mathrm{BH}}$	$M_{\mathrm{t0}}$	$R_{\mathrm{in},0}$	$R_{\mathrm{out},0}$	$R_0$	$s_0$	$Y_{\mathrm{e0}}$	$p_b/p_f$
$\left[M_{\odot }\right]$		$\left[M_{\odot }\right]$	$\left[M_{\mathrm{BH}}\right]$	$\left[M_{\mathrm{BH}}\right]$	$\left[\mathrm{km}\right]$	$\left[k_B/b\right]$
$3.00$	$0.8$	$0.03$	$4$	$24$	$30$	8	$0.1$

Table 2. Table 2: Initial torus mass after early relaxation and characteristics of the disk outflows as measured by tracer particles: mean electron fraction, specific entropy, expansion timescales, all at t = t 5 ​ G ​ K 𝑡 subscript 𝑡 5 G K t=t_{\mathrm{5GK}} and subdivided into equatorial ( 30 ∘ < θ < 150 ∘ superscript 30 𝜃 superscript 150 30^{\circ}<\theta<150^{\circ} ) and polar ( θ ≤ 30 ∘ 𝜃 superscript 30 \theta\leq 30^{\circ} and θ ≥ 150 ∘ 𝜃 superscript 150 \theta\geq 150^{\circ} ) outflow (the polar angle θ 𝜃 \theta being measured at the end of the simulation), as well as total integrated outflow mass (polar/equatorial and total). Corresponding values by Fernández  et al. ( 2015 ) (F15) and Just  et al. ( 2015 ) (J15) are also listed.

simulation	outflow type		equatorial outflow			polar outflow			total outflow
		$M_{t,\mathrm{in}}$	${\overline{Y}}_e$	$\overline{s}$	${\overline{t}}_{\exp }$	${\overline{Y}}_e$	$\overline{s}$	${\overline{t}}_{\exp }$	${\overline{Y}}_e$	$\overline{s}$	${\overline{t}}_{\exp }$	$M_{\mathrm{pol}}$	$M_{\mathrm{out}}$
		$\left[{10}^{-2}{M}_{\odot }\right]$		$\left[k_B/b\right]$	$\left[\mathrm{ms}\right]$		$\left[k_B/b\right]$	$\left[\mathrm{ms}\right]$		$\left[k_B/b\right]$	$\left[\mathrm{ms}\right]$	$\left[M_{\mathrm{eq}}\right]$	$\left[M_{t,\mathrm{in}}\right]$
this work	unbound	2.02	0.18	31	24	0.19	39	18	0.18	32	23	0.22	0.16
this work	total	2.02	0.17	28	26	0.19	43	18	0.17	30	25	0.15	0.23
F15 t-a80	total	3.00	0.22	21	35	0.31	38	9.4	-	-	-	0.01	0.17
J15 M3A8m03a2	total	3.00	-	-	-	-	-	-	0.27	30	-	-	0.23
J15 M3A8m03a5	total	3.00	-	-	-	-	-	-	0.25	33	-	-	0.24

Full text

*Introduction.—*Approximately half of the elements heavier than iron are synthesized by the capture of neutrons onto lighter seed nuclei in a dense neutron-rich environment in which the timescale for neutron capture is shorter than the $\beta-$decay timescale Burbidge et al. (1957); Cameron (1957). This ‘rapid neutron-capture process’ (r-process) occurs along a nuclear path far on the neutron-rich side of the valley of stable isotopes. Despite this realization 70 years ago, the identity of the astrophysical sites giving rise to the r-process remains an enduring mystery Qian and Wasserburg (2007); Arnould et al. (2007a); Thielemann et al. (2011).

Among the promising r-process sites are the mergers of compact binaries consisting of two neutron stars (NS-NS, BNS; Symbalisty and Schramm (1982)) or of a NS and stellar-mass black hole (NS-BH; Lattimer and Schramm (1974)). These violent events produce several sources of neutron-rich ejecta, which contribute to their total r-process yields Rosswog (2015); Fernández and Metzger (2016). Historically, most work has focused on matter ejected during the merger process itself, either by tidal forces or due to shock and compression-induced heating at the interface between merging bodies Rosswog (2005); Oechslin et al. (2007); Bauswein et al. (2013); Hotokezaka et al. (2013); Kyutoku et al. (2015); Kastaun and Galeazzi (2015); Radice et al. (2016). While there is broad agreement that a portion of this “dynamical ejecta” is sufficiently neutron-rich to create the heaviest r-process elements, its quantity is sensitive to the properties of the merging binary and the NS equation of state (EOS).

NS mergers are also accompanied by the formation of a massive accretion disk surrounding the central compact object (e.g., Ruffert et al. (1997); Shibata and Taniguchi (2006)). Soon after forming, the neutrino luminosity of the disk is high Popham et al. (1999), driving a small quantity of mass from the disk surface in a neutrino-driven wind Surman et al. (2008); Metzger et al. (2008); Dessart et al. (2009); Perego et al. (2014); Richers et al. (2015); Martin et al. (2015).

On longer timescales of hundreds of milliseconds, the disk expands radially due to the outwards transport of angular momentum. One-dimensional models of this spreading evolution using an $\alpha$-prescription for the effective turbulent viscosity Metzger et al. (2008, 2009) showed that, as the disk accretion rate drops, the midplane transitions from a neutrino-cooled state to a radiatively-inefficient one Lee et al. (2009); Beloborodov (2008). Powerful outflows were predicted following this transition, once heating from turbulent dissipation and nuclear recombination (chiefly alpha-particle formation) are no longer balanced by neutrino cooling.

These initial models were followed by two-dimensional hydrodynamical simulations of the disk evolution in a pseudo-Newtonian gravitational potential, which also adopted an $\alpha$-viscosity prescription. Fernández and Metzger (2013) and Fernández et al. (2015) employed an approximate leakage scheme to account for neutrino cooling, and a ‘light bulb’ irradiation model for the neutrino heating, while Just et al. (2015) used an energy-dependent two-moment closure scheme for the transport of electron neutrinos and antineutrinos. These works found unbound outflows with electron fractions in the range $Y_{\text{e}}\sim 0.2-0.4$ Fernández and Metzger (2013); Just et al. (2015), sufficient to produce the entire mass range of r-process elements Just et al. (2015); Wu et al. (2016); Lippuner et al. (2017). The total fraction of the original disk mass unbound in outflows ranged from $\sim\!5\%$ for a non-spinning BH, to $\sim\!30\%$ for high BH spin $\chi_{\rm BH}\simeq 0.95$ Just et al. (2015); Fernández et al. (2015).

Previous simulations of the remnant disk employ a parameterized hydrodynamical viscosity in place of a self-consistent physical mechanism for angular momentum transport as mediated by the magneto-rotational instability (MRI) Balbus and Hawley (1998). Shibata et al. (2007) performed two-dimensional general-relativistic magnetohydrodynamic (GRMHD) simulations of the disk evolution lasting 60 ms; however, they were not focused on nucleosynthesis and their restriction to 2D precluded a study of the saturated MRI due to the anti-dynamo theorem. In this Letter, we present the first fully three-dimensional GRMHD simulations of the remnant accretion disk evolution and its outflows over a timescale of $\approx\!400$ ms following the merger.

*Numerical setup and initial conditions.—*Simulations are performed in ideal GRMHD with a fixed background spacetime using the open-source EinsteinToolkit111http://einsteintoolkit.org Löffler et al. (2012) with the GRMHD code GRHydro Mösta et al. (2014). GRMHD is implemented using a finite-volume scheme with piecewise parabolic reconstruction Colella and Woodward (1984), the HLLE Riemann solver Harten (1983); Einfeldt (1988), and constrained transport Tóth (2000) for maintaining the magnetic field divergence-free. We have implemented a new framework for the recovery of primitive variables in GRHydro that provides support for any 3-parameter EOS, as well as a recovery scheme based on three-dimensional root finding according to Cerdá-Durán et al. (2008), which shows better and faster convergence properties than two-dimensional schemes; its ability to recover strongly magnetized regions is important for evolving low-density, magnetized disk winds.

Thermodynamic properties of matter are described by the Helmholtz EOS Timmes and Arnett (1999); Timmes and Swesty (2000), which includes contributions to the Helmholtz free energy from nuclei (treated as ideal gas) with Coulomb corrections, electrons and positrons with an arbitrary degree of relativity and degeneracy, and photons in local thermodynamic equilibrium. We consider free neutrons, protons, and alpha particles, whose abundances are calculated assuming nuclear statistical equilibrium (NSE). We add dissociation energies to the Helmholtz EOS as in Fernández and Metzger (2013) to account for the energy release from alpha-particle formation, as well as the additional terms to the thermodynamic derivatives arising from compositional changes.

Neutrino cooling is described by a leakage scheme newly implemented into GRHydro. Leakage schemes are widely used in both core-collapse supernovae and compact-binary merger simulations (e.g., van Riper and Lattimer (1981); Ruffert et al. (1996); Rosswog and Liebendörfer (2003); Sekiguchi et al. (2011); Ott et al. (2013); Perego et al. (2016)). Our implementation follows closely Radice et al. (2016), which is based on Galeazzi et al. (2013) and employs the formalism by Ruffert et al. (1996). We calculate optical depths following the procedure by Neilsen et al. (2014), which is well suited for the aspherical geometry of an accretion disk. We neglect neutrino absorption, which is expected to appreciably change the outflow dynamics only for significantly more massive accretion disks (Fernández and Metzger (2013); see also Fig. 1).

Initial data consists of an equilibrium torus of constant specific angular momentum and specific entropy around a Kerr BH Stergioulas (2011); Friedman and Stergioulas (2013) (Tab. 1). We compute a torus solution in horizon-penetrating Kerr-Schild coordinates, which we use in our simulation. The BH mass and spin correspond to a typical NS merger remnant. BHs promptly formed in BNS mergers show spins $\chi_{\mathrm{BH}}\approx 0.8$ Kiuchi et al. (2009); Rezzolla et al. (2010); Bernuzzi et al. (2014), and are unlikely to be significantly larger Kastaun et al. (2013), whereas BHs formed by delayed collapse typically show spins $\chi_{\mathrm{BH}}\lesssim 0.7$ Sekiguchi et al. (2016). Furthermore, $\chi_{\mathrm{BH}}\sim 0.8$ also represents a typical BH spin for BH–NS mergers required to tidally disrupt the NS and form a sufficiently massive torus Foucart (2012). The initial torus mass also corresponds to typical NS merger scenarios (e.g., Hotokezaka et al. (2011); Foucart et al. (2016)). We set up a weak initial magnetic seed field inside the torus with vector potential components $A^{r}=A^{\theta}=0$ and $A^{\phi}=A_{b}\,\mathrm{max}\{p-p_{\mathrm{cut}},0\}$. Here, $p$ denotes the fluid pressure, $p_{\mathrm{cut}}=1.3\times 10^{-2}p_{\mathrm{max}}$, where $p_{\mathrm{max}}$ is the pressure at maximum density in the torus, and $A_{b}$ sets the initial field strength, which we adjust such that the maximum magnetic-to-fluid pressure ratio in the torus is $<5\times 10^{-3}$; this results in a maximum magnetic field strength of $\approx\!3.3\times 10^{14}\,\mathrm{G}$. The torus is initially embedded in a uniform, tenuous atmosphere with $\rho\approx 37\,\text{g}\,\text{cm}^{-3}$, $T=10^{5}\,\text{K}$, and $Y_{\mathrm{e}}=1$. This density and temperature are sufficiently low to neither impact the dynamics nor the composition of the disk outflows.

Simulations are performed in full 3D without employing symmetries. The grid setup consists of a Cartesian grid hierarchy of 8 refinement levels, extending from the center of the BH to $1.53\times 10^{9}\,\mathrm{cm}$ in every coordinate direction. The finest refinement level corresponds to a spatial domain with a resolution of $856\,\mathrm{m}$ and a diameter of $240\,\mathrm{km}$, which entirely contains the initial accretion torus.

*Results.—*The initial torus is evolved from $t=0$ to $t=381\,\text{ms}$, after which an appreciable fraction of the initial torus mass has been unbound in powerful outflows. After an initial transient phase of about $20\,\text{ms}$ due to the onset of turbulence created by the MRI, the disk settles into a quasi-stationary state for the rest of the simulation. During this early relaxation, $\approx\!33\%$ of the initial torus mass is either accreted onto the BH or ejected via outflows, leaving an effective initial torus of $\approx\!0.02\,M_{\odot}$ (Tab. 2). We exclude matter ejected or accreted during this phase from all further analysis.

Figure 1 shows snapshots of several quantities at the beginning of the quasi-stationary evolution phase. Until the end of the simulation, the disk and outflows remain qualitatively similar as depicted here. In particular, the disk remains optically thin with respect to neutrinos, which have typical energies of a few MeV (Fig. 1, left, upper panel). Neutrino cooling mainly acts in regions closely to the disk midplane, as neutrino emission tracks density. Matter in the disk is heated as it gradually falls into the BH potential due to angular momentum transport via MHD turbulence mediated by the MRI. We check that the MRI is well resolved by monitoring the wavelength of the fastest-growing MRI mode, $\lambda_{\text{MRI}}$, which is typically resolved by 10 or more grid points (cf. Fig. 1, left); $\lambda_{\mathrm{MRI}}$ is estimated by $\lambda_{\mathrm{MRI}}=(2\pi/\Omega)(b/\sqrt{4\pi\rho h+b^{2}})$ Kiuchi et al. (2015), where $\Omega$ is the angular frequency, $\rho$ the rest-mass density, $h$ the specific enthalpy, and $b=\sqrt{b^{\mu}b_{\mu}}$ the comoving magnetic field strength. Very close to the BH resolving the MRI becomes challenging with current computational resources and $\lambda_{\mathrm{MRI}}$ is not resolved by $>10$ grid points at all times and spatial points. At the beginning of the simulation, after initial amplification by magnetic winding, the onset of the MRI further amplifies the weak initial seed magnetic field in the disk over a few rotational periods (resulting in a total amplification of roughly two orders of magnitude for the maximum field strength), before the disk settles into a saturated MRI state. Triggering the MRI both in the poloidal and toroidal components entirely without magnetic winding (for the same initial seed field strength) would require higher resolution and would thus be challenging with current computational resources; this simulation only represents a first attempt in this direction. We note that the resulting typical magnetic field strengths of up to $\sim\!10^{15}$ G close to the BH and the midplane, and typical magnetic-to-fluid pressure ratios of $\sim\!10^{-3}-10^{-1}$, are similar to values found in early BNS post-merger accretion systems Ciolfi et al. (2017); Kiuchi et al. (2015). Figure 3 shows the evolution of the (density-averaged) ratio of electromagnetic energy to internal energy of the disk, $\langle e_{\mathrm{EM}}/e_{\mathrm{int}}\rangle_{D}=\langle n_{\mu}n_{\nu}T^{\mu\nu}_{\mathrm{EM}}/\epsilon\rho W\rangle_{D}$, where $n^{\mu}$ denotes the unit normal to the spatial hypersurfaces of the spacetime foliation, $T^{\mu\nu}_{\mathrm{EM}}$ the stress-energy tensor of the electromagnetic field, $\epsilon$ the specific internal energy, $W$ the Lorentz factor, and $\langle\cdot\rangle_{D}\equiv\int\cdot D\mathrm{d}^{3}x/\int D\mathrm{d}^{3}x$, with $D=\sqrt{\gamma}\rho W$ the conserved rest-mass density and $\gamma$ the determinant of the spatial metric. This ratio stays roughly constant for $t>20\,\mathrm{ms}$ (in a time-averaged sense) and thus indicates that a steady turbulent state of the disk is indeed achieved.

Optically thin neutrino cooling in the midplane is balanced by MHD-driven heating, and the disk regulates itself to a mildly degenerate state with low $Y_{\text{e}}$ Beloborodov (2003). The latter results from a negative feedback process: higher electron degeneracy $\mu_{\text{e}}/k_{\text{B}}T$ results in less electrons (lower $Y_{\text{e}}$) and positrons, causing less neutrino emission, i.e., a lower cooling rate, therefore higher temperatures, and thus lower degeneracy; the resulting state is independent of the initial conditions. Figure 2 shows the disk once it has fully self-regulated itself into this mildly degenerate state ($\mu_{\text{e}}/k_{\text{B}}T\sim 1$). The inner disk remains neutron rich ($Y_{\text{e}}\approx 0.1$) over the course of the simulation up to radii $r\lesssim\!60\,\text{km}$ ($\lesssim\!14$ gravitational radii), consistent with previous one-dimensional models of neutrino-cooled disks Chen and Beloborodov (2007); Metzger et al. (2009).

Above the disk midplane powerful thermal outflows are generated. These are the result of a heating-cooling imbalance: in regions of lower density, viscous heating from MHD driven turbulence and energy release from recombination of free nucleons into alpha particles exceeds cooling by neutrino emission, and the weak interactions essentially ‘freeze-out’ (although further mixing can still change $Y_{\text{e}}$). In the polar funnel these outflows possess high-$Y_{\text{e}}$ ($>0.2$) and high specific-entropy ($s\gtrsim\!100\,k_{\text{B}}/\text{b}$), while the denser equatorially-directed outflows have lower specific entropy ($\sim\!10\,k_{\text{B}}/\text{b}$) and lower $Y_{\text{e}}$.

Thermodynamic properties of the outflow are recorded by $10^{4}$ passive tracer particles that are advected with the fluid. We place these tracer particles of equal mass in the initial setup with a probability proportional to the conserved rest-mass density $D=\sqrt{\gamma}\rho W$. Tab. 2 and Fig. 4 characterize the outflow properties relevant to the r-process, including $Y_{\text{e}}$, $s$, and the expansion timescale $t_{\text{exp}}=r/v$, where $v$ denotes the three-velocity (e.g., Lippuner and Roberts (2015)). These quantities are evaluated for each tracer particle at the last time $t=t_{5\text{GK}}$ when the temperature of the particle drops below 5 GK. At 5 GK, NSE breaks down and full nuclear reaction network calculations are required to track nuclear abundances. We distinguish between total outflow, defined as all tracer particles that have reached $r\geq 10^{3}\,\text{km}$ by the end of the simulation, and unbound outflow, defined as those that are additionally unbound according to the Bernoulli criterion $-hu_{t}>1$, where $u_{t}$ is the time-component of the four-velocity.

By the end of the simulation, $\approx\!(16-23)\%$ of the initial disk mass has been ejected into unbound outflows with $v\approx(0.03-0.1)c$. With the disk still launching outflows by the end of the simulation, our GRMHD setup potentially unbinds significantly more mass compared to two-dimensional, non-MHD, Newtonian simulations with similar disk parameters (Tab. 2; Fernández et al. (2015); Just et al. (2015)). Polar outflows show higher $Y_{\text{e}}$ and specific entropy, and smaller $t_{\text{exp}}$ than equatorial outflows, consistent with Fernández et al. (2015), while we find a factor $\sim\!20$ higher mass in polar outflows. Our $Y_{\text{e}}$ distribution shows a smaller mean and does not extend as high as in Fernández et al. (2015); Just et al. (2015). This may indicate that neutrino absorption (not included here) plays a dominant role in setting the high-$Y_{\text{e}}$ tail of the distribution. Indeed, a preliminary re-analysis including effects of neutrino absorption as in Roberts et al. (2017) shows the ejecta achieves a high-$Y_{\text{e}}$ tail extending up to $\gtrsim\!0.3$; however, our finding of a sizable quantity of low-$Y_{\text{e}}$ ejecta, capable of a full (2nd and 3rd peak) r-process, remains robust. Alternatively, previously employed pseudo-Newtonian potentials and $\alpha$-disks may not accurately capture the heating/cooling interplay which also controls the evolution of $Y_{\text{e}}$.

Full nuclear reaction network calculations with SkyNet Lippuner and Roberts (2015) were performed on the tracer particles in a post-processing step, starting at $t=t_{\text{10GK}}$. Figure 5 shows that the solar abundances Arnould et al. (2007b) are well reproduced throughout the mass number ($A$) range from the 2nd r-process peak ($A\sim 130$) to the rare-earth peak ($A\sim 165$) to the 3rd r-process peak ($A\sim 195$). There is also excellent agreement with observed abundances in metal-poor stars Sneden et al. (2003); Roederer and Lawler (2012); Roederer et al. (2012). We find an overproduction at $A=132$ as observed in Wu et al. (2016); Lippuner et al. (2017). Below the 2nd r-process peak, we recover the trends of the observed solar abundance pattern, but overall underproduce these nuclei, which is consistent with the absence of a significant high-$Y_{\text{e}}$ tail $Y_{\text{e}}>\!0.25$ (Fig. 4); however, a preliminary re-analysis including effects of neutrino absorption as in Roberts et al. (2017) indicates that the entire range of r-process nuclides can be obtained.

*Conclusion.—*We have shown that neutrino-cooled accretion disks in 3D GRMHD quickly self-regulate themselves into a state of moderate electron degeneracy (low $Y_{\text{e}}$) where heating from MRI-driven turbulence is balanced by neutrino cooling. The outflows launched self-consistently as a result of this state tend to unbind more mass with a lower average $Y_{\text{e}}$ than previous axisymmetric Newtonian simulations employing an $\alpha$-viscosity. The nucleosynthesis yields show that these outflows are sufficiently neutron rich to trigger a strong r-process and are well able to reproduce observed solar abundances and observed r-process abundances in metal poor stars from the 2nd to the 3rd r-process peak. Significant contributions to abundances below the 2nd r-process peak can also come from BNS mergers leading to an accretion disk around a metastable hot neutron star, which, due to its strong neutrino emission, may raise $Y_{\text{e}}$ in part of the outflow material Metzger and Fernández (2014); Wu et al. (2016); Lippuner et al. (2017).

Note added.—Following the submission of this paper, a BNS merger was detected by Advanced LIGO and Virgo The LIGO Scientific Collaboration and The Virgo Collaboration (2017). The properties of the infrared kilonova emission observed from this event (inferred total ejecta mass $\approx 0.05M_{\odot}$ and mean velocity $v\approx 0.1\,c$; e.g. Cowperthwaite et al. (2017)) are consistent with the lanthanide-rich matter predicted in this work from disk outflows from a torus of initial mass $\approx 0.1M_{\odot}$.