Emergence of Microphysical Bulk Viscosity in Binary Neutron Star Postmerger Dynamics

In nuclear matter in isolated neutron stars, the flavor content (e.g., proton fraction) is subject to weak interactions, establishing flavor (β-)equilibrium. However, there can be deviations from this equilibrium during the merger of two neutron stars. We study the resulting out-of-equilibrium dynamics during the collision by incorporating direct and modified Urca processes (in the neutrino-transparent regime) into general-relativistic hydrodynamics simulations with a simplified neutrino transport scheme. We demonstrate how weak-interaction-driven bulk viscosity in postmerger simulations can emerge and assess the bulk viscous dynamics of the resulting flow. We further place limits on the impact of the postmerger gravitational-wave strain. Our results show that weak-interaction-driven bulk viscosity can potentially lead to a phase shift of the postmerger gravitational-wave spectrum, although the effect is currently on the same level as the numerical errors of our simulation.

There are various uncertainties in merger simulations, such as the precise form of the cold EoS, finite-temperature effects (Bauswein et al. 2010;Figura et al. 2020;Raithel et al. 2021), and the appearance of exotic degrees of freedom or phase transitions (Oechslin et al. 2004;Most et al. 2019a;Bauswein et al. 2019;Most et al. 2020a;Weih et al. 2020;Chatziioannou & Han 2020;Prakash et al. 2021;Tootle et al. 2022;Tan et al. 2022).In this work, we focus on a dynamical effect, the equilibration of the proton fraction via weak interactions (Alford et al. 2018;Radice et al. 2022), and its potential to alter the post-merger dynamics and systematically bias our ability to infer dense matter properties (Most et al. 2021).This is possible because the dynamical timescales of a merger can be fast enough to drive the system out of β-equilibrium (Hammond et al. 2021;Most & Raithel 2021).One expected consequence of the resultant re-equilibration via weak interac-tions is bulk viscosity (Alford & Harris 2019;Gavassino et al. 2021).
Previous work has largely been devoted to understanding the properties (Alford & Harris 2019;Alford et al. 2019) and the potential appearance of bulk viscosity (Alford et al. 2018;Most et al. 2021;Hammond et al. 2021) in the remnant, or to study the idealized limit of perfect β−equilibrium (Ardevol-Pulpillo et al. 2019;Hammond et al. 2022).While effects of microphysical viscosity in neutron-star merger simulations have largely been ignored, they have been well-studied in the related context of dense matter formed heavy-ion collisions (Monnai & Hirano 2009;Song & Heinz 2010;Bozek 2010;Dusling & Schäfer 2012;Noronha-Hostler et al. 2013;Ryu et al. 2015Ryu et al. , 2018)), see also Romatschke & Romatschke (2019) for a review.
Going beyond all of the above, in this paper, we take essential steps towards fully self-consistent numerical simulations of the out-of-equilibrium dynamics of the flavor content of the matter in the post-merger phase of a neutron star collision.We begin by briefly reviewing the relevant processes (Urca processes) active in the core of neutron stars (Yakovlev & Levenfish 1995).We then provide a detailed discussion of the dynamics of relaxation to β−equilibrium after the merger.For the first time, we directly demonstrate the emergence of bulk viscosity (Gavassino & Noronha 2023) in post-merger simulations and provide a first quantification of its effects.We conclude by discussing how this may alter the post-merger gravitational wave emission.

