Jets from neutron-star merger remnants and massive blue kilonovae

Introduction.—The astrophysical origin of about half of the elements heavier than iron, created via rapid neutron capture (the r-process), remains an open question [cowan2021Origin, siegel2022r]. The first observed binary neutron-star merger, detected via gravitational waves [abbott2017GW170817], referred to as GW170817, was followed by quasi-thermal emission [abbott2017Multimessenger, coulter2017Swope, soares-santos2017Electromagnetic] consistent with radioactive heating from the nucleosynthesis of r-process elements in the merger debris—a kilonova [metzger2010electromagnetic]. Based on inferred event rates of BNS mergers, the ejecta yield as well as the photometric and spectroscopic properties [villar2017Combined, kasen2017Origin, smartt2017Kilonova, pian2017Spectroscopic, chornock2017Electromagnetic, tanaka2020Systematic, watson2019identification, kasliwal2022spitzer], this transient not only provided strong evidence for the production of both light (atomic mass number $A\lesssim 135$) and heavy (lanthanide-bearing; $A\gtrsim 136$) r-process elements in this particular event, but also for BNS mergers being a potentially dominant production site of r-process elements in the Universe.

Whereas the origin of the lanthanide-bearing red and infrared emission of the GW170817 kilonova peaking on timescales of a week is naturally explained by MHD-driven outflows from a massive, self-regulated, neutrino-cooled accretion disk around a final remnant black hole [siegel2017ThreeDimensional, kasen2017Origin, fernandez2019Longterm, christie2019role, li2021Neutrino], the origin of the early ($\sim$day) blue and ultraviolet (UV) emission remained more elusive. Both the blue and red GW170817 kilonova emission are difficult to explain by dynamical debris from the collision itself [siegel2019GW170817, metzger2020Kilonovae, radice2020Dynamics] (but see Ref. [kawaguchi2018Radiative] for a corner case); therefore, several alternative mechanisms for the origin of the blue emission post-merger have been considered, including a combination of magnetically driven and neutrino-driven winds from a remnant neutron star [metzger2018Magnetar, ciolfi2020Magnetically, mosta2020Magnetar], additional turbulent viscosity in the neutron-star remnant [fujibayashi2020postmerger, radice2018ViscousDynamical], spiral waves driven into a post-merger accretion disk by non-axisymmetric modes of a remnant neutron star [nedora2019Spiralwave], and outflows from a post-merger accretion disk both around a remnant neutron star [fahlman2018Hypermassive] and a black hole [miller2019Full].

Here, we demonstrate by means of self-consistent ab-initio simulations of the merger and post-merger phases that the inferred properties of the blue GW170817 kilonova emission can arise naturally from mass ejection within only $\lesssim\!50$ ms post-merger due to MHD-driven winds from an accretion disk, aided by non-linear hydrodynamic effects. The emergence of a jet from the remnant neutron star generates a UV/blue precursor signal and, upon collapse of the remnant, ‘seeds’ an ultra-relativistic jet to generate a short gamma-ray burst.

Computational setup.—We solve Einstein’s equations coupled to the general-relativistic magneto-hydrodynamic (GRMHD) equations with weak interactions using an enhanced version [siegel2018Threedimensional, combi2023grmhd] of the flux-conservative code GRHydro [baiotti2003New, mosta2013GRHydro], which is part of the EinsteinToolkit ¹¹1http://einsteintoolkit.org open-source code framework [goodale2003Cactus, schnetter2004Evolutions, thornburg2004Fast, loffler2012Einstein, babiuc2019einstein]. We use the numerical setup presented in detail in Ref. [combi2023grmhd] (henceforth CS23); in particular, we implement the recovery of primitive variables of Ref. [siegel2018Recovery], which provides support for tabulated nuclear equations of state (EOS), weak interactions, and neutrino radiation transport via a one-moment approximation of the general-relativistic Boltzmann equation [radice2016Dynamical, radice2018Binary]. A fixed Cartesian grid hierarchy composed of six nested refinement boxes is used. The finest mesh covers $\simeq\!76$ km in diameter with a resolution of $\Delta x=180$ m. The largest box has an extend of 3000 km. Reflection symmetry across the orbital plane is employed for computational efficiency. A comparison run without imposing the symmetry at slightly lower resolution shows that our conclusions are not affected by this setup (see CS23). The GRMHD scheme uses an atmosphere floor of $\rho_{\rm atm}\sim 5\times 10^{2}\,{\rm g}\,{\rm cm}^{-3}$.

