Secular Outflows from Long-Lived Neutron Star Merger Remnants

We study mass ejection from a binary neutron star merger producing a long-lived massive neutron star remnant with general-relativistic neutrino-radiation hydrodynamics simulations. In addition to outflows generated by shocks and tidal torques during and shortly after the merger, we observe the appearance of a wind driven by spiral density waves in the disk. This spiral-wave-driven outflow is predominantly located close to the disk orbital plane and have a broad distribution of electron fractions. At higher latitudes, a high electron-fraction wind is driven by neutrino radiation. The combined nucleosynthesis yields from all the ejecta components is in good agreement with Solar abundance measurements.


Introduction
Neutron star merger outflows are thought to be one of the main astrophysical sites of production of heavy r-process elements [1].Strong evidence linking mergers and r-process nucleosynthesis is provided by AT2017gfo, the kilonova associated with the neutron star merger GW170817.In particular, models of the light-curve and spectra of AT2017gfo require the presence of rare earth elements in the ejecta [2,3,4,5,6,7].More recently, kilonova-like signals were observed following the long gamma-ray bursts 211211A [8,9,10] and 230307A [11,12,13], which were also interpreted as binary neutron star mergers [14].However, with the possible exception of strontium [15,16] and tellurium [17], it has not been possible to identify specific heavy r-process elements in AT2017gfo, or other kilonovae.As such, nucleosynthesis yields need to be estimated from merger simulations.
A major difficulty for theoretical models is that the outflows in neutron star mergers are generated through a variety of mechanisms operating over different timescales [18].Some of the ejecta is generated by tidal torques and shocks over a timescale of milliseconds during and shortly after the merger [19,20,21,22,23,24,25,26,27,28,29,30,31] the so-called dynamical ejecta.Winds from the remnant operate on a timescale of seconds, the so-called secular ejecta [32,33,34,35,36,37,38,39,40,41,42,43,44].Most studies have focused on 3.0 Figure 1.Baryonic rest mass of the remnant (ρ > 10 13 g cm −3 ), of the disk (ρ < 10 13 g cm −3 ), and of the ejecta (−hu t ≥ h min ), and mass outflow rate as a function of time from merger.The disk is initially formed of material expelled from the remnant on a timescale of ∼20 ms.It then evolves on a secular timescale under the effect of the spiral-wave instability.The disk mass evolution curve has been smoothed using a square window with width of 1 ms.
either the secular, or the dynamical ejecta.A notable exception are the recent works by Kiuchi and collaborators [45], considering the dynamical and secular ejecta from a binary neutron star system producing a black hole shortly after merger, and of Hayashi et al. [46], considering the long-term evolution of a neutron-star black hole merger remnant.However, such simulations are extremely challenging because, while fully general-relativistic codes are required to model the merger, the large computational costs of these simulations make them unsuitable for long-term integration.On the other hand, various groups have successfully modeled the secular ejecta, particularly from black-hole torus systems, while neglecting the back-reaction of the metric, or in Newtonian gravity.This approach is viable for remnants forming black-holes, but not for mergers producing long-lived massive neutron star remnants.In the latter case, a major outflow component is the spiral-wave wind, which is driven by hydrodynamic instabilities inside the remnant.Capturing such instabilities requires the same-kind of full-3D simulations with dynamical spacetime that are used to model the merger phase.Understanding such mergers is important, since it is believed that a substantial fraction of mergers produce long-lived remnants.GW170817 is also likely to have resulted in the formation of a meta-stable remnant [47,48,49,50,51,52,53,44,54,55].In Ref. [56], we have reported on the first ab-initio, general-relativistic neutrino-radiation hydrodynamics simulations of a long-lived neutron star merger remnant to secular timescales.We reported on the structure of the remnant and of the disk and on their evolution in the first 100 ms after the merger.Here, we study the mass ejection and nucleosynthesis yield from our models.We find that a massive outflow is launched by the spiral-wave wind mechanism, in agreement with previous studies with more approximate neutrino treatment.We discuss the geometry and composition of the outflows.Finally, we show that our model predicts an r-process nucleosynthesis pattern in good agreement with Solar abundance measurements.

