End-to-end Kilonova Models of Neutron Star Mergers with Delayed Black Hole Formation

The recent first multimessenger observation of a binary neutron star (NS) merger, GW170817/AT 2017gfo (e.g., Abbott et al. 2017; Villar et al. 2017; Watson et al. 2019), lends strong support to the idea (Lattimer et al. 1977) that NS mergers are indeed significant sites of rapid neutron-capture (r-) process nucleosynthesis (Arnould et al. 2007; Cowan et al. 2021). Simulations of NS mergers and their aftermath predict that r-process viable outflows can be produced in each of the following three phases: during and right after the merger (called dynamical ejecta; e.g., Goriely et al. 2011; Korobkin et al. 2012; Wanajo et al. 2014), from an NS torus remnant in the case that it is formed (e.g., Metzger & Fernández 2014; Perego et al. 2014; Fujibayashi et al. 2018; Mösta et al. 2020), and from a black hole (BH) torus remnant formed promptly or after collapse of an NS remnant (e.g., Fernández & Metzger 2013; Just et al. 2015a; Siegel & Metzger 2018; Fujibayashi et al. 2020a). Depending on the masses of the initial two NSs and the nuclear equation of state (EOS), those components can make different contributions to and thus have different relative importance for the nucleosynthesis yields and the electromagnetic kilonova (KN) counterpart.

The role of each component is not well constrained so far, either theoretically or observationally (i.e., based on AT 2017gfo). One reason is that most previous theoretical studies treat each component individually or separately, and so far, only a few studies discuss models with a consistent inclusion of all components. Fujibayashi et al. (2020b) and Shibata et al. (2021), in combination with the corresponding KN studies of Kawaguchi et al. (2021, 2022), reported models with a very long-lived NS remnant, which do not produce ejecta from a BH torus system. On the other hand, Fujibayashi et al. (2023) considered systems in which matter ejection from the NS remnant is terminated early on because the remnant collapses shortly after the merger. Recently, Kiuchi et al. (2022) reported a magnetohydrodynamic (MHD) simulation covering the first second of evolution of a similarly short-lived case, which confirmed the basic results by Fujibayashi et al. (2023) obtained using a more approximate α-viscosity scheme (Shakura & Sunyaev 1973).

In this Letter, we present the first end-to-end models of NS mergers with intermediate remnant lifetimes (between ∼0.1 and 1 s). These systems are distinguished from the aforementioned scenarios, as they yield roughly comparable amounts of all three types of ejecta. Different from previous long-term evolution models of NS remnants (e.g., Metzger & Fernández 2014; Perego et al. 2014; Fujibayashi et al. 2020a), our simulations adopt an energy-dependent neutrino transport scheme, as well as an improved α-viscosity scheme guided by MHD results.

These first viscous neutrino transport models of mergers with significantly delayed BH formation lead to several new insights. (1) The lifetime of the NS remnant in such systems is shorter for asymmetric than for symmetric binaries, and it depends sensitively on the viscosity inside the NS. (2) Ejecta launched during the BH torus phase are less neutron-rich than predicted by models using manually constructed initial conditions. (3) In the considered systems of intermediate lifetimes, the synthesis of lanthanides and heavier elements is not efficient enough to explain the solar pattern. (4) The combination of all ejecta components is significantly more anisotropic than just the dynamical ejecta because of a massive, dominantly polar, neutrino-driven outflow from the NS remnant. (5) The KN produced by the combined ejecta can (may not) shine bright enough to explain AT 2017gfo at late (early) times. (6) For a given viscosity, both the summed mass of all ejecta components and their individual contributions are systematically higher for asymmetric than for equal-mass binaries.

After outlining our model setup in Section 2, we will report on the aforementioned findings in Section 3 and discuss some implications in Section 4. The Appendices provide additional information regarding selected properties of our models.