The initial data consists of two cold, $\beta$-equilibrated, equal-mass neutron stars of radius 11.6 km and mass $1.35M_{\odot}$ in quasi-circular orbit at a separation of 45 km. We build initial data with the elliptical solver LORENE [gourgoulhon2001Quasiequilibrium], employing the APR EOS [akmal1998Equation] in finite-temperature, tabulated form [schneider2019AkmalPandharipandeRavenhall]. This EOS generates cold non-rotating neutron stars with a maximum mass of $2.2M_{\odot}$ and a radius of 11.6 km for a $1.4M_{\odot}$ star, well within the ballpark of currently allowed values [landry2020nonparametric, cromartie2020relativistic]. While similar post-merger phenomenology (see below) is observed for the other, stiffer EOS configurations of CS23, we focus here on the APR configuration, which results in the longest-lived remnant neutron star. After setting the hydrodynamical variables, we initialize a weak poloidal magnetic seed field confined to the interior of each star with maximum strength of $B_{\rm max}=3\times 10^{15}$ G at the center and total energy of $10^{48}$ erg. The initial magnetic field is energetically and dynamically insignificant.

Magnetic field evolution and jet formation.—During inspiral, the magnetic seed field remains buried inside the stars. Starting at the merger (referred to as $t=0$ ms), the field is amplified exponentially by the Kelvin-Helmholtz instability (KHI) for $\approx\!2$ ms, as evidenced by the toroidal field in Fig. 1 [kiuchi2014high, kiuchi2015Efficient, kiuchi2018Global, price2006Producing]. Since the KHI has no characteristic spatial scale, the finite resolution of our simulation is only able to capture partial amplification of the average toroidal magnetic field from essentially zero G to $\approx\!5\times 10^{15}$ G. In the first $15$ ms post-merger, the generation of large-scale eddies in the remnant’s interior by quasi-radial core bounces further amplifies the toroidal field by an order of magnitude and roughly $\propto\!t^{2}$. These currents are associated with the formation of vortices around the inner core [kastaun2021numerical, kastaun2015Propertiesa, ciolfi2017General] which dissipate kinetic into magnetic energy.

The toroidal field eventually continues to grow at the expected linear growth rate of magnetic winding due to differential rotation inside the star. Amplification proceeds mainly in the slowly rotating core with positive angular velocity gradient interior to $r\approx 8$ km. At larger radii, a nearly rotationally supported ‘Keplerian envelope’ is established that gradually transitions into an accretion disk formed by merger debris, a well-established quasi-universal configuration [kastaun2015Properties, ciolfi2017General, kastaun2021numerical, fujibayashi2018Mass]. \rvwA larger resolution would lead to earlier saturation of the turbulent magnetic field via the KHI [kiuchi2018Global, aguilera2023role, chabanov2023crustal], only making the subsequent amplification processes obsolete, but unlikely altering the qualitative evolution of the system.

At $\approx\!25$ ms after merger, thanks to magnetic winding providing an inverse turbulent cascade, the small-scale, turbulent field has been wound up into a large-scale toroidal structure (Fig. 1). Owing to their magnetic buoyancy, toroidal fields eventually rise to the stellar surface in the polar regions, break out of the star and form a magnetic tower [Lynden-Bell1996] (Figs. 1, 2).

Outflow properties.—Prior to the emergence of the magnetic tower structure, strong neutrino radiation ($L_{\nu_{\rm e}}\approx 1.7\times 10^{52}\,{\rm erg}\,{\rm s}^{-1}$) from the hot merger remnant drives a neutrino-driven wind [dessart2009Neutrino, perego2014Neutrinodriven, desai2022three] of unbound material in polar regions with typical velocities $v\lesssim 0.1\,c$ and electron fraction $Y_{e}\approx 0.5$ via absorption of neutrinos in a gain layer above the stellar surface that extends out to $\lesssim\!50$ km (Fig. 2). In this gain layer, which is similar to that above proto-NSs in core-collapse supernovae, the net absorbed energy per unit time as seen by the Eulerian observer $\dot{Q}_{\rm net}=\int_{\theta<30^{\circ}}\dot{q}_{\rm net}\rho\Gamma\sqrt{\gamma}\,d^{3}x$ in a polar volume (polar angle $\theta<30^{\circ})$ roughly equals the kinetic power $\dot{E}_{\rm k}=\int_{\theta<30^{\circ}}\sqrt{\gamma}h\rho\Gamma(\alpha v^{i}-\beta^{i})\,dA_{i}$ of the wind leaving that volume along the polar axis. Here, $\dot{q}_{\rm net}$ is the net specific neutrino heating rate, $h$ the specific enthalpy, $\Gamma=-n_{\mu}u^{\mu}$ is the Lorentz factor of the plasma with three-velocity $v^{i}$ and four-velocity $u^{\mu}$, $dA_{i}$ is the surface element, and $\alpha$, $\mathbf{\beta}$, and $\gamma$ denote the lapse, shift, and determinant of the three-metric of the adopted 3+1 foliation of spacetime with normal vector $n^{\nu}$, which defines the Eulerian observer. A steady-state wind profile $\rho\propto r^{-2}$ emerges in polar regions, as expected from mass conservation, $\dot{M}=\Delta\Omega\,r^{2}\rho\Gamma v$, for a wind opening solid angle $\Delta\Omega$, constant mass-loss rate $\dot{M}$, and four-velocity $\Gamma v$ set by neutrino absorption at the base of the wind.