METHODS
We model the global dynamics of the system using general-relativistic (magneto-)hydrodynamics (Duez et al. 2005) in a dynamically evolved spacetime in the Z4c formulation (Hilditch et al. 2013).We evolve the baryon current J µ B = n B u µ and the perfect-fluid energy-momentum tensor T µν (Denicol & Rischke 2021) according to where Q β ν is an effective energy "sink" term associated with neutrino energy-losses (Sekiguchi et al. 2012;Galeazzi et al. 2013).Here, n B and u µ are the baryon number density and fluid four-velocity, respectively.
Assuming npe matter, the electron fraction, Y e = n e /n B , can only change by weak-interaction decays of neutrons (n) and protons (p).We incorporate these as follows.First, we include a standard set of weak interactions without Urca contributions, R νe , commonly used in leakage scheme approaches to neutron star mergers (Ruffert et al. 1996;Rosswog & Liebendörfer 2003).These include approximate direct Urca reactions n + ν e ↔ p + e − (though not n → p + e − + νe which is relevant at low neutrino density) which we replace as outlined in the following), pair annihilation (e + + e − → νe + ν e ), plasmon decay γ → νe + ν e .Those rates are taken from Ruffert et al. (1996) and will, at low neutrino densities, not drive the system to flavor equilibrium.Secondly, we consistently include direct, Γ dUrca , and -for the first time-modified, Γ mUrca , Urca (net) rates (Yakovlev et al. 2001;Alford et al. 2021b) that operate in the dense nuclear matter core of the merger; see our discussion of the neutrino transparent regime below.To track the relative number of protons, we evolve the relative fraction of electrons (e) We solve the coupled system of equations using the Frankfurt/IllinoisGRMHD code (Most et al. 2019b;Etienne et al. 2015), which utilizes the EinsteinToolkit infrastructure (Loffler et al. 2012), where the Urca processes are included via a fully implicit operator-split approach.The FIL code solves the equations of general-relativistic magnetohydrodynamics in dynamical spacetimes Duez et al. (2005).Since magnetic fields do not affect chemical equilibration, we set the magnetic field to zero throughout the evolution, making our simulations fully hydrodynamic.The system is integrated numerically using a fourth-order accurate version of the ECHO scheme (Del Zanna et al. 2007).The spacetime is evolved in the Z4c formulation of Einstein's equations (Bernuzzi & Hilditch 2010;Hilditch et al. 2013), using a fourth-order accurate finite-difference discretization (Zlochower et al. 2005).We solve the equation on a set of eight nested grids, with a finest resolution of 262 m and an outer domain size of 3022 km.
In addition, we also perform simulations at 150 m resolution when quantifying the numerical error of our results.More details on the numerical setup and the choice of initial data can be found in (Most & Raithel 2021).
This first study considers a system with fixed binary parameters loosely modeled after the GW170817 event (Abbott et al. 2017).We fix the total mass M = 2.74 M ⊙ with a mass ratio of q = 0.85.This results in a hot, (meta-)stable neutron star remnant after the merger, which does not immediately collapse to a black hole (Bauswein et al. 2013a;Köppel et al. 2019;Bauswein et al. 2020;Tootle et al. 2021;Kashyap et al. 2022;Kölsch et al. 2021;Perego et al. 2021).A detailed description of the (post-)merger dynamics and gravitational wave emission of this system for different EoSs can be found in (Most & Raithel 2021;Raithel & Most 2022).This work focuses on the impact of departures from and relaxation back to beta equilibrium (Hammond et al. 2021;Most & Raithel 2021).
To approximate the complex dynamics of neutrino populations, one can divide the merger into two regions: the neutrino-trapped regime and the neutrino-transparent regime.In the former, local opacities are so high that the neutrinos cannot escape, forming a gas.This leads to a correction in the fluid element's overall pressure, energy density, and lepton number.Our work mainly focuses on the merger in the neutrino transparent regime, but in the Appendix we provide an assessment of neutrino trapping which will establish chemical equilibrium much faster than in the free streaming regime (Alford et al. 2019(Alford et al. , 2021a;;Perego et al. 2019).Accordingly, we treat neutrinos numerically in a simplified Leakage manner (Galeazzi et al. 2013), assuming the free-streaming loss of neutrinos from each fluid element.While this will not affect most of our qualitative conclusions on how microphysical dissipation emerges in the post-merger dynamics, it will change their quantitative impact, especially in macroscopic observables such as the gravitational wave emission.We estimate the effect of neutrino transparency using additional simulations in the fully neutrino trapped regime-see the Appendix for details.
In the neutrino transparent regime beta equilibration occurs via the direct Urca (dUrca), n → p + e − + νe , p + e − → n + ν e , and modified Urca (mUrca) processes, n + X → p + X + e − + νe , p + e − + X → n + X + ν e , where X is a spectator nucleon (Yakovlev & Levenfish 1995;Yakovlev et al. 2001).For this first investigation of β-equilibration in merger dynamics, we use simple analytic expressions for the Urca rates obtained via the Fermi surface approximation in which all participating nucleons and electrons are assumed to be close to their Fermi surfaces.Realistically, there are finite temperature corrections that are neglected here because their main effect is to shift the temperature at which the resonant bulk viscosity reaches its maximum value by a few MeV (Alford & Harris 2018;Alford et al. 2021b).
The neutrino-transparent regime of npe matter can then be divided into two sub-regimes, according to whether the proton fraction Y e is above or below the direct Urca threshold value of 0.11 (Yakovlev et al. 2001;Lattimer et al. 1991).Above this threshold, the faster direct Urca process dominates over modified Urca.Below the threshold, direct Urca is suppressed and beta equilibrium occurs via the slower modified Urca process.To probe the influence of the dUrca threshold, we consider two EoSs, TMA (Hempel & Schaffner-Bielich 2010) and SFHo (Steiner et al. 2013;Alford et al. 2021b).At β-equilibrium, TMA has a direct Urca threshold just below 2n sat (where n sat ≈ 0.16 fm −3 is nuclear saturation density) while in SFHo the proton fraction never reaches the direct Urca threshold (in both cases referring to β-equilibrated matter at T = 0).More information can be found in the Appendix.
In the Fermi surface approximation, one can characterize the degree to which neutrino-transparent material is out of beta equilibrium by the chemical potential imbalance δµ = µ n − µ p − µ e , where µ a is the chemical potential for particle a.In beta equilibrium δµ = 0 and Y e = Y eq (n, ε), where ε is the total fluid energy density.Given that δµ = δµ(ε, n B , Y e ), one can use ( 1) and ( 2) to find that δµ evolves according to the following equation is a thermodynamic quantity, with P = P (ε, n B , Y e ) being the total pressure.Also, denotes a source term coming from neutrino radiative energy loss.Equation ( 3) is a complicated, highly nonlinear1 relaxation-type equation that describes how δµ varies given the local fluid expansion rate and the source.We shall show in the following that during a neutron star merger, the full relaxation equation ( 3) experiences different regimes that are nothing but different manifestations of the physics of bulk viscosity.
For slight deviations from beta equilibrium, 3) enters the linear response regime where δµ relaxes to zero on a timescale τ 0 , and the relaxation dynamics can be understood in terms of bulk viscosity (Gavassino et al. 2021;Camelio et al. 2022a;Celora et al. 2022;Most et al. 2021).Indeed, one can define the bulk scalar pressure correction, Π, by splitting the total pressure, P , into where P eq = P (n B , ε, Y e = Y eq e ) is the pressure in beta equilibrium.In the linear response regime, one can approximate Π ∼ I 1 δµ (with ) to find (Gavassino et al. 2021) where ζ 0 = B 0 I 1 τ 0 is the static bulk viscosity coefficient induced by the weak-interactions (Sawyer 1989;Sa'd & Schaffner-Bielich 2009;Harris 2020).For periodic density oscillations with frequency ω, δµ and Π also oscillate and ( 5) defines an AC bulk viscosity et al. 2023), which reaches a resonant maximum when ω = 1/τ 0 (see, e.g., Harris (2020)).Eq. ( 5) is a typical Israel-Stewart-like (Israel & Stewart 1979) equation of motion for the bulk scalar (apart from the source term) (Gavassino et al. 2021), which is commonly used in the description of the quark-gluon plasma formed in heavy ion collisions ( τ 0 are determined by the strong, not the weak interactions).
Applications of semi-analytic expressions for the bulk viscosity (Alford & Harris 2019) to the background of a nonviscous neutron star merger calculation have projected that large viscous corrections could arise after the merger (Most et al. 2021).Here, we go beyond such estimates, as well as first-order bulk-viscous approximations (Gavassino et al. 2021;Celora et al. 2022;Camelio et al. 2022a), and instead compute the full dynamics of the Urca process (in the Fermi Surface approximation -see Appendix for an extended discussion) and its influence on dense matter beyond leading order in δµ by directly extracting Π from the total pressure (4) in our simulations (see also Ref. Camelio et al. (2022b,a) for similar simulations of oscillating neutron stars in spherical symmetry).We show that neutron star mergers can reach rapid equilibration and resonant bulk-viscous regimes.

RESULTS
While weak interactions have sufficient time to establish equilibrium for all but the late stages of inspiral (Arras & Weinberg 2019), the violent dynamics of the merger can drive the electron fraction away from its pre-merger equilibrium value (Perego et al. 2019;Hammond et al. 2021).This is displayed in Fig. 1, which shows that the distribution of Y e over baryon density 5 ms after merger is quite different when Urca reactions, which restore equilibrium on a time scale that depends on the density and temperature, are included.EoS, this affects the structure of the star since, in the presence of Urca processes, the matter reaches higher densities, n B > 2.5 n sat .In the case of SFHo, the merger produces some regions with an excess of neutrons and some with an excess of protons; for TMA there is mostly an excess of neutrons (see Fig. 1).If weak interactions are correctly included, this leads to a rapid onset of neutron ↔ proton conversion in the star, as shown by the difference between equilibration and non-equilibration Y e profiles in Fig. 1.In the case of SFHo, this will be driven mainly by mUrca, as most fluid elements have Y e values below the dUrca threshold (Fig. 1), whereas for TMA it is a combination of dUrca and mUrca.We can also anticipate that this change in structure and composition might lead to a change in the (early) mass ejection (e.g., Rosswog et al. 1999;Bauswein et al. 2013b;Sekiguchi et al. 2015Sekiguchi et al. , 2016;;Radice et al. 2016;Lehner et al. 2016;Bovard et al. 2017;Radice et al. 2018a;Nedora et al. 2019), but leave the details to future work.
We now move on to the main part of our work, which is to quantify the influence on the global dynamics of beta equilibration via Urca processes.
Using our approach, we can provide a first quantification of the out-of-equilibrium dynamics in terms of the relative bulk pressure contribution Π/P eq as shown in Fig. 2 for a fixed time.We can see that bulk viscous corrections in terms of Π appear in the remnant.We will quantify them in the following and analyze their dynamical impact.
We will also utilize the concept of a dynamic inverse Reynolds number (Denicol et al. 2012;Most et al. 2021), χ = Π/ (ε eq + P eq ), with ε eq being the total energy density of the fluid in β−equilibrium.We will, in particular, rely on a density-weighted average value ⟨χ⟩.In Fig. 3 we show the the departure from β-equilibrium δµ relative to the tem- perature: δµ/T ≪ 1 is the linear response regime.In the very late inspiral, right before the collision, tidal forces drive the material inside the stars out of β-equilibrium.While we only show the beginning of the collision (Fig. 3, left panel), this could potentially set in even earlier (Arras & Weinberg 2019).Overall, we find that δµ/T ≳ 0.1 in most of the stellar material (Fig. 3, left panel).The subsequent evolution now strongly depends on whether weak-interaction effects are included.If Urca processes are not included, the matter does not re-equilibrate on the dynamical time scale of the merger and the out-of-equilibrium pressure contribution, χ, remains constant on average for both EoSs (Fig. 4, dashed lines).If Urca processes are included the system begins to re-equilibrate at a rate that depends on density and temperature.Indeed, we can see (Fig. 4, solid lines) that with the inclusion of Urca processes the system relaxes to chemical equilibrium with an average relaxation time comparable to the millisecond timescale of the merger dynamics.If we included neutrino-trapped regions, these would equilibrate much faster (Perego et al. 2019).
During the initial phase of largest out-of-equilibrium dynamics, we see in Fig. 4 (left panel) that both systems reach χ ≃ 0.01, which is comparable to the effective viscous damping of density oscillations caused by post-merger grav-itational wave emission alone (Most et al. 2021).In other words, corrections to the macroscopic dynamics can be sizeable in the neutrino free-streaming limit.Furthermore, in both cases |δµ|/T ≫ 1, meaning that the system is in a nonlinear far-from-equilibrium regime2 where there is effectively an amplitude-dependent ("suprathermal") bulk viscosity (Madsen 1992;Reisenegger 1995;Celora et al. 2022;Camelio et al. 2022b,a).
However, post-merger oscillations (Stergioulas et al. 2011) continue to locally drive the system out of β−equilibrium.It is especially in this nontrivial far-from-equilibrium regime characterized by the approximate plateau in Fig. 4, u µ ∇ µ δµ ≈ 0 at |δµ|/T ≫ 1, where we expect bulk viscous damping to be most pronounced.As a consequence of Eq. ( 5), this will also drive an approximate steady state in χ with fluctuations, which is consistent with Fig. 4.
With this in mind, we can now correlate chemical equilibration and its hydrodynamic feedback in Fig. 4. The equili-bration behavior for the two models is qualitatively different: For TMA, which is expected to equilibrate faster via dUrca as well as mUrca, we see that ⟨χ⟩ briefly jumps up after the merger, following the no-Urca simulation, but then follows an approximately exponential decay with a lifetime of about 1 ms as Urca processes establish β−equilibrium.For SFHo, which is expected to equilibrate more slowly via mUrca processes, the behavior is more complicated: (a) ⟨χ⟩ almost immediately reaches a plateau at a value about 10 times higher than in the no-Urca simulations; (b) the plateau persists for about 4 ms before giving way to exponential decay with lifetime ∼ 1.7 ms.
As a direct consequence, the plateau can only be sustained when matter is driven out of equilibrium by bulk fluid motion at a rate comparable to the flavor relaxation time scale, τ , in this region.From our simulations using the SFHo EoS, we can directly correlate the appearance of these outof-equilibrium regions with the periodic compressive motion of the stellar cores (see also the regions with n B > 3n sat in Fig. 3, center panel).
The plateau regions end at t ≈ 5 ms, which is when the two cores have fused into a single core, and out-of-equilibrium matter is now found in lower density layers (n B = 2n sat , right panel Fig. 3).At t ≳ 5 ms, in the absence of global compressive and expansive motion from core bounces, δµ shows exponential equilibration akin to the TMA case.
For both EoSs, the flavor relaxation times are shorter than the timescale (≈ 20 − 30 ms) of post-merger gravitational wave emission.This implies that for the short-term evolution of the merger remnant (e.g., Fujibayashi et al. 2017Fujibayashi et al. , 2018)), the matter can be treated as being in chemical equilibrium.As the system approaches chemical equilibrium, δµ/T ≪ 1, at times t > 5 ms, post-merger oscillations of the stellar remnant only mildly drive the system out of equilibrium, so ⟨χ⟩ ≃ 10 −5 and δµ/T ≃ 10 −2 .Furthermore, we can qualitatively understand why SFHo matter shows a sustained far-from-equilibrium behavior compared to TMA.The merger drives both systems out of β-equilibrium, such that there are regions at early times where |δµ|/T ≥ 1.According to Eq. ( 3), a non-equilibrium steady-state for δµ/T can only appear in the regions where δµ+Bτ θ ≈ 0 (assuming nearly constant T ).This means that the appearance of such a state (and its duration) depends on whether Bτ can be locally sufficiently large.Given that the expansion rate θ of both systems is comparable, a reasonable explanation for the approximate plateau in SFHo is that it has a larger τ because slower modified Urca processes dominate the rates, while TMA reaches the direct Urca threshold, leading to faster equilibration.Consequently, this approximate cancellation effect that drives the non-equilibrium steady state is less prominent in TMA, and Eq.(3) predicts the nearly exponential behavior for this EoS seen in Fig. 4 in contrast to the approximate plateau seen in SFHo for a few milliseconds.
In the light of the relation between chemical reactions and bulk-viscous processes discussed in (Gavassino & Noronha 2023), one may say that SFHo matter in this regime will be in a resummed (i.e., amplitude-dependent) Navier-Stokeslike limit (Yang et al. 2023) allowing us to extract effective bulk viscosities with magnitude ζ eff ≲ 10 28 g/(cm s) directly from the simulation by computing |ζ eff | ≈ |Π/θ| ≈ |Π/∇ i u i |.We point out that these viscosities can also be calculated precisely in the far-from-equilibrium limit using the same microphysical reactions employed in this work (Yang et al. 2023).This appearance of dissipation also leads to the production of entropy, increasing it on average by ∆s ∼ 0.1k B /baryon compared to a merger without Urca processes (see also Most et al. (2022)).
Ideally, one would like to identify a signature of dissipative re-equilibration in the gravitational wave signal emitted after the collision.To this end, in Fig. 5 we compare the gravitational emission for all four simulations.For both EoSs, we find good agreement between the gravitational wave strains with and without Urca processes in the late inspiral and the early merger phase.This is consistent with the Urca processes not enforcing β-equilibrium immediately, as shown in Fig. 4.However, after the characteristic re-equilibration time ≃ 1 − 2 ms , a phase difference begins to build up and continues to grow until the end of the simulation t ≃ 10 ms.
A quantitative estimate of the exact magnitude of how chemical equilibration affects the gravitational wave emission from our simulation is subject to both effects of numerical resolution and our assumption of neutrino-transparency of the remnant.The latter stems from the fact that flavor equilibration is much slower in the free-streaming regime so there are larger departures from chemical equilibrium.
In order to provide a complete and faithful quantification of these errors, we focus on the TMA system, which featured the strongest out-of-equilibrium effects in the gravitational wave signal.We have performed two additional simulations at significantly increased numerical resolution (150 m), one of which imposed the neutrino trapped regime above temperatures of T > 1 MeV, allowing us to quantify both sources of error (see Appendix).Since the neutrino mean free path depends on the neutrino energy, a clear separation of the neutrino free streaming and neutrino trapped regime by a single transition temperature is not possible.Here, we err on the side of caution and choose a very low cutoff temperature to compute an upper bound on the effect of neutrino trapping.Based on these simulations, we now comment on potential phase differences between Urca, trapped and frozen composition models.
First, viscous processes are expected to damp out gravitational wave emission.Using the waveforms shown in the supplemental material, we observe a slight reduction in power for the Urca cases compared to neutrino-trapped and regular evolution cases, providing an indirect consistent indication for a bulk viscous damping of gravitational wave emission.Second, dissipation is expected to shift the dominant frequency of post-merger gravitational wave emission.Our simulations in the neutrino transparent regime (Urca vs. frozen composition) feature characteristic shifts of ∆f ≃ 40 Hz in both cases.These are roughly comparable to the effects of finite temperature (Raithel & Most 2022;Fields et al. 2023).However, we caution that these are smaller than the resolution error, which we quantify at 70 Hz and thus larger than the shift we measure.Third, the shift in gravitational wave frequency can also depend on the assumptions on neutrino transparency.When including neutrino trapping above T > 1 MeV we find an overall shift of about 50 Hz (again smaller than our numerical error budget): see also Hammond et al. (2022) for the use of infinitely fast equilibration rates.Additionally, numerical errors can move the shift ∆f upwards or downwards in these simulations.Ultimately, simulations at higher resolutions and with varying microphysics will be needed to provide a reliable quantification of this effect on the post-merger gravitational wave signal.

