Ab-initio General-relativistic Neutrino-radiation Hydrodynamics Simulations of Long-lived Neutron Star Merger Remnants to Neutrino Cooling Timescales

We perform the first 3D ab-initio general-relativistic neutrino-radiation hydrodynamics of a long-lived neutron star merger remnant spanning a fraction of its cooling timescale. We find that neutrino cooling becomes the dominant energy loss mechanism after the gravitational-wave dominated phase (∼20 ms postmerger). Electron flavor antineutrino luminosity dominates over electron flavor neutrino luminosity at early times, resulting in a secular increase of the electron fraction in the outer layers of the remnant. However, the two luminosities become comparable ∼20–40 ms postmerger. A dense gas of electron antineutrinos is formed in the outer core of the remnant at densities ∼1014.5 g cm−3, corresponding to temperature hot spots. The neutrinos account for ∼10% of the lepton number in this region. Despite the negative radial temperature gradient, the radial entropy gradient remains positive, and the remnant is stably stratified according to the Ledoux criterion for convection. A massive accretion disk is formed from the material squeezed out of the collisional interface between the stars. The disk carries a large fraction of the angular momentum of the system, allowing the remnant massive neutron star to settle to a quasi-steady equilibrium within the region of possible, stable, rigidly rotating configurations. The remnant is differentially rotating, but it is stable against the magnetorotational instability. Other MHD mechanisms operating on longer timescales are likely responsible for the removal of the differential rotation. Our results indicate the remnant massive neutron star is thus qualitatively different from a protoneutron stars formed in core-collapse supernovae.