Each model consists of three successive hydrodynamics simulations and two postprocessing steps that provide the nucleosynthesis yields and KN light curve. The hydrodynamical evolution of the merger is followed by a 3D general relativistic (GR) smoothed particle hydrodynamics (SPH) code (Oechslin et al. 2007; Bauswein et al. 2010) that employs a modern leakage-plus-absorption scheme (ILEAS; Ardevol-Pulpillo et al. 2019) to describe neutrino cooling and heating, including electron neutrinos (ν_e ), electron antineutrinos (${\bar{\nu }}_{e}$), and a third species (ν_x ) representative of all heavy-lepton neutrinos. At a postmerger time, t_pm, of t_pm = t_map = 10 ms, we azimuthally average the SPH configuration, map it to a spherical polar grid, and, assuming axisymmetry, continue the postmerger evolution using the special relativistic code ALCAR-AENUS (Obergaulinger 2008; Just et al. 2015b), which adopts an energy-dependent M1 neutrino transport scheme. We employ the same GR corrections in the transport equations and the same neutrino interaction rates and formulations that were used in Just et al. (2018; and are based on Bruenn 1985; Hannestad & Raffelt 1998; Pons et al. 2000; Horowitz 2002), except that in the present study, we neglect inelastic neutrino-electron scattering and use the approximation by O'Connor (2015) to describe pair processes. For the transition from the SPH simulations (which do not evolve local neutrino energy and flux densities), the neutrino energies are initially (i.e., at t_pm = t_map) set to Fermi distributions corresponding to the local thermodynamic state above densities of 10⁹ g cm⁻³ and vanish everywhere else, and the fluxes vanish everywhere.

The numerical settings adopted in the SPH simulations are the same as detailed in Ardevol-Pulpillo et al. (2019) and Kullmann et al. (2021). The total number of SPH particles is 3 × 10⁵, and the neutrino source terms are computed on a uniform Cartesian grid having 305 cells of size 738 m along each direction (see Appendix A for a test of the neutrino grid resolution; uncertainties in the dynamical ejecta masses due to the particle resolution are estimated to be several tens of percent (Bauswein et al. 2013), comparable to the estimated error bars of grid-based merger simulations (e.g., Radice et al. 2018)). The postmerger simulations are conducted using a radial grid with a constant cell size of Δr = 100 m within radii of r < 20 km and afterward increasing by ≈2.3% per cell, as well as a uniform polar grid with a resolution of 225. The neutrino energy range between zero and 400 MeV is discretized using 15 bins, the size of which increases by 40% per bin. In order to prevent extremely small time integration steps, we assume a uniformly rotating core with spherically symmetric thermodynamic properties at radii below 1.5 km. We verified that postmerger models with a higher resolution and smaller 1D core produce essentially the same results.

In order to reduce the inconsistency between the curved-spacetime merger models and flat-spacetime postmerger models, we map the primitive variables such that the radial volume element of the SPH model ${\psi }^{6}d({r}_{\mathrm{SPH}}^{3}/3)$ (with conformal factor ψ; see Oechslin et al. 2007) equals the postmerger volume element d(r³/3) along each radial direction, thus approximately preserving the volume integrals of the conserved variables (baryonic rest mass, etc.), and we define the three-velocities in the postmerger model as functions of the corresponding SPH velocities (see Oechslin et al. 2007) as ${v}^{i}={v}_{\mathrm{SPH}}^{i}{\psi }^{2}$. Gravitation is treated by solving a Poisson equation augmented with relativistic corrections (Müller et al. 2008), in which the monopole contribution is replaced either by an effective Tolman–Oppenheimer–Volkof (TOV) potential (Marek et al. 2006; for times t_pm earlier than the time of BH formation, t_BH) or by the pseudo-Newtonian BH potential of Artemova et al. (1996; for t_pm > t_BH). Once the NS remnant becomes gravitationally unstable (i.e., at t_pm = t_BH), we replace the innermost region with an outflow boundary mimicking the central BH while consistently updating its size, mass, and angular momentum through time integration of the boundary fluxes.

For describing turbulent viscosity driven by the magnetorotational instability (MRI), we extend the classical α-viscosity scheme of Shakura & Sunyaev (1973) such that it can capture MRI-related viscosity in both the rotation-supported regime (i.e., the accretion torus) and the pressure-supported regime (i.e., the NS remnant). In a pressure-supported object with a subsonic shear velocity, the MRI-driven viscosity is indeed not expected to scale with the sound speed but rather to behave in a quasi-incompressible manner (Reboul-Salze et al. 2021, 2022). Our formulation therefore expresses the kinematic viscosity as

$$\begin{eqnarray}&&{\nu }_{\mathrm{vis}}={\alpha }_{\mathrm{vis}}{H}_{\mathrm{vis}}^{2}| {\rm{\Omega }}| {\tilde{q}}^{{n}_{\mathrm{vis}}},\end{eqnarray} \tag{ 1 }$$