Methods
We consider a binary with component masses 1.35 M ⊙ and 1.35 M ⊙ .We adopt the HS(DD2) equation of state [57,58], which predicts a maximum nonrotating neutron-star mass of 2.42 M ⊙ and the radius of a nonrotating 1.4 M ⊙ neutron star to be R 1.4 = 13.2 km.We construct irrotational initial data with the Lorene pseudo-spectral code [59].The initial separation is of 50 km.Evolutions are carried out using the gray moment-based (M1) general-relativistic neutrino-radiation hydrodynamics code THC M1 [60,61,62,63,64], which is based on the Einstein Toolkit [65,66].For the simulations discussed here, we use the Carpet adaptive mesh-refinement (AMR) driver [67,68], which implements the Berger-Oilger scheme with refluxing [69,70], and evolve the spacetime geometry using the CTGamma code [71,72], which solves the Z4c formulation of Einstein's equations [73,74].The simulations reported here neglect magnetic fields.On the one hand, it is known that if the stars are originally endowed with magnetar-level fields, the magnetic stresses can qualitatively affect the outflow dynamics [44].On the other hand, simulations starting with realistic initial fields and employing a sophisticated treatment to capture unresolved turbulence and dynamo action show that the magnetic fields are unlikely to be dynamically significant over the timescales considered here [75,76].In particular, the spiral-wave wind is expected to dominate the angular momentum transport in the disk, for the values of turbulent viscosity inferred by those simulations [51].The results of additional simulations, including with more equations of state and with the inclusion of angular momentum transport due to turbulent viscosity will be reported elsewhere.
The evolution grid employs 7 levels of AMR.Our simulation is performed at two resolutions with the finest grid having finest grid spacing of h = 0.167 GM ⊙ /c 2 ≃ 247 m and h = 0.125 GM ⊙ /c 2 ≃ 185 m, respectively denoted as LR and SR setups.Here, we report the results from the higher resolution simulations.There lower-resolution data is in good quantitative and qualitative agreement.The curves have been smoothed using a square windows with width of 1 ms.

Results
The binary considered in this study results in the formation of a remnant that is dynamically stable, as documented in, e.g.[77].Figure 1 shows the evolution of the baryon mass density inside the remnant massive neutron star and in the disk.Following a standard convention, we flag material as belonging to the remnant massive neutron star if its density is larger than 10 13 g cm −3 .Material with lower density and situated within a cubical box of half-diameter 160 GM ⊙ /c 2 ≃ 236 km is consider as belonging to the disk.The latter is formed of matter that is squeezed out of the collisional interface between the neutron stars during and shortly after merger [26,78].After ∼20−30 ms mass ejection from the central remnant terminates.The subsequent disk mass evolution is primarily driven by mass loss to outflows with negligible accretion onto the central object.
The ejecta are extracted on a sphere of coordinate radius R = 300 GM ⊙ /c 2 ≃ 443 km, only material with outgoing radial velocity and that is unbound according to the Bernoulli condition1 is considered in our analysis.Our simulations show two distinct mass ejection components: a first pulse of material (∼3 × 10 −3 M ⊙ ) component ejected during and shortly after the merger, followed by a steady wind with Ṁ ≃ 0.05 M ⊙ s −1 .The first ejecta component comprises the tidal tails of the stars and the shocked ejecta launched shortly after merger [26].The second component is dominated by the spiral-wave wind [51].However, we also observe the emergence of a neutrino-driven wind at high latitudes [28,64].
The spiral-wave wind dominates the overall mass outflow rate in the postmerger.This wind is the result of hydrodynamic and gravitational torques exerted by the remnant on the disk.The former is deformed as a result of the so-called low-T /|W | instability [80,81,82,83,84,85].
To diagnose this we show in Fig. 2 the density contrast on the orbital plane of the disk: where ρ is the fluid rest-mass density and ⟨•⟩ denotes average in the azimuthal direction.The m = 2 pattern in the remnant is clearly visible at the center, as are the spiral waves.These waves carry angular momentum outwards and drive the mass outflow, as anticipated by previous longterm post-merger simulations with simplified neutrino physics [51,28].This spiral wave pattern forms shortly after merger and persists for the entire duration of our simulations.The spiralwave wind is predominantly emitted close to the disk orbital plane and has a broad distribution of electron fractions (see Fig. 4

below).
Figure 3 shows the evolution of the density modes on the orbital plane, as a function of time.The odd-m modes are initially zero, because equal mass binaries have a πsymmetry2 .However, the symmetry is broken by turbulence during merger and subsequently the odd-m modes are excited by the low-T /|W | instability [84,86].The even m modes are present from the beginning, particularly the m = 2 mode.In absence of the low-T /|W | instability these modes would be efficiently damped by gravitational radiation [87].However, in our simulation they are continuously driven by the instability and saturate at a roughly constant amplitude or decay slowly.The one-arm (m = 1 mode) and the triaxial (m = 3 mode) instabilities dominate the remnant dynamics from ∼ 20 ms to ∼60 ms postmerger.Afterwards, the m = 1, 2, 3 modes equally contribute to the remnant dynamics.This has been observed in both isolated neutron star [81,88] and binary mergers [84] simulations, although the details of the evolution are dependent on the equation of state, the compactness and possibly the microphysics [86].
The gravitational waves associated to these nonaxisymmetric modes carry the imprinting of the equation of state but, even under optimistic assumption, an observation will likely require third-generation detectors and a source at distance of ∼30 Mpc [84].We caution the reader that, because our code uses a Cartesian grid, the m = 4 mode might be artificially enhanced.