INTRODUCTION
Binary neutron star (NS) mergers can result in a variety of outcomes ranging from prompt black hole (BH) formation to the formation of absolutely stable remnants, depending on the component masses and spins and on the poorly known equation of state (EOS) of NSs (Shibata et al. 2005;Hotokezaka et al. 2011;Bauswein et al. 2013;Agathos et al. 2020;Köppel et al. 2019;Bauswein et al. 2020Bauswein et al. , 2021;;Perego et al. 2022;Kashyap et al. 2022;Kölsch et al. 2022;C ¸okluk et al. 2023).The evolution of remnant massive neutron star (RMNS) that do not collapse promptly to BH is governed by david.radice@psu.edu* Alfred P. Sloan fellow gravitational wave (GW) and neutrino energy and angular momentum losses and by angular momentum redistribution due to turbulence (Radice et al. 2020;Bernuzzi 2020).GW losses are dominant in the first ∼10−20 ms of the merger (Bernuzzi et al. 2016;Zappa et al. 2018).Considerable effort has been dedicated to quantifying the GW signal from merger remnants (Shibata 2005;Hotokezaka et al. 2011;Bauswein et al. 2013;Takami et al. 2014;Bernuzzi et al. 2015;Bose et al. 2018;Most et al. 2019;Breschi et al. 2022b;Wijngaarden et al. 2022;Breschi et al. 2022a;Espino et al. 2023;Fields et al. 2023;Dhani et al. 2023), because this is one of the key targets for next-generation GW experiments Cosmic Explorer (Reitze et al. 2019), Einstein Telescope (Punturo et al. 2010), and NEMO (Ackley et al. 2020).In some cases, the loss of angular momentum support due to GW emission can drive the RMNS to collapse during this phase.We refer to this outcome as a short-lived remnant (Baumgarte et al. 2000;Rosswog & Davies 2002;Shibata & Taniguchi 2006;Baiotti et al. 2008;Sekiguchi et al. 2011;Palenzuela et al. 2015;Kiuchi et al. 2022).If this does not happen, i.e., for long-lived remnants, then the dynamics transitions to being dominated by neutrino losses and turbulence (Hotokezaka et al. 2013;Kiuchi et al. 2018;Radice 2017;Radice et al. 2018a;Palenzuela et al. 2022a,b).
It is expected that a significant fraction of binary NS mergers will result in the formation of long-lived RMNSs (Piro et al. 2017;Radice et al. 2018a).Such remnants have been invoked to explain late time X-ray tails observed in a sizable fraction of short gamma-ray bursts (SGRBs; Dai & Lu 1998a,b;Zhang & Meszaros 2001;Dai et al. 2006;Metzger et al. 2008;Coward et al. 2012;Bucciantini et al. 2012;Rowlinson et al. 2013;Metzger & Piro 2014;Gao et al. 2017;Geng et al. 2016;Murase et al. 2018;Sarin et al. 2022;Dimple et al. 2023;Yuan et al. 2023), however searches of radio flares powered by long-term energy injection of a central engine in the ejecta following SGRB have not yet identified a candidate (Schroeder et al. 2020;Ghosh et al. 2022;Eddins et al. 2023); see also Lu & Quataert (2022) for an alternative explanation.Models of the "blue" kilonova in GW170817 also suggest that the remnant was stable for a timescale of at least few tens of milliseconds (Rezzolla et al. 2018;Li et al. 2018;Metzger et al. 2018;Ai et al. 2018;Nedora et al. 2019;Ciolfi & Kalinani 2020;Mösta et al. 2020;Combi & Siegel 2023;Curtis et al. 2023;Kawaguchi et al. 2023).
Here, we present the first ab-initio neutrino-radiation hydrodynamics simulations of a long-lived NS merger remnant extending to more than 100 ms after the merger.We quantify for the first time the evolution of angular momentum in the RMNS and in the disk, we discuss the secular evolution of GW and neutrino luminosities, and the internal structure of the RMNS.A brief description of the simulation and analysis methodology is given in Sec. 2 and our results are presented in detail in Sec. 3. Finally, we summarize and discuss our findings in Sec. 4. The appendices present additional technical details and analysis of the simulations.
We perform simulations using the general-relativistic large-eddy simulation (GRLES) formalism (Radice 2017(Radice , 2020) ) to account for angular momentum transport due to magnetohydrodynamics (MHD) effects in the remnant.Within this formalism, turbulent viscosity is parametrized in terms of a characteristic length scale ℓ mix .We perform simulations with ℓ mix = 0 (no viscosity) and with a new density dependent prescription for ℓ mix (ρ) calibrated using data from high-resolution GRMHD simulations of Kiuchi et al. (2018Kiuchi et al. ( , 2022)).More details are given in Appendix A. The effective viscosity used in our simulations is calibrated to high resolution simulations that span up to two seconds after merger, so it is expected to be realistic.However.we cannot exclude that GRMHD effects not captured by the viscous prescription could also influence the dynamics.That said, at least in the case of postmerger accretion disks around BHs there are indications suggesting that the viscous prescription is qualitatively and quantitatively accurate (Fernández et al. 2019).The key conclusions of this study appear to be independent on ℓ mix , so we discuss primarily the ℓ mix = 0 case in the main text, leaving a detailed comparison of the two ℓ mix prescriptions to Appendix D.
The evolution grid employs 7 levels of AMR.Each model is simulated 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.In total, these 4 simulations required more than 10 million CPU-core hours.

Timescales
Our simulations cover the last ∼6 orbits prior to merger and extend to more than ∼100 milliseconds after merger1 We define the time of merger t mrg as the time at which the amplitude of the GW (ℓ, m) = (2, 2) mode peaks, e.g.(Bernuzzi et al. 2014).Up to merger and shortly afterwards, the dynamics is driven by GW losses, while neutrinos become dominant on longer timescales.To quantify this, we estimate the GW and neutrino timescales as where L GW and L ν are the total (all modes and all flavors) GW and neutrino luminosities, respectively.We have chosen 0.1 M ⊙ as the typical energy scale, since this is a typical value of the energy radiated in GWs and neutrinos after the merger (Zappa et al. 2018;Radice et al. 2018a).An alternative way to estimate the GW luminosity is through the angular momentum as J/ JGW (Radice et al. 2018a).This alternative method estimates values that are broadly consistent with τ GW .The timescales defined in Eq. ( 1) are shown in Fig. 1.Shortly after merger τ GW ≃ 0.01 s, while τ ν ≃ 1 s.However, τ GW grows rapidly, because the remnant relaxes to a nearly axisymmetric configuration.The neutrino luminosities decay more slowly, so 10−20 ms after merger neutrinos they become the dominant mechanism through which energy is lost by the RMNS.At this moment, both GW and neutrino timescales significantly exceed the dynamical timescale of the remnant τ dyn ∼ 1 ms.By ∼50 ms after merger we can consider the dynamical phase of evolution of the remnant to be concluded (see also Appendix B).