CONCLUSIONS
This work presents the first study of weak-interactiondriven out-of-equilibrium dynamics in a binary neutron star merger.We have shown how Urca processes acting on a millisecond timescale restore the departures from β−equilibrium that arise during the initial stages of the merger.We demonstrate the emergence of a microphysical bulk viscosity and explicitly demonstrate the dissipation induced by Urca processes, and the establishment of far-fromequilibrium conditions (Gavassino & Noronha 2023).We finally assess the impact on the gravitational wave signal and place upper bounds on the impact of viscosity and chemical equilibration (see also (Hammond et al. 2022)).This might change (upwards or downwards) if significantly finer grids were used (e.g., Breschi et al. 2019;Kiuchi et al. 2020).
Further work will be required to fully map out the impact of out-of-equilibrium dynamics on the post-merger evolution.Firstly, not all equations of state lead to strong violation of β−equilibrium at merger (Most & Raithel 2021;Hammond et al. 2021).Those that do will be subject to uncertainties in the finite-temperature part (Raithel et al. 2021), critically affecting what fraction of the merger will be in the neutrino trapped regime (Perego et al. 2019;Alford et al. 2021b;Perego et al. 2021), in which re-equilibration happens almost instantaneously (Alford et al. 2019(Alford et al. , 2021a)).Similarly, trapped neutrino pressure contributions might be of similar strength, but will affect hot instead of cold regions (Perego et al. 2019).Furthermore, it remains to be seen what impact the change in remnant structure has on the mass ejection and lifetime of the post-merger system (Gill et al. 2019), see also Zappa et al. (2022).Effects that are not captured by the Fermi Surface approximation, such as finite-temperature blurring of the direct Urca threshold and the resultant modification of the β−equilibrium condition (Alford et al. 2021b;Alford & Harris 2018), should be explored.Different phases of matter, such as hyperonic and quark matter, have different channels of chemical equilibration giving rise to bulk viscosities with different dependencies on temperature and density (Alford & Haber 2021;Schmitt & Shternin 2018).The presence of a bulk-viscous phase in the merger opens a new window towards connecting dense matter in neutron star mergers to that in heavy-ion collisions (Most et al. 2022), where out-ofequilibrium viscous corrections are relevant (Romatschke & Romatschke 2019).ACKNOWLEDGMENTS ERM thanks F. Foucart, J. Noronha-Hostler, A. Pandya, F. Pretorius, C. Raithel and J. Ripley for insightful discussions related to this work.All authors are grateful to M. Antonelli, S. Bernuzzi, G. Camelio, L. Gavassino, and B. Haskell for very helpful comments on the manuscript, and to L. Gavassino for pointing out the question of hydrodynamic frame choice associated with Eq. (4).ERM gratefully acknowledges support from a joint fellowship at the Princeton Center for Theoretical Science, the Princeton Gravity Initiative and the Institute for Advanced Study.ERM acknowledges support for compute time allocations on the NSF Frontera supercomputer under grants AST21006.This work used the Extreme Science and Engineering Discovery Environment (XSEDE) (Towns et al. 2014) through Expanse at SDSC and Bridges-2 at PSC through allocations PHY210053 and PHY210074.The simulations were also in part performed on computational resources managed and supported by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology's High Performance Computing Center and Visualization Laboratory at Princeton University.ERM also acknowledges the use of high-performance computing at the Institute for Advanced Study.ERM acknowledges partial support from the National Science Foundation through PHY- In this initial work, we use the Fermi surface approximation to calculate the Urca rates.In this assumption, the nuclear matter is treated as strongly degenerate Alford et al. (2021b); Alford & Harris (2018).The net rate of direct Urca is given by (see also Camelio et al. (2022b,a)) where G ≡ G F cos θ c , G F = 1.166 × 10 −11 MeV −2 is the Fermi constant, θ c = 13.04 • is the Cabbibo angle, g A = 1.26, and N are the nucleon Fermi energies, and δµ ≡ µ n − µ p − µ e quantifies the departure of the system from chemical equilibrium.For simplicity, we fix the effective mass m * N ≈ 0.7 m n relative to the neutron mass m n , although in reality the effective mass drops with increasing density Alford et al. (2021b).In the Fermi surface approximation, the direct Urca process operates only above a threshold density, where (2) We neglect finite temperature corrections to the Fermi surface approximation which would blur this threshold and lead to a nonzero δµ eq Alford & Harris (2018); Alford et al. (2021b).Above temperatures of T ≈ 1 − 10 MeV we expect neutrinos to be trapped Roberts & Reddy (2017); Alford & Harris (2018).Since the associated Urca timescales are significantly shorter than the time steps of our simulations Alford et al. (2021a), we capture this regime approximately by enforcing cold β-equilibrium.Due to the implicit time-stepping, all rates acting on time scales shorter than the numerical time step, will lead to an effective instantaneous re-equilibration.Since this happens automatically for large T we use the same rates also in hot matter.

