In this work we present a set of binary neutron star BNS merger simulations including the net muon fraction as an additional degreeoffreedom in the equation of state EoS and hydrodynamics evolution using the numericalrelativity code BAM. Neutrino cooling is modeled via a neutrinos leakage scheme including inmedium corrections to the opacities and emission rates of semileptonic chargedcurrent reactions. We show that for our particular choice of baseline baryonic EoS the presence of muons delays the gravitational collapse of the remnant compared to the case where muons are neglected. Furthermore when muons and muondriven neutrino reactions are considered no gravitational collapse occurs within our simulation timespan and muons are confined in the densest portions of the remnant while the disk is effectively colder less protonized and demuonized. Accordingly ejecta properties are affected e.g. ejecta masses are systematically smaller for the muonic setups and exhibit a larger fraction of neutronrich small velocity material. Overall our results suggest that the inclusion of muons and muonflavored neutrino reactions in the context of BNS merger simulations should not be neglected thus representing an important step towards more realistic modelling of such systems.
['General Relativity and Quantum Cosmology', 'Astrophysics - High Energy Astrophysical Phenomena', 'High Energy Physics - Phenomenology']
On the Role of Muons in Binary Neutron Star Mergers First Simulations
{'On the Role of Muons in Binary Neutron Star Mergers: First Simulations': 'Henrique Gieg 1 , Federico Schianchi 1 , Maximiliano Ujevic 2 , and Tim Dietrich 1 , 3  \n1 Institut f¨ur Physik und Astronomie, Universit¨at Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, 14476, Potsdam, Germany 2 Centro de Ciˆencias Naturais e Humanas, Universidade Federal do ABC, 09210-170, Santo Andr´e, S˜ao Paulo, Brazil Max Planck Institute for Gravitational Physics (Albert Einstein Institute),  \n3 Am M¨uhlenberg 1, Potsdam 14476, Germany  \n(Dated: September 9, 2024)  \nIn this work we present a set of binary neutron star (BNS) merger simulations including the net muon fraction as an additional degree-of-freedom in the equation of state (EoS) and hydrodynamics evolution using the numerical-relativity code BAM. Neutrino cooling is modeled via a neutrinos leakage scheme, including in-medium corrections to the opacities and emission rates of semi-leptonic charged-current reactions. We show that, for our particular choice of baseline baryonic EoS, the presence of muons delays the gravitational collapse of the remnant compared to the case where muons are neglected. Furthermore, when muons and muon-driven neutrino reactions are considered, no gravitational collapse occurs within our simulation timespan and muons are confined in the densest portions of the remnant, while the disk is effectively colder, less protonized and de-muonized. Accordingly, ejecta properties are affected, e.g., ejecta masses are systematically smaller for the muonic setups and exhibit a larger fraction of neutron-rich, small velocity material. Overall, our results suggest that the inclusion of muons and muon-flavored neutrino reactions in the context of BNS merger simulations should not be neglected, thus representing an important step towards more realistic modelling of such systems.', 'I. INTRODUCTION': 'Merging binary neutron stars (BNS) are among the most prominent sources of multimessenger signals, which combine gravitational-wave (GW) measurements [1-3] with their associated electromagnetic (EM) counterparts, such as gamma-ray bursts [4-10] and kilonovae, e.g. AT2017gfo [11-16]. Such observations allow, for instance, to establish constraints on the uncertain Equation of State (EoS) governing neutron-star matter at supranuclear densities [17-32], assess cosmological properties of the Universe [28, 33-38] and to study the formation of heavy elements [39-41].  \nOver the last years, many efforts have been made to interpret the available multimessenger data. One particularly successful approach consists in systematically connecting astrophysical observables to theoretical predictions produced by numerical-relativity simulations, which are needed given the strongly relativistic nature of such systems. In the context of BNS merger simulations, state-of-art comprises attempts to capture, as realistically as possible, a myriad of phenomena that is expected to take place under the extreme conditions encountered throughout the inspiral, merger and postmerger stages. Some examples include accurate modelling of the matter and hydrodynamics [42-52], the role of magnetic fields [53-58] and the inclusion of neutrinosdriven mechanisms [59-72].  \nOne key assumption regarding the modelling of matter is that the EoS contains only electrons e -and positrons e + as representative leptonic species. Hence, matter  \nproperties are described by EoSs that are represented as three-dimensional functions of the baryon rest-mass density ρ , temperature T and net electron fraction Y e . However, Core-Collapse Supernovae simulations [73-75] show that muons µ -and antimuons µ + are produced in nonnegligible amounts via neutrino-driven reactions during the formation of a neutron star. Over larger timescales, when the matter cools and reaches the neutrinoless β -equilibrium, muons are expected to be present wherever the electron chemical potential µ e exceeds the muon restmass m µ , i.e., µ e ≥ m µ c 2 ≈ 106 MeV [76-79]. Hence, the description of matter encompassed by the EoS should include an additional degree-of-freedom, accounting for the net muon fraction Y µ .  \nBesides, it has been shown that the presence of muons in the interior of neutron stars leads to important microphysical consequences, e.g. the arise of bulk viscosity due to leptonic reactions [80, 81], the modification of the direct Urca threshold [82], the occurrence of muon-flavored neutrinos-driven reactions and an overall increased proton fraction in locally neutral, β -stable matter when compared to EoSs that include only e -and e + .  \nInterestingly, up to our knowledge, there is only one work that addresses the possible impacts of considering muons within the remnant of a BNS [83], although their method consist in post-processing data from BNS merger simulations that were produced neglecting the presence of muons. Such a procedure provides important insights about the role of muons during the post-merger stage with respect to the hydrodynamics of matter and the behavior of neutrinos in the trapped regime. How-  \never, the method is not able to capture in details the complete dynamics and post-merger evolution of a BNS that is simulated ab-initio including muons and treating the neutrinos-driven interactions on-the-fly. Therefore, in this work we intend to surpass this shortcoming and present numerical-relativity simulations of BNS mergers carried out with the inclusion of muons in the EoS and hydrodynamics, and the use of a Neutrinos Leakage Scheme (NLS) to model the cooling of matter in response to the production of neutrinos, in particular accounting for muon-flavored neutrinos. The structure of this paper is as follows: in Sec. II we describe our procedures for the construction of EoSs including muons and the modifications of the general-relativistic hydrodynamics (GRHD) equations. In Sec. III we present details about the NLS implementation and the computation of emissivities and opacities for the neutrinos-driven reactions. In Sec. IV we state our methods and BNS setups considered in this worked, which were evolved with the BAM code [47, 49, 58, 71, 84-86]. In Sec. V we present a qualitative discussion about the merger and post-merger dynamics. In Sec. VI we make an analysis of the ejecta properties. Finally, in Sec. VII we state our concluding remarks. Throughout this work we adopt units in which the gravitational constant G , the speed of light in vacuum c , the solar mass M ⊙ and the Boltzmann constant k B are equal to one. Greek letters represent spacetime indices ranging from 0 to 3, while Latin letters are used for spacelike tensor and range from 1 to 3. The spacetime metric g µν has signature ( -, + , + , +).', 'II. EQUATION OF STATE AND HYDRODYNAMICS': "Generally, a muonic EoS may be constructed by 'dressing' a baseline baryonic EoS 1 , parameterized by the baryon number density n b , temperature T and proton fraction Y p , with a leptonic EoS. In the following, we consider the leptons and the anti-leptons l = { e -, µ -, e + , µ + } as relativistic ideal Fermi gases. Hence, the lepton number density n l ∓ , energy density ε l ∓ and pressure p l ∓ read [87, 88]  \nn l ∓ = K l β 3 / 2 l [ F 1 / 2 ( η 0 l ∓ , β l ) + β l F 3 / 2 ( η 0 l ∓ , β l ) ] ,  \nε l ∓ = K l m l c 2 β 5 / 2 l [ F 3 / 2 ( η 0 l ∓ , β l ) + β l F 5 / 2 ( η 0 l ∓ , β l ) ] + m l c 2 n l ∓ ,  \n(1) (2)  \np l ∓ = K l m l c 2 3 β 5 / 2 l [ 2 F 3 / 2 ( η 0 l ∓ , β l ) + β l F 5 / 2 ( η 0 l ∓ , β l ) ] ,  \n(3)  \nwhere  \nε γ = 8 π 5 15 T 4 ( hc ) 3 , p γ = ε γ / 3 , (11)  \ncorrespond to the EoS of photons at zero chemical potential and in thermal equilibrium with the matter. Note that we neglect the pressure exerted by trapped neutrinos and their contribution to the matter energy density. See, e.g., [83, 90] for investigations concerning the role of trapped neutrinos.  \nFor the purposes of this work, there are two special cases of interest for the muonic EoS: the first one is when Y µ = 0, which trivially sets η 0 µ -= η 0 µ + = -1 /β µ . The second one, relevant for the construction of initial data, comes by imposing the neutrinoless β -equilibrium condition for the reactions  \nn + ν l ↔ l -+ p. (12)  \nwhere m l is the lepton rest-mass, β l = T/m l c 2 is the relativity parameter, K l is a constant  \nK l = 8 π √ 2( m l c 2 /hc ) 3 , (4)  \nη 0 l ∓ are the non-relativistic degeneracy parameters  \nη 0 l -= µ l --m l c 2 T , (5)  \nη 0 l + = -( η 0 l -+ 2 β l ) , (6)  \nand µ l -is the relativistic lepton chemical potential (including rest-mass). Note that Eq. (5) is a definition, while Eq. (6) arises from the equilibrium between particles, antiparticles and photons (with zero chemical potential). Finally, F k ( η 0 l ∓ , β l ) is the generalized Fermi integral of order k , whose evaluation is performed numerically following [89].  \nFrom Eqs. (1), (5) and (6), it is straightforward to define the (net) lepton fraction Y l as  \nY l ( n b , T, η 0 l -) = [ n l -( η 0 l -, T ) -n l + ( η 0 l + , T )] /n b . (7)  \nSince our goal is to produce a tabulated EoS, our scheme begins by fixing the range of tabulated muon fractions Y µ . Then, for each ( n b , T, Y µ ), we numerically solve Eq. (7) for η 0 µ -. Next, for each ( n b , T, Y p , Y µ ), we define the electron fraction Y e by imposing local charge neutrality, i.e.  \nY e = Y p -Y µ , (8)  \nwhich is plugged in the LHS of Eq. (7) and numerically solved for η 0 e -.  \nOnce the lepton degeneracies are known for all ( n b , T, Y p , Y µ ) by means of the procedure above, contributions from leptons [Eqs. (2), (3)] and photons γ are added to the baseline baryonic ( b ) EoS, giving  \np = p b + p e -+ p e + + p µ -+ p µ + + p γ , (9)  \nε = ε b + ε e -+ ε e + + ε µ -+ ε µ + + ε γ , (10)  \nThus, the neutrino chemical potential vanishes and the lepton chemical potentials are given by  \nµ e ( n b , T, Y p ) = µ µ ( n b , T, Y p ) = µ n ( n b , T, Y p ) -µ p ( n b , T, Y p ) . (13)  \nSetting a constant temperature T = T 0 , e.g., the lowest tabulated temperature of the baryonic EoS, Eq. (8) is solved for each n b along with Eq. (13) for Y p . Then, adding the leptons and photons contributions according to Eqs. (2), (3), a one-dimensional cold neutrinoless β -equilibrated EoS is produced.  \nTo illustrate the changes introduced by the presence of muons in the composition of the matter, we depict in Fig. 1 the proton and muon fraction for a cold T = 0 . 1 MeV, β -equilibrated muonic EoS adopting the SFHo baryonic EoS [91]. For small baryon densities n b ≲ 0 . 125 fm -3 , the muonic (thick black line) and electronic (dashed black line) EoSs have the same proton fraction. Once µ e ≥ m µ c 2 (which is represented by the red line), muons are present within the matter ( Y µ > 0) and, accordingly, due to the local charge neutrality condition, the proton fraction becomes larger for the muonic EoS by a factor of at most 31% at n b = 1 . 5 fm -3 . It should be noted that this leads to a slight decrease of internal energy compared to when muons are absent because of (i) the conversion of electrons into muons and the subsequent loss of electron degeneracy energy by deoccupation of energy levels, which is consistent with the (small) reduction of the electron chemical potential when muons are present and (ii) the larger proton fraction leads to a loss of the neutron rest-mass contribution to the internal energy. But since muons only appear at moderately high densities, where the baryonic pressure and energy density dominate over the leptonic contributions, the impact of muons in macroscopic properties of a cold NS (e.g., mass, radius and tidal deformability) is negligible.  \nIn the following we summarize relevant information regarding the EoS used in this work. For comparisons purpose with results from the literature, we adopt the SFHo baseline EoSs. The baryon mass constant is chosen as m b = 1 . 659 × 10 -24 g and the rest-mass density is given by ρ = m b n b . The range of validity for the EoS parameters are ρ = [1 . 695 × 10 3 , 2 . 489 × 10 15 ] g / cm 3 , equispaced in log-scale with 30 points per decade, T = [0 . 1 , 120] MeV, equispaced in log-scale with 30 points per decade, Y p = [0 . 01 , 0 . 60] equispaced in linear scale with stride 0 . 01 and Y µ = [1 × 10 -4 , 1 × 10 -1 ], equispaced in log-scale with 20 points per decade plus Y µ = 0, for a total of 62 points in the Y µ dimension. All necessary thermodynamical information is obtained by means of a quadrilinear interpolation over the aforementioned EoS validity region. For such a parameterization choice, the 3+1 form of the GRHD equations of [92] are the same as in [93, 94] for the conserved rest-mass density D , internal energy density τ and momentum density S i , but here we  \nFIG. 1. Proton fraction for a cold, β -equilibrated slice of the SFHo EoS with muons (black thick line), its correspondent muon fraction (blue thick line), and the proton fraction without muons (black dashed line) as a function of the baryonic number density n b . The red vertical line represents the fractional difference ( µ e -m µ c 2 ) / ( m µ c 2 ). For µ e ≥ m µ c 2 , muons are present.  \n<!-- image -->  \nevolve Y p and Y µ according to  \n∂ 0 ( √ γDY p ) + ∂ i [ √ γDY p ( αv i -β i )] = α √ γS Y p , (14) ∂ 0 ( √ γDY µ ) + ∂ i [ √ γDY µ ( αv i -β i )] = α √ γS Y µ , (15)  \nwhere γ is the determinant of the spatial metric, α is the lapse function, β i is the shift vector, v i is the 3-velocity measured in the Eulerian frame, S Y p and S Y µ are source terms that we introduce in the next section. The extension of our previous high-resolution shock-capturing (HRSC) scheme [86] to the present case is straightforward.", 'A. A Brief Overview': 'The NLS [59-61, 68, 93, 95-104] is a simplified method to account for the radiative transport of neutrinos, devised to produce order-of-magnitude estimates of the role of neutrinos in astrophysical scenarios without resorting to approximate solutions of the Boltzmann equation, such as in lattice-Boltzmann transport schemes [105], moments-based schemes [62, 66, 67, 70, 106] or Monte Carlo schemes [107-109]. Differently than in previous BNS merger studies, instead of three neutrino species { ν e , ¯ ν e , ν x } , we consider five neutrino species { ν e , ¯ ν e , ν µ , ¯ ν µ , ν x } , where the heavy lepton neutrinos ν x = { ν τ , ¯ ν τ } are grouped as a single species with statistical weight 2. In the following, neutrinos are assumed  \nto be governed by the ultrarelativistic Fermi-Dirac distribution in local thermal and β -equilibrium with the matter [60], i.e., the degeneracy parameters read  \nη ν e = ( µ p + µ e -µ n ) /T, η ¯ ν e = -η ν e , (16)  \nη ν µ = ( µ p + µ µ -µ n ) /T, η ¯ ν µ = -η ν µ , (17)  \nη ν x = 0 , (18)  \nwhere the above chemical potentials include rest-mass.  \nThe NLS prescription adopted in this work assumes that the energy-momentum conservation applied to a matter element is modified according to  \n∇ ν T µν matter = -Q u µ , (19)  \nwhere ∇ ν is the covariant derivative compatible with the spacetime metric g µν , T µν matter is the stress-energy tensor of matter, considered here as an ideal fluid, u µ is the four-velocity of the matter element and Q is the total effective energy production rate, given by the sum of effective energy production rates from each neutrino species Q eff I , I = { ν e , ¯ ν e , ν µ , ¯ ν µ , ν x }  \nQ = ∑ I Q eff I . (20)  \nAs stated, the standard NLS only accounts for the cooling of the matter by emission of neutrinos. On the other hand, since neutrinos are leptons, the set of considered neutrinos-driven reactions consistently modify the lepton family number conservation laws as  \n∇ ν ( ρY e u ν ) = S Y e , (21)  \n∇ ν ( ρY µ u ν ) = S Y µ , (22)  \nwhile the baryon number conservation law reads  \n∇ ν ( ρu ν ) = 0 . (23)  \nIn face of the above equation, Eqs. (21), (22) become, respectively  \nS Y e = ρu ν ∇ ν ( Y e ) = ρ dY e dτ ≡ m b ( R eff ¯ ν e -R eff ν e ) , (24)  \nS Y µ = ρu ν ∇ ν ( Y µ ) = ρ dY µ dτ ≡ m b ( R eff ¯ ν µ -R eff ν µ ) , (25)  \nwhere d/dτ is the derivative with respect to the proper time of a fluid element. Hence R eff I is interpreted as the effective particle production rate of I in the fluid restframe. Finally, applying the local charge neutrality condition Eq. (8), the source term for Y p in Eq. (14) reads  \nS Y p = S Y e + S Y µ = m b ( R eff ¯ ν e -R eff ν e + R eff ¯ ν µ -R eff ν µ ) . (26)  \nFollowing [59, 69, 102, 104], the effective energy and particle production rates are computed, respectively, according to  \nQ eff I = Q I ( 1 + t diff I, 1 /t prod I, 1 ) -1 , (27)  \nR eff I = R I ( 1 + t diff I, 0 /t prod I, 0 ) -1 , (28)  \nwhere Q I , R I are the free energy and particle production rates, the production timescales are  \nt prod I, 1 = B I, 1 /Q I , t prod I, 0 = B I, 0 /R I , (29)  \nwith the neutrino energy density B I, 1  \nB I, 1 = g I 4 π ( hc ) 3 T 4 F 3 ( η I ) , (30)  \nthe neutrino number density B I, 0  \nB I, 0 = g I 4 π ( hc ) 3 T 3 F 2 ( η I ) , (31)  \nF k ( η I ) the ultrarelativistic Fermi integral of order k and the degeneracy factors g ν e = g ¯ ν e = g ν µ = g ¯ ν µ = 1, g ν x = 2. Note that given a set of reactions that produce neutrinos, all the aforementioned quantities may be estimated within our approach by direct interpolation from the EoS since thermal and chemical equilibrium is assumed.  \nHowever, the estimation of the diffusion timescale t diff I, 0 ( t diff I, 1 ) is more involved, since it depends on the local number-averaged opacitiy κ I, 0 (energy-averaged opacity κ I, 1 ), and on the non-local optical depth τ I, 0 ( τ I, 1 ) according to  \nt diff I,j = D τ 2 I,j cκ I,j , j = { 0 , 1 } , (32)  \nwith D = 6 chosen following [95]. For future convenience, we define here the I neutrino-sphere as the surface where τ I,j = 1, which represents the location outside of which neutrinos are effectively decoupled from matter [110].  \nThe optical depths are estimated following the iterative procedure of Ref. [111], i.e., during the initial timestep, the optical depths are iterated until convergence for all grid points. During the evolution, optical depths are recomputed at each point once per timestep by a single iteration using the optical depths from the previous timestep. An alternative approach, based on the solution of the Eikonal equation for the optical depths may be found in Ref. [112].  \nIn the following, we present in details the methods employed for the computation of opacities and emission rates for the processes considered in this work, which are summarized in Table I.', 'B. Opacities Computation': 'A crucial part of modelling neutrinos-driven processes consist in the evaluation of opacities associated with scattering and absorption reactions. As will become clear, a few differences are found between our opacities estimates and the widely adopted prescription for BNS studies, originally due to Ref. [59]. Instead, we closely follow Ref. [75].  \nTABLE I. Neutrino-driven reactions considered in this work. All the charged current processes are computed within the elastic approximation. Note that for pair processes and elastic scatterings, neutrinos of all species may participate.  \nWe begin by considering that the charged-current (CC) absorption processes of Table I may be generically represented as  \nν + N 1 → l + N 2 , (33)  \nwhich corresponds to the absorption of a neutrino ν by the nucleon N 1 , yielding the lepton l and the nucleon N 2 , where N 1 , N 2 = { n, p } .  \nFor simplicity, we restrict to model such reactions by means of the elastic approximation, i.e., neglecting the momentum transferred to nucleons by neutrinos. In this case, the absorption opacity is given by [75, 113, 115, 116]  \nκ abs I ( ϵ ) = σ 0 ( m e c 2 ) 2 (1 + 3 g 2 A ) 4 ( ϵ + Q ) 2 √ 1 -( m l c 2 ϵ + Q ) 2 [1 -f l ( ϵ + Q )] η 12 , (34)  \nwhere I = { ν e , ¯ ν e , ν µ , ¯ ν µ , ν x } , κ abs ν x ( ϵ ) = 0, ϵ is the incoming neutrino energy, σ 0 ≈ 1 . 705 × 10 -44 cm 2 , the axial coupling constant is g A ≈ 1 . 23 [114], E l = ϵ + Q is the energy of the lepton l and the medium-modified Q value is  \nQ = m ∗ 1 c 2 + U 1 -m ∗ 2 c 2 -U 2 , (35)  \nwhere m ∗ 1 / 2 is the effective mass and U 1 / 2 is the singleparticle vector-interaction potential of N 1 / 2 , generally provided by the EoS. Otherwise, estimates of U may be obtained following the procedure of [117]. The lepton distribution function f l is the Fermi-Dirac function  \nf l ( ϵ + Q ) = 1 1 + exp[( ϵ + Q ) /T -η l ] , (36)  \nand the nucleon phase-space blocking factor η 12 is  \nη 12 = n 2 -n 1 exp[( µ 2 -µ 1 + Q ) /T ] -1 , (37)  \nwhere n 1 / 2 is the nucleon number density and µ 1 / 2 is the nucleon chemical potential (including rest-mass).  \nTo avoid unphysical behavior of η 12 in the nondegenerate regime, we follow the prescription found  \nin [118] and set  \nη np = n n , (38)  \nη pn = n p , (39)  \nif µ n -µ p -Q < 0 . 01 MeV. Note that, although no other corrections are consided, we include in-medium effects in the limited kinematics of the absorption reaction Eq. (33) by means of the medium-modified Q factor.  \nThe next step, common to all energy-independent schemes, is to consider the spectral-average of the absorption opacity Eq. (34). To do so, the usual procedure consists in dropping the square root term in Eq. (34), which is equivalent to state that the energy of the outcoming lepton E l ≫ m l c 2 , i.e., that the produced leptons are ultrarelativistic. This is reasonable for electrons, since in general Q > m e c 2 , but for the case of muons, due to their substantially larger rest-mass, such an approximation is not adequate.  \nInstead, we keep the square-root term, but follow [59, 103] and approximate the lepton phase-space blocking through averaging the energy ¯ E l of the produced lepton via reaction Eq. (33)  \n[1 -f l ( ϵ + Q )] ≈ ⟨ 1 -f l ( ¯ E l ) ⟩ = { 1 + exp [ -( ¯ E l /T -η l )] } -1 . (40)  \nHence, the spectrally-averaged absorption opacity reads  \nκ abs I,j = 1 B I,j ⟨ 1 -f l ( ¯ E l ) ⟩ σ 0 ( m e c 2 ) 2 (1 + 3 g 2 A ) 4 η 12 I I,j , (41)  \nwhich is written in terms of the integral  \nI I,j ( m l , Q, T, η I ) = 4 π ( hc ) 3 T 5+ j ∫ ∞ x min ( x + Q/T ) 2 √ 1 -( m l c 2 xT + Q ) 2 x 2+ j f I ( x ) dx. (42)  \nIn Eq. (42) the lower integration limit is x min = max[0 , ( m l c 2 -Q ) /T ], which ensures that (i) the squareroot term is real and (ii) that only neutrinos with energies larger than m l c 2 -Q > 0 are absorbed. Naturally, f I ( x ) is the ultrarelativistic Fermi Dirac distribution function describing neutrinos, i.e.,  \nf I ( x ) = 1 1 + exp( x -η I ) . (43)  \nBefore proceeding to the methods employed to perform the integral Eq. (42), we are in position of defining the average energy of the produced lepton ¯ E l by noting that the average energy of the absorbed neutrinos may be estimated as  \n¯ E I = κ abs I, 1 /κ abs I, 0 = T ( I I, 1 / I I, 0 ) . (44)  \nThus, energy conservation implies  \n¯ E l = T I I, 1 I I, 0 + Q. (45)  \nWe verified that such a prescription, along with the lower integration bound x min defined later ensures that ¯ E l ≥ m l c 2 . It is straightforward to verify that when the square-root term of the integral Eq. (42) is neglected, one recovers Eq. (B13) of Ref. [103] from Eq. (41). Furthermore, neglecting Q , one recovers the widely adopted estimate of Ref. [59]  \n¯ E l = T F 5 ( η I ) F 4 ( η I ) .  \nSo far we have restated the problem of computing opacities as that of evaluating the integral Eq. (42). The current, widely adopted procedure of neglecting the squareroot term and the Q factor have the clear advantage of reducing the problem to the evaluation of the ultrarelativistic Fermi-Dirac integral, which is easily computed along a simulation given the pair ( T, η I ) by means of, e.g.,  \nthe sufficiently accurate formulas of Ref. [119]. However, as said, such an approach is not justified when applied to muon-driven reactions.  \nOn the other hand, the numerical integration of Eq. (42) on-the-fly is computationally intensive, since it may take up to hundreds of function evaluations per integral. Therefore, we resort to a pre-computation of the integrals as to produce, from the EoS as input, a table of spectrally-averaged opacities and emission rates parameterized by ( ρ, T, Y p , Y µ ), which are then used to compute opacities and emission rates along our simulations by means of quadrilinear interpolations.  \nTherefore, for each EoS point ( ρ, T, Y p , Y µ ), the integration of Eq. (42) is performed by adaptative quadratures up to desired accuracy with the Double Exponential method [120] as implemented in Refs. [121, 122]. Naturally, modifications concerning the kernel of the integral and the lower integration boundary are in order, since the original method is devised to integrate moments of the Fermi-Dirac distribution from x min = 0, which can be handled by simple variable transformations.  \nMore specifically, we first distinguish two cases:  \n- (i) If m l c 2 -Q ≤ 0, we integrate Eq. (42) with x min = 0, since neutrinos with all energies may participate of the reaction. The I degeneracy ¯ η I = η I and distribution function ¯ f I ( x ) = f I ( x ) remain unchanged,\n- (ii) If m l c 2 -Q > 0, we make ϵ + Q = E + m l c 2 and x = E/T . The I degeneracy is, thus, re-scaled by the transformation such that  \n˜ η I = η I -( m l c 2 -Q ) T , (46)  \n˜ f I ( x ) = 1 1 + exp ( x -˜ η I ) , (47)  \nand the integral Eq. (42) becomes  \n˜ I I,j = 4 π ( hc ) 3 T 5+ j ∫ ∞ 0 ( x + m l c 2 /T ) 2 √ 1 -( m l c 2 xT + m l c 2 ) 2 ( x + m l c 2 -Q T ) 2+ j ˜ f I ( x ) dx, (48)  \nwhere we omitted the dependencies of ˜ I I,j for a shorter notation.  \nOne last limit has to be considered when computing Eq. (41): the black-body function B I,j may be evaluated to zero, which generally occurs for very negative neutrino  \ndegeneracies, although the ratio I I,j /B I,j is finite. Thus, in order to circumvent such a possible issue, we first compute the black-body functions B I, 0 , B I, 1 . If one of those  \nfunctions evaluate to zero, we make  \nB I,j = 4 π ( hc ) 3 T 3+ j exp( η I )(2 + j )! , (49)  \nwhich comes from expanding f I ( x ) ≈ exp( η I ) exp( -x ) for exp( η I ) ≪ 1.  \nOn the other hand, when η I ≤ -100 (for m l c 2 -Q ≤ 0) or ¯ η I ≤ -100 (for m l c 2 -Q > 0), we proceed to similar expansions  \nf I ( x ) ≈ exp( η I ) exp( -x ) , (50)  \n¯ f I ( x ) ≈ exp(¯ η I ) exp( -x ) , (51)  \nand carry out the integrations Eq. (42) or Eq. (48) with a 64 points Gauss-Laguerre quadrature. By doing so, we have explicitly factored out the exp( η I ) term that may drive B I,j → 0, thus allowing the computation of finite ratios I I,j /B I,j .  \nFinally, for the elastic scattering of neutrinos ν in free nucleons N and heavy nuclei A  \nν + N → ν + N, (52)  \nν + A → ν + A, (53)  \nwe compute the respective scattering opacities κ scatt I,j ( N ) according to [93] and κ scatt I,j ( A ) according to [103]. Hence, the opacities used in Eq. (32) and for the computation of optical depths is simply the sum of the opacities over all processes, i.e.,  \nκ I,j = κ abs I,j + κ scatt I,j ( n ) + κ scatt I,j ( p ) + κ scatt I,j ( A ) , (54) with κ abs ν x ,j = 0.', 'C. Emission Rates Computation': 'For the emission via charged-current processes we consider the inverse of the reaction presented in Eq. (33). Following [115], detailed-balance sets the spectrallyaveraged emission rates as  \nQ CC I,j = ⟨ 1 -f I ( ¯ E I ) ⟩ σ 0 c ( m e c 2 ) 2 (1 + 3 g 2 A ) 4 η 21 I ∗ I,j , (55)  \nwhere the integral I ∗ I,j reads  \nI ∗ I,j ( m l , Q, T, η l ) = 4 π ( hc ) 3 T 5+ j ∫ ∞ x min ( x + Q/T ) 2 √ 1 -( m l c 2 xT + Q ) 2 x 2+ j f l ( x + Q/T ) dx. (56)  \nFor easy of notation and consistency with the text, we note that Q CC I, 0 = R CC I , Q CC I, 1 = Q CC I and Q CC ν x , 1 = Q CC ν x , 0 = 0.  \nIn this case the neutrinos produce the phase-space blocking, thus in complete analogy to Eq. (40) we define  \n⟨ 1 -f I ( ¯ E I ) ⟩ = { 1 + exp [ -( ¯ E I /T -η I )] } -1 , (57)  \nsuch that the average energy of the produced neutrino ¯ E I is given by  \n¯ E I = max [ 0 , T I ∗ I, 1 I ∗ I, 0 -Q ] . (58)  \nSimilarly to the calculation of absorption opacities, we distinguish two cases:  \n- (i) If m l c 2 -Q ≤ 0, we set x min = 0 and compute the integral Eq. (56) with the modified lepton distribution function f ∗ l ( x )  \nf ∗ l ( x ) = f l ( x + Q/T ) = 1 1 + exp ( x -η ∗ l ) , (59)  \nη ∗ l = η l -Q T , (60)  \nwhich is the same integral as Eq. (42), up to a substitution f I ( x ) → f ∗ l ( x ).  \n- (ii) If m l c 2 -Q > 0, we make again ϵ + Q = E + m l c 2 , x = E/T , which transforms the lepton distribution function and the lepton chemical potential, respectively, as  \n˜ f ∗ l ( x ) = 1 1 + exp ( x -¯ η ∗ l ) , (61)  \n¯ η ∗ l = η l -m l c 2 T . (62)  \nThe resulting integral, then, is the same as Eq. (48), up to a substitution ˜ f I ( x ) → ˜ f ∗ l ( x ).  \nAs per the pair processes, we follow the expressions of Ref. [59] for the electron-positron pair annihilation ( e -e + ) and transversal plasmon decay ( γ ) with a few adaptations, namely:  \n- (i) For ν x we divide their Eqs. (B10), (B12) by 2 to account for our statistical weight 2 instead of 4.\n- (ii) For ν µ and ¯ ν µ produced via e -e + , we use their Eq. (B8), changing the term ( C 1 + C 2 ) ν e ¯ ν e → ( C 1 + C 2 ) ν x ¯ ν x and employ the degeneracies Eq. (17) to  \ncompute the corresponding blocking factors in their Eq. (B9).  \n- (iii) Analogous adaptions were made in their Eqs. (B11), (B13) for the production of ν µ and ¯ ν µ via γ .  \nThen the free production rates in Eqs. (27), (28) are given by the sum of production rates over the chargedcurrent and pair processes.  \nWe end this section with a brief discussion about the possible contributions of Q , henceforth specified as  \nQ = m ∗ n c 2 + U n -m ∗ p c 2 -U p , (63)  \nto the opacities and emission rates for the CC processes. In the context of BNS simulations, such a correction has been either neglected, following Ref. [59] or partially considered, as in Ref. [103], where a constant Q = m n c 2 -m p c 2 = 1 . 2935 MeV is assumed. In the later, the authors show that CC opacities and emissivities involving ultrarelativistic leptons receive corrections that may (i) re-scale the degeneracy factors for leptons and neutrinos as Q/T and (ii) introduce additional terms proportional to Q , Q 2 and Q 3 , which are in general mild for Q = 1 . 2935 MeV, but become important for larger Q , in particular around its maximum.  \nTo illustrate the behavior of the medium-modified Q factor, we depict in the upper panels of Fig. 2 diagrams of Q (color-coded) as a function of ρ and T for three representative Y p = { 0 . 06 , 0 . 20 , 0 . 40 } .  \nThere we observe that Q has a very weak temperature dependency, only deviating from 1 . 2935 MeV at moderate densities ρ ≥ 10 12 g / cm 3 and reaching its maximum at 10 14 g / cm 3 < ρ < 10 15 g / cm 3 . Hence, we observe that the in-medium modifications set by Q are sizeable only at high densities. For concreteness, we show in the lower panel of Fig. 2 that in the case of an isolated, cold, β -equilibrated NS, the Q factor (blue line) is substantially larger than 1 . 2935 MeV within most of the NS interior, as ρ ≤ 10 12 g / cm 3 (black line) close to the edge of the star. Since Q is roughly independent of T and provided that β -equilibrium is reached sufficiently fast, we argue that this panel also qualitatively represents the Q factor profile in the densest portions of a BNS remnant.', 'IV. METHODS AND SETUPS': "In this work we performed four equal-mass, irrotational, BNS merger simulations. For our simulations we employed the SFHo baseline EoS. For comparison purposes, we ran one simulation with electrons and positrons only, which corresponds to a three-dimensional EoS, incorporating the usual three neutrino species { ν e , ¯ ν e , ν x } , ν x = ν µ , ¯ ν µ , ν τ , ¯ ν τ , with g ν x = 4. This setup is referred to as SFHo 3D. The other three setups were simulated with the full four-dimensional EoS. In order to single out the role of muons-driven reactions, one  \nof the four-dimensional EoS setup was simulated with the same aforementioned three neutrino species, hence named SFHo 4D 3, while the remaining two were simulated with five neutrino species { ν e , ¯ ν e , ν µ , ¯ ν µ , ν x } , ν x = ν τ , ¯ ν τ , g ν x = 2, identified as SFHo 4D 5, with same spatial grid resolution as the previously described runs, and SFHo 4D 5 High, with higher spatial grid resolution. All systems have total gravitational mass M = 2 . 70 M ⊙ and are initially governed by cold T = 0 . 1 MeV, neutrinoless β -equilibrated EoSs. The initial data was produced with the SGRID code [123-125], adapted to the use of one-dimensional tabulated EoSs as input. One caveat is that finite-temperature EoSs hardly reproduce the limit ( p, ε ) → 0 for n b → 0, which is needed for the proper imposition of boundary conditions. Instead, there is typically a (small) critical density n ∗ b for which de/dn b | n ∗ b = 0 and such that de/dn b > 0 for n b < n ∗ b [126]. To ensure the desired behavior at densities below the critical one, we assume that in this region the EoS is described by a cold polytrope.  \nThe dynamical evolution of the spacetime and matter was performed with the BAM code, which received the updates described in Sec. II and in Sec. III. The numerical domain consists in a hierarchy of 7 nested Cartesian grids (referred to as levels and indexed by L ≥ 0) with ∆ x L / ∆ x L +1 = 2 grid spacing refinement. For L ≥ 3 the levels move following the motion of the stars.  \nThe finest level ( L = 6) has grid spacing ∆ x 6 = 178 m and ∆ x 6 = 142 m for the high resolution run. The space is discretized with fourth-order finite differencing, while all fields evolve in time using the method of lines with an explicit fourth-order Runge-Kutta integrator and the Berger-Oliger time stepper [127]. Geometry quantities are evolved with the Z4c formulation of Einstein's equations [128, 129] along with the moving punctures method [84, 130]. Gauge fields are evolved adopting the 1 + log slicing [131] and the Gamma-driver conditions [132]. The HRSC scheme adopted for the evolution of matter fields employs the WENOZ reconstruction [133] and the HLL Riemann solver [134, 135] for the computation of inter-cell fluxes. Finally, a conservative adaptative mesh refinement strategy is used to ensure mass conservation across refinement levels [47]. The low density regions outside of the stars are treated with an artificial atmosphere prescription according to which matter elements are static and the thermodynamical properties follow from a cold and β -equilibrated slice of the EoS.", 'A. Matter Evolution': 'All simulations begin with a coordinate distance between stars ≈ 41 . 4 km, merging after ∼ 4 orbits. As expected, the presence of muons and muon-driven neutrino reactions does not affect the orbital dynamics during the inspiral.  \nFIG. 2. Upper panels: Medium-modified Q factor for the SFHo EoS at three representative proton fractions Y p = 0 . 06 (left panel), Y p = 0 . 20 (middle panel) and Y p = 0 . 40 (right panel). Note that the scales are different for each plot. Lower panel: Radial dependency of the Q factor (blue line) and log 10 ρ (black line) for a cold, β -equilibrated, M = 1 . 35 M ⊙ isolated star with SFHo EoS.  \n<!-- image -->  \nIn Fig. 3 we depict the time evolution of the maximum rest-mass density ρ max for our simulated setups. We begin noting that the muonic EoSs exhibit a slightly larger maximum rest-mass density before the merger, which is a consequence of the slightly smaller internal energy at same density introduced by the inclusion of muons in the cold, β -equilibrated EoSs employed for the construction of the ID.  \nIn good agreement with the results from Ref. [69], our SFHo 3D simulation collapsed to a black-hole in ∼ 12 ms after the merger which, despite important differences between the hydrodynamics and NLS implementations, is reassuring. Next, we note that the SFHo 4D 3 run has its collapse delayed by ∼ 7 ms compared to SFHo 3D, which suggests a stabilizing role of the muons in the densest portions of the remnant with a noticeable damping of ρ max oscillations until the collapse. This result is in line with the pressure increase in the densest portions of the remnant for the SFHo setup of Ref. [83] due to the presence of muons inherited from the cold coalescing NS cores. In fact, an increased pressure may postpone the collapse by preventing matter from contracting. Furthermore, the remnant of the 5-species run SFHo 4D 5 experiences a stronger damping of the oscillations and does not collapse within our simulation timespan. It is worth pointing out that gravitational collapses are rather sensitive to grid  \nresolution, thus the observation of a longer-lived remnant in SFHo 4D 5 High, albeit exhibiting weaker damping as the oscillations are sustained for longer, suggests that the stabilization is robust.  \nNow we proceed with an analysis of the remnant and disk structure. In Fig. 4 we depict the matter state on the x -y plane 10 ms after the merger. First we note that the SFHo 3D setup develops the most compact disk (upper row, first column), with a hot coredisk interface, pronounced shocked-tidal arms (middle row, first column) and a highly protonized disk, with Y p ≳ 0 . 25 up to 40 km from the origin (lower row, first column). The SFHo 4D 3 run (second column), exhibits more pronounced tidally-shocked arms, although with overall smaller temperatures throughout the remnant core and disk, leading to a less protonized disk. It is worth noticing that the higher proton fraction in the densest portions of the muonic runs is reminiscent from the more protonized initial data (see Fig. 1).  \nThe 5-species runs SFHo 4D 5 and SFHo 4D 5 High (third and fourth columns) develop an extended ρ ≥ 10 13 g / cm 3 region, with a noticeable suppression of the formation of tidally shocked arms. Accordingly, the whole remnant and disk are cooler than for SFHo 4D 3 and SFHo 3D and the proton fractions are smaller (lower row, third and fourth columns).  \nFIG. 3. Maximum rest-mass density evolution for our simulated setups. The merger instant t mrg is defined as the instant in which the amplitude of the dominant (2 , 2) mode of the GW strain is maximum.  \n<!-- image -->  \n-  \nThe same qualitative features extend to the remnant and disk structures as seen in the x -z plane slice presented in Fig. 5. The most noticeable difference concerns Y p , which is higher in the polar cap z ≥ 20 km for the 5-species runs.', 'B. Neutrino Emission': "In Fig. 6 we present the evolution of the neutrinos source luminosities. The pronounced peaks in the upper panels are associated to the gravitational collapse. Overall we note that the emissions peak around ∼ 2 ms for all setups. The higher total luminosity peak of SFHo 4D 3 (upper right panel) when compared to the SFHo 3D (upper left panel) at this instant follows from the higher ρ max reached by the former (see Fig. 3). In this case, the compression of matter increases the reaction rates across all neutrinos species.  \nAlong the post-merger stage, ¯ ν e dominates up to 6 ms, followed by ν x for the 3-species runs SFHo 3D and SFHo 4D 3. After that, L ν x is otherwise comparable to L ¯ ν e because of the high temperatures reached by the remnants and the strong temperature scaling of pair processes yielding ν x . It is worth pointing out that during the formation of a black hole, we don't adopt any particular excision strategy. Instead, we let the rest-mass density evolve and linearly extrapolate thermodynamical quantities for densities above the maximum tabulated one. Such a procedure is somewhat arbitrary, but since this only happens within the apparent horizon (hence causally disconnected from the remaining of the grid),  \nno significant effect is observed in the matter properties outside of the apparent horizon. However, since the optical depths computation depends on neighboring points, which might include points within the apparent horizon, unphysical optical depths may develop as a consequence of the linear extrapolation of opacities in those regimes. This is what we observe after the collapse for SFHo 4D 3: the maximum rest-mass density reaches between two and three times the maximum tabulated value, which produces unphysically high opacities and, consequently, very high optical depths within the apparent horizon. The optical depths at neighboring points then increase in response, leading to very small effective emission rates Eqs. (27), (28) and source luminosities. Contrary, for the SFHo 3D run, the maximum rest-mass density within the apparent horizon is ∼ 50% larger than the maximum tabulated value, thus the opacities don't reach as high of values and the optical depths are not contaminated by the region within the apparent horizon. Therefore, we are able to capture the fading luminosity emitted by the disk. For the 5-species runs SFHo 4D 5 (lower left panel) and SFHo 4D 5 High (lower right panel), we observe an overall larger total luminosity, which we associate to the additional cooling channels provided by the CC muonic reactions. The smaller luminosities reached by ν x compared to the 3-species simulations is due to the statistical weight 2 of the former instead of 4 of the later (see Section III C). The early post-merger neutrinos burst is such that the luminosities for ν e , ν µ and ¯ ν µ are comparable up to ∼ 5 ms. Once the emissions stabilize (from around 10 ms on), L ν µ < L ¯ ν µ because the hot and dense remnant is essentially optically thick and neutrinos mostly diffuse with average energies ⟨ E ¯ ν µ ⟩ > ⟨ E ν µ ⟩ , as suggested by the average neutrino energies presented in Tab. II. We note that here the average neutrino energy is estimated as the ratio between energy and particle source luminosities, defined as in Ref. [93], thus based on volume integrals over a grid level and without gravitational redshift corrections. Naturally, our reported average energies are higher by ∼ 50% for all species when compared to the literature (e.g. Refs. [67, 70, 71, 110]), which is an expected feature for leakage schemes given that neutrino luminosity and energy estimates include trapped neutrinos in the hot and dense remnant. On the other hand, more advanced treatments include absorption and neutrino properties are extracted far from the remnant, hence in the freely-streaming regime. Therefore, our estimates should be regarded as semi-quantitative.  \nNevertheless, for all cases we recover the usual hierarchy ⟨ E ν x ⟩ > ⟨ E ¯ ν e ⟩ > ⟨ E ν e ⟩ . Irrespective to the presence of muons and muon-driven reactions, ⟨ E ν e ⟩ , ⟨ E ¯ ν e ⟩ and ⟨ E ν x ⟩ agree within ∼ 5%. Interestingly, for the 5species runs we note that ∼ 10 ms after the merger, the neutrino-spheres for ¯ ν µ , ν x and ν µ are, respectively, located at increasing radii from the remnant center, although somewhat close, and are found deeper within the remnant, hence at higher matter temperatures, than the ¯ ν e and ν e neutrino-spheres. Thus, our results are in qual-  \nFIG. 4. Matter properties in the x -y plane for the SFHo runs at t -t mrg = 10 ms. Columns are referred from left to right as first to fourth. First column: SFHo 3D. Second column: SFHo 4D 3. Third column: SFHo 4D 5. Fourth column: SFHo 4D 5 High. Upper row: rest-mass density. Middle row: temperature. Lower row: proton fraction.  \n<!-- image -->  \nFIG. 5. Same as Fig. 4, but in the x -z plane.  \n<!-- image -->  \n[110] that  \nthe observed energy hierarchy is related to the higher  \nTABLE II. Average neutrino energy per species, in MeV, 10 ms after the merger for our simulated setups. The columns read simulation name, average electron neutrino energy, average anti-electron neutrino energy, average muon neutrino energy, average anti-muon neutrino energy and average heavy lepton-neutrino energy.  \ntemperature of the medium from which ν x , ν µ and ¯ ν µ decouple, with a caveat that there the authors employ an M0 scheme for the transport of energy and particle number. Besides, since the emission rates for pair processes are the same for ν µ and ¯ ν µ , their difference in average energies come from the CC reactions. In fact, by enforcing the lower bound on the energy of neutrinos that may be emitted via CC processes, we have E ¯ ν µ , min = m µ c 2 + Q ≈ 107 -170 MeV, while E ν µ , min = m µ c 2 -Q ≈ 46 -104 MeV. With our results we add to the currently described energy hierarchy the observation that ⟨ E ¯ ν µ ⟩ > ⟨ E ν x ⟩ > ⟨ E ν µ ⟩ > ⟨ E ¯ ν e ⟩ > ⟨ E ν e ⟩ , although further simulations employing more EoSs and an improved neutrino treatment are desirable in order to draw firmer conclusions.", 'C. Muon Content': 'In order to describe the evolution of the muonic content we introduce the conserved muon number, defined as  \nm b N µ = ∫ V DY µ √ γd 3 x, (64)  \nwhere D = Wρ is the rest-mass density in the Eulerian frame, W is the usual Lorentz factor and the integration volume V is a grid level. In the absence of muonneutrinos reactions, the balance-law Eq. (22) implies the approximate constancy of N µ along the evolution, which is verified for SFHo 4D 3 up to the collapse in the upper panel of Fig. 7. Conversely, when including the CC muonic reactions (that could alter the net muon fraction of a fluid element), as in the SFHo 4D 5( High) run, we observe that the conserved muon number decreases by as much as 8% (4%) with respect to the initial condition within our simulation timespan. In fact, we observe an early de-muonization, prompted by artificial shockheating between the stars surface and the atmosphere. However, such an effect is diminished both before and after the merger with increased grid resolution.  \nThe de-muonization is better visualized in the lower panels of Fig 7, where we show Y µ in the x -y plane 10 ms after the merger. For the SFHo 4D 3 run (lower left panel), we note that a substantial Y µ > 0 . 01 is  \ndistributed throughout the disk, which is a hydrodynamical effect attributed to the sourceless advection of muons coming from the remnant core. Contrary, for the SFHo 4D 5 (lower middle panel) and SFHo 4D 5 High (lower right panel) muons are found only within the remnant, where ρ ≥ 10 14 g cm -3 (see third and fourth, left panels of Fig. 4), rapidly decreasing outside of this region, as indicated by the Y µ = 10 -6 red dashed contour line.  \nHere we make some remarks about our findings. First, the muon fraction distributions along the x -y plane at t -t mrg = 5 ms is similar to Fig. 3 of Ref. [83], i.e., around the final instants of the neutrino bursts, we found Y µ = 10 -4 -10 -2 within 25 km of the recently formed remnant for the 5-species runs, which is due to the early redistribution of Y µ from the merging cores. Thus the de-muonization reported by us operates on longer timescales, mostly affecting the outermost disk regions.  \nNext, it is worth pointing out that, contrary to the CCSNe simulations of Refs. [73, 75], we do not observe the muonization of high-density matter, which is explainable by the following: first, their simulations start without muons within the matter. Then a substantial Y µ builds up only shortly before the core bounce and after it, most importantly via a two-stages process comprised of production of high-energy muon-(anti)neutrinos through pair processes that may then participate in muonic absorption reactions. Such a mechanism cannot be properly modeled by a neutrinos leakage scheme, because absorption is not realistically captured by the treatment. Second, the protonization observed in NLS simulations is based on an excess emission of ¯ ν e with respect to ν e , which occurs in the intermediate region between spatially separated ¯ ν e and ν e neutrino-spheres such that τ ¯ ν e , 0 = 1 is located closer to the remnant than τ ν e , 0 = 1. This is not the case for ν µ and ¯ ν µ neutrino-spheres. Instead, what we observe is that the ¯ ν µ neutrino-sphere is slightly wider than the ν µ neutrinosphere, as depicted in the lower panels of Fig. 7. This happens because the spectrally-averaged opacities employed in this work are heavily dominated by scattering processes in the muon(anti)-neutrino neutrino-spheres. Besides, there the neutrino degeneracies are relatively small η ν µ = -η ¯ ν µ ≈ -2, such that the number-averaged opacity is typically a few tens of percent larger for ¯ ν µ than for ν µ . On the other hand, the maximum values of the optical depths found during the post-merger are around one order of magnitude smaller for ¯ ν µ than for ν µ , which is expected given that the energy-dependent CC opacities are one to two orders of magnitude smaller for ¯ ν µ than for ν µ in remnant and disk conditions [75].  \nFinally, compared to late stages of the post-bounce reported by Ref. [75], we found Y µ in excess of about one order of magnitude. This is because most of the Y µ found in the remnant cores come from the initially cold, neutrinoless, β -equilibrated NSs, while in the former matter is still hot ( T ≈ 15 MeV) and far from β -equilibrium, i.e., there µ n -µ p < µ µ < µ e . On the other hand, our Y µ  \nFIG. 6. Source neutrino luminosity L ν evolution. Upper left panel: SFHo 3D, where the peak around 12 ms is related to the gravitational collapse. Upper right panel: SFHo 4D 3, where we followed the simulation up to 1 ms after the collapse for the purpose of comparing the remnant evolution with the remaining setups. Lower left panel: SFHo 4D 5. Lower right panel: SFHo 4D 5 High. In all setups we note a burst of neutrinos shortly after the merger, prompted by the heating that follows the compression of matter elements. Before the merger, neutrinos are produced due to artificial heating produced by shocks in the interface between the stars and the atmosphere, although to a lesser degree with increased grid resolution.  \n<!-- image -->  \n-  \n-  \nvalues are comparable to Ref. [73] in similar thermodynamic conditions, which also agrees with the reported in Ref. [83].', 'VI. EJECTA ANALYSIS': 'In this Section we present an analysis of the ejecta properties for our simulations. The relevant quantities are summarized in Table III. We note that the geodesic  \ncriterion [136] was adopted and the reported averages were extracted at a fixed sphere of coordinate radius r = 300 M ⊙ by the procedure outlined in Ref. [71].  \nWe begin comparing our SFHo 3D results with works that employ similar physical setups and neutrinos treatment. Our ejecta masses extracted at r = 200 M ⊙ are, respectively, ≈ 30%( ≈ 23%) smaller than those in Ref. [65] (Ref. [64]). When comparing ejecta masses extracted at r = 300 M ⊙ with those of Ref. [69], we have ≈ 37% less. Such differences may be partly explained by  \n<!-- image -->  \nFIG. 7. Evolution of the muon fraction for the muonic SFHo runs. Upper panel: Evolution of the conserved number of muons in the grid level L = 1. The 3-species run SFHo 4D 3 conserves m b N µ along the evolution, while the 5-species exhibits demuonization of matter. Here we observe that such effect takes place earlier due to the higher temperatures induced by artificial shock-heating, but to a lesser extent with increased resolution. Lower panels: Y µ on the x -y plane 10 ms after the merger. Lower left panel: SFHo 4D 3 run, where the region with Y µ > 0 . 01 extends up to 60 km from the origin. Lower middle panel: SFHo 4D 5 run, where the disk is found de-muonized, while Y µ ≳ 0 . 015 is present in regions where ρ ≥ 10 14 g cm -3 . Lower right panel: SFHo 4D 5 High run. The red dashed line marks Y µ = 10 -6 and the thick contour lines mark the neutrino-spheres of ν µ (cyan) and ¯ ν µ (gray), where τ ν µ , 0 = τ ¯ ν µ , 0 = 1.  \n<!-- image -->  \nTABLE III. Ejecta properties for our simulated setups. Columns show the simulation, ejecta mass extracted at r = 200 M ⊙ , ejecta mass extracted at r = 300 M ⊙ , average proton fraction, average entropy per baryon and average velocity measured by an observer at infinity. All the average quantities were extracted at r = 300 M ⊙ . For the muonic runs, ⟨ Y µ ⟩ ≲ 10 -3 , thus ⟨ Y p ⟩ ≈ ⟨ Y e ⟩ .  \ndifferences in the hydrodynamics implementations, e.g., Refs. [65, 69] employ a positivity-preserving limiter with the MP5 reconstruction and a different prescription for the atmosphere [137]. Additionally, we point out our approach for the computation of opacities and emissivities as a plausible source of discrepancy, given that ejecta  \nproperties are importantly impacted by the treatment of neutrinos.  \nRegarding average properties, good agreement is found for ⟨ Y p ⟩ , although our ⟨ s ⟩ is smaller by ∼ 24 % and v ∞ is smaller by ∼ 27 % compared to Ref. [69].  \nComparing our SFHo runs, the SFHo 4D 3 setup has  \na more protonized ⟨ Y p ⟩ = 0 . 20 and slightly more entropic ⟨ s ⟩ = 12 . 0 ejecta than SFHo 3D, although less massive by a factor 1 . 5 -1 . 8. The smaller amount of ejecta is consistent with the larger total luminosity of the former compared to the later (see Fig. 6), by means of which energetic matter elements become gravitationally bound due to neutrinos emission. The same rationale applies to the SFHo 4D 5 and SFHo 4D 5 High setups, which exhibit even higher total luminosities because of the additional muonic charged-current cooling channels. It should be noted, however, that our interpretation does not rule out the possibility that pressure changes due to the presence of muons may also affect the ejecta mass, as pointed out in Ref. [83].  \nAt same grid resolution, SFHo 3D and SFHo 4D 5 have very similar ⟨ Y p ⟩ and ⟨ s ⟩ , with differences mainly in v ∞ by a factor of ∼ 1 . 2 and in ejecta masses by a factor of 1 . 9 -2 . 2. With increasing resolution we note that SFHo 4D 5 and SFHo 4D 5 High vary by less than 10% in the ejecta masses and in less than 20% in ⟨ Y p ⟩ , ⟨ s ⟩ and v ∞ , allowing us to estimate numerical uncertainties of at least ∼ 20%.  \nIn Fig. 8, we present the distributions of ejected mass fractions with respect to the proton fraction Y p (left panel), entropy per baryon s (middle panel) and asymptotic velocity v ∞ (right panel) for our simulations, extracted at r = 300 M ⊙ . Since at this position the muon fraction Y µ for the muonic runs are much smaller than Y p , the distributions are identical with respect to Y e and our results may be compared to others from the literature. Comparing our SFHo 3D result (left panel) to similar runs of Refs. [65, 69], our Y p distribution is flatter from Y p = 0 . 06 up to the peak at Y p = 0 . 20, followed by a similar fall-off for Y p ≥ 0 . 3. For SFHo 4D 3 the distribution is more clearly peaked at Y p = 0 . 20, with considerably smaller fraction of neutron-rich material and a tail of neutron-poor material Y p ≤ 0 . 30. We verified that, in both simulations, the (neutron-rich) tidal component of the ejecta is rapidly reached by the (neutron-poor) shock-driven component of the ejecta, such that the material is reprocessed and the distribution is shifted towards higher Y p values. In the case of SFHo 4D 5 and SFHo 4D 5 High runs, the additional energy and momentum losses associated to the muonic reactions yield higher fractions of small velocity material (see right panel of Fig. 8). Hence, the reprocessing mechanism is inhibited and we note a pronounced neutronrich secondary peak along with the expected dominant neutron-poor peak. The entropy per baryon distribution (middle panel) is very similar across all setups at same grid resolution peaking at s ∼ 8 k B and followed by a rapid decay such that a negligible fraction of the ejecta is found with s ≥ 32 k B . For the SFHo 4D 5 High setup, the peak is shifted towards s ∼ 14 k B and a tail is found up to s ∼ 35 k B . The asymptotic velocity (right panel) follows a similar pattern for all simulations with a trend that increased amounts of small velocity ejecta are found in runs with higher total neutrinos luminosities. Besides,  \nwe did not find in our simulations a fast ejecta component as reported in Ref. [69].', 'VII. CONCLUSIONS': "In this work we presented, for the first time, a set of binary neutron star merger simulations that include muons and muon-driven neutrino reactions. To do so, we introduced a scheme to produce 4-dimensonal EoS tables parameterized by ( ρ, T, Y p , Y µ ) by 'dressing' a 3dimensional baryonic baseline EoS with a leptonic EoS modeling electrons, positrons, muons and antimuons as relativistic, ideal Fermi gases.  \nNext, we introduced a scheme for the tabulation of neutrinos opacities and emission rates that, differently from previous works [59, 103] in which most of the BNS studies are based on, here we included in-medium corrections embodied by the medium-modified Q factor for the charged-current absorption and emission reactions, as in Refs. [74, 75, 116, 117], but we restricted our treatment to the elastic approximation. In future works it would be important to consider the full kinematics approach of Ref. [75], since it implies important corrections to muonneutrinos opacities. Our particular interest for including Q comes from considering that Q = 30 -60 MeV leads to a noticeably smaller (higher) minimum energy that ν µ (¯ ν µ ) must have to participate in CC reactions. Thus, for consistency, we extended the same methods and formalism for ν e and ¯ ν e . We showed that indeed Q may differ substantially from the rest-mass difference Q = 1 . 29 MeV for most of the interior of an isolated, cold, β -equilibrated NS, suggesting that this correction may affect neutrinos properties in BNS remnant conditons. However, a detailed study of the impact of the use of the medium-modified Q factor would require a larger set of simulations specifically devised to investigate its role, which we reserve for future works. We ran a set of BNS simulations adopting the SFHo baseline baryonic EoS, modeling neutrinos as per a leakage scheme incorporating the aforementioned updated sets of neutrinos opacities and emission rates. For comparison purposes, our SFHo 3D setup was simulated with a 3dimensional EoS including electrons and positrons and the usual 3-neutrino species { ν e , ¯ ν e , ν x } . For this particular setup we observed gravitational collapse, in good agreement with Ref. [69], with a systematic underestimate of ejecta mass of a few tens of percent compared to Refs. [64, 65, 69]. At this stage, it is impossible to point out to which extent those differences arise due to different hydrodynamics implementations and to neutrinos treatment. We also ran three simulations including muons, namely SFHo 4D 3, SFHo 4D 5 and SFHo 4D 5 High, the first one with 3-neutrino species and the remaining with 5-neutrino species { ν e , ¯ ν e , ν µ , ¯ ν µ , ν x } , explicitly separating the muon-flavored neutrinos ν µ , ¯ ν µ from the heavy-lepton neutrinos ν x and adopting chargedcurrent muonic reactions. The last setup has increased  \nFIG. 8. Distribution of ejecta mass detected at r = 300 M ⊙ with respect to the proton fraction Y p (left panel), entropy per baryon (middle panel) and asymptotic velocity (right panel).  \n<!-- image -->  \ngrid resolution with respect to the remaining ones. For SFHo 4D 3 the collapse was delayed by ∼ 7 ms, while for SFHo 4D 5 there is no collapse within our simulation timespan (up to 20 ms after the merger) and we note a significant damping of the central rest-mass density oscillations. At increased resolution, SFHo 4D 5 High also does not collapse within 20 ms, suggesting that such stabilization feature is robust. Therefore, we conclude that the inclusion of muons indeed delays the collapse, in line with the prediction of Ref. [83], which describes an overall pressure increase for the SFHo EoS in the densest portions of the remnant due to the presence of muons. Regarding the evolution of the muon content along the postmerger stage we note that when muonic CC reactions are neglected, Y µ is advected from the remnant core, distributing Y µ > 0 . 01 relatively far within the disk. When muonic CC reactions are included, the disk is found effectively de-muonized and Y µ > 0 . 01 is restricted to dense portions of the remnant with ρ ≥ 10 14 g / cm 3 . For the muonic simulations, we found systematically smaller amounts of ejecta, but also higher total neutrino luminosities compared to SFHo 3D. This suggests that the ejecta production is affected both by the pressure increase in the remnant core and by the additional energy loss due to neutrino emission. Structure-wise, the muonic simulations exhibit less compact and cooler disks, with higher proton fraction in the polar cap. In particular, we noted that the inclusion of muon-driven CC reactions lead to the suppression of the formation of shocked arms along the disk. Furthermore, despite the initially higher proton fraction in the interior of NSs contain-  \ns, which is retained throughout the post-merger stage, the disk is less protonized than in the non-muonic counterpart. Overall our results suggest that the inclusion of muons and muon-driven reactions lead to significant consequences regarding the outcomes of BNS merger simulations, mostly affecting the post-merger evolution, the thermal and compositional structure of the remnant and impacting the ejecta properties. In future works we plan on extending the M1 scheme of Ref. [71] to account for muon-flavored neutrinos interactions and simulate a number of different baseline baryonic EoSs in order to assess the impacts of a more realistic neutrinos treatment. It would also be important to include the full kinematics approach of Ref. [74, 75] for the computation of semileptonic CC opacities and emission rates.", 'VIII. ACKNOWLEDGEMENTS': "HG, FS and TD acknowledge funding from the EU Horizon under ERC Starting Grant, no. SMArt101076369. The simulations were performed on HPE Apollo Hawk at the High Performance Computing (HPC) Center Stuttgart (HLRS) under the grant number GWanalysis/44189, on the GCS Supercomputer SuperMUC-NG at the Leibniz Supercomputing Centre (LRZ) [project pn29ba], and on the HPC systems Lise/Emmy of the North German Supercomputing Alliance (HLRN) [project bbp00049].  \n- [3] F. Acernese et al. (VIRGO), Class. Quant. Grav. 32 , 024001 (2015), arXiv:1408.3978 [gr-qc].  \n- [5] A. Hajela et al. , Astrophys. J. Lett. 927 , L17 (2022), arXiv:2104.02070 [astro-ph.HE].\n- [6] A. Balasubramanian, A. Corsi, K. P. Mooley, K. Hotokezaka, D. L. Kaplan, D. A. Frail, G. Hallinan, D. Lazzati, and E. J. Murphy, Astrophys. J. 938 , 12 (2022), arXiv:2205.14788 [astro-ph.HE].\n- [7] V. Nedora, T. Dietrich, M. Shibata, M. Pohl, and L. C. Menegazzi, (2022), 10.1093/mnras/stad175, arXiv:2208.01558 [astro-ph.HE].\n- [8] H. Wang, R. G. Dastidar, D. Giannios, and P. C. Duffell, (2024), arXiv:2402.19359 [astro-ph.HE].\n- [9] V. Savchenko et al. , Astrophys. J. Lett. 848 , L15 (2017), arXiv:1710.05449 [astro-ph.HE].\n- [10] B. P. Abbott et al. (LIGO Scientific, Virgo, FermiGBM, INTEGRAL), Astrophys. J. Lett. 848 , L13 (2017), arXiv:1710.05834 [astro-ph.HE].\n- [11] B. P. Abbott et al. (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Astrophys. J. Lett. 848 , L12 (2017), arXiv:1710.05833 [astro-ph.HE].\n- [12] I. Arcavi et al. , Nature 551 , 64 (2017), arXiv:1710.05843 [astro-ph.HE].\n- [13] D. A. Coulter et al. , Science (2017), 10.1126/science.aap9811, [Science358,1556(2017)], arXiv:1710.05452 [astro-ph.HE].\n- [14] V. M. Lipunov et al. , Astrophys. J. Lett. 850 , L1 (2017), arXiv:1710.05461 [astro-ph.HE].\n- [15] N. R. Tanvir et al. , Astrophys. J. 848 , L27 (2017), arXiv:1710.05455 [astro-ph.HE].\n- [16] S. Valenti, D. J. Sand, S. Yang, E. Cappellaro, L. Tartaglia, A. Corsi, S. W. Jha, D. E. Reichart, J. Haislip, and V. Kouprianov, Astrophys. J. 848 , L24 (2017), arXiv:1710.05854 [astro-ph.HE].\n- [17] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120 , 172703 (2018), arXiv:1711.02644 [astro-ph.HE].\n- [18] A. Bauswein, O. Just, H.-T. Janka, and N. Stergioulas, Astrophys. J. 850 , L34 (2017), arXiv:1710.06843 [astroph.HE].\n- [19] F. J. Fattoyev, J. Piekarewicz, and C. J. Horowitz, Phys. Rev. Lett. 120 , 172702 (2018), arXiv:1711.06615 [nucl-th].\n- [20] M. Ruiz, S. L. Shapiro, and A. Tsokaros, Phys. Rev. D 97 , 021501 (2018), arXiv:1711.00473 [astro-ph.HE].\n- [21] M. Shibata, S. Fujibayashi, K. Hotokezaka, K. Kiuchi, K. Kyutoku, Y. Sekiguchi, and M. Tanaka, Phys. Rev. D96 , 123012 (2017), arXiv:1710.07579 [astro-ph.HE].\n- [22] D. Radice, A. Perego, F. Zappa, and S. Bernuzzi, Astrophys. J. 852 , L29 (2018), arXiv:1711.03647 [astro-  \nph.HE].  \n- [23] E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, Phys. Rev. Lett. 120 , 261103 (2018), arXiv:1803.00549 [gr-qc].\n- [24] I. Tews, J. Margueron, and S. Reddy, Phys. Rev. C 98 , 045804 (2018), arXiv:1804.02783 [nucl-th].\n- [25] M. W. Coughlin et al. , Mon. Not. Roy. Astron. Soc. 480 , 3871 (2018), arXiv:1805.09371 [astro-ph.HE].\n- [26] M. W. Coughlin, T. Dietrich, B. Margalit, and B. D. Metzger, Mon. Not. Roy. Astron. Soc. 489 , L91 (2019), arXiv:1812.04803 [astro-ph.HE].\n- [27] C. D. Capano, I. Tews, S. M. Brown, B. Margalit, S. De, S. Kumar, D. A. Brown, B. Krishnan, and S. Reddy, Nature Astron. 4 , 625 (2020), arXiv:1908.10352 [astroph.HE].\n- [28] T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Science 370 , 1450 (2020), arXiv:2002.11355 [astro-ph.HE].\n- [29] V. Nedora, D. Radice, S. Bernuzzi, A. Perego, B. Daszuta, A. Endrizzi, A. Prakash, and F. Schianchi, Mon. Not. Roy. Astron. Soc. 506 , 5908 (2021), arXiv:2104.04537 [astro-ph.HE].\n- [30] S. Huth et al. , (2021), arXiv:2107.06229 [nucl-th].\n- [31] H. Koehn, T. Wouters, H. Rose, P. T. H. Pang, R. Somasundaram, I. Tews, and T. Dietrich, (2024), arXiv:2407.07837 [astro-ph.HE].\n- [32] H. Koehn et al. , (2024), arXiv:2402.04172 [astroph.HE].\n- [33] C. Guidorzi et al. , Astrophys. J. Lett. 851 , L36 (2017), arXiv:1710.06426 [astro-ph.CO].\n- [34] K. Hotokezaka, E. Nakar, O. Gottlieb, S. Nissanke, K. Masuda, G. Hallinan, K. P. Mooley, and A. T. Deller, Nature Astron. 3 , 940 (2019), arXiv:1806.10596 [astro-ph.CO].\n- [35] M. W. Coughlin, T. Dietrich, J. Heinzel, N. Khetan, S. Antier, M. Bulla, N. Christensen, D. A. Coulter, and R. J. Foley, Phys. Rev. Res. 2 , 022006 (2020), arXiv:1908.00889 [astro-ph.HE].\n- [36] M. A. P'erez-Garc'ıa et al. , (2022), arXiv:2204.00022 [astro-ph.CO].\n- [37] H. Wang and D. Giannios, Astrophys. J. 908 , 200 (2021), arXiv:2009.04427 [astro-ph.HE].\n- [38] M. Bulla, M. W. Coughlin, S. Dhawan, and T. Dietrich, (2022), arXiv:2205.09145 [astro-ph.HE].\n- [39] J. M. Lattimer and D. N. Schramm, Astrophys. J. 192 , L145 (1974).\n- [40] S. Rosswog, M. Liebendoerfer, F. Thielemann, M. Davies, W. Benz, et al. , Astron.Astrophys. 341 , 499 (1999), arXiv:astro-ph/9811367 [astro-ph].\n- [41] S. Rosswog, Astrophys. J. 634 , 1202 (2005), arXiv:astro-ph/0508138 [astro-ph].\n- [42] M. Shibata, K. Taniguchi, and K. Uryu, Phys. Rev. D68 , 084020 (2003), arXiv:gr-qc/0310030.\n- [43] M. Shibata, K. Taniguchi, and K. Uryu, Phys. Rev. D71 , 084021 (2005), arXiv:gr-qc/0503119.\n- [44] L. Rezzolla, L. Baiotti, B. Giacomazzo, D. Link, and J. A. Font, Class. Quant. Grav. 27 , 114105 (2010), arXiv:1001.3074 [gr-qc].\n- [45] A. Bauswein, S. Goriely, and H.-T. Janka, Astrophys.J. 773 , 78 (2013), arXiv:1302.6530 [astro-ph.SR].\n- [46] K. Hotokezaka, K. Kiuchi, K. Kyutoku, T. Muranushi, Y.-i. Sekiguchi, M. Shibata, and K. Taniguchi, Phys. Rev. D 88 , 044026 (2013), arXiv:1307.5888 [astroph.HE].\n- [47] T. Dietrich, S. Bernuzzi, M. Ujevic, and B. Brugmann, Phys. Rev. D91 , 124041 (2015), arXiv:1504.01266 [grqc].\n- [48] T. Dietrich, N. Moldenhauer, N. K. JohnsonMcDaniel, S. Bernuzzi, C. M. Markakis, B. Brugmann, and W. Tichy, Phys. Rev. D92 , 124007 (2015), arXiv:1507.07100 [gr-qc].\n- [49] S. Bernuzzi and T. Dietrich, Phys. Rev. D94 , 064062 (2016), arXiv:1604.07999 [gr-qc].\n- [50] D. Radice, L. Rezzolla, and F. Galeazzi, Mon.Not.Roy.Astron.Soc. 437 , L46 (2014), arXiv:1306.6052 [gr-qc].\n- [51] D. Radice, L. Rezzolla, and F. Galeazzi, Proceedings, Numerical Modeling of Space Plasma Flows (ASTRONUM-2014): Long Beach, CA, USA, June 23-27, 2014 , ASP Conf. Ser. 498 , 121 (2015), arXiv:1502.00551 [gr-qc].\n- [52] D. Radice, Astrophys. J. 838 , L2 (2017), arXiv:1703.02046 [astro-ph.HE].\n- [53] B. Giacomazzo, J. Zrake, P. Duffell, A. I. MacFadyen, and R. Perna, Astrophys. J. 809 , 39 (2015), arXiv:1410.0013 [astro-ph.HE].\n- [54] D. M. Siegel and R. Ciolfi, Springer Proc. Phys. 170 , 119 (2016), arXiv:1401.5275 [gr-qc].\n- [55] K. Kiuchi, P. Cerd'a-Dur'an, K. Kyutoku, Y. Sekiguchi, and M. Shibata, Phys. Rev. D 92 , 124034 (2015), arXiv:1509.09205 [astro-ph.HE].\n- [56] R. Ciolfi, Gen. Rel. Grav. 52 , 59 (2020), arXiv:2003.07572 [astro-ph.HE].\n- [57] C. Palenzuela, R. Aguilera-Miret, F. Carrasco, R. Ciolfi, J. V. Kalinani, W. Kastaun, B. Mi˜nano, and D. Vigan'o, Phys. Rev. D 106 , 023013 (2022), arXiv:2112.08413 [grqc].\n- [58] A. Neuweiler, T. Dietrich, B. Brugmann, E. Giangrandi, K. Kiuchi, F. Schianchi, P. Mosta, S. Shankar, B. Giacomazzo, and M. Shibata, (2024), arXiv:2407.20946 [gr-qc].\n- [59] M. H. Ruffert, H. T. Janka, and G. Schafer, Astron. Astrophys. 311 , 532 (1996), arXiv:astro-ph/9509006.\n- [60] S. Rosswog and M. Liebendoerfer, Mon.Not.Roy.Astron.Soc. 342 , 673 (2003), arXiv:astroph/0302301 [astro-ph].\n- [61] Y. Sekiguchi, Class.Quant.Grav. 27 , 114107 (2010), arXiv:1009.3358 [astro-ph.HE].\n- [62] M. Shibata, K. Kiuchi, Y.-i. Sekiguchi, and Y. Suwa, Prog.Theor.Phys. 125 , 1255 (2011), arXiv:1104.3937 [astro-ph.HE].\n- [63] F. Foucart, R. Haas, M. D. Duez, E. O'Connor, C. D. Ott, L. Roberts, L. E. Kidder, J. Lippuner, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D93 , 044019 (2016), arXiv:1510.06398 [astro-ph.HE].\n- [64] L. Lehner, S. L. Liebling, C. Palenzuela, O. L. Caballero, E. O'Connor, M. Anderson, and D. Neilsen, Class. Quant. Grav. 33 , 184002 (2016), arXiv:1603.00501 [gr-qc].\n- [65] L. Bovard, D. Martin, F. Guercilena, A. Arcones, L. Rezzolla, and O. Korobkin, Phys. Rev. D96 , 124005 (2017), arXiv:1709.09630 [gr-qc].\n- [66] F. Foucart, E. O'Connor, L. Roberts, M. D. Duez, R. Haas, L. E. Kidder, C. D. Ott, H. P. Pfeiffer, M. A. Scheel, and B. Szilagyi, Phys. Rev. D91 , 124021 (2015), arXiv:1502.04146 [astro-ph.HE].\n- [67] F. Foucart, E. O'Connor, L. Roberts, L. E. Kidder, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D94 , 123016  \n- (2016), arXiv:1607.07450 [astro-ph.HE].\n- [68] D. Radice, F. Galeazzi, J. Lippuner, L. F. Roberts, C. D. Ott, and L. Rezzolla, Mon. Not. Roy. Astron. Soc. 460 , 3255 (2016), arXiv:1601.02426 [astro-ph.HE].\n- [69] D. Radice, A. Perego, K. Hotokezaka, S. A. Fromm, S. Bernuzzi, and L. F. Roberts, Astrophys. J. 869 , 130 (2018), arXiv:1809.11161 [astro-ph.HE].\n- [70] D. Radice, S. Bernuzzi, A. Perego, and R. Haas, Mon. Not. Roy. Astron. Soc. 512 , 1499 (2022), arXiv:2111.14858 [astro-ph.HE].\n- [71] F. Schianchi, H. Gieg, V. Nedora, A. Neuweiler, M. Ujevic, M. Bulla, and T. Dietrich, Phys. Rev. D 109 , 044012 (2024), arXiv:2307.04572 [gr-qc].\n- [72] D. Radice and S. Bernuzzi, (2023), arXiv:2306.13709 [astro-ph.HE].\n- [73] R. Bollig, H. T. Janka, A. Lohs, G. Martinez-Pinedo, C. J. Horowitz, and T. Melson, Phys. Rev. Lett. 119 , 242702 (2017), arXiv:1706.04630 [astro-ph.HE].\n- [74] G. Guo, G. Mart'ınez-Pinedo, A. Lohs, and T. Fischer, Phys. Rev. D 102 , 023037 (2020), arXiv:2006.12051 [hep-ph].\n- [75] T. Fischer, G. Guo, G. Mart'ınez-Pinedo, M. Liebendorfer, and A. Mezzacappa, Phys. Rev. D 102 , 123001 (2020), arXiv:2008.13628 [astro-ph.HE].\n- [76] S. L. Shapiro and S. A. Teukolsky, Black holes, white dwarfs, and neutron stars: The physics of compact objects (Wiley, New York, USA, 1983).\n- [77] N. K. Glendenning, Compact stars: Nuclear physics, particle physics, and general relativity (1997).\n- [78] P. Haensel, K. Levenfish, and D. Yakovlev, Astronomy & Astrophysics 357 , 1157 (2000).\n- [79] Haensel, P., Levenfish, K. P., and Yakovlev, D. G., Astronomy & Astrophysics 372 , 130 (2001).\n- [80] M. G. Alford and G. Good, Phys. Rev. C 82 , 055805 (2010), arXiv:1003.1093 [nucl-th].\n- [81] M. Alford, A. Harutyunyan, and A. Sedrakian, Phys. Rev. D 104 , 103027 (2021), arXiv:2108.07523 [astroph.HE].\n- [82] J. M. Lattimer, M. Prakash, C. J. Pethick, and P. Haensel, Phys. Rev. Lett. 66 , 2701 (1991).\n- [83] E. Loffredo, A. Perego, D. Logoteta, and M. Branchesi, Astron. Astrophys. 672 , A124 (2023), arXiv:2209.04458 [astro-ph.HE].\n- [84] Brugmann, Bernd and Gonzalez, Jose A. and Hannam, Mark and Husa, Sascha and Sperhake, Ulrich and Tichy, Wolfgang, Physical Review D 77 , 024027 (2008), arXiv:gr-qc/0610128 [gr-qc].\n- [85] Thierfelder, Marcus and Bernuzzi, Sebastiano and Brugmann, Bernd, Physical Review D 84 , 044012 (2011), arXiv:1104.4751 [gr-qc].\n- [86] H. Gieg, F. Schianchi, T. Dietrich, and M. Ujevic, Universe 8 , 370 (2022), arXiv:2206.01337 [gr-qc].\n- [87] S. A. Bludman and K. A. van Riper, Astrophys. J. 212 , 859 (1977).\n- [88] F. X. Timmes and D. Arnett, The Astrophysical Journal Supplement Series 125 , 277 (1999).\n- [89] J. M. Aparicio, The Astrophysical Journal Supplement Series 117 , 627 (1998).\n- [90] A. Perego, S. Bernuzzi, and D. Radice, Eur. Phys. J. A55 , 124 (2019), arXiv:1903.07898 [gr-qc].\n- [91] A. W. Steiner, M. Hempel, and T. Fischer, Astrophys. J. 774 , 17 (2013), arXiv:1207.2184 [astro-ph.SR].\n- [92] J. A. Font, Living Rev. Rel. 11 , 7 (2008).\n- [93] F. Galeazzi, W. Kastaun, L. Rezzolla, and J. A. Font, Phys. Rev. D 88 , 064009 (2013), arXiv:1306.4953 [grqc].\n- [94] L. Rezzolla and O. Zanotti, Relativistic Hydrodynamics , 1st ed., Mathematics (Oxford University Press, Oxford, 2013).\n- [95] E. O'Connor and C. D. Ott, Class. Quant. Grav. 27 , 114103 (2010), arXiv:0912.2393 [astro-ph.HE].\n- [96] E. O'Connor, Astrophys. J. Suppl. 219 , 24 (2015), arXiv:1411.7058 [astro-ph.HE].\n- [97] M. B. Deaton, M. D. Duez, F. Foucart, E. O'Connor, C. D. Ott, L. E. Kidder, C. D. Muhlberger, M. A. Scheel, and B. Szilagyi, Astrophys. J. 776 , 47 (2013), arXiv:1304.3384 [astro-ph.HE].\n- [98] F. Foucart, M. B. Deaton, M. D. Duez, E. O'Connor, C. D. Ott, R. Haas, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and B. Szilagyi, Phys. Rev. D90 , 024026 (2014), arXiv:1405.1121 [astro-ph.HE].\n- [99] A. Perego, E. Gafton, R. Cabez'on, S. Rosswog, and M. Liebendorfer, Astron. Astrophys. 568 , A11 (2014), arXiv:1403.1297 [astro-ph.HE].\n- [100] A. Perego, R. Cabezon, and R. Kaeppeli, Astrophys. J. Suppl. 223 , 22 (2016), arXiv:1511.08519 [astro-ph.IM].\n- [101] C. Palenzuela, S. L. Liebling, D. Neilsen, L. Lehner, O. L. Caballero, E. O'Connor, and M. Anderson, Phys. Rev. D92 , 044045 (2015), arXiv:1505.01607 [gr-qc].\n- [102] D. M. Siegel and B. D. Metzger, Astrophys. J. 858 , 52 (2018), arXiv:1711.00868 [astro-ph.HE].\n- [103] R. Ardevol-Pulpillo, H. T. Janka, O. Just, and A. Bauswein, Mon. Not. Roy. Astron. Soc. 485 , 4754 (2019), arXiv:1808.00006 [astro-ph.HE].\n- [104] A. Murguia-Berthier et al. , Astrophys. J. 919 , 95 (2021), arXiv:2106.05356 [astro-ph.HE].\n- [105] L. R. Weih, A. Gabbana, D. Simeoni, L. Rezzolla, S. Succi, and R. Tripiccione, Mon. Not. Roy. Astron. Soc. 498 , 3374 (2020), arXiv:2007.05718 [physics.compph].\n- [106] P. Anninos and P. C. Fragile, The Astrophysical Journal 900 , 71 (2020).\n- [107] F. Foucart, M. D. Duez, F. Hebert, L. E. Kidder, P. Kovarik, H. P. Pfeiffer, and M. A. Scheel, Astrophys. J. 920 , 82 (2021), arXiv:2103.16588 [astro-ph.HE].\n- [108] K. Kawaguchi, S. Fujibayashi, and M. Shibata, Phys. Rev. D 107 , 023026 (2023), arXiv:2209.12472 [astroph.HE].\n- [109] F. Foucart, M. D. Duez, R. Haas, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and E. Spira-Savett, Phys. Rev. D 107 , 103055 (2023), arXiv:2210.05670 [astro-ph.HE].\n- [110] M. Cusinato, F. M. Guercilena, A. Perego, D. Logoteta, D. Radice, S. Bernuzzi, and S. Ansoldi, (2021), 10.1140/epja/s10050-022-00743-5, arXiv:2111.13005 [astro-ph.HE].\n- [111] D. Neilsen, S. L. Liebling, M. Anderson, L. Lehner, E. O'Connor, et al. , Phys.Rev. D89 , 104029 (2014), arXiv:1403.3680 [gr-qc].\n- [112] C. Palenzuela, S. Liebling, and B. Mi˜nano, Phys. Rev. D 105 , 103020 (2022), arXiv:2204.02721 [gr-qc].\n- [113] M. Rampp and H. T. Janka, Astron. Astrophys. 396 , 361 (2002), arXiv:astro-ph/0203101.  \n- [114] A. Burrows, S. Reddy, and T. A. Thompson, Nuclear Physics A 777 , 356 (2006).\n- [115] S. W. Bruenn, apjs 58 , 771 (1985).\n- [116] H. H.-Y. Ng, P. C.-K. Cheong, A. T.-L. Lam, and T. G. F. Li, (2023), arXiv:2309.03526 [astro-ph.HE].\n- [117] G. Martinez-Pinedo, T. Fischer, A. Lohs, and L. Huther, Phys. Rev. Lett. 109 , 251104 (2012), arXiv:1205.2793 [astro-ph.HE].\n- [118] T. Kuroda, T. Takiwaki, and K. Kotake, The Astrophysical Journal Supplement Series 222 , 20 (2016).\n- [119] K. Takahashi, M. F. El Eid, and W. Hillebrandt, AAP 67 , 185 (1978).\n- [120] H. Takahasi and M. Mori, Publications of the Research Institute for Mathematical Sciences 9 , 721 - 741 (1974).\n- [121] N. Mohankumar, T. Kannan, and S. Kanmani, Computer Physics Communications 168 , 71 (2005).\n- [122] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes 3rd Edition: The Art of Scientific Computing , 3rd ed. (Cambridge University Press, USA, 2007).\n- [123] W. Tichy, Classical and Quantum Gravity 26 , 175018 (2009), arXiv:0908.0620 [gr-qc].\n- [124] W. Tichy, Phys. Rev. D 86 , 064024 (2012), arXiv:1209.5336 [gr-qc].\n- [125] W. Tichy, A. Rashti, T. Dietrich, R. Dudi, and B. Brugmann, Phys. Rev. D 100 , 124046 (2019), arXiv:1910.09690 [gr-qc].\n- [126] M. G. Alford, L. Brodie, A. Haber, and I. Tews, Phys. Rev. C 106 , 055804 (2022), arXiv:2205.10283 [nucl-th].\n- [127] M. J. Berger and J. Oliger, Journal of Computational Physics 53 , 484 (1984).\n- [128] S. Bernuzzi and D. Hilditch, Phys. Rev. D81 , 084003 (2010), arXiv:0912.2920 [gr-qc].\n- [129] D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao, W. Tichy, and B. Brugmann, Phys. Rev. D88 , 084057 (2013), arXiv:1212.2901 [gr-qc].\n- [130] M. Campanelli, C. Lousto, P. Marronetti, and Y. Zlochower, Phys. Rev. Lett. 96 , 111101 (2006), arXiv:grqc/0511048.\n- [131] C. Bona, J. Masso, E. Seidel, and J. Stela, Phys. Rev. Lett. 75 , 600 (1995), arXiv:gr-qc/9412071.\n- [132] M. Alcubierre, B. Brugmann, P. Diener, M. Koppitz, D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D 67 , 084023 (2003), arXiv:gr-qc/0206072.\n- [133] R. Borges, M. Carmona, B. Costa, and W. S. Don, Journal of Computational Physics 227 , 3191 (2008).\n- [134] A. Harten, P. D. Lax, and B. v. Leer, SIAM Review 25 , 35 (1983).\n- [135] E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics (Springer-Verlag, 1999).\n- [136] K. Hotokezaka, K. Kiuchi, K. Kyutoku, H. Okawa, Y.i. Sekiguchi, M. Shibata, and K. Taniguchi, Physical Review D 87.2 , 024001 (2013).\n- [137] D. Radice, L. Rezzolla, and F. Galeazzi, Class.Quant.Grav. 31 , 075012 (2014), arXiv:1312.5004 [gr-qc]."}