Neutrino and Gravitational-Wave Luminosities
Figure 2 shows the GW luminosity for the three dominant modes, (ℓ, m) = (2, 1), (2, 2), and (3, 3), to- Gravitational wave and neutrino cooling timescales.The data is smoothed using a running average with window width ∆t = 0.5 ms.Neutrino radiation becomes the dominant mechanism for the evolution of the remnant ∼10 ms after merger.
gether with neutrino luminosities and average energies.The GW luminosity peaks at merger and then decays rapidly, particularly for the simulations with ℓ mix (ρ).The (ℓ, m) = (2, 2) mode remains the largest throughout the evolution, even though it is initially decaying more rapidly than the other two.At late times, its luminosity becomes comparable to those of the (ℓ, m) = (3, 3) and (ℓ, m) = (2, 1) modes.
The νe luminosity is the largest in the first ∼20 ms of the merger.This is due to the fact that the betaequilibrium shifts to higher Y e at the extreme densities attained in the postmerger (Cusinato et al. 2021).As a result, we observe the electron fraction close to the surface of the RMNS, which we define to be ρ = 10 13 g cm −3 , to grow from ∼0.05 to ∼0.07 within the first ∼100 ms of postmerger evolution (Fig. 3).The Y e grows to 0.09 at the inner edge of the accretion disk.Eventually, at t − t mrg ≃ 20−40 ms, the ν e luminosity becomes comparable to the νe luminosity.The average neutrino energies always satisfy ⟨ϵ νe ⟩ < ⟨ϵ νe ⟩ < ⟨ϵ νµ ⟩, as a result of the fact that the ν e 's decouple from matter at lower densities and temperatures compared to the other two species (Endrizzi et al. 2020).

Internal Stratification of the RMNS
Figure 3 shows azimuthally-averaged space and time diagrams (as a function of time and cylindrical radius) for various thermodynamic quantities on the equatorial plane for the ℓ mix = 0 SR simulation.Appendix C discusses a similar analysis performed along the z-axis, while Appendix D discusses the results for the ℓ mix (ρ) simulation.We find that the material heated up to high temperatures during the merger settles in a hot annular  (2, 1), (2, 2), and (3, 3).Lν µ is the luminosity of one of the heavy-lepton neutrino species.We do not track separately νµ, ντ and their antiparticles, so the overall neutrino luminosity of the four species combined is 4 × Lν µ .
region with ρ ≃ 10 14.5 g cm −3 , as also previously found by Bernuzzi et al. (2016) and Kastaun et al. (2016).This region is rich of νe , which contribute ∼10% of the overall lepton number in this region.The inner core and the mantle of the RMNS remain at lower temperatures, creating a negative radial temperature gradient.Nevertheless, the radial entropy gradient is positive at all times, stabilizing the star against convection and we find no evidences of radial overturn in the RMNS.
To quantify the stability to convection we compute where ρ is the rest-mass density, p is the pressure, s is the entropy, Y e is the proton fraction, and ϱ is the cylindrical radius.The first term in parenthesis in Eq. ( 2) is ob-tained from the EOS.The factor ρ 0 = 2.7×10 1 4 g cm −3 is introduced so that χ has dimension of an inverse length scale.If χ > 0, then fluid elements displaced adiabatically upwards are more dense than their environment and will not become buoyant.In other words, the fluid is stable against convection if χ > 0, this is the so-called Ledoux condition for instability2 .χ, also shown in Fig. 3, is always positive, with the exception of a few regions around the time of merger (highlighted in red).However, we remark that the Ledoux instability condition only applies when the background flow is stationary, so it is not applicable during the most dynamical phase of evolution around t ≃ t mrg .The angular velocity increases radially within the star and peaks at the surface, as also found in Kastaun & Galeazzi (2015) and Hanauske et al. (2017).This implies that the remnant is not only stable against convection, but also against the magnetorotational instability (MRI).