A. MODIFIED URCA RATES
The modified Urca rate is also calculated using the Fermi surface approximation.The individual rates are presented in Sec. 5 of Alford et al. (2021b).The net rates of modified Urca, with a neutron or a proton spectator, are × δµ(1835π 6 T 6 + 945π 4 δµ 2 T 4 + 105π 2 δµ 4 T 2 + 3δµ 6 ).
Here we have introduced ) and The total modified Urca rate is given by the sum B. TOTAL URCA RATE The total net Urca rate in the Fermi Surface approximation is Since bulk viscosity is typically defined in the linear (also sometimes called "subthermal" Haensel et al. (2002); Alford et al. (2010)) regime where δµ ≪ T , it is natural to write this as Γ = Γ bulk + higher order in δµ/T (B9) where Γ bulk ≡ dΓ dδµ δµ=0 δµ. (B10) C. COMPARISON OF TIMESCALES In the following, we briefly compare the time scales provided by the Urca reaction rates in the previous section.In order to demonstrate that these rates can indeed affect the dynamics of the post-merger, we compare them to two relevant time scales present in the system.These effectively govern the propagation of sound waves through the remnant, as well as the characteristic time scale for post-merger gravitational wave emission.In order to bulk viciously damp sound waves, the damping time must be comparable or larger than the sound crossing time, t sound = d/c s .The latter is defined as the time it would take a sound wave to cross the diameter,d, of the remnant at its local sound speed c s .We can see in Fig. 6, that this is indeed always the case in cold regions, where the Urca rates act on timescales of about ≃ t sound .In hot regions, where the rates equilibrate instantaneously, bulk viscous damping cannot happen.These regions are clearly shown in dark blue, but constitute only a small fraction of the remnant.Similarly, we can show that the Urca rates operate on timescales much faster than the gravitational wave damping timescale t GW ≃ 30 ms Most et al. (2021).This means that bulk viscous damping and re-equilibration can efficiently affect the post-merger gravitational wave dynamics.Finally, we compare the reaction timescale to that of post-merger oscillations in the gravitational wave signal, as quantified by the timescale t f2 ≃ 1/f 2, associated with the dominant peak frequency of post-merger oscillations Bauswein et al. (2012); Takami et al. (2015).We can see that this timescale is similarly comparable across the merger remnant.