As the buoyant toroidal magnetic field structures break out of the star and neutrino absorption helps to form a magnetic tower along the rotational axis, a strongly magnetically dominated outflow with magnetic-to-fluid pressure ratio of up to $p_{\rm B}/p_{\rm f}\approx 10^{2}$ is established (Fig. 2). A steady-state wind profile $\rho\propto r^{-2}$ quickly emerges (Fig. 2), with a steady-state mass-loss rate $\dot{M}$ enhanced by roughly one order of magnitude—consistent with the wind solutions of Ref. [metzger2018Magnetar]. Within the same time window of 5–10 ms, the kinetic power of the outflow increases by more than an order of magnitude and comes into equipartition with neutrino heating and the (dominant) Poynting luminosity of the magnetic structure, $\dot{E}_{\rm k}\approx\dot{Q}_{\rm net}+L_{\rm EM}\approx L_{\rm EM}$. The magnetization, $\sigma=L_{\rm EM}/\dot{M}$, in the polar region ($\theta\lesssim 30^{\circ}$ or specific entropy $s>25k_{\rm B}$ baryon^-1), rapidly increases to $\sigma\approx 0.1$ (Fig. 4). Through strong neutrino heating at the base, shock-heating at the jet head, and in internal shocks, the plasma outflow attains high specific entropy $\gtrsim\!(25-100)k_{\rm B}$ baryon^-1 (Fig. 3). Magneto-centrifugal forces accelerate the plasma in the polar funnel from average velocities of $v\lesssim 0.1c$ (neutrino-driven wind) to $v\approx(0.25-0.6)c$. For acceleration of the wind primarily by magnetic fields, the flow should acquire a four-velocity $u=v^{r}\Gamma\approx c\sigma^{1/3}$ when it reaches the fast magnetosonic surface [michel1969relativistic]. This limit with the density-averaged value of $\langle u\rangle\approx 0.35c$ is indeed reached in the polar region (Fig. 4).

A magnetized ($\sigma\sim 5-10$) jet structure emerges in the polar funnel consisting of exclusively dynamically unbound material (geodesic criterion; $-u_{0}<1$) with an half-opening angle of $\approx 20^{\circ}$ at the base, and high entropy, $s\gtrsim 50k_{\rm B}$ baryon^-1. The jet head propagates with $v\gtrsim 0.6c$ through and breaks out from the envelope of merger debris (Fig. 3). The jet is stabilized by subdominant toroidal magnetic fields in the core, which reduce instabilities at the jet-envelope boundary layer [gottlieb2020structure] while avoiding global (kink) instabilities that can develop in strongly magnetized jets ($\sigma\gg 1$) [bromberg2016relativistic]. The terminal velocity of the jet outflow can be further boosted by neutrino pair-annihilation (expected to be subdominant; not included here) and dissipation of magnetic and thermal energy into kinetic energy at larger spatial scales, up to $\Gamma\lesssim-u_{0}(h/h_{\infty}+b^{2}/\rho)\approx 5-10$. Here, $h_{\infty}$ denotes the EOS-specific asymptotic value of $h$ and $b$ is the co-moving magnetic field strength.