namely, the product of the generalized characteristic length scale

$$\begin{eqnarray}&&{H}_{\mathrm{vis}}=\min \{\left|\rho /{\rm{\nabla }}\rho \right|,r,{c}_{i}/{{\rm{\Omega }}}_{{\rm{K}}}\}\end{eqnarray} \tag{ 2 }$$

(with density ρ, spherical radius r, isothermal sound speed ${c}_{i}=\sqrt{P/\rho }$, gas pressure P, and Keplerian angular velocity Ω_K) and the characteristic velocity scale H_visΩ (with angular velocity Ω). The additional quenching factor

$$\begin{eqnarray}{\tilde{q}}^{{n}_{\mathrm{vis}}}=\left\{\begin{array}{ll}0, & \mathrm{if}\ d{\rm{\Omega }}/{dR}\gt 0,\\ 1, & \mathrm{else}\ \mathrm{if}\ \left|d\mathrm{ln}{\rm{\Omega }}/d\mathrm{ln}R\right|\gt {q}_{0},\\ {\left(\displaystyle \frac{1}{{q}_{0}}\left|\displaystyle \frac{d\mathrm{ln}{\rm{\Omega }}}{d\mathrm{ln}R}\right|\right)}^{{n}_{\mathrm{vis}}}, & \mathrm{else}\end{array}\right.\end{eqnarray} \tag{ 3 }$$

(with cylindrical radius R) accounts for the tendency (e.g., Pessah et al. 2008) of the MRI to be reduced in regions where the shear is sub-Keplerian, i.e., $q=\left|d\mathrm{ln}{\rm{\Omega }}/d\mathrm{ln}R\right|\lt {q}_{0}\sim 1.5$. ⁹ The parameter n_vis thus varies the strength of the turbulent viscosity in the NS remnant (where q < q₀) relatively independently of that in the surrounding disk (where q ≃ q₀). This allows one to parametrically explore the sensitivity to the viscosity inside the NS remnant (which is poorly constrained so far by existing simulations; Kiuchi et al. 2018; Margalit et al. 2022; Palenzuela et al. 2022) while keeping the viscosity in the disk (which is known to be fairly well reproduced by a conventional α-viscosity scheme; Fernández et al. 2019; Hayashi et al. 2022; Just et al. 2022b) unchanged. Significant uncertainty also comes from the dependence on the diffusive processes through the magnetic Prandtl number (Guilet et al. 2022; Held & Mamatsashvili 2022), which justifies exploring different parameter values.

At t_pm = 10 s, the inner (outer) radial boundary is moved to a radius of 10⁴ km (4 × 10⁷ km), and a third simulation is conducted to follow the expansion of just the ejected material until t_pm = 100 s. The ejecta configuration at t_pm = 100 s is assumed to be homologous, with r(t_pm) = v t_pm, and equatorially symmetric and gets sampled in the northern hemisphere up to velocities of 0.7c by 1000–2000 tracer particles per model, the time evolution of which is obtained by path integration backward in time using the available simulation outputs. The sampling of the dynamical ejecta takes into account the entire evolution; i.e., it also utilizes data from the SPH simulations by splitting the trajectories constructed from postmerger simulation data at t_pm = t_map and associating them with a number of SPH particles. Since this step, by which the effective number of tracers is increased to 4000–5000 per model, is nontrivial, we provide further details on this procedure in Appendix D. The tracers are input to a nuclear network solver that predicts the nucleosynthesis yields. We use two independent solvers here (hereafter networks A and B), allowing us to cross-validate the yields and heating rates and isolate uncertainties related to the network code and its physics input from other modeling uncertainties. Network A (used previously in, e.g., Goriely et al. 2011; Just et al. 2015a; Kullmann et al. 2023) takes nuclear ingredients from experiments where available and, where not, from theoretical models, namely, nuclear masses from the BSkG2 mass model (Ryssens et al. 2022); β-decay rates from Marketin et al. (2016); reaction rates from TALYS estimates with microscopic inputs (Goriely et al. 2018), including BSkG2 masses; and fission probabilities and fragment distributions from Lemaître et al. (2021). Network B (employed previously in, e.g., Wu et al. 2016; Collins et al. 2023) uses the reaction rates for neutron captures, photodissociation, and fission based on the HFB21 mass model (Goriely et al. 2010) as described in Mendoza-Temis et al. (2015) and β-decay rates from Marketin et al. (2016).

Finally, for assessing the KN light curve, the tracers, including their composition and radioactive heating rates, are used as input for an approximate photon transport scheme to estimate the KN light curve in the same way as detailed in Just et al. (2022).

Table 1 summarizes the parameters for all investigated models. We consider both a symmetric (ratio of gravitational masses of M₁/M₂ = 1) and an asymmetric (M₁/M₂ = 0.75) progenitor configuration (with M₁ + M₂ = 2.75 M_⊙), and for both cases, we vary the viscosity parameters n_vis ∈ {0.5, 1, 10} and α_vis ∈ {0.03, 0.06}. The SFHo EOS (Steiner et al. 2013), extended to low densities with a four-species EOS (e.g., Just et al. 2015a), is adopted. ¹⁰

Table 1. Model Properties and Results

Model Name	Mass Ratio	nvis	αvis	tBH	${m}_{\mathrm{tor}}^{\mathrm{BH}}$	${Y}_{e,\mathrm{tor}}^{\mathrm{BH}}$	${m}_{\mathrm{ej}}^{\mathrm{total}}$	${m}_{\mathrm{ej}}^{\mathrm{dyn}/\mathrm{NS}/\mathrm{BH}}$	${Y}_{e,\mathrm{ej}}^{\mathrm{dyn}/\mathrm{NS}/\mathrm{BH}}$	${v}_{\mathrm{ej}}^{\mathrm{dyn}/\mathrm{NS}/\mathrm{BH}}$	${X}_{\mathrm{LA}}^{\mathrm{dyn}/\mathrm{NS}/\mathrm{BH}}$
				(ms)	(10−2 M⊙)		(10−3 M⊙)	(10−3 M⊙)		(10−2 c)	(10−3)
sym-n1-a6	1	1	0.06	122	12.5	0.257	74	6/20/47	0.24/0.41/0.30	22/18/5.5	142/0.00/2.85
sym-n05-a3	1	0.5	0.03	186	13.4	0.214	57	6/21/31	0.23/0.42/0.30	22/16/4.3	135/0.00/8.56
sym-n05-a6	1	0.5	0.06	104	14.1	0.255	76	6/18/52	0.24/0.42/0.31	22/20/5.8	143/0.00/2.37
sym-n10-a3	1	10	0.03	915	1.78	0.317	33	6/21/5	0.24/0.39/0.32	21/11/3.8	141/0.36/0.73
sym-n10-a6	1	10	0.06	815	1.87	0.318	37	6/29/2	0.23/0.37/0.33	21/10/5.1	147/0.52/0.00
asy-n1-a6	0.75	1	0.06	96	16.2	0.250	86	11/21/55	0.25/0.41/0.30	22/20/5.9	125/0.06/4.90
asy-n05-a3	0.75	0.5	0.03	148	16.3	0.224	71	11/25/35	0.25/0.41/0.31	20/17/4.8	118/0.21/7.03
asy-n05-a6	0.75	0.5	0.06	88	17.5	0.250	87	11/21/56	0.25/0.41/0.31	21/20/6.1	121/0.16/6.90
asy-n10-a3	0.75	10	0.03	680	5.57	0.252	61	11/30/20	0.25/0.39/0.31	19/13/3.5	126/0.08/6.78
asy-n10-a6	0.75	10	0.06	520	5.77	0.283	76	13/39/24	0.24/0.37/0.29	18/12/5.0	131/4.40/3.97

Note. From left to right: model name, mass ratio, the two parameters entering the viscosity scheme, time of BH formation, mass and average electron fraction of the torus measured at t_BH + 10 ms, total ejecta mass, as well as for each ejecta component (dynamical/NS torus/BH torus ejecta) the mass, average electron fraction (measured at temperature of 5 GK), average velocity, and mass fraction of lanthanides plus actinides. Values for X_LA and q_heat were obtained using nuclear network A.

Download table as:  ASCIITypeset image

The following sections address the collapse behavior, torus properties, nucleosynthesis yields, ejecta geometry, and KN signal. Figure 1 illustrates snapshots at different times for model sym-n1-a6, and Figure 2 shows the time evolution of global properties for several models. The Y_e distribution, nucleosynthesis yields, and radioactive heating rate are depicted in Figure 3, and KN observables are provided in Figure 4. We adopt the (somewhat ambiguous) criterion ρ > 10¹² g cm⁻³ to discriminate material located in the NS remnant from that in the surrounding torus.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Lifetime Dependence on Mass Ratio and Viscosity

3.2. Torus Properties at BH Formation

3.3. Ejecta Interaction and Geometry

3.4. Nucleosynthesis Yields

3.5. Comparison with AT 2017gfo

A major unknown in merger models containing metastable NS remnants is the effective viscosity produced by MHD effects inside the NS. Compared to purely hydrodynamic systems, viscous merger remnants are able to tap the large reservoir of rotational and gravitational energy in the system and partially convert it to thermal energy. Since neutrino emission rates grow with high powers of the temperature, viscous merger remnants are stronger sources of neutrinos than nonviscous remnants, implying that their NDWs are expected to be more powerful compared to nonviscous rotating NSs or nonrotating proto-NSs (of which the NDWs have been extensively studied in the past; e.g., Fischer et al. 2010; Hüdepohl et al. 2010). Using for the first time an energy-dependent neutrino transport scheme, our simulations confirm this expectation ¹⁵ and demonstrate that NDWs in viscous NS remnants can be as massive as a few percent of M_⊙ and exhibit velocity distributions reaching up to v/c ∼ 0.5c. Judging from this result, mechanisms invoking genuinely magnetically driven outflows (Metzger et al. 2018; Ciolfi & Kalinani 2020; Mösta et al. 2020; Shibata et al. 2021) may not be necessary in order to explain the early blue component of AT 2017gfo.

The adopted viscosity scheme is a parameterization of actually more complex physics. Our models aim at bracketing this uncertainty by generalizing the original Shakura–Sunyaev viscosity and introducing the viscosity parameter n_vis that effectively regulates the viscosity just inside the NS. Despite a significant impact on the BH formation time, t_BH, of about 1 order of magnitude, the total ejecta masses only vary by a factor of ∼2 (and even less in asymmetric merger models). Moreover, we find that the nucleosynthesis pattern is relatively robust with variations of n_vis, and the KN light curve varies only moderately. Although the absolute values of the total ejecta mass depend on the viscosity, we find for our (admittedly small) set of models that the total ejecta mass is systematically increased for asymmetric binaries, independent of the chosen set of viscosity parameters. This implies that for two observations with similar chirp mass, the mass ratios can be related to each other, possibly allowing mass-ratio constraints based on the KN properties.

If we assume that our models are representative, the systematic underproduction found for A > 140 elements would imply that mergers with intermediate (and probably also long, as suggested by Fujibayashi et al. 2020b) NS remnant lifetimes could not be main r-process sites but would likely be dominated by events that do not overproduce light relative to heavy elements, such as prompt or briefly delayed collapse scenarios (Just et al. 2015b; Fujibayashi et al. 2023).

As for the KN, our results suggest the considered delayed-collapse scenario to be a viable progenitor for GW17087/AT 2017gfo based on the overall good agreement of the trends seen in the bolometric light curve, photospheric temperature, and photospheric velocity at times t_pm ≳ 1–3 days. We find three features in which our models agree less well with observational analyses of the early (t_pm ≲ few days) data from AT 2017gfo, namely, too-faint emission, a nonspherical photosphere, and partial enrichment of the photospheric region by lanthanides (see Section 3.5). We stress, however, that since radiative transfer calculations using detailed atomic data–based opacities are not yet available for delayed-collapse hydro models, our approximate KN modeling remains a nonnegligible source of uncertainty.

Since our set of models is not exhaustive, other combinations of binary masses, nuclear EOSs, and viscosity parameters may yield better agreement with AT 2017gfo. For instance, all three of the aforementioned inconsistencies could possibly be mitigated in cases where only a small amount of dynamical ejecta is launched such that the hot NDW is able to expand in a nearly spherical manner above and below the equatorial plane. Alternatively, it may also be possible that the anisotropic outflow structure seen in our present models is generic for delayed-collapse scenarios with significant NDW components. In this case, one would expect that systems with shorter NS remnant lifetimes, i.e., smaller binary masses, would generally yield more spherical ejecta distributions (see panels (a) and (b) of Figure 4) dominated by dynamical ejecta (for velocities v/c ≳ 0.1). The two competing scenarios for the relative role of NDWs (i.e., increase versus decrease of sphericity with NS remnant lifetime) could be distinguished by future KN observations from systems with different binary masses relative to AT 2017gfo using the P Cygni method developed by Sneppen et al. (2023). At any rate, since the degree of sphericity in either of the two aforementioned scenarios depends on the lifetime, the P Cygni method—in combination with end-to-end hydrodynamical models that capture all ejecta components—could be a powerful new tool for constraining binary properties and the nuclear EOS from KN observations.

We finally point out that our models, despite featuring a coherent end-to-end modeling strategy, still carry nonnegligible uncertainties connected to, e.g., the approximate treatment of GR, turbulent viscosity, and neutrino transport; the omission of magnetic fields; and the simplified KN physics. Moreover, our models are missing some physics ingredients that are potentially relevant to the KN problem, such as jets (e.g., Nativi et al. 2021), neutrino flavor conversion (e.g., Just et al. 2022a), and nonaxisymmetric NS oscillation modes that could produce additional ejecta components (Nedora et al. 2019).

We are grateful for inspiring discussions with Ninoy Rahman and Brian Metzger and to the anonymous referee for comments that improved the manuscript. O.J., V.V., and A.B. acknowledge support by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 759253. V.V. and A.B. acknowledge support by the Deutsche Forschungsgemeinschaft (DFG; German Research Foundation), project ID 138713538—Sonderforschungsbereich (Collaborative Research Center) SFB 881 ("The Milky Way System," subproject A10). Z.X. and G.M.P. acknowledge support by the ERC under the European Union's Horizon 2020 research and innovation program (ERC advanced grant KILONOVA No. 885281). O.J., Z.X., A.B., and G.M.P. acknowledge support by the DFG, project ID 279384907—SFB 1245. O.J., Z.X., A.B., G.M.P., and T.S. acknowledge support by the State of Hesse within the Cluster Project ELEMENTS. S.G. is an F.R.S.-FRNS research associate. This work has been supported by the Fonds de la Recherche Scientifique (FNRS; Belgium) and the Research Foundation Flanders (FWO; Belgium) under EOS project Nos. O022818F and O000422. J.G. acknowledges support from the European Research Council (MagBURST, grant agreement No. 715368). H.T.J. is grateful for support by the DFG through SFB-1258–283604770 "Neutrinos and Dark Matter in Astro- and Particle Physics (NDM)" and under Germany's Excellence Strategy through Cluster of Excellence ORIGINS (EXC-2094)-390783311. O.J. and Z.X. acknowledge computational support by the VIRGO cluster at GSI, and O.J. by the HOKUSAI computing facility at RIKEN. The nucleosynthesis calculations benefited from computational resources made available on the Tier-1 supercomputer of the Fédération Wallonie-Bruxelles, infrastructure funded by the Walloon Region under grant agreement No. 1117545. The publication is funded by the DFG—491382106, and by the Open Access Publishing Fund of GSI Helmholtzzentrum fuer Schwerionenforschung.

The ILEAS scheme coupled to the SPH solver, which is used to simulate the merger phase until mapping to the ALCAR code at t_map = 10 ms, computes the source terms related to neutrino emission and absorption on a uniform Cartesian grid. In this appendix, we briefly test the dependence of basic neutrino quantities on the grid resolution, to which end we pick a snapshot of our symmetric merger model at t_pm ≈ 5 ms and run only the ILEAS scheme on it, keeping all hydrodynamic quantities fixed. Figure 5 depicts the resulting polar angle–dependent isotropic-equivalent luminosities, L = 4π r² F_r (with radial neutrino flux F_r ), and mean energy, 〈〉, measured at r = 100 km for various cases of the grid resolution Δx. The resolution dependence turns out to be mild, suggesting that the grid width of Δx = 0.738 km adopted in our dynamical models is small enough to ensure that the grid-related discretization errors are subdominant (at least considering the early postmerger phase until t_pm = 10 ms, when the surface of the NS remnant is still relatively hot and the density gradient is rather shallow).

Zoom In Zoom Out Reset image size

As pointed out in Section 3.1, we suspect the shorter NS lifetimes seen for our asymmetric models to be a consequence of the smaller amount of angular momentum carried by the high-density NS compared to the symmetric case. Here we show that this tendency—i.e., smaller angular momentum carried by NS remnants of more asymmetric mergers—is not only fulfilled for the two SPH simulations in the main part of this study but also supported by another set of SPH simulations adopting a different total mass (M₁ + M₂ = 3 M_⊙); different mass ratios, q = M₁/M₂; and different nuclear EOS (MPA1; Müther et al. 1987), as well as by a similar set of models simulated with an entirely independent GR hydro code, namely, the Einstein Toolkit (ET; Etienne et al. 2021). In contrast to the SPH simulations, the ET does not assume the conformal flatness condition (CFC) to approximate GR but solves the full GR equations. The numerical setup is the same as that described in Soultanis et al. (2022).

In Figure 6, we compare for different values of q the angular momentum, J, as well as the specific angular momentum, j = J/M_b (with baryonic mass M_b ), each normalized to the corresponding instantaneous value of the q = 1 configuration. The angular momenta and masses are computed as (ignoring noncompact spacetime contributions to J for this qualitative analysis)

$$\begin{eqnarray}&&{J}_{\mathrm{NS}/\mathrm{tot}}=\int {d}^{3}x\,{\rho }_{* }\left(x{\hat{u}}_{y}-y{\hat{u}}_{x}\right),\end{eqnarray} \tag{ B1a }$$

$$\begin{eqnarray}&&{M}_{b,\mathrm{NS}/\mathrm{tot}}=\int {d}^{3}x\,{\rho }_{* },\end{eqnarray} \tag{ B1b }$$

where ${\rho }_{* }=\sqrt{\gamma }\rho W$ with the determinant of the spatial metric γ, Lorentz factor W, and ${\hat{u}}_{i}={{hu}}_{i}$ with the specific enthalpy h and the spatial component of the covariant fluid four-velocity u_i . As can be seen from the solid lines in Figure 6, if only the angular momentum carried by the NS remnant, J_NS, is considered (by restricting the integration in Equation (B1a) to regions with ρ > 10¹² g cm⁻³), both sets of simulations support the tendency of a faster reduction of J_NS and J_NS/M_b,NS for lower q, although in a somewhat less pronounced manner in the ET models than in the SPH models. On the other hand, if one considers the total angular momentum, J_tot (obtained from integration over the entire computational domain in Equation (B1a); see dashed lines in Figure 6), one rather obtains the opposite tendency, namely, a more slowly declining J_tot for lower q. In these purely hydrodynamic models, the total angular momentum can change only due to the emission of gravitational waves (assuming that the numerical nonconservation errors are small). Thus, the J_tot evolution likely reflects the fact that for a given total mass, the postmerger gravitational-wave emission in equal-mass systems is stronger, and therefore gravitational waves more efficiently reduce the J_tot than in asymmetric mergers (e.g., Kiuchi et al. 2020). Note, however, that very asymmetric configurations initially contain slightly smaller J_tot at the time of merging.

Zoom In Zoom Out Reset image size

Thus, the q-dependence of the angular momentum remaining in the NS remnant is mainly shaped by two counteracting effects, namely, gravitational-wave emission (which is less efficient for low q values than for q = 1) and redistribution of angular momentum into the low-density torus (which is more efficient for low q). Our results suggest that the second effect is stronger than the first one, ultimately leading to lower values of J_NS for q < 1 at t_pm ∼ 10–15 ms. We stress, however, that our analysis is rather tentative at this point, and that the above characteristics, as well as their consequences for the lifetime of the NS remnant, have to be tested and explored more systematically using models evolved for a longer time and with more realistic physics ingredients (i.e., neutrinos and magnetic fields), as well as including a detailed analysis of the impact of numerical discretization errors. In fact, we suspect that numerical effects (and to a lesser extent the different treatment of GR) are the main reason for the differences between the SPH and ET models in the current test. The relatively high numerical viscosity of the SPH models may dampen postmerger oscillations, and therefore reduce angular momentum losses by gravitational-wave emission, more strongly than in the ET models, which would explain the more pronounced reduction of J_NS with lower q in the SPH models. This is suggested by simulations performed with our recently developed moving mesh code employing CFC and the same gravitational-wave backreaction scheme (Lioutas et al. 2022), where we find the angular momentum loss to be quantitatively more comparable to full GR static mesh simulations, corroborating that CFC is not the main reason for the quantitative differences between the SPH and ET results in Figure 6.

The radial profiles of the angular velocity, measured at the time of mapping from the 3D merger models to the 2D postmerger models, t_pm = t_map = 10 ms, as well as for various later times, are shown in Figure 7 for three models. The inner core of the merger remnant is already rotating nearly uniformly at the time of mapping, while the profile in the surrounding disk corresponds to Keplerian rotation (∝r^−3/2). These characteristics remain essentially unchanged throughout the entire evolution of the NS remnant. We note that other results in the literature exist that report somewhat different behavior shortly after the merger, namely, a combination of a slowly rotating inner and fast-rotating outer core (e.g., Hanauske et al. 2017), or a double-core structure, surviving for a significantly longer time (Lioutas et al. 2022). We suspect that these differences are connected to the different numerical discretization schemes adopted by the aforementioned models (SPH versus Cartesian grid versus moving mesh, respectively); however, a detailed understanding of these differences has yet to be obtained.

Zoom In Zoom Out Reset image size

When constructing the outflow trajectories via backward time integration, special care must be taken to ensure that the distribution of thermodynamic properties (most importantly, Y_e ) in the dynamical ejecta remains consistent with that of the original merger simulations. This is because variations of fluid properties on small spatial scales or along the azimuthal direction get averaged out by the mapping from the 3D SPH configuration to the 2D grid at t_map = 10 ms. In order to approximately retain the Y_e pattern of the SPH simulations, we construct the dynamical ejecta trajectories as follows. After backward integration from t_pm = 100 s to t_map, all trajectories fulfilling r(t_map) > 250 km are split into five trajectories that differ only in their mass and Y_e . The mass and Y_e values for these five trajectories are taken from the five SPH particles of the corresponding merger model with the closest locations to the original postmerger trajectory at that time. The masses of these SPH particles are normalized such that their sum equals the mass of the postmerger trajectory. In the case that for these trajectories, the temperature already dropped below 10 GK at t_map, their expansion history at earlier times is taken directly and consistently from the SPH simulation. The resulting Y_e histograms and abundance yields are compared with those of the original SPH data in Figure 8. We find overall good consistency between both data sets but also noticeable differences for postmerger models based on the same SPH model, particularly at low Y_e values (e.g., the Y_e ≈ 0.05 peak in the asymmetric models). However, these differences are not necessarily connected only to sampling errors (i.e., errors related to the aforementioned mapping at t_map, as well as to the finite number of postmerger tracers) but could also be caused to some extent by different late-time (t_pm > t_map) behavior. While the postmerger simulations capture the hydrodynamic evolution of the ejecta far beyond t_map—including effects such as fallback or interaction with other ejecta components, which can all be sensitive to the viscosity—the plotted SPH data assume that spherical adiabatic expansion holds subsequent to t_map. Nevertheless, despite the approximate mapping and different assumptions at late times, the global pattern and most relevant features of the nucleosynthesis yields agree very well.

Zoom In Zoom Out Reset image size

In the main text, we argue that the polar outflow observed before BH formation is driven by neutrino heating; however, we do not explicitly back this statement. We also do not discuss the role of neutrino pair annihilation in driving this outflow. In order to briefly address these aspects, we ran two additional simulations similar to model sym-n1-a6: one in which only heating due to neutrino pair annihilation is ignored and one in which heating due to neutrino-nucleon absorption is also neglected. As for the technical implementation of these modifications, at each integration step, we first compute all source terms as usual but then set to zero all source terms corresponding to the aforementioned neutrino interactions (in both the hydro and moment equations) in regions where the density is smaller than 10¹¹ g cm⁻³ and neutrino interactions would otherwise heat up the fluid. The results are shown in Figure 9, which depicts the mass of all material with velocities v/c > 0.1 ejected beyond r = 3000 km, as well as the mass–velocity distribution of the same material at t_pm = 0.4 s, late enough to capture all fast ejecta from the NS remnant that collapses at t_BH = 122 ms. Without any neutrino heating (purple lines), the mass ejected within the first few hundred milliseconds is about five times smaller than in the unmodified model sym-n1-a6 (black lines) and corresponds to just about the mass of the dynamical ejecta with v/c > 0.1, demonstrating that neutrino heating is indeed the main driver of the fast polar outflow. The relative impact of pair annihilation can be assessed when comparing with the orange lines, which reveal that the ejecta exhibit a slightly less extended high-velocity tail, reaching only up to v/c ≈ 0.45 instead of 0.6, when pair annihilation is not taken into account. However, the total ejecta mass is only reduced by about 4 × 10⁻³ M_⊙ (corresponding to ≈5% of the total ejecta mass for this model), suggesting that pair annihilation has only a small impact on the r-process and KN-related features of the NS torus ejecta in our models.

Zoom In Zoom Out Reset image size