Angular Momentum Balance
Figure 4 shows the evolution of the ℓ mix = 0 SR binary in the total angular momentum J and baryonic mass M 0 plane (see Appendix D for the same analysis in the ℓ mix (ρ) case).We estimate the angular momentum of the binary in two ways.First, as J ADM − J GW (green dashed line in the figure), where J ADM is the ADM angular momentum of the system (which is constant in time), and J GW is the angular momentum carried away by GWs.Second, as the volume integral of the stress energy tensor (red line), T µν n µ ϕ ν where n µ is the normal to the t = const hypersurface and ϕ µ is the generator of the rotations in the orbital plane, in a cylinder of (coordinate) radius R = 300 G/c 2 M ⊙ ≃ 443 km.This second estimate is strictly only valid when the spacetime is stationary and axisymmetric, as such the dynamics close to merger, when a kink appear in the angular momentum of the RMNS might not be physically meaningful.Nevertheless, the two estimates agree well even at merger.At later times, angular momentum evolution is determined by mass ejection and disk formation, while JGW ≃ 0. The shaded region in Fig. 4 shows the allowed range for rigidly rotating equilibrium configurations.We remark that R = 300 G/c 2 M ⊙ is the approximate location at which the energy due to nuclear recombination of the disk material is approximately equal to the gravitational potential energy, so matter can be considered to be unbound past this point (Radice et al. 2018a).The purple contour denotes ρ = 10 13 g • cm −3 , while the green contours are ρ = 10 13.5 , 10 14 , and 10 14.5 g • cm −3 .
3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 Every 10 ms RNS Disk + RMNS RMNS only . Baryonic mass and angular momentum at the end of the simulations for the SR binary.The grey shaded area shows the set of all rigidly-rotating equilibrium configurations.The green line shows the evolution of the angular momentum due to gravitational wave emission, while the red and blue line show the evolution of the baryonic mass and angular momentum for the disk and the remnant, respectively.The remnant settles to a region corresponding to a stable rigidly-rotating configuration, but it is still rotating differentially at the end of our simulations We see that binary moves on an horizontal line during the inspiral and through merger, as GW carry away angular momentum, but not rest mass.Angular momentum loss is dominated by GW emission for ∼20 ms after the merger (see also Fig. 1).Within this time, the two estimates for the angular momentum agree well.On longer timescales angular momentum and total mass are carried away by neutrino-and spiral-wave-driven outflows (Nedora et al. 2019(Nedora et al. , 2021)), so the GW estimate is no longer reliable and the red line in Fig. 4 provides the correct representation of the state of the system.
The blue line in Fig. 4 shows the angular momentum and rest-mass pertaining to the inner core of the RMNS (ρ > 10 13 g cm −3 ).The difference between the blue and the red line is the total mass and angular momentum of the disk.The latter is formed of material squeezed out of the collisional interface between the two cores of the NSs (Radice et al. 2018b;Zenati et al. 2023).The mass of the disk saturates within ∼20 ms of the merger at 0.25 M ⊙ .The disk contains a substantial fraction of the angular momentum of the remnant, allowing the RMNS to settle to a quasi-steady state within the range of allowed equilibria.We remark that the RMNS is still differentially rotating at this stage (Fig. 3), so the ultimate outcome of the evolution of this binary will depend on the interplay between accretion, driven by MHD torques, and dissipation of the differential rotation within the remnant (Radice et al. 2018a;Just et al. 2023).