During the first 50 ms post-merger, the system ejects $\gtrsim 2\times 10^{-2}M_{\odot}$ of matter at a time-averaged rate of $\dot{M}\sim 0.5M_{\odot}\,{\rm s}^{-1}$. Before the emergence of the jet structure ($t\lesssim 30$ ms), spiral waves propagating outward in the accretion disk [nedora2019Spiralwave, nedora2021Dynamical] dominate angular momentum transport and mass ejection of the system of $\lesssim\!1\times 10^{-2}M_{\odot}$ with mass-averaged velocity $\langle v\rangle\approx 0.1c$. These waves and associated mass ejection, generated mainly by $m=2$ and $m=1$ non-axisymmetric density modes of the remnant through hydrodynamical instabilities [lehner2016Instability, radice2016Onearmed, east2016Relativistic], are visible as concentric waves in the radial velocity (Fig. 3), paralleled by oscillations of the unbound mass flux (Fig. 4). As the vicinity of the remnant becomes strongly magnetized and the jet structure emerges ($t\approx 30$ ms), the unbound mass flux increases by an order of magnitude to $\dot{M}\approx 1M_{\odot}\,{\rm s}^{-1}$ (Fig. 4), quickly dominating the total mass unbound post-merger. Within only $\approx 15$ ms a total of $\gtrsim 1\times 10^{-2}M_{\odot}$ of material ($\gtrsim 50\%$ of the cumulative mass ejected post-merger) is unbound with mass-averaged velocity $\langle v\rangle\approx 0.1c$.

As accretion blocks outward radial mass flux in equatorial regions over timescales of interest, the neutrino and magnetically-driven wind from the remnant totaling $1\times 10^{-3}M_{\odot}$ escapes in polar directions ($\theta\lesssim\!30^{\circ}$) with only $0.2\times 10^{-3}M_{\odot}$ being launched within the jet core. The dominant contribution to ejected material, however, is launched as winds from the accretion disk. Disk winds intensify after $\approx\!30$ ms when angular momentum transport by spiral waves through the compact bound merger debris, magnetic stresses in the vicinity of the star, and the onset of magnetohydrodynamic turbulence driven by the magnetorotational instability have established and enlarged the accretion disk to a radius of $\gtrsim\!150$ km with approximate inflow–outflow equilibrium. The onset of strongly enhanced disk winds also coincides with the first cycles of an emerging dynamo as evident from a ‘butterfly diagram’ similar to that obtained in previous work [siegel2018Threedimensional]. Despite intense neutrino irradiation from the remnant, the disk then settles into a self-regulated state of moderate electron degeneracy $\mu_{e}/k_{\rm B}T\sim 1$ [siegel2017ThreeDimensional, siegel2018Threedimensional], which implies high neutron-richness of $Y_{e}\approx 0.1-0.15$ [beloborodov2003nuclear, chen2007Neutrinocooled, siegel2017ThreeDimensional] (Fig. 5). The mass averaged $Y_{e}$ of the disk indeed shifts from $\approx\!0.25$ ($t<30$ ms) to $\approx\!0.15$ as it approaches a quasi-stationary state with an accretion rate of $\gtrsim\!1M_{\odot}\,{\rm s}^{-1}$ and mass of $\approx\!0.19M_{\odot}$.

Nucleosynthesis & kilonova.—Figure 6 shows properties of unbound outflows at the onset of neutron capture reactions ($T\approx\!5$ GK) as sampled by multiple families of $\approx 2\times 10^{4}$ unbound passive tracer particles injected into the simulation domain (see CS23 for details of tracer placement). Fast outflow speeds $>0.2c$ are almost exclusively associated with polar outflows. As a result of neutrino absorption but high outflow speeds due to magnetic fields, material ejected from the highly neutron-rich degenerate surface layer of the star is protonized to asymptotic values of $Y_{\rm e}\approx 0.3-0.4$, much lower than $Y_{e}\approx 0.5$ as in purely neutrino-driven winds of hot proto-neutron stars, even in the presence of fast rotation [desai2022three]. Outflows from the self-regulated neutron-rich reservoir of the disk are protonized by absorption of intense neutrino radiation from the remnant (cf. Figs. 2, 5) to a mass averaged value of $\langle Y_{e}\rangle\approx 0.3$ at 5 GK.

Nucleosynthesis calculations based on the unbound tracer particles are conducted with the nuclear reaction network SkyNet [lippuner2017SkyNet] using 7843 nuclides and 140 000 nuclear reactions with the setup described in Refs. [li2021Neutrino] and CS23. They start in nuclear statistical equilibrium at a temperature of $T=7$ GK and take neutrino irradiation into account using neutrino fluxes directly extracted from our simulation as in Ref. [li2021Neutrino]. Final abundances at $t=10^{9}$ s are shown in Fig. 6. Elements beyond the 2nd r-process peak ($A\approx 130$) are entirely suppressed, similar to results of long-term 2D hydrodynamical simulations of similar systems [fujibayashi2020mass], but different from Ref. [curtis2021Process], which may result in part from differences in the treatment of neutrino absorption in the nuclear network calculation.