D. NEUTRINO TRAPPING
Under conditions probed in the merger neutrinos can become locally trapped Perego et al. (2019).In this case, local opacities are so high that the neutrinos cannot escape, effectively form a gas, constituting a correction to the overall pressure, energy density, and lepton number of a fluid element.In this Section, we will try to provide systematic estimates of the effect of neutrino trapping on our results.To this end, we will utilize a joint post-processing and direct simulation approach, where neutrino trapping is either estimated or directly included into our simulations following a similar approach to Kaplan et al. (2014).We begin by reviewing the basic description of a neutrino trapping.In a neutrino trapped fluid element, absorption and scattering will become as efficient as neutrino emission processes.This implies the formation of a trapped neutrino gas, which is in chemical equilibrium with the fluid, for non-vanishing neutrino chemical potential µ ν .While for non-trapped npe-matter, the conserved lepton Y l is given by the electron fraction Y e , in the case of trapping the total lepton number is given as the sum of electrons and trapped neutrinos.To this end, we introduce a trapped neutrino fraction Y ν = Y l − Y e = n ν /n B , where n ν is the trapped neutrino number density and n B the baryon density.For a degenerate Fermi gas in thermal equilibrium (i.e., having the same temperature) such as the trapped neutrino component, we may write Kaplan et al. (2014), Given an initial lepton fraction Y l , which we evolve in our simulations, we can then compute Y ν using the above expression, and obtain µ ν from Eq. (D11).The trapped neutrino component further exerts a pressure, P ν .Summed over all three neutrino species, we find that Kaplan et al. (2014) We note that due to their massless nature, the energy density e ν = 3P ν is trivially related to the pressure.