Relative abundances
. Nucleosynthesis yields of the ejecta (solid line) and Solar r-process abundances from Ref. [92] (points).The spiral-wave wind produces a full r-process yield.
While the structure of the wind close to the equator is similar to that observed in simulations with simplified M0 neutrino transport we have published in the past [51,28,89], the wind at high latitude is substantially enhanced in our new M1 simulations, as we have also documented for a lower resolution simulation in [64].The angular structure of the ejecta is shown in Fig. 4, upper panel.This figure should be contrasted with Fig. 1 of [6], where the same diagnostics is shown from our older simulations using M0 neutrino transport.In contrast to the M0 simulations, which had angular profiles well fitted by sin 2 (θ) -θ being the angle from the equatorial plane -the new M1 simulations show a shallower profile, closer to sin 3/2 (θ).The composition of the ejecta is also quite different: the M0 simulations predict ejecta with electron fraction Y e < 0.4 at all latitudes, while the newer M1 simulations predict a broader distribution with Y e ≲ 0.5.These differences arise from the fact that the M0 scheme artificially suppresses neutrino heating in the polar region of the remnant, suppressing the wind [64].
Figure 5 shows the combined nucleosynthesis yields for all ejecta components.These have been obtained using the SkyNet code [90] with the procedure discussed in Ref. [26].We find that our model well reproduces the Solar r-process pattern.The only significant discrepancy is in the region with A ∼ 140, however those can be ascribed to the particular nuclear-mass model used in our nucleosynthesis calculations [91].

Conclusions
We have performed long-term simulations of a binary neutron star merger producing a longlived remnant with a general-relativistic neutrino-radiation hydrodynamics code.We find that such events drive massive outflows through a combination of tidal torques and shocks, during and shortly after the merger, and spiral-wave and neutrino-driven winds after the merger.The dynamics of the spiral-wave wind is in good agreement with predictions from simulations with lower-fidelity neutrino transport.However, there are substantial differences in the dynamics and composition of the neutrino-driven wind.In particular, our previous studies underestimated both the mass and the electron fraction of this component of the ejecta.The combined nucleosynthesis yield of all ejecta components is in good agreement with the Solar r-process abundance pattern, indicating that neutron star mergers producing long-lived remnants provide a natural environment for the formation of these elements.

Figure 2 .
Figure 2. Density contrast in the disk 111.36 ms after merger.The spiral-wave in the disk is generated by the m = 2 deformation of the remnant (the so-called bar mode) and persists throughout our simulation.

6 Figure 3 .
Figure 3. Normalized density modes on the xy-plane (see Eq. 2).The m = 1, 2, 3 density modes are the dominant modes.They are driven by the low-T /|W | instability after merger.The curves have been smoothed using a square windows with width of 1 ms.

Figure 4 .
Figure 4. Angular distribution and electron fraction of the ejecta.Most of the ejecta is equatorial and has a broad range of electron fractions.The polar ejecta is less massive and characterized by significantly larger electron fraction Y e ≳ 0.4.