DISCUSSION
We have performed long-term, neutrino-radiation general-relativistic hydrodynamics (GRHD) simulations of a binary NS merger producing a long-lived RMNS.Simulations were performed at two grid resolutions and both with and without a treatment for subgrid scale angular momentum transport due to turbulent viscosity.While there are some quantitative differences between the results found in our simulations, the trends discussed here are robust.
During the inspiral and in the first ∼10−20 ms after merger, GW emission drives the dynamics of the system.However, neutrino radiation becomes the dominant energy-loss mechanism ∼20 ms after merger.We find no evidence for a revival of the GW signal due to convective instabilities.In contrast, we show that the remnant is stable against convection, in tension with the claims of De Pietri et al. (2020).The GW luminosity is sensitive to resolutions and viscosity, while neutrino radiation is robustly captured in our simulations, within the approximations inherent to the gray-M1 scheme we are using.This highlights the importance of neutrinos and resolution for quantitative predictions of the postmerger dynamics, as anticipated by Zappa et al. (2022).
Material that has been heated up by stirring during the merger forms a hot annular region in the outer core of the RMNS, at densities ρ ≃ 10 14.5 g cm −3 .Our simulations did not include MHD-effects explicitly, other than through the GRLES formalism.However, this is the region in which strong, small-scale magnetic fields generated by the Kelvin-Helmholtz instability (KHI) are expected to settle in the aftermath of the merger (Kiuchi et al. 2014;Combi & Siegel 2023).Along the rotational axis, this corresponds to a depth of 1 km below the surface.If RMNS are a viable central engine for SGRB, then the field need somehow to bubble out of the remnant and form large scale magnetic structures (Mösta et al. 2020).Recent GRMHD simulations suggest that it is possible for the field to break out of the star (Most & Quataert 2023;Combi & Siegel 2023).However, our simulations indicate that the RMNS is stably stratified, so it remains unclear how the magnetic fields can emerge from it.On the other hand, the simulations of Most & Quataert (2023) and Combi & Siegel (2023) used grid resolutions of 250 m and 180 m, respectively, which are likely insufficient to resolve the dynamics in this layer.Understanding whether or not, and on which timescale the magnetic fields amplified by the KHI can emerge from the RMNS will require high-resolution GRMHD simulations with realistic microphysics.
Trapped neutrinos also play an important role in the dynamics of material in the hot layers, since they ac-count for about 10% of the lepton number.The antineutrinos are formed as the result of the conversion of degenerate neutrons into protons.This process relieves some of the neutron degeneracy and results in an overall decrease of the pressure of matter (Perego et al. 2019).Such effect has been neglected in previous simulations that use leakage-based approaches, which assume nearly-frozen composition at high-optical depths.However, it might have implications for the stability of remnants closer to the BH formation threshold (Radice et al. 2018a;Zappa et al. 2022).
A massive accretion disk is formed by the ejection of material squeezed out of the collisional interface between the two stars forming a massive disk in the first ∼20 ms after merger.The disk carries a significant fraction of the angular momentum of the binary.As the disk forms, the inner part of the remnant settles to a quasi-steady state that corresponds to a stable uniformly rotating equilibrium configuration.However, differential rotation persists throughout our simulations, even when we include a treatment of subgrid angular momentum transport calibrated against very high-resolution GRMHD simulations.Given that the remnant is stable against both convection and the MRI, it is not clear what physical mechanism will operate to bring it to solid body rotation.This could occur as a result of magnetic winding (Baumgarte et al. 2000), possible doubly diffusive instabilities (Bruenn & Dineva 1996), or other GRMHD instabilities, such as the Tayler (Tayler 1973;Markey & Tayler 1973;Wright 1973;Lasky et al. 2011;Ciolfi et al. 2011;Sur et al. 2022), or Tayler-Spruit dynamo (Fuller et al. 2019;Margalit et al. 2022).Addressing this question will the object of future work.Within the GRLES formalism, angular momentum transport in the remnant is described as an effective turbulent subgrid stress tensor (Radice 2017(Radice , 2020)): where D i is the covariant derivative consistent with the spatial metric γ ij , e is the energy density, p, is the pressure, v i is the velocity, W is the Lorentz factor, and ν T is turbulent viscosity, which we write in terms of the mixing length ℓ mix and the sound speed c s as ν T = ℓ mix c s .
See Duez et al. (2020) for a discussion of possible alternative formulations.In the context of accretion disk theory, turbulent viscosity is typically parametrized in terms of a dimensionless constant α linked to ℓ mix through the relation where Ω is the angular velocity of the fluid (Shakura & Sunyaev 1973;Pringle 1981).Kiuchi et al. (2018Kiuchi et al. ( , 2022) ) with angular velocities and sound speed values of Radice (2020).The solid line is the fitting function adopted here.
Table 1.Values of the fitting coefficients used in Eq.A1 in units in which To estimate ℓ mix we take the values of α as a function of rest mass density ρ reported by Kiuchi et al. (2018Kiuchi et al. ( , 2022) and values of Ω and c s from Radice (2020).These values are reported in Fig. 5.We construct piecewise linear fit for ℓ mix in the form where ξ = log 10 ρ.The fitting coefficients are given in Tab. 1 and the results of the fit are shown in Fig. 5.Note that, because the piecewise linear fit would predict divergent ℓ mix at low densities, we need to limit the value of ℓ mix at low densities ρ < 10 ξ0 = 10 −10 c 6 G −3 M −2 ⊙ ≃ 6 × 10 7 g cm −3 , as shown in Eq. (A1).We remark that the simulations of Kiuchi et al. (2018) employed extremely high resolution of ∆x = 12.5 m and extended for over 30 ms.They also employed a very large initial magnetic field (∼10 15 G).As such they provide a realistic upper limit for the viscosity inside the remnant.The more recent simulations by Kiuchi et al. (2022) not only employ very high resolution, but also extend for over two seconds.As such, they provide realistic estimates for the viscosity also in the disk.Moreover, Fernández et al. (2019) found excellent agreement between GRMHD and alphaviscosity simulations of postmerger accretion disks around BHs.As such, we expect that realistic merger evolution should lay in between the predictions of our simulations with and without GRLES.

B. DIAGNOSTICS
Figure 6 shows the time evolution of the minimum of the lapse function, the maximum rest-mass density, and the maximum temperature in our simulations.The remnant is heated to temperatures of up to 70 MeV during merger.As material initially in the "hot spots" formed during merger is mixed with the rest of the NS matter the maximum temperature decreases (Bernuzzi et al. 2016;Kastaun et al. 2016).The RMNS also undergoes violent oscillations in the first milliseconds of its formation, seen in both the maximum density and central lapse values.However, the most dynamical phase terminates ∼50 ms after the merger, after which the oscillations have subsided and the temperature stops evolving on a dynamical timescale.
C. VERTICAL STRATIFICATION Figure 7 shows space and time profiles of various thermodynamic quantities along the rotational axis of the ℓ mix = 0 SR model.This figure should be contrasted with Fig. 3, which shows the same quantities in the equatorial plane.The main difference is that the density scale-height is smaller in the vertical direction, due to the absence of centrifugal support in the z-direction.As in the equatorial plane, the peak of the temperature is at ρ ≃ 10 14.5 g cm −3 , although the maximum temperature in the vertical direction (∼25 MeV) is somewhat smaller than the peak temperature on the xy-plane (∼30 MeV).This is also the region at which Y νe peaks, Y νe is about a factor two smaller than in the equatorial plane.The electron fraction Y e slowly increases with time in the regions close to the stellar surface, with Y e going from 0.06 to 0.08 at ρ ≃ 10 13 g cm −3 .Finally, we compute the Ledoux instability criterion χ using the data along the z-axis and find that the RMNS is stable against convection along the rotational axis.

D. EFFECT OF VISCOSITY
Figure 8 shows azimuthally averaged space and time profiles for several thermodynamic quantities in the equatorial for the ℓ mix (ρ) SR binary.This figure should be contrasted with Fig. 3, which shows the same quantities for the ℓ mix = 0 binary.The main effect of the viscosity is to suppress small scale oscillations of the remnant, as can be seen The purple contour denotes ρ = 10 13 g • cm −3 , while the green contours are ρ = 10 13.5 , 10 14 , and 10 14.5 g • cm −3 .This figure should be contrasted with Fig. 3, which shows the same quantities in the equatorial plane.
from the lack of oscillations in the ρ = 10 13 g cm −3 surface position (marked in purple in the figure), and to increase the entropy in the outer core of the RMNS.Despite these changes, the viscous model remains stable against both convection and the MRI.In fact, the increase in the radial entropy gradient due to viscosity has a stabilizing effect.
Figure 9 shows the evolution of the ℓ mix (ρ) SR binary in the J−M 0 plane.It should be contrasted to Fig. 4, which shows the same analysis for the ℓ mix = 0 SR binary.There are some quantitative differences between the two binaries, most notably the disk mass for the ℓ mix (ρ) binary is slightly smaller than that of the ℓ mix = 0 binary (0.2 M ⊙ vs 0.25 M ⊙ ).The mass loss rate from the disk is also slightly enhanced with viscosity, however this difference is within the numerical uncertainty we estimate for our simulations.That said, the evolution of the ℓ mix (ρ) binary is qualitatively identical to that of the ℓ mix = 0 binary.In both cases, the RMNS settles to a quasi-steady configuration within the phase-space region allowed for rigidly rotating equilibria, while disk and ejecta take the excess angular momentum. .Baryonic mass and angular momentum at the end of the simulations for the ℓmix(ρ) SR binary.The grey shaded area shows the set of all rigidly-rotating equilibrium configurations.The green line shows the evolution of the angular momentum due to gravitational wave emission, while the red and blue line show the evolution of the Baryonic mass and angular momentum for the disk and the remnant, respectively.This figure should be contrasted with Fig. 4, which shows the same quantities for the ℓmix = 0 binary.

Figure 2 .
Figure2.GW luminosity and neutrino luminosity and average energies.The data is smoothed using a running average with window width ∆t = 0.5 ms.The GW luminosity is shown for the three most dominant modes: (ℓ, m) = (2, 1), (2, 2), and (3, 3).Lν µ is the luminosity of one of the heavy-lepton neutrino species.We do not track separately νµ, ντ and their antiparticles, so the overall neutrino luminosity of the four species combined is 4 × Lν µ .

Figure 3 .
Figure 3. Angularly averaged profiles of entropy per baryon s, electron fraction Ye, convection instability criterion χ, neutrino fractions Yν e , Yν e , and Yν µ , and angular frequency Ω on the equatorial plane for the ℓmix = 0 SR binary.The profiles are shown as a function of cylindrical radius ϱ and time from merger t − tmrg.The purple contour denotes ρ = 10 13 g • cm −3 , while the green contours are ρ = 10 13.5 , 10 14 , and 10 14.5 g • cm −3 .

Figure 5 .
Figure5.Mixing length prescription used for the viscous simulations.The data points are obtained by combining the values of α reported byKiuchi et al. (2018Kiuchi et al. ( , 2022) ) with angular velocities and sound speed values ofRadice (2020).The solid line is the fitting function adopted here.

Figure 6 .
Figure 6.Time evolution of diagnostic quantities in the simulations: minimum of the lapse function, αmin, maximum rest-mass density, ρmax, and maximum temperature Tmax.The dynamical phase of evolution terminates ∼50 ms after merger.

Figure 7 .
Figure 7. Profiles of entropy per baryon s, electron fraction Ye, convection instability criterion χ, and neutrino fractions Yν e , Yν e , and Yν µ along the z-axis for the ℓmix = 0 SR binary.The profiles are shown as a function of time from merger t − tmrg.The purple contour denotes ρ = 10 13 g • cm −3 , while the green contours are ρ = 10 13.5 , 10 14 , and 10 14.5 g • cm −3 .This figure should be contrasted with Fig.3, which shows the same quantities in the equatorial plane.
Figure9.Baryonic mass and angular momentum at the end of the simulations for the ℓmix(ρ) SR binary.The grey shaded area shows the set of all rigidly-rotating equilibrium configurations.The green line shows the evolution of the angular momentum due to gravitational wave emission, while the red and blue line show the evolution of the Baryonic mass and angular momentum for the disk and the remnant, respectively.This figure should be contrasted with Fig.4, which shows the same quantities for the ℓmix = 0 binary.