Issue 
A&A
Volume 672, April 2023



Article Number  A124  
Number of page(s)  19  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/202244927  
Published online  10 April 2023 
Muons in the aftermath of neutron star mergers and their impact on trapped neutrinos
^{1}
Gran Sasso Science Institute, Viale Francesco Crispi 7, 67100 L’Aquila, Italy
email: eleonora.loffredo@gssi.it
^{2}
INFN – Laboratori Nazionali del Gran Sasso, Via G. Acitelli 22, 67100 Assergi L’Aquila, Italy
^{3}
INAF – Osservatorio Astronomico d’Abruzzo, Via M. Maggini snc, 64100 Teramo, Italy
^{4}
Dipartimento di Fisica, Universitá di Trento, Via Sommarive 14, 38123 Trento, Italy
^{5}
INFNTIFPA, Trento Institute for Fundamental Physics and Applications, Via Sommarive 14, 38123 Trento, Italy
^{6}
Dipartimento di Fisica, Universitá di Pisa, Largo B. Pontecorvo, 3, 56127 Pisa, Italy
^{7}
INFN, Sezione di Pisa, Largo B. Pontecorvo, 3, 56127 Pisa, Italy
Received:
8
September
2022
Accepted:
15
February
2023
Context. In the upcoming years, present and nextgeneration gravitational wave observatories will detect a larger number of binary neutron star (BNS) mergers with increasing accuracy. In this context, improving BNS merger numerical simulations is crucial to correctly interpret the data and constrain the equation of state (EOS) of neutron stars (NSs).
Aims. Stateoftheart simulations of BNS mergers do not include muons. However, muons are known to be relevant in the microphysics of cold NSs and are expected to have a significant role in mergers, where the typical thermodynamic conditions favour their production. Our work is aimed at investigating the impact of muons on the merger remnant.
Methods. We postprocess the outcome of four numerical relativity simulations of BNS mergers performed with three different baryonic EOSs and two mass ratios considering the first 15 milliseconds after merger. We compute the abundance of muons in the remnant and analyse how muons affect the trapped neutrino component and the fluid pressure.
Results. We find that depending on the baryonic EOS, the net fraction of muons is between 30% and 70% the net fraction of electrons. Muons change the flavour hierarchy of trapped (anti)neutrinos such that deep inside the remnant, muon antineutrinos are the most abundant, followed by electron antineutrinos. Finally, muons and trapped neutrinos modify the neutrontoproton ratio, affecting the remnant pressure by up to 7% when compared with calculations neglecting them.
Conclusions. This work demonstrates that muons have a nonnegligible effect on the outcome of BNS merger simulations, and they should be included to improve the accuracy of a simulation.
Key words: dense matter / equation of state / gravitational waves / neutrinos / methods: numerical
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
The detection of the gravitational wave (GW) event GW170817 from the inspiral phase of a binary neutron star (BNS) merger (Abbott et al. 2017a) opened a new window into the microphysics of microphysics of neutron stars (NSs). The measurement of the tidal deformability in addition to the constraint represented by the maximum measured NS mass has allowed researchers to put constraints on the equation of state (EOS) of cold NSs (Abbott et al. 2018; De et al. 2018). More stringent constraints have been obtained by combining information from the GW signal with the properties of the detected electromagnetic counterparts (Abbott et al. 2017b), specifically a short gammaray burst and a kilonova (e.g. Bauswein et al. 2017; Margalit & Metzger 2017; Radice et al. 2018b). The improved sensitivity planned in the next runs of observations by present GW detectors, such as the Laser Interferometer GravitationalWave Observatory (LIGO), the Virgo interferometer, and the Kamioka Gravitational Wave Detector (KAGRA) will increase the number of detected BNS mergers and the chance to observe their electromagnetic counterparts (Abbott et al. 2020; Patricelli et al. 2022; Colombo et al. 2022). The next generation of GW observatories, such as the Einstein Telescope (Punturo et al. 2010) and the Cosmic Explorer (Evans et al. 2021), are expected to significantly enlarge the horizon of detectable BNS mergers and to dramatically improve the estimate of the source parameters (Grimm & Harms 2020; Maggiore et al. 2020; Ronchini et al. 2022; Iacovelli et al. 2022). To date, postmerger GWs from BNS mergers have not been detected, but the Einstein Telescope and the Cosmic Explorer should reach the required sensitivity and largely increase the probability of detecting these signals. Such a discovery would strongly impact our knowledge of the nuclear EOS at high densities, even at a finite temperature (Sekiguchi et al. 2011b; Bauswein et al. 2012; Takami et al. 2014; Weih et al. 2020; Perego et al. 2022; Breschi et al. 2022), which is still quite uncertain. In this context, precise numerical simulations based on accurate modelling of the merger microphysics are mandatory. On the one hand, they are pivotal to interpret the collected data and to help constrain the nuclear EOS in both the zero and finite temperature regimes. On the other hand, they can provide reliable predictions of the properties of electromagnetic counterparts that can be used to optimise multimessenger observational campaigns.
The EOS provides the relation between matter density, temperature, and the thermodynamical variables characterising a certain system. Building generalpurpose EOSs for astrophysical simulations is extremely challenging because of the wide range of densities, temperatures, and charge fractions involved (Oertel et al. 2017). For this reason, it is important to assess the relevant degrees of freedom (d.o.f.) for a given astrophysical system. The thermodynamical conditions of BNS mergers are extreme, with temperatures reaching T ∼ 50 − 100 MeV and densities approaching n_{b} ∼ 3 − 6 n_{0}, where n_{b} and n_{0} are the baryon density and the nuclear saturation density, respectively (Sekiguchi et al. 2011a; Bernuzzi et al. 2016; Radice et al. 2018a; Perego et al. 2019). The d.o.f. usually included in a BNS merger simulation are nucleons, nuclei, electrons, positrons, and photons. Some simulations take into account hyperons, quarks, and/or trapped neutrinos (Sekiguchi et al. 2011a; Foucart et al. 2016a; Most et al. 2019; Radice et al. 2022). However, stateoftheart simulations do not include muons even though the typical merger temperatures and densities allow for their presence.
Muons are already known to play a relevant role in the thermodynamics of protoneutron stars (Prakash et al. 1997) and have a nonnegligible abundance inside cold NSs in neutrinoless βequilibrium (Cohen et al. 1970; Cameron 1970; Glendenning 1997; Haensel et al. 2000, 2001; Steiner et al. 2005; Alford & Good 2010). Their potential role in the context of BNS mergers was investigated by Fore & Reddy (2020) and Alford et al. (2021). For example, it was found that weak reactions involving muons and negatively charged pions at nuclear densities n_{b} ≲ n_{0} and temperatures T ∼ 30 MeV significantly reduce the mean free path of low energy muon (anti)neutrinos. These processes are expected to modify the energy transport inside the merger remnant and to affect its stability (Fore & Reddy 2020). Moreover, muonic Urca processes and (anti)neutrino scattering off muons produce a relevant contribution to the equilibration processes taking place in the postmerger remnant (Alford et al. 2021). Muons have already been included in some corecollapse supernova (CCSN) simulations (Bollig et al. 2017, 2020; Guo et al. 2020; Fischer et al. 2020). The appearance of a significant muonic fraction has been demonstrated after core bounce, and, albeit to a lesser extent, even prior to it (Fischer et al. 2020). In this context, muons induce a softening of the EOS, prompt a burst of muon neutrinos soon after the core bounce, and possibly facilitate the supernova explosion by enhancing the neutrino emission and changing the emitted neutrino spectra (Bollig et al. 2017; Guo et al. 2020; Fischer et al. 2020).
In addition to the nuclear EOS and the related d.o.f., other variables, such as thermal pressure, differential rotation, and rotation rate, influence the evolution of the merger remnant (Kaplan et al. 2014). In particular, thermal pressure affects the angular velocity threshold of mass shedding and the remnant evolution during the phase of secular instability, with nontrivial implications for the time of collapse (Kaplan et al. 2014). Another important contribution to the merger thermodynamics comes from trapped neutrinos, which are present deep inside the merger remnant. Indeed, they are relevant, especially for the implications on the stability of differentially rotating remnants. Their inclusion in BNS merger simulations is extremely challenging and requires a proper radiation transport scheme. To date, only a few BNS simulations have explicitly included them (e.g. Foucart et al. 2016a; Radice et al. 2022). These simulations showed electron antineutrinos are the dominant trapped neutrinos and their appearance results in an increase of the electron and proton fractions. These results are consistent with the outcome of Perego et al. (2019), where the relevance of the trapped neutrino component was quantified within a postprocessing approach. By changing the protontoneutron content, trapped neutrinos induce a pressure decrease of approximately 5%−10% inside the remnant (Perego et al. 2019), so they may speed up the collapse of massive binaries to black holes. Nevertheless, muonic and tauonic (anti)neutrinos have been treated at the same level as heavylepton (anti)neutrinos in previous works because muons were not included in the microphysics modelling. To the best of our knowledge, there are no BNS merger simulations taking into account the coupling between neutrinos and muons.
In the present work, we develop a postprocessing technique to quantify the fractions of muons and trapped neutrinos in the remnants of BNS mergers in order to assess their combined effect on fluid pressure. We apply our analysis to the outcome of four fully relativistic numerical simulations employing three different microphysical EOSs and two different mass ratios. We aim to show that the production of muons and the trapping of neutrinos modify the remnant pressure, which changes asymmetrically in space with respect to the case in which muons and neutrinos are neglected, depending on the assumed baryonic EOS and binary mass ratio.
This paper is structured as follows. In Sect. 2, we briefly discuss the mechanism of muon production and neutrino trapping in BNS mergers and estimate the typical timescales and energies involved. In Sect. 3, we introduce our method based on a postprocessing approach. From Sect. 4 to Sect. 6, we present our results, analysing the abundance of muons in the remnant (Sect. 4), the trapping of neutrinos (Sect. 5), and how muons and trapped neutrinos affect the remnant pressure (Sect. 6). In Sect. 7, we discuss our results, comparing them to stateoftheart simulations. Finally, we summarise our work and present our conclusions in Sect. 8.
2. Preliminary estimates
2.1. Muons in BNS mergers
The fraction of muons in a BNS remnant is determined both by the fraction of muons in the two cold NSs before merger and by the efficiency of (anti)muon production during and after merger. We consider as an example a fluid element taking part in a BNS merger. During the inspiral, muons are already present in the two cold NSs for sufficiently high densities when the electron chemical potential, μ_{e}, exceeds the muon rest mass, m_{μ}c^{2} ∼ 106 MeV (Cohen et al. 1970; Cameron 1970; Glendenning 1997; Haensel et al. 2000, 2001; Steiner et al. 2005; Alford & Good 2010). For relativistic degenerate electrons, the former can be estimated as:
where p_{F, e} is the Fermi impulse of electrons, while Y_{e−} = n_{e−}/n_{b} and n_{e−} are the electron fraction and the electron number density, respectively. Clearly, for neutrinoless βequilibrated nuclear matter, Y_{e−} ∼ 0.05 − 0.1, muons are expected to be present as nonrelativistic degenerate particles already close to saturation density. A simple estimate of their chemical potential, μ_{μ−}, is given by:
where p_{F, μ} is the Fermi impulse of muons, while Y_{μ−} = n_{μ−}/n_{b} and n_{μ−} are the fraction of muons and the muon number density, respectively.
We next investigate whether the typical thermodynamical conditions in BNS mergers allow for the creation of muons. During the merger, nuclear matter is heated and the average photon energy, E_{γ} ≈ 2.7 k_{B}T, overcomes the muon mass threshold when the temperature reaches k_{B}T ≳ 40 MeV. Therefore, in this temperature regime, thermal processes such as
drive the creation of μ^{±} pairs. Additionally, neutrino pairs of all flavours are produced via thermal processes (see e.g. Haft et al. 1994; Ruffert et al. 1996)
or nucleonnucleon bremsstrahlung (see e.g. Hannestad & Raffelt 1998)
where ν denotes any flavour neutrino, while N is a nucleon. For example, in the case of e^{±} annihilation, the average energy of the neutrino pairs is (Cooperstein 1988):
where η_{e−} = μ_{e−}/k_{B}T is the (relativistic) electron degeneracy parameter and F_{n}(η) is the Fermi function of order n. For very degenerate electrons (η_{e−} ≫ 1) and accounting for Eq. (1)
Even when muonic (anti)neutrinos thermalise, their average energy is
Thus, once high energy ν_{μ} and are formed, (anti)muons can be created via chargedcurrent (CC) semileptonic reactions
as in the case of CCSNe (Bollig et al. 2017; Guo et al. 2020; Fischer et al. 2020). In addition, neutron decay, , plays a relevant role in the high densitylow temperature regime (Alford et al. 2021) typical of the remnant core. Finally, purely CC leptonic processes, such as ν_{μ} + e^{−} → ν_{e} + μ^{−}, and inverse muon decay can enhance the abundance of muons at low average neutrino energies ≲50 MeV since the amount of electrons exceeds the amount of positrons (Guo et al. 2020; Fischer et al. 2020; Alford et al. 2021).
In general, matter in BNS mergers is not in weak and thermal equilibrium, and it is not obvious that weak reactions are fast enough to allow for muon creation at short enough timescales everywhere inside the remnant (see e.g. Hammond et al. 2022). We estimate the rate at which these reactions occur by considering the reaction ν_{μ} + n → p + μ^{−} as a representative example. In the zero momentum transfer limit, neglecting the momentum transferred to nucleons and assuming p_{N}≪m_{N}, where p_{N} and m_{N} are the nucleon momentum and mass, respectively, the reaction rate of ν_{μ} + n → p + μ^{−} is at leading order:
where G is the Fermi constant, g_{V} and g_{A} are the weak couplings, Q = m_{n} − m_{p} is the neutronproton mass difference, f_{i}(E_{i}) is the FermiDirac distribution function of particle i with energy E_{i}, and ω_{np} is given by
where n_{N} and indicate the number density and the nonrelativistic chemical potential of nucleons N = p, n, respectively. An order of magnitude estimate of the reaction rate can be easily obtained in the nondegenerate regime, ω_{np} → n_{n}, by neglecting the neutronproton mass difference^{1} and considering neutrino energies E_{νμ} ≃ μ_{μ−} ≳ m_{μ}c^{2} (see Eq. (2)). The inverse rate thus gives the muon production timescale:
We note that for high enough densities and neutrino energies where t_{dyn} ∼ 10^{−4} s is the merger dynamical timescale (see e.g. Radice et al. 2020). Because of the high degeneracy of neutrons, CC semileptonic reactions are expected to favour the production of μ^{−} over μ^{+}. Accordingly, the opacity of ν_{μ} to CC semileptonic reactions exceeds that of (at least in cases with high enough neutrino energy; see also the discussion in Guo et al. 2020; Fischer et al. 2020), and the net muon fraction becomes enhanced.
2.2. Neutrino trapping
In this section, we provide an estimate of the neutrino diffusion timescale from a BNS merger remnant to explore the conditions in which neutrinos can be considered trapped, that is, when their diffusion timescale, t_{diff}, exceeds the dynamical timescale, as well as the conditions in which trapped neutrinos are in thermal and weak equilibrium with matter. Since neutrino scattering off nucleons is among the largest sources of opacity for all neutrino flavours, we considered an approximate expression of its corresponding mean free path, λ_{νN}(E_{ν})∼1/(n_{b}σ_{νN}(E_{ν})), where
gives the ν − N scattering cross section and σ_{0} ≈ 1.76 × 10^{−44} cm^{2} (Shapiro & Teukolsky 1983). Using randomwalk arguments, we estimated t_{diff} as
where ℓ is the characteristic length scale of the diffusion process and τ is the optical depth, which we estimated as τ ∼ ℓ/λ. For example, in the case of thermal neutrinos diffusing from a remnant core of size R_{NS},
and R_{NS} and T_{NS} are the characteristic radius and temperature of the massive NS remnant, respectively (see e.g. Bernuzzi 2020). We can repeat the estimate for thermal neutrinos diffusing from the innermost part of the disc surrounding the massive central remnant:
where H_{disc} and T_{disc} are the characteristic disc height and temperature, respectively. We note that the average energy of thermal neutrinos at k_{B}T ∼ 5 MeV is E_{ν} ∼ 15 MeV, which corresponds to the typical energy of the neutrinos emitted at infinity by BNS merger remnants (e.g. Ruffert & Janka 1998; Rosswog & Liebendoerfer 2003; Sekiguchi et al. 2016; Foucart et al. 2016a,b; Cusinato et al. 2022). As implied by Eq. (12), the diffusion timescale of neutrinos with energy E_{ν} ∼ 15 MeV from the disc becomes comparable to the dynamical timescale around n_{b} ∼ 5 × 10^{−3} fm^{−3}, corresponding to a rest mass density ρ = m_{amu}n_{b} ∼ 10^{12} g cm^{−3}, where m_{amu} is the atomic mass unit. Other interactions involving neutrinos in addition to ν − N scattering provide further opacity. Thus, the above estimates can be understood as conservative upper limits for the trapping density.
Due to the presence of inelastic processes (such as absorption of neutrinos on free nucleons, inverse pair processes, and scattering off electrons and positrons), neutrinos are not only efficiently trapped, but they couple to matter and photons if the corresponding mean free path is short enough. The situation is however complicated by the fact that neutrinos have a spectral distribution and that cross sections have a nontrivial energy dependence. As a result, neutrinos of different flavours and energies decouple from matter at different locations. Rest mass density has been shown to be the most relevant matter property in determining the location of both the last scattering surface and the surfaces where weak and thermal equilibrium freezes out (Endrizzi et al. 2020). In particular, according to Endrizzi et al. (2020), we observe that for densities larger than 10^{12} g cm^{−3}, almost all relevant neutrinos are within the last scattering surface (Fig. 8), while for densities larger than 10^{13} g cm^{−3}, they are also in thermal and weak equilibrium (Fig. 10). In the case of ν_{e} and , these surfaces are even characterised by densities one order of magnitude smaller. Thus, we defined a limiting density ρ_{lim} such that neutrinos can be considered as a trapped gas if ρ ≳ ρ_{lim}. For electron (anti)neutrinos, we considered ρ_{lim, e} = 10^{11} g cm^{−3}, while for muonic and tauonic (anti)neutrinos, we set ρ_{lim, x} = 10^{12} g cm^{−3}. Moreover, at ρ ≳ 10 × ρ_{lim} neutrinos can be approximately modelled as a gas in equilibrium with other constituents of matter.
3. Method
We present our method to estimate the fraction of muons and their impact on the trapped neutrino component in BNS merger remnants by postprocessing the outcome of a significant sample of numerical relativity simulations. First, we review the properties of the simulations considered in this work (Sect. 3.1). Next, we describe how we model the nuclear EOS, including muons and trapped neutrinos (Sect. 3.2). Finally, we explain our postprocessing technique by arguing for the physical rationale and the limits of validity (Sect. 3.3).
3.1. The simulation sample
We considered the outcomes of four BNS merger simulations in numerical relativity targeted to GW170817 (see Table 1) published in Nedora et al. (2019, 2021), and Bernuzzi et al. (2020). The term M_{A, B} refers to the gravitational mass of each of the two NSs at infinity such that the binary mass ratio is defined as q = M_{B}/M_{A} ≥ 1. More details about the simulations and their numerical setup can be found in the references listed in Table 1 and in Radice et al. (2018a). In the following paragraphs, we provide a summary of the simulation properties necessary to understand our postprocessing procedure.
List of the BNS merger simulations considered in this work.
Our simulation sample contains three symmetric (q = 1) simulations. These simulations were performed with three different nuclear EOSs, namely, BLh (Bombaci & Logoteta 2018; Logoteta et al. 2021), HS (DD2; Typel et al. 2010; Hempel & SchaffnerBielich 2010), and SFHo (Steiner et al. 2013). In the following sections, we refer to the second EOS simply as DD2. Details about these nuclear EOSs are provided in Sect. 3.2.1. These simulations predicted quite different merger outcomes. In particular, for SFHo the remnant collapsed into a black hole within a few milliseconds after the merger, while it survived for at least 90 ms in the other two cases. In addition, we considered a simulation of a significantly asymmetric BNS merger (q = 1.43) using the BLh EOS. Within this sample, we could explore how the presence of muons and trapped neutrinos correlates with the properties of the nuclear EOS and with the mass asymmetry degree of the binary. Additionally, half of the simulations included the effect of physical viscosity of magnetic origin (Radice et al. 2016). The simulation sample is however too small to address the impact of viscosity on our results, and we leave such an investigation to future studies. In the following sections, we refer to each simulation by specifying the nuclear EOS and q according to the notation: (EOS, q).
All simulations modelled NS matter as made of neutrons, protons (both free and bound in nuclei), photons, electrons, and positrons. No muons and trapped neutrinos were considered. All the simulations included the effect of neutrino radiation through a leakage plus M0 scheme, whose details can be found in Radice et al. (2016, 2018a). Electron neutrinos and antineutrinos were considered separately, while there was no distinction between muonic and tauonic (anti)neutrinos, which were treated as a single species (heavylepton neutrinos). The leakage prescription was used to compute the net rate of change in the internal energy and in the lepton fraction for matter in optically thick conditions due to neutrino diffusion. Even though the particle and energy density of trapped neutrinos were estimated by the leakage algorithm in order to compute the neutrino diffusion rates, they were not dynamically evolved, and their contribution to the thermodynamical state of matter was not explicitly taken into account in the simulations. This also implies that the energy density and pressure of the neutrinos were not included in the stressenergy tensor of matter and radiation.
In our analysis, we considered snapshots of the computational domain from the postmerger phase in the time interval ∼5 − 15 ms postmerger where the merger time is defined by the maximum in the ℓ = m = 2 mode in the GW waveform. In this phase, the merger remnant is still in the GWdominated phase (e.g. Bernuzzi et al. 2016; Zappa et al. 2018), but the initial, highly dynamical transients that characterise the very first milliseconds after the merger have disappeared and a disc around the central remnant has formed. From each snapshot, we extracted the full 3D profile of the local baryon number density , net electron fraction , and energy density e_{sim}, as directly obtained by the simulation. The net electron fraction was defined as Y_{e} = (n_{e−} − n_{e+})/n_{b}, where n_{e−}(n_{e+}) is the electron (positron) number density. By using the same EOS used in the simulation, we could reconstruct the full thermodynamic states anywhere inside the computational domain, including the values of the temperature T_{sim}, chemical potentials , particle fractions , partial pressures , and total pressure P_{sim}, where the subscript j = {n, p, e^{±}, γ} indicates neutrons, protons, electrons, positrons, and photons. In Fig. 1, we show as an example the rest mass density , the temperature, and the electron fraction from the outcome of the simulation (BLh, 1.00) at approximately 5 ms after merger. The dense and cold core of the remnant can be recognised, characterised by ρ_{sim} ≳ 10^{15} g cm^{−3}, T_{sim} < 20 MeV, and . The core is surrounded by warm matter in the density regime ρ_{sim} ∼ 10^{14} − 10^{15} g cm^{−3}, with the temperature reaching up T_{sim} ∼ 50 MeV in correspondence to the hot spots and a slightly smaller electron fraction . At a lower density ρ_{sim} ≲ 10^{14} g cm^{−3}, cold and warm matter streams alternate. In the following sections, we explain how these data are used to postprocess the original simulations.
Fig. 1. Outcome of the numerical simulation (BLh, 1.00) at 4.6 ms after merger. The left, central, and right panels show the rest mass density, the temperature, and the electron fraction on the equatorial plane, respectively. The black lines mark isodensity contours corresponding to 10^{12} g cm^{−3} (solid line), 10^{13} g cm^{−3} (dashed line), 10^{14} g cm^{−3} (dasheddotted line), and 10^{15} g cm^{−3} (dotted line). 
3.2. EOS modelling
In our analysis, we considered baryons, massive leptons, photons, and neutrinos as the relevant d.o.f. to describe matter and radiation in BNS mergers. Baryons, massive leptons, and photons are assumed to be in thermal and nuclear statistical equilibrium everywhere inside the remnant. Nuclear statistical equilibrium is reached when strong and electromagnetic reactions are in equilibrium, and the nuclear composition is provided by the minimum of the free energy. For astrophysically relevant conditions, this is verified as long as k_{B}T ≳ 0.6 MeV (e.g. Hix & Thielemann 1999). In astrophysical plasma, the negative net charge of massive leptons is balanced by the positive charge of protons. In the presence of more than one massive lepton, charge neutrality implies:
where Y_{p} = n_{p}/n_{b} is the proton fraction; n_{p} is the proton number density; Y_{l} = (n_{l−} − n_{l+})/n_{b} is the net fraction of the massive lepton l; and e, μ, and τ stand for electron, muon, and tauon, respectively.
Every relevant species was characterised by its own EOS. These EOSs can be expressed in terms of any thermodynamical potential, such as the Helmholtz free energy F with the volume V, the number of particles N, and the temperature T, as independent variables. Neglecting the modifications of the thermodynamical potentials due to the interaction between different d.o.f., the total free energy was given by
where the subscript i = {b, l, γ, ν} indicates baryons, massive leptons, photons, and neutrinos, respectively. All the thermodynamic variables, such as the chemical potentials, μ_{i}, the energy densities, e_{i}, and the pressures, P_{i}, were derived from F_{i} at fixed temperature T, number of particles N_{i}, and volume V according to the standard rules of the canonical ensemble. Hence, the total energy density and pressure were simply given by
In the following sections, we provide a detailed description of the EOSs used in this work for each species.
3.2.1. Baryons
Our knowledge of the baryonic matter EOS is still affected by large uncertainties. Therefore, to bracket possible uncertainties in this work, the baryonic contribution was taken into account via three different nuclear EOSs in tabulated form corresponding to the ones used in the simulations presented in Sect. 3.1: the BLh (Bombaci & Logoteta 2018; Logoteta et al. 2021), the DD2 (Typel et al. 2010; Hempel & SchaffnerBielich 2010), and the SFHo (Steiner et al. 2013)^{2}. We note that we did not consider the possible formation of hyperons (Oertel et al. 2012; Fortin et al. 2018) or a phase transition to quark matter (Weissenborn et al. 2011; Klähn et al. 2013; Chatterjee & Vidaña 2016; Bombaci et al. 2016; Logoteta et al. 2019; Logoteta 2021). The BLh EOS is a microscopic EOS constructed following the framework of the BruecknerBetheGoldstone manybody theory extended at finite temperatures within the BruecknerHartreeFock approximation. In contrast, the modelling of nuclear interactions in DD2 and SFHo is based on a relativistic mean field theory. For each nuclear EOS, the independent variables are (n_{b}, T, Y_{p}). The range of (n_{b}, T, Y_{p}) for each EOS table is reported in Table 2. The three EOSs are in reasonable agreement with the experimental constraints on the properties of nuclear matter at n_{0} and T = 0 (for experimental constraints, see Shlomo et al. 2006; Danielewicz & Lee 2014; Oertel et al. 2017; Drischler et al. 2019). Moreover, all three EOSs satisfy to a good extent present astrophysical constraints on the NS maximum mass M_{max} > 2.01 M_{⊙} and NS radius (Antoniadis et al. 2013; Lattimer & Steiner 2014; Oertel et al. 2017; Riley 2019; Miller et al. 2019, 2021; Raaijmakers et al. 2021).^{3} In Table 2, we report for each EOS the values of the nuclear saturation density, the maximum mass M_{max} of a cold nonrotating NS, the corresponding radius R_{max}, the radius R_{1.4} of a 1.4 M_{⊙} NS, and the dimensionless tidal deformability (e.g. Hinderer 2008; Damour et al. 2012; Favata 2014) for GW170817 targeted binaries. We note that SFHo supports a smaller maximum NS mass with a smaller radius. More generally, SFHo produces more compact NSs compared to the other two EOSs. In the case of BLh, it predicts slightly larger M_{max} and R_{max}, but the values are still close to the ones from SFHo. Conversely, DD2 exhibits significantly larger values for M_{max} and R_{max} and a smaller compactness. Despite the fact that the maximum NS properties of SFHo and BLh are relatively close, it is worth noting that for other properties at densities equal to several times n_{0}, SFHo and DD2, regarding the symmetry energy, and DD2 and BLh, for incompressibility, are more similar. With this selection of baryonic EOSs, we could explore very different microphysical descriptions and a wide range of thermodynamical conditions and macroscopic properties of merging BNSs.
Properties of nuclear EOS tables used in this work and of cold, spherically symmetric NSs for each nuclear EOS.
3.2.2. Massive leptons and photons
We treated massive leptons as ideal Fermi gases in thermal equilibrium. In particular, we incorporated electrons and positrons e^{±} as well as muons and antimuons μ^{±}, while we neglected tauons and antitauons τ^{±} due to their large mass compared to the typical temperatures and Fermi energies in the system. The EOS of a noninteracting gas^{4} of fermions at a finite temperature is well known in the literature, and it provides the expressions for the number density n_{l±}, the specific energy density e_{l±}, and the pressure P_{l±} (Bludman & van Riper 1977):
where m_{l} is the lepton mass and K_{l} is a constant
θ_{l} and are the relativity and the nonrelativistic degeneracy parameters, respectively,
while are the generalised Fermi functions of order k. At thermal equilibrium, the degeneracy parameter of antiparticles is related to that of particles (Timmes & Arnett 1999):
By inverting Eq. (16), the term can be expressed as a function of (n_{b}, Y_{l}, T), as was suggested in Timmes & Arnett (1999).
Photons form an ideal Bose gas in thermal equilibrium with matter and with zero chemical potential, μ_{γ} = 0. The resulting EOS of photons depends only on T:
and P_{γ} = e_{γ}/3.
3.2.3. Trapped neutrinos
Based on the discussion in Sect. 2 and as typical neutrino energies inside the remnant are much larger than any neutrino mass, we modelled the trapped neutrino component as a massless Fermi gas in weak and thermal equilibrium with matter and radiation. Accordingly, the particle number density n_{ν}, the energy density e_{ν}, and the pressure P_{ν} of any flavour neutrino ν were given by
where the exponential factor e^{−(ρlim/ρ)} was introduced in order to model the fading of the trapped component at ρ ≲ ρ_{lim} (see also Kaplan et al. 2014 and Perego et al. 2019 for a similar choice). The degeneracy parameter of antineutrinos was fixed by because of thermal equilibrium. Furthermore, η_{ντ} = 0 since the fraction of tauons is negligible, while weak equilibrium defined the degeneracy parameters of electronic and muonic neutrinos:
where and are the nonrelativistic degeneracy parameters of neutrons and protons, respectively. Accordingly, the EOS of electron (muon) neutrinos depends on T, Y_{p}, and Y_{e} (Y_{μ}), while the EOS of tauon neutrinos depends only on T and does not distinguish between ν_{τ} and . The fraction of neutrinos and their mean energy were defined as Y_{ν} = n_{ν}/n_{b} and E_{ν} = e_{ν}/n_{ν}, respectively.
From the previous discussion, if all the species are assumed to be in thermal and reaction equilibrium, the set of variables needed to describe the complete EOS is given by (n_{b}, T, Y_{p}, Y_{e}, Y_{μ}). However, the charge neutrality of NS matter constrains the proton fraction Y_{p} = Y_{e} + Y_{μ}, reducing the set of independent variables to (n_{b}, T, Y_{e}, Y_{μ}).
3.3. The postprocessing technique
Based on the simulation outcome discussed in Sect. 3.1, we estimated the amount and the impact of muons in postprocessing as well as their influence on the trapped neutrino properties. In particular, we were interested in the implications for the remnant on a timescale Δt ∼ 5 − 15 ms after the merger and in the density region ρ ≳ 10^{13} g cm^{−3} where the approximation of weak and thermal equilibrium better applies. The absence of muons inside the original simulations posed a significant issue and required some modelling (shown in Sect. 3.3.1) before the postprocessing technique could be applied (discussed in Sect. 3.3.2).
3.3.1. Modelling
We modelled the NS fluid elements before the merger as being composed of neutrons, protons, electrons, and muons in cold neutrinoless weak equilibrium:
In Fig. 2, we show the fractions of electrons, muons, and protons as a function of n_{b}, obtained by solving Eq. (27). At n_{b} ≈ 0.3fm^{−3} ∼ 2 n_{0}, the muon fraction is approximately 40% of the electron fraction for all the considered EOSs, and it increases monotonically with n_{b}. For comparison, we also show the electron fraction calculated by assuming βequilibrium without muons (see dashed black line in Fig. 2), which corresponds to the initial conditions used in the simulations presented in Sect. 3.1. A common feature of all three models is that when muons are included, the electron fraction becomes smaller, while the proton fraction tends to increase, that is, Y_{p} = Y_{e} + Y_{μ}. In particular, the system lowers the Fermi energy of electrons by converting electrons into muons up to the point where the chemical potentials of electrons and muons are equal, balancing the difference between the neutron and proton chemical potentials. Moreover, when muons are included, neutrons can turn into protons plus muons via Urca processes (Leinson 2002; Yakovlev & Pethick 2004; Potekhin et al. 2015). The three EOSs show pronounced differences in the density behaviour of the electron fraction. In particular, BLh predicts a higher Y_{p} at large densities compared to DD2 and SFHo. Such differences are mainly due to the different density behaviours of the symmetry energy in the three models. Symmetry energy is indeed the principal source that determines matter composition at a given density and temperature. In the following paragraphs, we refer to the particle fractions of the cold NSs obtained by solving Eq. (27) as and which are represented as solid lines in Fig. 2, and to the net electron fraction computed in cold weak equilibrium without muons (see dashed black line in Fig. 2) as .
Fig. 2. Particle fractions computed at T = 0 in neutrinoless βequilibrium for three different EOS: BLh (left), DD2 (centre), and SFHo (right). Black lines show the electron fraction computed in two different cases: equilibrium with muons (solid lines) or without muons (dashed lines). The green line corresponds to the fraction of protons at equilibrium with muons, while the red line corresponds to the muon fraction. We note that when muons are not included, the proton fraction coincides with the electron one (dashed black line). 
During a merger, the density, temperature, and composition of the fluid elements change because of matter compression and decompression, and neutrinos are produced. If a fluid element is in good approximation within the last scattering neutrino surface (i.e. ρ ≳ ρ_{lim}), a trapped neutrino gas is formed in weak and thermal equilibrium with matter. In these conditions, it is convenient to introduce for each fluid element the electron, muon, and tauon lepton numbers
and the total fluid energy
The central part of the remnant, characterised by ρ ≳ 10^{14} g cm^{−3}, is formed by the fusion of the two NSs cores, and it is characterised by fluid elements whose density is significantly in excess of the neutrino surface density during their evolution. Such fluid elements are not affected by shocks and significant neutrino emission. In particular, as discussed in Sect. 2, the neutrino diffusion timescale, t_{diff}, is in the order of seconds; thus t_{diff} ≫ Δt ∼ 5 − 15 ms (see Eq. (11)). So these fluid elements conserve their electron and muon lepton fractions over Δt, which are ultimately equal to the electron and muon lepton fractions set by the initial cold neutrinoless equilibrium, Eq. (27)^{5} A close inspection of the thermodynamic conditions experienced by matter in a BNS merger event and of their corresponding evolution (Perego et al. 2019) revealed that the conservation of the lepton numbers is well realised for matter originally inside the outer core of the two cold NSs, that is, at ρ ≳ 10^{14} g cm^{−3} also during the inspiral. In contrast, most of the matter at a lower density inside the merger remnant, ρ < 10^{14}g cm^{−3}, originates from the decompression of matter originally around ρ ∼ 10^{14} g cm^{−3} inside the two merging NSs. In this case, matter is significantly heated by shocks and compression such that the electron and muon lepton numbers of each fluid element become considerably altered. In this density domain, the electron fraction evolved through the numerical relativity simulations is a better proxy of Y_{l, e}: . The muon lepton fraction is instead approximately given by Y_{l, μ} ∼ 0.01, which corresponds to the typical muon fraction in the original NSs at the edge of the outer core, that is, around and just below saturation density (see Fig. 2).
Concerning the relativistic internal energy, no significant energy leaks out in the form of neutrinos from the fluid elements until t_{diff} ≫ Δt in neutrino trapping conditions. Thus, the evolution of e_{sim} inside the simulation is also expected to reproduce in good approximation the evolution (i.e. the variation) of the total internal energy once muons and neutrinos are present. However, the presence of physical muons inside the two cold NSs alters the absolute value of the relativistic internal energy due to their rest mass contribution. Nevertheless, this contribution remains approximately constant in time as long as the bulk of the muons present in the remnant comes from the two cold NSs.
3.3.2. Algorithm
Based on the above considerations, we designed the following postprocessing algorithm. We considered the time snapshots of the simulations listed in Table 1 in the time interval ∼5 − 15 ms. From each snapshot, we read the full 3D profile of the baryon number density , the net electron fraction , and the internal energy density e_{sim} on a Cartesian grid (x, y, z) with −30 km ≤ x, y, z ≤ +30 km. Given on the Cartesian grid and for each snapshot, we performed our analysis according to the following steps.
Given , we computed the baryon density such that (See Sect. 3.3.1 and dashed black line in Fig. 2.) We note that corresponds to the original density of the fluid element inside the NSs during the inspiral, and it is, in general, different from .
We computed the cold equilibrium electron and muon net fractions corresponding to (i.e. and (See Sect. 3.3.1 and solid black and red lines in Fig. 2.)
Given the simulation baryon number density , we computed . If ρ_{sim} > 10^{14} g cm^{−3}, then the lepton numbers were separately conserved (see discussion in Sect. 3.3.1). Therefore, we identified the electron and muon lepton fractions with
If instead ρ_{sim} < 10^{14} g cm^{−3}, we imposed
In accordance with our assumptions, Y_{l, τ} = 0 everywhere and at any time. We note that the proton fraction (which in the simulation is equal to ) is now larger than since .
Given the simulation energy e_{sim}, if ρ_{sim} > 10^{14} g cm^{−3}, we computed the total internal energy density of each fluid element as
while if ρ_{sim} < 10^{14} g cm^{−3}, we computed
The terms proportional to correspond to the rest mass energy density of the muons already present in the two cold NSs.
We imposed , and we solved the system
with respect to (T, Y_{e}, Y_{μ}), where and j = {e, μ, τ}, to find the equilibrium configuration in the presence of (anti)muons and trapped neutrinos. The (anti)neutrino fractions, Y_{ν}, and energies, e_{ν}, are defined in Sect. 3.2.3. We note that the condition on the tauon lepton number, Y_{l, τ} = 0, is automatically satisfied in the absence of physical tauons and by the resulting η_{ντ} = 0 assumption. Nevertheless, the contribution of tauon (anti)neutrinos was taken into account in the calculation of the total energy and pressure.
We computed all the other thermodynamic variables, such as fractions of trapped neutrinos, chemical potentials, and pressures as functions of (n_{b}, T, Y_{e}, Y_{μ}), as prescribed in Sect. 3.2.
We note that in our procedure, we did not prescribe any time evolution. Rather, we postprocessed each available time snapshot of the original simulation in the interval ∼5 − 15 ms after merger.
3.3.3. Limits of validity
Our postprocessing analysis allowed us to compute, a posteriori, the thermodynamics of the merger aftermath taking into account the presence of muons and trapped neutrinos. We note that since there are no muons nor any explicit neutrinos in the original simulations, the definition of Y_{l, e}, Y_{l, μ}, and e_{tot} based on the simulation outcome cannot be done in a unique and fully consistent way. The choice described in the previous paragraph has the advantage of including the contribution given by the muons already present in the cold NSs. However, it relies on the assumption that the merger dynamics and e_{sim}, once corrected for the muons rest mass energy, are not significantly altered by them. To bracket the uncertainties due to this choice, in Appendix A, we discuss our results for a different choice of Y_{l, e}, Y_{l, μ}, and e_{tot} where the initial muon fraction from the cold NSs is neglected.
In our modelling of neutrinos, we assumed that above ρ_{lim, e}, there is a full FermiDirac distribution of trapped neutrinos. However, neutrinos with energy smaller than 10 − 20 MeV were not expected to be in equilibrium with matter in this density regime. Moreover, from Eq. (12) it follows that the diffusion timescale of neutrinos with energy E_{ν} ∼ 15 MeV becomes comparable to Δt ∼ 10 ms around ρ ∼ 10^{12} g cm^{−3}. For this reason, even though we applied our postprocessing procedure in the density domain ρ > ρ_{lim, e}, we restricted our analysis of the results to ρ > 10^{13} g cm^{−3}, where the assumption of neutrino trapping and weak equilibrium is robust, being t_{diff} ≫ t_{dyn}, and the postprocessing approach is more reliable, being t_{diff} ≫ Δt. We also note that our modelling of the trapped neutrino component closely follows the one used by Perego et al. (2019). Recent numerical results obtained by including neutrino transport through a grey M1 moment scheme in BNS merger simulations (Zappa et al. 2023) are very consistent with the outcome of the postprocessing analysis of Perego et al. (2019).
Finally, we highlight that in our procedure, the population of antimuons is given by a FermiDirac distribution in thermal equilibrium. Therefore, antimuons arise as a thermal tail at a high enough temperature.
4. The fraction of muons in the remnant
We started our postprocessing analysis from the simulation (BLh, 1.00). In Fig. 3, we show the fractions of μ^{±}, ν_{μ}, and as well as the muon chemical potential approximately 5 ms after merger. The fraction of muons, Y_{μ−}, is ∼0.01 − 0.07 in the highdensity region characterised by ρ > 10^{14} g cm^{−3}, while it decreases down to ∼10^{−3} − 10^{−2} at the lower density ρ ∼ 10^{13} − 10^{14} g cm^{−3}. By computing the difference (upperright panel of Fig. 3 ), we deduced that the bulk of muons comes from the two cold NSs. This result is consistent with the hypothesis of our postprocessing procedure (see Sect. 3.3.1). However, a fraction of the muons was created during the merger, corresponding to ∼10^{−3} − 10^{−2}. Muon creation became enhanced in correspondence to the hot spots where k_{B}T_{sim} ≳ 50 MeV (see the second panel of Fig. 1) and at high density, that is, above 10^{15} g cm^{−3}. By inspecting the change in temperature and particle fractions
Fig. 3. Particle fractions and chemical potential in the muon flavour sector computed by postprocessing the simulation (BLh, 1.00) at 4.6 ms after merger. In the upper panels, we show the fraction of μ^{−} (left), the fraction of μ^{+} (centre), and the difference between the muon fraction and the one inherited from the cold NSs, that is, (right). In the lower panels, we show the fraction of ν_{μ} (left), the fraction of (centre), and the muon chemical potential (rest mass subtracted) (right). As in Fig. 1, the black lines mark isodensity contours. 
where T and Y_{i} result from postprocessing calculations, while T_{sim} and are read from the simulation, we could guess which were the processes more likely to be involved in muon creation. In the hot spots, both muons and antimuons were produced, and the fraction of newly created muons Y_{μ−} ∼ 10^{−2} exceeded the fraction of antimuons Y_{μ+} ∼ 10^{−3} (see Fig. 3, uppercentre and right panels). At the same time, we observed ΔT < 0, ΔY_{n} < 0 and ΔY_{e−} > 0. We concluded that, on the one hand, thermal processes drive the creation of μ^{±} pairs at the expense of the system’s internal energy, and on the other hand, the reduction in the degeneracy of muons due to the high temperature (see Fig. 3, lowerright panel) favours the conversion of n into p + μ^{−}, further enhancing the production of μ^{−}. In addition, we found an enhancement in the creation of e^{±} pairs and e^{−}, but the amount of newly created μ^{−} systematically exceeded that of the newly created e^{−}. In the core of the remnant, the density of fluid elements was enhanced by matter compression, and μ^{−} as well as e^{−} were created at the expense of the neutrons’ degeneracy (ΔY_{n} < 0). Finally, at a lower density, ρ < 10^{14} g cm^{−3}, the initial muon fraction from cold, decompressing NSs matter was almost entirely converted into muon neutrinos. We note that the values of Y_{μ−} are stable in the time domain of our analysis.
Next, we compare the results obtained from the simulation (BLh, 1.00) with the ones from (DD2, 1.00) and (SFHo, 1.00). In the latter cases, muons were also present in the fulldensity region, ρ > 10^{13} g cm^{−3}. However, the muon fraction Y_{μ−} ∼ 10^{−3} − 0.05 presented a slightly smaller maximum with respect to (BLh, 1.00). This is expected since, as discussed in Sect. 3.3.1, DD2 and SFHo predict a smaller muon fraction for cold nuclear matter in weak equilibrium (see also Fig. 2). The amount of muons created during the merger (i.e. is comparable for the three simulations and correlates in the same way with ΔT and ΔY_{i}, defined in Eq. (35). The amount of μ^{−} and produced in the centre of the remnant, however, was smaller for (DD2, 1.00) and (SFHo, 1.00) with respect to (BLh, 1.00), with a difference that can reach one order of magnitude if we compare (DD2, 1.00) to (BLh, 1.00). As in the case of the cold, neutrinoless weak equilibrium conditions presented in Fig. 2, such a difference between the BLh EOS and the DD2 EOS is mostly due to the larger slope of the symmetry energy in the former case. Additionally, the fraction of μ^{−} was on average slightly larger in (SFHo, 1.00) than in (DD2, 1.00) because the softer EOS SFHo exhibits larger temperatures and densities than the stiffer DD2 EOS (see e.g. Foucart et al. 2016a; Sekiguchi et al. 2016; Radice et al. 2018a; Perego et al. 2019).
To assess the relevance of muons in the merger remnants, we compared the equilibrium net muon fraction obtained in our postprocessing procedure with the equilibrium net electron fraction by analysing the ratio Y_{μ}/Y_{e}, as shown in Fig. 4. We note that for all the simulations, the net muon fraction was at least 30%−40% of the net electron fraction. The ratio Y_{μ}/Y_{e} was the largest for (BLh, 1.00), where it can reach a maximum of approximately 0.7, while it reached up to 0.6 and 0.5 for (SFHo, 1.00) and (DD2, 1.00), respectively. We stress, however, that due to the large mass difference between electrons and muons, the amount of both electrons and positrons was largely increased where the temperature was high such that Y_{e−} ≳ Y_{e+} and Y_{e} = Y_{e−} − Y_{e+} ∼ 0.1, while Y_{μ−} ≫ Y_{μ+} and Y_{μ} ≈ Y_{μ−}.
Fig. 4. Ratio of the net muon fraction over the net electron fraction, Y_{μ}/Y_{e}, on the equatorial plane obtained by postprocessing simulations (BLh, 1.00) (upper left), (BLh, 1.43) (upper right), (DD2, 1.00) (lower left), and (SFHo, 1.00) (lower right) at 5 − 6 ms after merger. As in Fig. 1, the black lines mark isodensity contours. 
Finally, we postprocessed data from the simulation (BLh, 1.43), which had a significantly larger mass asymmetry. We did not observe notable differences while comparing the maximum value of Y_{μ}/Y_{e} in (BLh, 1.00) and (BLh, 1.43). However, (BLh, 1.43) exhibited a more extended spatial region where Y_{μ}/Y_{e} ∼ 0.5 − 0.6 (see Fig. 4). This region is characterised by a large temperature, ≳40 MeV, because it embeds the core of the secondary NS, which was broadened and strongly heated by compression and shocks during the merger.
5. Trapped neutrino properties
In Figs. 3 and 5, we show all flavour (anti)neutrino abundances as well as the chemical potentials of baryons and massive leptons computed by postprocessing the simulation (BLh, 1.00) at approximately 5 ms after merger. Trapped gases of (anti)neutrinos were present in the fulldensity regime of our analysis, ρ > 10^{13} g cm^{−3}, but the abundances and the species hierarchy varied according to the system density and temperature.
Fig. 5. Fractions of ν_{e}, , ν_{τ}, and chemical potentials of electrons, protons, and neutrons computed in postprocessing for (BLh, 1.00) at 4.6 ms after merger. In the upper panels, we show the fraction of ν_{e} (left), the fraction of (centre), and the fraction of ν_{τ} (right). In the lower panels, we show the electron chemical potential (rest mass subtracted) μ_{e} − m_{e} (left), the proton chemical potential (rest mass subtracted) μ_{p} − m_{p} (centre), and the neutron chemical potential (rest mass subtracted) μ_{n} − m_{n} (right). We note that the colour bar scales are different for the various chemical potentials. As in Fig. 1, the black lines mark the isodensity contours. 
In the veryhighdensity regime, ρ ≳ 7 × 10^{14} g cm^{−3}, ν_{e} and ν_{μ} were suppressed (Y_{νe, μ} < 10^{−4}), while and had typical fractions , with being slightly more abundant than . We note that the excess of and exactly compensates for the amount of newly created electrons and muons (see the discussion in the previous section). Among the remaining neutrino species, ν_{τ} and were the most abundant, with a fraction ∼10^{−4} − 10^{−3}.
When ρ ∼ 2 − 6 × 10^{14} g cm^{−3}, the fractions of neutrinos increased up to ∼10^{−3}, but electron and muon antineutrinos still remained the most abundant neutrino species, with . All the (anti)neutrino abundances peaked in the hot spots. However, muon production dominated in this density and temperature regime (see the discussion in the previous section) so that the species was the most abundant, followed by and then ν_{τ}, ν_{e}, and ν_{μ}, in that order.
At the lower density, ρ ∼ 10^{13} − 10^{14} g cm^{−3}, all neutrinos and antineutrinos were present in warm matter streams characterised by k_{B}T_{sim} ∼ 25 MeV, with order of magnitude fractions ∼10^{−2}, and the species dominated, followed in order by , ν_{τ}, ν_{μ}, and ν_{e}. On the contrary, in cold matter streams with k_{B}T_{sim} ∼ 10 − 15 MeV, antineutrinos were suppressed, and the ν_{μ} species dominated, with a maximum fraction of ∼0.01, followed by ν_{e}, with Y_{νe} ≳ 10^{−3}, and ν_{τ}, with Y_{ντ} ∼ 10^{−4}.
The different abundances and the hierarchy reflect the behaviour of the neutrino chemical potentials at equilibrium (see Eq. (26)), which ultimately depend on the EOS both directly, through the nuclear interaction, and indirectly, through the different thermodynamic conditions realised in the remnant. The regions where antineutrinos dominated over neutrinos were characterised by a higher density and a higher value of the neutron chemical potential μ_{n}, which exceeded μ_{p} + μ_{l−} with l = {e, μ}, as we deduced from Figs. 3 and 5. In such regions, the density was increased by matter compression, and the fraction of protons increased as well since the lepton number was fixed. This coincided with the enhancement of electrons and muon production at the expense of free neutrons (see discussion in Sect. 4), and the equilibrium condition Eq. (26) enforced the creation of the corresponding antineutrinos. In contrast, neutrinos dominated in relatively cold matter streams in the outer layers, where the proton chemical potential increased (see Fig. 5) so that μ_{p} + μ_{l−} > μ_{n}. The enhancement of μ_{p} is indeed correlated with the decrease of both density and temperature experienced by the expanding matter in such regions. This is illustrated in Fig. 6, where we plotted μ_{p} and μ_{n} in the ρ − T plane for Y_{p} fixed to a relevant value for the analysis (i.e. Y_{p} = 0.08). We note that μ_{p} started increasing when ρ ≲ 10^{14} g cm^{−3} and k_{B}T ≲ 15 MeV, while μ_{n} decreased in the same ρ − T region. This is because when ρ ≲ 10^{14} g cm^{−3} and k_{B}T ≲ 15 MeV, nucleons start clustering in nuclei so that the chemical potential of free protons is enhanced and the amount of free protons drops (see e.g. Fig. 8 of Hempel & SchaffnerBielich 2010). At a fixed lepton fraction and at equilibrium, this results in a conversion of electrons and muons into the corresponding neutrinos.
Fig. 6. Proton and neutron chemical potentials (rest mass subtracted) of BLh in the densitytemperature plane at fixed Y_{p} = 0.08. The proton fraction Y_{p} = 0.08 is the equilibrium proton fraction at 4.6 ms after merger around the point x = 0.64 km, y = 12.76 km, z = 0 km, where ρ ∼ 5 × 10^{13} g cm^{−3}, T_{sim} ∼ 10 MeV, and neutrinos dominate, while antineutrinos are suppressed (see also Figs. 1, 3 and 5). We notice that the colour bars are scaled differently. 
The properties of the neutrino gases are better characterised in terms of their degeneracy parameters. In Fig. 7, we show the degeneracy parameter of electron and muon neutrinos computed for three snapshots of (BLh, 1.00): ∼5 ms, ∼7 ms, and ∼12 ms after merger. The spatial distribution of η_{νe} and η_{νμ} clearly indicates the presence of trapped degenerate gases of the and species in the remnant core and of the ν_{e} and ν_{μ} species in the outer layers. In this latter region, the trapped gas of electron neutrinos is characterised by a degeneracy parameter η_{νe} ranging between three and six, while in the remnant core degenerate electron antineutrinos have . The degenerate gases of muonic neutrinos (in the outer layers) and muonic antineutrinos (in the core) exhibited the same qualitative behaviour as the ν_{e} and species, with a degeneracy parameter η_{νμ} ranging between 4 and 9 so that these muon (anti)neutrino gases were slightly more degenerate than the electron ones.
Fig. 7. Degeneracy parameter of ν_{e} (first row) and ν_{μ} (lower row) on the equatorial plane obtained by postprocessing the simulation (BLh, 1.00). Three different time slices are shown corresponding from left to right to 4.6 ms, 6.7 ms, and 11.8 ms after merger. As in Fig. 1, the black lines mark the isodensity contours as in Fig. 1. 
The spatial distribution of the abundances as well as the neutrino hierarchy in (DD2, 1.00) and (SFHo, 1.00) are fully analogous to the ones of (BLh, 1.00). Accordingly, we can analyse the differences among the three EOSs by inspecting the degeneracy parameters. Similarly to (BLh, 1.00), and formed degenerate trapped gases in the remnant core, while ν_{e} and ν_{μ} formed degenerate trapped gases in the cold matter streams of the outer layers. The degeneracy parameter of electron neutrinos, η_{νe}, ranged between −2 and 6 for both (DD2, 1.00) and (SFHo, 1.00), while in the case of muonic neutrinos, we found η_{νμ} ranging between −2 and 11 for (DD2, 1.00) and between −2 and 12 for (SFHo, 1.00). In general, the gases of and are more degenerate for (BLh, 1.00) than for (DD2, 1.00) and (SFHo, 1.00). The main cause of this resides in the differences among the neutron and the proton chemical potentials, which depend on the selected EOS and on the values of ρ, T, and Y_{p} in the remnant core. In particular, the average values of (rest masses subtracted) at a radius of approximately 3 km from the centre of the remnants correspond to ∼200 MeV for DD2, ∼245 MeV for SFHo, and ∼325 MeV for BLh. Hence, the degeneracy of the trapped antineutrinos is comparable for DD2 and SFHo but quite different for BLh (see Eq. (26)). In this sense, antineutrino trapping in the core of the remnant shows a significant and nontrivial dependence on the EOS through the nucleons’ chemical potentials. If the merger dynamics are significantly affected, this characteristic could become a probe for investigating matter properties in the highdensity regime.
Finally, by postprocessing the results of the simulation (BLh, 1.43), we observed that the degeneracy of the trapped component of antineutrinos was the same as in (BLh, 1.00). However, the trapped neutrino gases at ρ ∼ 10^{13} − 10^{14} g cm^{−3} were less degenerate for q = 1.43, displaying η_{νe} ∼ 4 and η_{νμ} ∼ 7. This feature depends mostly on the temperature in the density region 10^{13} − 10^{14} g cm^{−3}, which is systematically larger in the case of q = 1.43 compared to q = 1.00.
To test the validity of the approximations used in Sect. 2 and to further characterise the properties of neutrinos trapped inside the remnant in the presence of muons for the first time, we computed the neutrino mean energies, E_{ν}, in postprocessing. In Fig. 8, we show the spatial distribution of E_{ν} for the remnant of simulation (BLh, 1.00) at approximately 5 ms after merger. In the density region ρ > 10^{13} g cm^{−3}, we found E_{ν} > 30 MeV for all flavour neutrinos, unless a certain flavour was suppressed in a portion of space. Electron and muon antineutrinos were the most energetic, with MeV. Then, for ν_{τ} we found E_{ντ} ≲ 170 MeV, while ν_{e} and ν_{μ} were the least energetic, with E_{νe, μ} ≲ 160 MeV. In the case of (DD2, 1.00), the values of E_{ν} closely followed the ones discussed for (BLh, 1.00). However, (SFHo, 1.00) exhibited larger energies for all flavour neutrinos, with E_{νe, μ, τ} ≲ 180 MeV and MeV, due to larger remnant temperatures.
Fig. 8. Mean energy E_{ν} of trapped neutrinos on the equatorial plane computed by postprocessing simulation (BLh, 1.00) at 4.6 ms after merger. As in Fig. 1., the black lines mark isodensity contours. 
6. Changes in the remnant pressure
In this section, we investigate how the result of the postprocessing, which consistently included muons and trapped neutrinos, modifies the pressure of the remnant with respect to the values extracted from the original simulation. In Fig. 9, we show the ratio between the pressure computed in postprocessing and the simulation pressure, P/P_{sim}, for (BLh, 1.00), (DD2, 1.00), and (SFHo, 1.00). Similar to Sect. 4 and Sect. 5, we start our discussion with the (BLh, 1.00) simulation. In the time domain of our analysis, the equilibrium pressure computed in postprocessing decreased by 6 − 7% in the region characterised by ρ ∼ 5 × 10^{14} g cm^{−3}, while it increased by a maximum of 3% around 10 km from the centre at ρ ∼ 10^{14} g cm^{−3} compared to the pressure from the original simulation.
Fig. 9. Pressure ratio P/P_{sim} for (BLh, 1.00) (left), (DD2, 1.00) (centre), and (SFHo, 1.00) (right) where P is obtained in postprocessing, while P_{sim} is the simulation pressure. For (BLh, 1.00), we explicitly mark the maximum P/P_{sim} ∼ 1.03 (white circle) and the minimum P/P_{sim} ∼ 0.94 (white cross) at ρ > 10^{14} g cm^{−3}. As in Fig. 1., the black lines mark isodensity contours. 
To understand which processes are responsible for the changes in pressure, we analysed the role of each particle species in detail. In Fig. 10, we have plotted the difference in pressures and particle fractions after muons and neutrinos were introduced
Fig. 10. Change of particle pressure and fractions obtained in correspondence to the spacetime coordinates explicitly marked in the left panel of Fig. 9 and corresponding to P/P_{sim} = 0.94 (upper row) and P/P_{sim} = 1.03 (lower row). The bar plots in blue (left) give the difference between the equilibrium pressure computed in postprocessing and in the simulation, , where , and P_{tot} refers to the total pressure. The bar plots in red (centre) show the difference between the particle fractions computed in postprocessing and the ones from the simulation, , where , and e (μ) refers to the net electron (muon) fractions. The colourcoded plots (right) show P_{b} of BLh as a function of temperature and proton fraction for the same points. The density ρ is fixed to the value obtained in the simulation, that is, ρ = 4.3 × 10^{14} g cm^{−3} (upper row) and ρ = 1.55 ⋅ 10^{14} g cm^{−3} (lower row). The black star and the diamond mark are the values from the simulation (before postprocessing) and at equilibrium (after postprocessing), respectively. 
for all the particle species at two representative points inside the simulations, explicitly marked in Fig. 9 and corresponding to P/P_{sim} ∼ 0.94 and P/P_{sim} ∼ 1.03, respectively. The drop in pressure was localised inside a hot spot at k_{B}T_{sim} ∼ 55 MeV, and it was mainly due to a decrease of baryonic pressure. Once muons were included in the microphysics, the system reached a new equilibrium characterised by a smaller neutron fraction; larger proton, electron, and muon fractions; and a smaller temperature (see also discussion in Sect. 4). The decrease in both the neutron fraction and of the temperature lowered the baryonic pressure, and the pressure provided by muons and trapped antineutrinos was not large enough to compensate. Accordingly, the total pressure also decreased. This effect is strongly correlated with the enhancement in the production of massive leptons and (anti)neutrinos at high temperatures. An increase of pressure occurred at a lower initial temperature, k_{B}T_{sim} ∼ 15 MeV, and it was correlated with the inclusion of muons from the cold NSs. As shown in Fig. 2, including muons in cold neutrinoless βequilibrium implies lowering the initial electron and neutron fraction in favour of the muon and proton fractions. Therefore, the pressure of the system increased because some relativistic particles (electrons) were replaced with nonrelativistic particles (muons). An additional contribution came from the increase of the baryon pressure , where P_{b} is computed in postprocessing while is read from the simulation outcome. In general, we would expect ΔP_{b} < 0 since the neutron fraction was reduced. However, we observed a small increase of temperature ΔT = 1.3 MeV, which was large enough to result in ΔP_{b} > 0, as shown in the colourcoded plot in the second row of Fig. 10. We note that the increase in temperature was expected because it is the only way the system could simultaneously reduce the number of neutrons and electrons, which are highly degenerate at k_{B}T_{sim} ∼ 15 MeV and at a fixed internal energy. Interestingly, adding new d.o.f. without supplementary repulsive forces did not simply result in a softening of the EOS, as we would have expected, because of the role played by the equilibrium temperature.
The causes of the changes in pressure in (DD2, 1.00) and (SFHo, 1.00) are the same as in (BLh, 1.00). From our analysis up to 15 ms after merger, we found that P/P_{sim} ranges from 0.94 to 1.05 for (DD2, 1.00) and from 0.93 to 1.04 for (SFHo, 1.00). As shown in Fig. 9, (DD2, 1.00) was subjected mostly to a pressure increase, while (SFHo, 1.00) exhibited a pressure increase in the region closer to the core and a pressure decrease in the outer shells. The amount of pressure decrease is comparable to the one found in Perego et al. (2019), while the pressure increase is a new prediction.
Pressure modifications correlate in a complex way with other thermodynamical variables, such as density, electron and muon fractions, temperature, and chemical potentials. However, a more important role is played by temperature. A close inspection of the distribution of P/P_{sim} in the space revealed that for all the EOSs, the pressure decrease was more frequent at larger T_{sim}. This is because electrons and muons are less degenerate, so the of conversion neutrons into protons is favoured and the baryonic pressure reduces to a greater extent. Regarding the pressure increase, it was frequent at smaller k_{B}T_{sim} < 20 MeV because of the dominant role played by the muons coming from the NSs in colder temperature regimes. We concluded that the remnant formed in (DD2, 1.00) exhibits mostly a pressure increase because its temperature is, on average, smaller compared to (BLh, 1.00) and (SFHo, 1.00).
We conclude this section by analysing how the pressure is modified in the case of an asymmetric binary mass ratio. If q = 1.43, then the ratio P/P_{sim} ranges in the interval 0.93 − 1.05, showing a pressure decrease compatible with the case where q = 1.00 but has a pressure increase that is significantly larger. This behaviour depends mostly on the difference between the electron fractions in the two simulations. In the case of q = 1.43, the regions with the highest P/P_{sim} also have a larger compared to q = 1.00 and as a consequence a larger (see Fig. 2). Since the pressure increase is mainly driven by the muon fraction coming from the cold NSs, a larger results in a larger P/P_{sim}. In Fig. 11, we show the comparison between the pressure ratios for q = 1.00 and q = 1.43 at approximately 7 ms after merger. We found that the spatial distribution of P/P_{sim} for q = 1.43 is quite asymmetric compared to q = 1.00 in the fulltime domain of our analysis. In this case, the pressure decrease is localised in the core of the secondary NS, which is strongly heated by the shocks developing at the contact surface. The pressure increases in the cold external layers of the secondary NS, which are stripped away during the dynamical evolution without being significantly heated.
Fig. 11. Pressure ratio P/P_{sim} computed by postprocessing (BLh, 1.00) (top) and (BLh, 1.43) (bottom) at 6.7 ms after merger. As in Fig. 1, the black lines mark the isodensity contours. 
7. Discussion
Based on our results, in this section we speculate about the possible consequences for the evolution of the remnant and the emitted neutrinos. The fraction of muons is enhanced at high temperatures. Therefore, muons are expected to play a more important role in BNS mergers characterised by larger masses and in the case of soft EOSs producing more compact NSs. Such high temperatures develop mostly at the contact interface between the two fusing cores during the first milliseconds after a merger. In addition, CC reactions involving μ^{±}, ν_{μ}, and provide a new source of opacity for muonic (anti)neutrinos, as these reactions influence their diffusion. The neutrino surfaces of and ν_{μ} are expected to split in a manner similar to ν_{e} and (Endrizzi et al. 2020), with the surface of located at a larger rest mass density compared to ν_{μ}. Hence, the luminosity and the spectrum of the emitted will deviate from that of ν_{μ}, with possible implications for the emitted spectrum and for the neutrino oscillations (see Fischer et al. 2020 for luminosity and spectra of ν_{μ} and in CCSN). Furthermore, the formation of trapped degenerate gases of ν_{e} and ν_{μ} in the outer layers can favour the development of neutrinobursts during the dynamical evolution, similar to the case of CCSN (Fischer et al. 2020), when dense, hot matter expands in the forming disc and becomes transparent. The change in pressure induced by muons and trapped neutrinos can change the time of the collapse of the remnant and affect the dynamical mass ejection. In particular, in (DD2, 1.00) and (SFHo, 1.00), we observed a diffuse pressure increase near the core (see discussion in Sect. 6), which can possibly postpone the collapse of the remnant. For (BLh, 1.00), the drop in pressure at high density induces a softening that could speed up the collapse of more massive binaries, while the larger pressure at smaller densities could affect the lowdensity dynamics, including the ejecta production, the disc formation, and the stability of the remnant produced in low mass BNS mergers. In Sect. 6, we established a correlation between P/P_{sim} and T_{sim}, pointing out that P/P_{sim} > 1 (P/P_{sim} < 1) when k_{B}T_{sim} is below (above) ∼20 MeV. This feature could induce positive feedback that amplifies pressure variations once muons and trapped neutrinos are included at the dynamical stage of the simulation. When the pressure decreases, the temperature of the more compressed system will increase, possibly triggering processes that lower the pressure even more. Nevertheless, when the pressure increases, the temperature of the expanded system will decrease, possibly favouring further pressure enhancement. However, feedback effects can be assessed only by including muons and trapped neutrinos in the dynamical evolution of the system. Moreover, the effect of pressure variations is, in general, nontrivial to predict. Indeed, pressure is one of the source terms in the Einstein equations and its significant increase could result in a faster collapse.
It is interesting to compare our results to the findings from CCSN simulations that include muons (Bollig et al. 2017; Fischer et al. 2020). Our estimate of Y_{μ} exceeds the muon fraction in the proximity of the core bounce by one order of magnitude (Fischer et al. 2020), while it is compatible with Y_{μ} in the late postbounce phase (Bollig et al. 2017). This is because in CCSNe, matter density and temperature increase after the core bounce, while Y_{e} decreases, thus approaching the regime of our analysis. Since temperatures are significantly larger in BNS mergers than in CCSNe, thermal processes could contribute more significantly to the creation of μ^{±}, while they are negligible in CCSNe (Bollig et al. 2017; Guo et al. 2020; Fischer et al. 2020). However, the enhancement of the net muon fraction, Y_{μ}, with the temperature is in agreement with the results of Bollig et al. (2017), Fischer et al. (2020).
We note that while the inclusion of muons always implies a softening of the EOS in the postbounce phase in CCSNe (Bollig et al. 2017; Fischer et al. 2020), in our case, the pressure can either increase or decrease with respect to the case where muons are neglected from the beginning. This difference is expected since in BNS mergers, the pressure increase obtained in our analysis is driven by the muons already present in the two cold NSs (see Sect. 6). Notably, muons are not present in stellar iron cores before a collapse.
Next, we relate our results to stateoftheart BNS merger simulations that include trapped neutrinos either in postprocessing (Perego et al. 2019) or through an M1 scheme (e.g. Foucart et al. 2016a; Radice et al. 2022; Zappa et al. 2023). The neutrino hierarchy resulting from our analysis is in disagreement with these works. In these previous papers, matter is dominated by and followed by heavylepton neutrinos and then ν_{e}, whereas we found that are the most abundant, followed in order by , , ν_{e}, and ν_{μ}. According to Perego et al. (2019), the fluid pressure decreases after the inclusion of trapped neutrinos, but we show that it can also increase in some regions of the remnant. This is an indirect confirmation that the enhancement of pressure is mainly due to muons. Nonetheless, our estimate of the pressure decrease is compatible with the one computed in Perego et al. (2019).
8. Conclusions
In this paper, we evaluated the contribution of muons in the postmerger remnant formed during BNS mergers as well as their impact on the trapped neutrino component during the first milliseconds after the merger. In particular, we considered a sample of numerical relativity simulations targeted at GW170817 that were performed using three different nuclear EOSs (BLh, DD2, and SFHo). Based on the remnant conditions, we estimated the abundances and the properties of muons and trapped neutrinos in postprocessing, assuming weak and thermal equilibrium for rest mass densities in excess of 10^{13} g cm^{−3}. We found that a nonnegligible amount of muons is present in the remnant, and they significantly alter the neutrino hierarchy, favouring the production of muon antineutrinos. Moreover, the new equilibrium conditions modify the pressure inside the remnant in a nontrivial way with respect to the case in which muons and trapped neutrinos are neglected.
More specifically, we found that a significant amount of muons is present at baryon densities ρ > 10^{13} g cm^{−3} with a fraction Y_{μ−} ∼ 10^{−3} − 0.07 for the BLh EOS and Y_{μ−} ∼ 10^{−3} − 0.05 for the DD2 and SFHo EOSs. The net fraction of muons is between 30% and 50–70% of the net electron fraction, with a maximum depending on the nuclear EOS. The bulk of muons comes from the two cold NSs, while a fraction corresponding to 10^{−3} − 10^{−2} is created during the merger via conversion of neutrons into protons plus muons and thermal processes. Accordingly, muon production is enhanced in the highdensity fusing cores and in the hot spots characterising the remnant structure during the first milliseconds after the merger.
The presence of muons modifies the flavour hierarchy of the trapped neutrino component. In the core, matter is dominated by antineutrinos, with being the most abundant species. Immediately outside the core and in correspondence with the hot spots, gases of trapped neutrinos and antineutrinos coexist, but the species still dominates, followed in order by , , ν_{e}, and ν_{μ}. In contrast, in the cold matter streams at ρ ∼ 10^{13} − 10^{14} g cm^{−3}, antineutrinos are suppressed, while neutrinos dominate and the ν_{μ} species is the most abundant. Antineutrinos and neutrinos are degenerate in the core and in the cold matter streams, respectively. The level of degeneracy depends on the thermodynamical conditions of the remnant and on μ_{n} − μ_{p}. In particular, typical values of the degeneracy parameters of and are , with larger (smaller) values associated with larger (smaller) μ_{n} − μ_{p}. Therefore, the properties of trapped (anti)neutrinos in BNS mergers have a nontrivial dependence on the properties of the EOS at densities well above n_{0}.
The presence of muons and trapped neutrinos affect the remnant pressure, whose value can decrease by 7% or increase by up to 5%, depending on the nuclear EOS, compared with the case in which they are neglected. The pressure drops mainly in the hot spots, where the processes driving muon production reduce the temperature and number of neutrons. The greater pressure mostly observed at low temperatures and high density is due to the muons coming from the cold NSs. This new result shows how adding a d.o.f. in the system microphysics does not always imply a decrease in pressure, as we could naively expect. Finally, we demonstrated that in asymmetric binaries, the pressure decrease is localised in the shocked core of the secondary NS, while the pressure increase is found in the cold outer layers of the secondary NS, which are stripped away by the primary NS.
Since the procedure to postprocess the simulation data is not unique, in Appendix A we discuss the robustness of our results with respect to the arbitrary choices in our procedure. For example, we show that when neglecting muons inside the cold NSs before merger, the major results are the same as when the presence of muons in the inspiraling NSs is considered. The comparison with the first postprocessing procedure emphasises that a significant muon contribution originates from the muons already present in the cold NSs.
Our work is nevertheless limited by the postprocessing approach. To understand how muons really affect the merger and postmerger dynamics, it is necessary to implement a transport scheme that includes muons as an independent d.o.f.. Also, adding the rates of the reactions involved in muon production (see Sect. 2) would be important for estimating the effect of muons on the neutrino opacity and luminosity, as done in CCSN simulations (Guo et al. 2020; Fischer et al. 2020). Additionally, the role of pions coupling to muons and muonic (anti)neutrinos should be considered while computing the mean free path of low energy ν_{μ} and (Fore & Reddy 2020). We also highlight that in the description of NS matter, we have neglected the possible formation of hyperons (Oertel et al. 2012; Fortin et al. 2018) and/or that a deconfinement phase transition to quark matter may take place (Weissenborn et al. 2011; Klähn et al. 2013; Chatterjee & Vidaña 2016; Bombaci et al. 2016; Logoteta et al. 2019; Logoteta 2021). The inclusion of these additional d.o.f. will be the focus of future investigations.
This work shows, for the first time, that muons have a discernible impact on the trapped neutrino component and possibly on the thermodynamics and dynamics of BNS mergers. Our results emphasise the need of including muons in future simulation modelling of the postmerger remnant of BNS mergers.
The corresponding EOS tables are publicly available from the CompOSE repository https://compose.obspm.fr/.
We note, however, that the dimensionless tidal deformability predicted by DD2 for GW170817like events is > 800, which is in possible tension with constraints derived from GW170817 (e.g. Abbott et al. 2017a, 2019; Breschi et al. 2021).
We note that the electromagnetic interaction of negatively charged leptons with positively charged ions is however taken into account by applying a correction to the lepton chemical potential, as prescribed in Hempel & SchaffnerBielich (2010).
At high density, neutrino oscillations are suppressed (see e.g. Richers et al. 2019) so that electron and muon lepton fractions are separately conserved.
Acknowledgments
The authors thank David Radice, Sebastiano Bernuzzi, Micaela Oertel, and Francesco Pederiva for useful discussions. They also thank the anonymous Referee for their valuable comments and suggestions. They acknowledge the INFN for the usage of computing and storage resources through the tullio cluster in Turin, and the Computational Relativity (CoRe) collaboration for providing access to the simulation data used in the work. MB and EL acknowledge financial support from MIUR (PRIN 2020 grant 2020KB33TP). MB acknowledges financial support from MIUR (PRIN 2017 grant 20179ZF5KS).
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, Phys. Rev. Lett., 119, 161101 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, ApJ, 848, L12 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018, Phys. Rev. Lett., 121, 161101 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. X, 9, 011001 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020, Liv. Rev. Relativ., 23, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Alford, M., Harutyunyan, A., & Sedrakian, A. 2021, Phys. Rev. D, 104, 103027 [NASA ADS] [CrossRef] [Google Scholar]
 Alford, M. G., & Good, G. 2010, Phys. Rev. C, 82, 055805 [NASA ADS] [CrossRef] [Google Scholar]
 Antoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340, 448 [Google Scholar]
 Bauswein, A., Janka, H. T., Hebeler, K., & Schwenk, A. 2012, Phys. Rev. D, 86, 063001 [Google Scholar]
 Bauswein, A., Just, O., Janka, H.T., & Stergioulas, N. 2017, ApJ, 850, L34 [Google Scholar]
 Bernuzzi, S. 2020, General Relativ. Gravitation, 52, 108 [CrossRef] [Google Scholar]
 Bernuzzi, S., Radice, D., Ott, C. D., et al. 2016, Phys. Rev. D, 94, 024023 [Google Scholar]
 Bernuzzi, S., Breschi, M., Daszuta, B., et al. 2020, MNRAS, 497, 1488 [Google Scholar]
 Bludman, S. A., & van Riper, K. A. 1977, ApJ, 212, 859 [NASA ADS] [CrossRef] [Google Scholar]
 Bollig, R., DeRocco, W., Graham, P. W., & Janka, H.T. 2020, Phys. Rev. Lett., 125, 051104 [NASA ADS] [CrossRef] [Google Scholar]
 Bollig, R., Janka, H. T., Lohs, A., et al. 2017, Phys. Rev. Lett., 119, 242702 [NASA ADS] [CrossRef] [Google Scholar]
 Bombaci, I., & Logoteta, D. 2018, A&A, 609, A128 [EDP Sciences] [Google Scholar]
 Bombaci, I., Logoteta, D., Vidaña, I., & Providência, C. 2016, Eur. Phys. J. A, 52, 58 [NASA ADS] [CrossRef] [Google Scholar]
 Breschi, M., Gamba, R., & Bernuzzi, S. 2021, Phys. Rev. D, 104, 042001 [NASA ADS] [CrossRef] [Google Scholar]
 Breschi, M., Bernuzzi, S., Godzieba, D., Perego, A., & Radice, D. 2022, Phys. Rev. Lett., 128, 161102 [NASA ADS] [CrossRef] [Google Scholar]
 Cameron, A. G. W. 1970, ARA&A, 8, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Chatterjee, D., & Vidaña, I. 2016, Eur. Phys. J. A, 52, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Cohen, J. M., Langer, W. D., Rosen, L. C., & Cameron, A. G. W. 1970, Ap&SS, 6, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Colombo, A., Salafia, O. S., Gabrielli, F., et al. 2022, ApJ, 937, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Cooperstein, J. 1988, Phys. Rept., 163, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Cusinato, M., Guercilena, F. M., Perego, A., et al. 2022, EPJA, 58, 99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Damour, T., Nagar, A., & Villain, L. 2012, Phys. Rev. D, 85, 123007 [NASA ADS] [CrossRef] [Google Scholar]
 Danielewicz, P., & Lee, J. 2014, Nucl. Phys. A, 922, 1 [NASA ADS] [CrossRef] [Google Scholar]
 De, S., Finstad, D., Lattimer, J. M., et al. 2018, Phys. Rev. Lett., 121, 091102; Erratum: 121, 259902 [Google Scholar]
 Drischler, C., Hebeler, K., & Schwenk, A. 2019, Phys. Rev. Lett., 122, 042501 [Google Scholar]
 Endrizzi, A., Perego, A., Fabbri, F. M., et al. 2020, Eur. Phys. J. A, 56, 15 [EDP Sciences] [Google Scholar]
 Evans, M., Adhikari, R. X., Afle, C., et al. 2021, ApJ, submitted [arXiv:2109.09882] [Google Scholar]
 Favata, M. 2014, Phys. Rev. Lett., 112, 101101 [NASA ADS] [CrossRef] [Google Scholar]
 Fischer, T., Guo, G., MartínezPinedo, G., Liebendörfer, M., & Mezzacappa, A. 2020, Phys. Rev. D, 102, 123001 [NASA ADS] [CrossRef] [Google Scholar]
 Fore, B., & Reddy, S. 2020, Phys. Rev. C, 101, 035809 [NASA ADS] [CrossRef] [Google Scholar]
 Fortin, M., Oertel, M., & Providência, C. 2018, PASA, 35, e044 [NASA ADS] [CrossRef] [Google Scholar]
 Foucart, F., Haas, R., Duez, M. D., et al. 2016a, Phys. Rev. D, 93, 044019 [NASA ADS] [CrossRef] [Google Scholar]
 Foucart, F., O’Connor, E., Roberts, L., et al. 2016b, Phys. Rev. D, 94, 123016 [NASA ADS] [CrossRef] [Google Scholar]
 Glendenning, N. K. 1997, Compact stars. Nuclear Physics, Particle Physics, and General Relativity [Google Scholar]
 Grimm, S., & Harms, J. 2020, Phys. Rev. D, 102, 022007 [NASA ADS] [CrossRef] [Google Scholar]
 Guo, G., MartínezPinedo, G., Lohs, A., & Fischer, T. 2020, Phys. Rev. D, 102, 023037 [NASA ADS] [CrossRef] [Google Scholar]
 Haensel, P., Levenfish, K. P., & Yakovlev, D. G. 2000, A&A, 357, 1157 [NASA ADS] [Google Scholar]
 Haensel, P., Levenfish, K. P., & Yakovlev, D. G. 2001, A&A, 372, 130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haft, M., Raffelt, G., & Weiss, A. 1994, ApJ, 425, 222; Erratum: 1995, 438, 1017 [NASA ADS] [CrossRef] [Google Scholar]
 Hammond, P., Hawke, I., & Andersson, N. 2022, ApJ, submitted [arXiv:2205.11377] [Google Scholar]
 Hannestad, S., & Raffelt, G. 1998, ApJ, 507, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Hempel, M., & SchaffnerBielich, J. 2010, Nucl. Phys. A, 837, 210 [Google Scholar]
 Hinderer, T. 2008, ApJ, 677, 1216 [NASA ADS] [CrossRef] [Google Scholar]
 Hix, W. R., & Thielemann, F.K. 1999, ApJ, 511, 862 [NASA ADS] [CrossRef] [Google Scholar]
 Iacovelli, F., Mancarella, M., Foffa, S., & Maggiore, M. 2022, ApJ, 941, 208 [NASA ADS] [CrossRef] [Google Scholar]
 Kaplan, J. D., Ott, C. D., O’Connor, E. P., et al. 2014, ApJ, 790, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Klähn, T., Łastowiecki, R., & Blaschke, D. B. 2013, Phys. Rev. D, 88, 085001 [CrossRef] [Google Scholar]
 Lattimer, J. M., & Steiner, A. W. 2014, ApJ, 784, 123 [Google Scholar]
 Leinson, L. B. 2002, Nucl. Phys. A, 707, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Logoteta, D. 2021, Universe, 7, 408 [NASA ADS] [CrossRef] [Google Scholar]
 Logoteta, D., Vidana, I., & Bombaci, I. 2019, Eur. Phys. J. A, 55, 207 [CrossRef] [Google Scholar]
 Logoteta, D., Perego, A., & Bombaci, I. 2021, A&A, 646, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Maggiore, M., Van Den Broeck, C., Bartolo, N., et al. 2020, J. Cosmol. Astropart. Phys., 2020, 050 [CrossRef] [Google Scholar]
 Margalit, B., & Metzger, B. D. 2017, ApJ, 850, L19 [Google Scholar]
 Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, ApJ, 887, L24 [Google Scholar]
 Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2021, ApJ, 918, L28 [NASA ADS] [CrossRef] [Google Scholar]
 Most, E. R., Papenfort, L. J., Dexheimer, V., et al. 2019, Phys. Rev. Lett., 122, 061101 [NASA ADS] [CrossRef] [Google Scholar]
 Nedora, V., Bernuzzi, S., Radice, D., et al. 2019, ApJ, 886, L30 [NASA ADS] [CrossRef] [Google Scholar]
 Nedora, V., Bernuzzi, S., Radice, D., et al. 2021, ApJ, 906, 98 [NASA ADS] [CrossRef] [Google Scholar]
 Oertel, M., Fantina, A. F., & Novak, J. 2012, Phys. Rev. C, 85, 055806 [NASA ADS] [CrossRef] [Google Scholar]
 Oertel, M., Hempel, M., Klähn, T., & Typel, S. 2017, Rev. Modern Phys., 89, 015007 [NASA ADS] [CrossRef] [Google Scholar]
 Patricelli, B., Bernardini, M. G., Mapelli, M., et al. 2022, MNRAS, 513, 4159 [NASA ADS] [CrossRef] [Google Scholar]
 Perego, A., Bernuzzi, S., & Radice, D. 2019, Eur. Phys. J. A, 55, 124 [EDP Sciences] [Google Scholar]
 Perego, A., Logoteta, D., Radice, D., et al. 2022, Phys. Rev. Lett., 129, 032701 [NASA ADS] [CrossRef] [Google Scholar]
 Potekhin, A. Y., Pons, J. A., & Page, D. 2015, Space Sci. Rev., 191, 239 [Google Scholar]
 Prakash, M., Bombaci, I., Prakash, M., et al. 1997, Phys. Rep., 280, 1 [Google Scholar]
 Punturo, M., Abernathy, M., Acernese, F., et al. 2010, CQG, 27, 194002 [NASA ADS] [CrossRef] [Google Scholar]
 Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, ApJ, 918, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Radice, D., Galeazzi, F., Lippuner, J., et al. 2016, MNRAS, 460, 3255 [NASA ADS] [CrossRef] [Google Scholar]
 Radice, D., Perego, A., Hotokezaka, K., et al. 2018a, ApJ, 869, 130 [Google Scholar]
 Radice, D., Perego, A., Zappa, F., & Bernuzzi, S. 2018b, ApJ, 852, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Radice, D., Bernuzzi, S., & Perego, A. 2020, Ann. Rev. Nucl. Part. Sci., 70, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Radice, D., Bernuzzi, S., Perego, A., & Haas, R. 2022, MNRAS, 512, 1499 [NASA ADS] [CrossRef] [Google Scholar]
 Richers, S. A., McLaughlin, G. C., Kneller, J. P., & Vlasenko, A. 2019, Phys. Rev. D, 99, 123014 [NASA ADS] [CrossRef] [Google Scholar]
 Riley, T. E., et al. 2019, ApJ, 887, L21 [CrossRef] [Google Scholar]
 Ronchini, S., Branchesi, M., Oganesyan, G., et al. 2022, A&A, 665, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rosswog, S., & Liebendoerfer, M. 2003, MNRAS, 342, 673 [NASA ADS] [CrossRef] [Google Scholar]
 Ruffert, M., & Janka, H. T. 1998, A&A, 338, 535 [NASA ADS] [Google Scholar]
 Ruffert, M., Janka, H. T., & Schaefer, G. 1996, A&A, 311, 532 [NASA ADS] [Google Scholar]
 Sekiguchi, Y., Kiuchi, K., Kyutoku, K., & Shibata, M. 2011a, Phys. Rev. Lett., 107, 211101 [NASA ADS] [CrossRef] [Google Scholar]
 Sekiguchi, Y., Kiuchi, K., Kyutoku, K., & Shibata, M. 2011b, Phys. Rev. Lett., 107, 051102 [NASA ADS] [CrossRef] [Google Scholar]
 Sekiguchi, Y., Kiuchi, K., Kyutoku, K., Shibata, M., & Taniguchi, K. 2016, Phys. Rev. D, 93, 124046 [NASA ADS] [CrossRef] [Google Scholar]
 Shapiro, S. L., & Teukolsky, S. A. 1983, Black Holes, White Dwarfs, and Neutron Stars : the Physics of Compact Objects [Google Scholar]
 Shlomo, S., Kolomietz, V. M., & Colò, G. 2006, Eur. Phys. J. A, 30, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Steiner, A. W., Prakash, M., Lattimer, J. M., & Ellis, P. J. 2005, Phys. Rep., 411, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Steiner, A. W., Hempel, M., & Fischer, T. 2013, ApJ, 774, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Takami, K., Rezzolla, L., & Baiotti, L. 2014, Phys. Rev. Lett., 113, 091104 [NASA ADS] [CrossRef] [Google Scholar]
 Timmes, F. X., & Arnett, D. 1999, ApJS, 125, 277 [Google Scholar]
 Typel, S., Röpke, G., Klähn, T., Blaschke, D., & Wolter, H. H. 2010, Phys. Rev. C, 81, 015803 [NASA ADS] [CrossRef] [Google Scholar]
 Weih, L. R., Hanauske, M., & Rezzolla, L. 2020, Phys. Rev. Lett., 124, 171103 [NASA ADS] [CrossRef] [Google Scholar]
 Weissenborn, S., Sagert, I., Pagliara, G., Hempel, M., & SchaffnerBielich, J. 2011, ApJ, 740, L14 [NASA ADS] [CrossRef] [Google Scholar]
 Yakovlev, D. G., & Pethick, C. J. 2004, ARA&A, 42, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Zappa, F., Bernuzzi, S., Radice, D., Perego, A., & Dietrich, T. 2018, Phys. Rev. Lett., 120, 111101 [NASA ADS] [CrossRef] [Google Scholar]
 Zappa, F., Bernuzzi, S., Radice, D., & Perego, A. 2023, MNRAS, 520, 1481 [Google Scholar]
Appendix A: An alternative assumption on Y_{l, μ}
As discussed in Sect. 3.3, there is not a unique way to assign Y_{l, μ}, Y_{l, e}, and e_{tot} appearing on the lefthand side of Eq. (34) given the outcome of the simulations presented in Sect. 3.1. Therefore, in this Appendix, we consider a limiting case opposite to the one discussed in Sect. 3.3 as a different possible choice. In particular, we assume that muons are negligible in the initial conditions, that is, in the two cold NSs, and we impose Y_{l, μ} = 0. Then, we identify the simulation electron fraction, , and energy density, e_{sim}, with Y_{l, e} and e_{tot}, respectively:
(see from Eq. (30) to Eq. (33)). This choice, despite being inconsistent with what we expect, has the advantage of being fully coherent with the simulation assumptions, and it provides an easier initialisation of our postprocessing calculation. In the following paragraphs, we discuss the results obtained with this choice, and we compare them with the ones discussed in the previous sections. From this point on, we refer to this particular choice as Y_{l, μ} = 0 and to the one discussed in Sect. 3.3 as .
A.1. The fraction of muons
If Y_{l, μ} = 0, all the muons in the final state are produced during the merger. For all the simulations in Table 1, the equilibrium muon fraction at density ρ > 10^{14} g cm^{−3} is Y_{μ−} ∼ 10^{−3} − 3 × 10^{−2}. At a lower density, ρ ∼ 10^{13} − 10^{14} g cm^{−3}, the muon fraction is Y_{μ−} ∼ 10^{−3} if T_{sim} ≳ 10 MeV; otherwise, muons are not present. We note that the amount of muons produced only during the merger is slightly smaller when (see Sect. 4) because of the Pauli blocking coming from the muons already present in the cold NSs.
For all the simulations, Y_{μ}/Y_{e} ∼ 0.3 − 0.4, but the maximum Y_{μ}/Y_{e} is larger for DD2 and SFHo than for BLh because the latter exhibits a larger Y_{e} (see also Fig. 2). The creation of muons is enhanced in correspondence to the hot spots for all the EOSs. In the case of (BLh, 1.00), muon production also increases in the density region ρ > 10^{15} g cm^{−3}. This effect depends on the significantly larger electron fraction obtained with the BLh EOS since the more degenerate electrons block the conversion of neutrons into protons plus electrons, favouring instead the conversion into protons plus muons.
A.2. The trapping of neutrinos
The initial assumption on Y_{l, μ} does not significantly affect the behaviour of electron (anti)neutrinos, but it does imply discernible consequences on the muon neutrinos. If Y_{l, μ} = 0, muon neutrinos appear only in correspondence to the hot spots and the warm matter streams, with a maximum fraction Y_{νμ} ∼ 0.01. All ν_{μ} are produced together with since there is not an initial muon fraction that can be converted into ν_{μ}. In contrast, the amount of at equilibrium is strongly enhanced to guarantee the conservation of Y_{l, μ}. Therefore, for all the simulations considered, and constitute the most abundant species amongst all flavour (anti)neutrinos in the regions that are not dominated by thermal production.
The differences in the equilibrium fractions of neutrinos depending on the chosen Y_{l, μ} reflect the differences in the degeneracy parameters. For all the simulations η_{νμ} ≤ 0 and the production of ν_{μ}, only thermal and trapped degenerate ν_{μ} are suppressed. In addition, the gases of and trapped in the remnant core are more degenerate when Y_{l, μ} = 0 compared to when . This follows in part from the enhanced production of when Y_{l, μ} = 0, but it also depends on the equilibrium temperature T, which is systematically larger when . Actually, in the latter case, the system lowers the fraction of electrons and neutrons to take into account the fraction of muons from the cold NSs, which practically means reducing the system degeneracy by increasing the temperature (see the discussions in Sect. 3.3 and Sect. 6).
We stress that the presence of a trapped degenerate component of and in the core and of ν_{e} in the outer layers of the remnant is a robust prediction independent of the initial assumption on Y_{l, μ}. Instead, the presence of a degenerate gas of ν_{μ} in the cold matter streams depends on the initial choice of Y_{l, μ}.
A.3. The changes in pressure
In conclusion, we analysed the ratio P/P_{sim} where P is computed in postprocessing under the assumption Y_{l, μ} = 0. In this case, we did not observe a significant pressure increase after muons were introduced (i.e. P/P_{sim} > 1) for any of the simulations in Table 1. This result is expected since, within this assumption, we neglected the muons from the cold NSs, which are indeed responsible for the pressure increase when , as discussed in Sect. 6. In fact, the amount and the spatial distribution of the pressure decrease are common to all the simulations irrespective of the initial assumption on Y_{l, μ}. In particular, the pressure decrease is still favoured in correspondence to the hot spots because, in these regions, the final equilibrium fractions of muons and (anti)neutrinos are comparable for the two choices of Y_{l, μ}.
All Tables
Properties of nuclear EOS tables used in this work and of cold, spherically symmetric NSs for each nuclear EOS.
All Figures
Fig. 1. Outcome of the numerical simulation (BLh, 1.00) at 4.6 ms after merger. The left, central, and right panels show the rest mass density, the temperature, and the electron fraction on the equatorial plane, respectively. The black lines mark isodensity contours corresponding to 10^{12} g cm^{−3} (solid line), 10^{13} g cm^{−3} (dashed line), 10^{14} g cm^{−3} (dasheddotted line), and 10^{15} g cm^{−3} (dotted line). 

In the text 
Fig. 2. Particle fractions computed at T = 0 in neutrinoless βequilibrium for three different EOS: BLh (left), DD2 (centre), and SFHo (right). Black lines show the electron fraction computed in two different cases: equilibrium with muons (solid lines) or without muons (dashed lines). The green line corresponds to the fraction of protons at equilibrium with muons, while the red line corresponds to the muon fraction. We note that when muons are not included, the proton fraction coincides with the electron one (dashed black line). 

In the text 
Fig. 3. Particle fractions and chemical potential in the muon flavour sector computed by postprocessing the simulation (BLh, 1.00) at 4.6 ms after merger. In the upper panels, we show the fraction of μ^{−} (left), the fraction of μ^{+} (centre), and the difference between the muon fraction and the one inherited from the cold NSs, that is, (right). In the lower panels, we show the fraction of ν_{μ} (left), the fraction of (centre), and the muon chemical potential (rest mass subtracted) (right). As in Fig. 1, the black lines mark isodensity contours. 

In the text 
Fig. 4. Ratio of the net muon fraction over the net electron fraction, Y_{μ}/Y_{e}, on the equatorial plane obtained by postprocessing simulations (BLh, 1.00) (upper left), (BLh, 1.43) (upper right), (DD2, 1.00) (lower left), and (SFHo, 1.00) (lower right) at 5 − 6 ms after merger. As in Fig. 1, the black lines mark isodensity contours. 

In the text 
Fig. 5. Fractions of ν_{e}, , ν_{τ}, and chemical potentials of electrons, protons, and neutrons computed in postprocessing for (BLh, 1.00) at 4.6 ms after merger. In the upper panels, we show the fraction of ν_{e} (left), the fraction of (centre), and the fraction of ν_{τ} (right). In the lower panels, we show the electron chemical potential (rest mass subtracted) μ_{e} − m_{e} (left), the proton chemical potential (rest mass subtracted) μ_{p} − m_{p} (centre), and the neutron chemical potential (rest mass subtracted) μ_{n} − m_{n} (right). We note that the colour bar scales are different for the various chemical potentials. As in Fig. 1, the black lines mark the isodensity contours. 

In the text 
Fig. 6. Proton and neutron chemical potentials (rest mass subtracted) of BLh in the densitytemperature plane at fixed Y_{p} = 0.08. The proton fraction Y_{p} = 0.08 is the equilibrium proton fraction at 4.6 ms after merger around the point x = 0.64 km, y = 12.76 km, z = 0 km, where ρ ∼ 5 × 10^{13} g cm^{−3}, T_{sim} ∼ 10 MeV, and neutrinos dominate, while antineutrinos are suppressed (see also Figs. 1, 3 and 5). We notice that the colour bars are scaled differently. 

In the text 
Fig. 7. Degeneracy parameter of ν_{e} (first row) and ν_{μ} (lower row) on the equatorial plane obtained by postprocessing the simulation (BLh, 1.00). Three different time slices are shown corresponding from left to right to 4.6 ms, 6.7 ms, and 11.8 ms after merger. As in Fig. 1, the black lines mark the isodensity contours as in Fig. 1. 

In the text 
Fig. 8. Mean energy E_{ν} of trapped neutrinos on the equatorial plane computed by postprocessing simulation (BLh, 1.00) at 4.6 ms after merger. As in Fig. 1., the black lines mark isodensity contours. 

In the text 
Fig. 9. Pressure ratio P/P_{sim} for (BLh, 1.00) (left), (DD2, 1.00) (centre), and (SFHo, 1.00) (right) where P is obtained in postprocessing, while P_{sim} is the simulation pressure. For (BLh, 1.00), we explicitly mark the maximum P/P_{sim} ∼ 1.03 (white circle) and the minimum P/P_{sim} ∼ 0.94 (white cross) at ρ > 10^{14} g cm^{−3}. As in Fig. 1., the black lines mark isodensity contours. 

In the text 
Fig. 10. Change of particle pressure and fractions obtained in correspondence to the spacetime coordinates explicitly marked in the left panel of Fig. 9 and corresponding to P/P_{sim} = 0.94 (upper row) and P/P_{sim} = 1.03 (lower row). The bar plots in blue (left) give the difference between the equilibrium pressure computed in postprocessing and in the simulation, , where , and P_{tot} refers to the total pressure. The bar plots in red (centre) show the difference between the particle fractions computed in postprocessing and the ones from the simulation, , where , and e (μ) refers to the net electron (muon) fractions. The colourcoded plots (right) show P_{b} of BLh as a function of temperature and proton fraction for the same points. The density ρ is fixed to the value obtained in the simulation, that is, ρ = 4.3 × 10^{14} g cm^{−3} (upper row) and ρ = 1.55 ⋅ 10^{14} g cm^{−3} (lower row). The black star and the diamond mark are the values from the simulation (before postprocessing) and at equilibrium (after postprocessing), respectively. 

In the text 
Fig. 11. Pressure ratio P/P_{sim} computed by postprocessing (BLh, 1.00) (top) and (BLh, 1.43) (bottom) at 6.7 ms after merger. As in Fig. 1, the black lines mark the isodensity contours. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.