We compute kilonova light curves based on angular-dependent ejecta mass profiles extracted from the simulation, using the method presented in CS23. The resulting kilonova signal from post-merger ejecta is consistent with observations of GW170817 in the UV and blue bands up to several days (Fig. 7). Underestimation on timescales $\gtrsim\!5$ d can be explained by additional ‘redder’ (lanthanide bearing) components [villar2017Combined] not included here, which can be generated by neutron-richer accretion disk winds upon collapse of the remnant into a BH [siegel2017ThreeDimensional, siegel2018Threedimensional, kasen2017Origin, li2021Neutrino]. The $\sim\!{\rm day}$ kilonova is determined by the disk outflows. Fast material from the jet region carries most of the kinetic energy but only 10% of the total ejected mass; the latter dominates the signal on a few-hours timescale and can boost the $\sim\!{\rm hr}$ UV/blue signal depending on the line of sight. In direction of the jet, the signal is enhanced by up to 1.5 magnitudes, reaching similar luminosities to the kilonova precursor signal from free neutron decay in fast dynamical ejecta [metzger2015Neutronpowered] when the boost due to relativistic effects is taken into account (CS23).

Conclusion.—These results provide strong evidence for massive ($\gtrsim 10^{-2}M_{\odot}$) kilonovae such as GW170817 with early ($\sim\!{\rm day}$) blue and late-time ($\sim\!{\rm week}$) red emission being dominated by post-merger disk outflows. This provides additional support to the conjecture of Ref. [siegel2022r] that outflows from accretion disks are the main astrophysical site of the Galactic r-process. In particular, we show that binaries consistent with GW170817 and typical of galactic BNS require a remnant lifetime of only $\approx\!50$ ms to generate a lanthanide-poor blue kilonova component of $\gtrsim\!2\times 10^{-2}\,M_{\odot}$ with expansion velocity $v\approx 0.1c$ and lightcurves consistent with GW170817. While we find elements of previously proposed mechanisms such as magnetized winds from the remnant [metzger2018Magnetar] and spiral waves [nedora2019Spiralwave], bulk mass ejection here is due to a combination of a magnetic jet ($\sigma\sim 5-10$) that emerges $\lesssim 30$ ms post-merger, associated global magnetic stresses, and the onset of magnetohydrodynamic turbulence, which reconfigure the accretion disk, enhancing outflows that quickly dominate the cumulative post-merger mass ejected. At $\gtrsim\!50$ ms post-merger the accretion disk with mass of $\approx\!0.19M_{\odot}$ and accretion rate of $\gtrsim\!1M_{\odot}\,{\rm s}^{-1}$ is in a self-regulated neutrino-cooled state with properties in good agreement with initial conditions of previous work [siegel2017ThreeDimensional, siegel2018Threedimensional, li2021Neutrino]. We conclude that upon collapse of the remnant and its neutrino irradiation, lanthanide-bearing outflows of $\gtrsim 0.05\,M_{\odot}$ ($\sim\!30\%$ of the remaining disk mass [siegel2017ThreeDimensional, siegel2018Threedimensional, de2021Igniting, fujibayashi2020postmerger, fernandez2019Longterm, christie2019role]) consistent with the red kilonova emission of GW170817 are generated over the subsequent few hundred milliseconds [siegel2017ThreeDimensional, siegel2018Threedimensional, li2021Neutrino, fernandez2019Longterm].

The rapid and self-consistent emergence of a weakly magnetized ($\sigma\sim 0.1$), mildly relativistic ($v\lesssim\!0.6c$) wind from a merger remnant reported here leads to a $\sim\!{\rm hr}$ UV/blue kilonova signal that can be degenerate with the kilonova precursor signal from free neutron decay in the fast tail of dynamical merger ejecta. This novel precursor provides an additional discriminant to distinguish between BNS and NS–black-hole mergers and highlights the importance of early UV and optical follow-up observations of future merger events. The successful break-out of the jet from the surrounding merger debris found here may have additional non-thermal emission signatures.

The results also suggest a novel formation mechanism for the central engine of short gamma-ray bursts. Upon collapse of the remnant, the magnetic jet ‘seeds’ the black hole with magnetic flux and forms a strongly magnetized ($\sigma=L_{\rm EM}/\dot{M}\gg 1$), highly-relativistic jet, owing to the absence of baryon loading from the stellar wind. The emergence of a jet structure from small-scale turbulent fields in the final BH–disk configuration of binary systems with long-lived remnants has so far remained elusive in numerical simulations. The results here suggest that jets from the remnant could ‘seed’ gamma-ray burst jets for remnant lifetimes of $\gtrsim\,30$ ms.