E. ERROR BUDGET OF THE GRAVITATIONAL WAVEFORMS
In the following, we will discuss the error budget of our gravitational wave signals.Since we are interested in comparing signals with different assumptions on chemical equilibration, our error budget is two-fold.First, all numerical relativity simulations are subject to finite resolution error.For the code used in this study, this error has been assessed and quantified in Most et al. (2019b).In particular, it has been demonstrated that the code is capable of achieving third-order convergence in the post-merger gravitational wave signal, even when using tabulated finite-temperature equations of state.This is particularly important when assessing dynamical aspect of the post-merger phase, such as those considered here.Second, our result are naturally subject to assumptions on neutrino transparencies, with the free-streaming assumption used in this work leading to enhanced chemical equilibration effects compared to a fully neutrino trapped regime.
In order to quantify these errors we have performed a new set of simulations for the system using the TMA EoS, which had displayed the largest differences in the gravitational wave signal.For all these simulations, we have adopted a resolution of 150 m, which is significantly higher than our standard 256 m resolution used previously.One of these new simulations applies a neutrino trapping correction (as outlined above) at densities n > 0.5 n sat and temperatures T > 1 MeV, roughly consistent with previous estimates in the literature Perego et al. (2019).This allows us to assess both errors independently.
The resulting gravitational waveforms are shown in Fig. 7.We can see a few milliseconds after merger, chemical equilibration leads to a dephasing of these waveforms, initially aligned at merger.For each of these waveforms we compute the power spectrum and identify the dominant peak in post-merger gravitational wave emission Bauswein et al. (2012).We find differences in those frequencies of ∆f ≈ (90,40,110) Hz, between a) trapping and no equilibration at high resolution, b) free-streaming and no-trapping at high resolution, and c) free-streaming and no equilibration at standard resolutions.For comparison, at standard resolutions, we find ∆f ≃ 40 Hz for the SFHo simulations.Furthermore, we also show the power spectrum for the highresolution TMA waveforms, which feature suppressed power in the f 2 mode for the free-streaming Urca case (red curve), which is a clear sign of bulk viscous damping, which can only be active in this case since it requires a dynamical relaxation of the chemical composition, rather than a fixed equilibrium.

Figure 1 .
Figure 1.Distribution over density of electron fractions, Ye, present in the system at time t = 5 ms after merger.Shown are the results for the simulations using the SFHo (top) and TMA (bottom) equations of state, with and without Urca processes.

Figure 2 .
Figure2.Bulk viscous pressure correction Π relative to the equilibrium pressure in the post-merger phase when using the TMA EoS with Urca processes.The green contours denote the rest-mass density of the merger remnant in units of nuclear saturation, at 3.6 ms relative to the time of merger tmer.

Figure 3 .Figure 4 .
Figure3.Chemical potential imbalance δµ (normalized by the temperature T ) as a measure of out-of-β-equilibrium effects during the merger and post-merger phase of two neutron stars when using the SFHo EoS with Urca processes (the orbital plane is shown).The green contours denote the rest-mass density of the merger remnant in units of nuclear saturation.All times, t, are stated relative to the time of merger tmer.

Figure 5 .
Figure 5. Normalized gravitational wave strain, h 22 + for the l = m = 2 component.Differences due to weak interactions in the post-merger are clearly visible.All times, t, are stated relative to the merger time and extraction radius r.The vertical lines indicate times of black hole formation.
2309210.MGA, AH, and ZZ are partly supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Award No. #DE-FG02-05ER41375.SPH is supported by the U. S. Department of Energy grant DE-FG02-00ER41132 as well as the National Science Foundation grant No. PHY-1430152 (JINA Center for the Evolution of the Elements).JN is partially supported by the U.S. Department of Energy, Office of Science, Office for Nuclear Physics under Award No. DE-SC0023861.ZZ was supported in part by the National Science Foundation (NSF) within the framework of the MUSES collaboration, under grant number OAC-2103680.APPENDIX DIRECT URCA RATES

Figure 6 .
Figure 6.Comparison of the Urca reaction timescale tUrca = 1/Γ to the sound crossing time t sound , the gravitational wave damping time tGW, and the timescale t f 2 of post-merger gravitational waves, about 6 ms after merger.

Figure 7 .
Figure 7. Gravitational wave strain h 22 + (l = m = 2 component) for the merger performed using the TMA equation of state.Different curves refer to various degrees of neutrino trapping.The bottom panel shows the frequency domain waveform around the dominant mode.A bulk viscous suppression in power for the Urca simulation is visible.