Thermodynamics conditions of matter in the neutrino decoupling region during neutron star mergers

In this work we investigate the thermodynamics conditions at which neutrinos decouple from matter in neutron star merger remnants by post-processing results of merger simulations. We find that the matter density and the neutrino energies are the most relevant quantities in determining the decoupling surface location. For mean energy neutrinos (∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 9, 15 and 24 MeV for νe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _e$$\end{document}, ν¯e\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\nu }_e$$\end{document} and νμ,τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _{\mu ,\tau }$$\end{document}, respectively) the transition between diffusion and free-streaming conditions occurs around 1011gcm-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{11}\mathrm{g}~\mathrm{cm}^{-3}$$\end{document} for all neutrino species. Weak and thermal equilibrium freeze-out occurs deeper (several 1012gcm-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{12}\mathrm{g}~\mathrm{cm}^{-3}$$\end{document}) for heavy-flavor neutrinos than for ν¯e\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\nu }_e$$\end{document} and νe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu _e$$\end{document} (≳1011gcm-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gtrsim 10^{11}\mathrm{g}~\mathrm{cm}^{-3}$$\end{document}). Decoupling temperatures are broadly in agreement with the average neutrino energies, with softer equations of state characterized by ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document}1 MeV larger decoupling temperatures. Neutrinos streaming at infinity with different energies come from different remnant parts. While low-energy neutrinos (∼3MeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sim 3~\mathrm{MeV}$$\end{document}) decouple at ρ∼1013gcm-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \rho \sim 10^{13}\mathrm{g}~\mathrm{cm}^{-3}$$\end{document}, T∼10MeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T \sim 10~\mathrm{MeV}$$\end{document} and Ye≲0.1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_e \lesssim 0.1$$\end{document} close to weak equilibrium, high-energy ones (∼50MeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sim 50~\mathrm{MeV}$$\end{document}) decouple from the disk at ρ∼109gcm-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho \sim 10^{9}\mathrm{g}~\mathrm{cm}^{-3}$$\end{document}, T∼2MeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T \sim 2~\mathrm{MeV}$$\end{document} and Ye≳0.25\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_e \gtrsim 0.25$$\end{document}. The presence of a massive NS or a BH influences the neutrino thermalization. While in the former case decoupling surfaces are present for all relevant energies, the lower maximum density (≲1012gcm-3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim 10^{12}\mathrm{g}~\mathrm{cm}^{-3}$$\end{document}) in BH-torus systems does not allow softer neutrinos to thermalize and diffuse.


Introduction
The merger of two neutron stars (NSs) results in the emission of a burst of gravitational waves (GWs), in the formation of a massive remnant (possibly collapsing to a black hole, BH) and in the ejection of matter in the interstellar medium; see [1][2][3] for some recent reviews. The remnant is usually formed by a central object surrounded by a massive disk. The properties of the disk, as well as the timescale for BH collapse of the central object (if gravitationally unstable), depend mainly on the NS masses and on the still uncertain nuclear equation of state (EOS); see e.g. [4][5][6][7][8][9][10]. Due to its high neutron content, the ejecta expelled in the interstellar medium can immediately synthesize heavy elements via the so-called rapid neutron capture process (r -process) nucleosynthesis [11][12][13]. The radioactive decay of these freshly synthetized elements powers a peculiar UV/optical/NIR electromagnetic transient called kilonova on a timescale ranging from a few hours up to months after the merger [14][15][16]. Moreover, the merger remnant has been long thought to be the central engine of short-hard gamma-ray bursts, one of the most energetic and elusive extragalactic events detected up to cosmological distances [17,18].
This picture was recently confirmed by the detection of GW170817, the first GW signal from a binary neutron star (BNS) merger reported by LIGO and Virgo collaborations [19]. Indeed, this detection was followed by the observation of a short gamma-ray burst, GRB170817A (although a very weak one, due to the off-axis relative orientation of the jet with respect to the Earth, [20,21]) and of a kilo-nova (AT2017gfo) produced by the merger [22,23]. During the first day, the kilonova was characterized by an intense quasi-thermal emission peaking in the UV and visible part of the electromagnetic spectrum, while on longer timescales it became redder and peaked at NIR frequences. Modeling of the associated light curves revealed the presence of more than one component in the ejecta; see e.g. [24][25][26][27][28][29]. The components differ by their physical properties and possibly by their non-isotropic angular distributions. The low photon opacity required to explain the early kilonova emission suggests the presence of matter that was significantly irradiated by neutrinos. In this case, if the electron fraction increases above ∼ 0.25 for a significant fraction of the ejecta then the production of r -process elements between the second and the third r -process peak is strongly suppressed. This observation has thus revealed the potential impact of neutrinos and weak interactions on the outcome and on the observables associated with BNS mergers.
During the coalescence the production of neutrinos of all flavors is boosted by the large temperatures (up to 100 MeV, [30]) reached by dense matter, once a significant fraction of the kinetic energy of the two NSs has been converted into internal energy; see e.g. [18,[31][32][33]. Deep inside the forming hot remnant, for densities in excess of 10 12 g cm −3 , the neutrino mean free path for absorption and scattering processes becomes smaller than any relevant length scale over which thermodynamics quantities change significantly. Thus neutrinos form a trapped gas, in equilibrium with the plasma, that diffuses out on the diffusion timescale. When absorption reactions and inelastic scattering processes become less effective, weak and thermal equilibrium freezes out. Neutrinos can still diffuse, due to the relatively large contribution of quasi-elastic scattering off free baryons to the total opacity. Only when lower densities (∼ 10 11 g cm −3 ) have been reached, radiation can finally stream out freely. The region where the decoupling occurs is often referred as neutrino surface, or neutrinosphere in analogy with the photosphere of approximately spherically symmetric systems emitting photons, e.g. stars. Even in free-streaming conditions, neutrinos can still irradiate matter in the low-density part of the remnant and in the expanding ejecta, potentially altering its neutronto-proton content through charged-current reactions involving electron neutrinos and antineutrinos, see e.g. [34][35][36][37][38][39].
Neutrino cross-sections depend significantly on the incoming neutrino energy: for energies in the range 1-100 MeV the dependence is approximately quadratic [40][41][42]. The relevance of this dependence is twofold. On the one hand, neutrinos of different energies decouple at different locations inside the remnant. Thus neutrino spectra recorded far from the merger site result from matter-neutrino interaction over a wide range of thermodynamics conditions. On the other hand, the properties of the expanding ejecta, and ultimately of the resulting kilonova, depends significantly on the electron (anti)neutrino spectra emerging from the remnant through the effect of the free-streaming neutrino absorption on the ejecta composition; see e.g. [43][44][45].
The multidimensional nature of compact binary mergers has represented a challenge in the production of quantitative models. Only in the last three decades multidimensional numerical simulations have shed light on the merger dynamics and on its many observables. The inclusion of weak interactions in the context of BNS mergers and of their aftermath has been done at different levels of sophistications, including light bulbs, leakage schemes, moment schemes, and Monte Carlo approaches; see e.g. [8,37,[46][47][48][49][50][51][52][53][54]. In most of the cases, gray (i.e. energy-integrated) treatments have been adopted, while spectral treatments are still less common due to their significantly higher computational costs. One of the key ingredients of any treatment of weak interactions are the neutrino opacities and their relation with the thermodynamical properties provided by the nuclear EOS. While in the context of core-collapse supernovae (CCSNe) and proto-neutron stars (PNSs) many sensitivities studies on the relevant neutrino physics have been performed, e.g. [55][56][57][58][59][60][61][62][63][64], very little have been said in the case of compact binary mergers. In the latter scenario the different physical conditions, as well as the different geometry of the system, can provide qualitatively different answers compared to the former. First studies about the properties of the neutrino surfaces in BNS merger remnants were presented in [31,32]. In this work, the location of gray neutrino surfaces was determined by evaluating the optical depth along a finite number of directions, usually assuming cylindrical symmetry. A similar approach was adopted also by [35], but in this case both the locations of dynamical and the thermal decoupling were computed for different neutrino energies. More sophisticated algorithms to compute neutrino optical depths in multidimensional context have been developed, e.g. [65], but never applied to detailed studies. We recall that the evaluation of the optical depth in radiation hydrodynamical simulations is also required in a large number of approximated neutrino treatments [46,48,49,66].
In this work, we present a first quantitative study of the thermodynamics conditions of matter at the neutrinosphere in BNS remnants. Neutrinos opacities and optical depths are calculated by post-processing numerical relativity simulations employing different nuclear EOSs and producing either NS or BH remnants. We identify the regions where the decoupling occurs and compute there the properties of the matter. In our analysis, we distinguish between the neutrino surfaces that denote the transition from the diffusion to the free-streaming regime, and the ones where weak and thermal equilibrium freezes out.
The paper is structured as follows: In Sect. 2 we introduce the concepts of optical depth and the neutrino opacities used in this work. Then we present our algorithm for the optical depth evaluation and the analysis procedures followed in this work. In Sect. 3 we present the results obtained in terms of the distribution of thermodynamic quantities in the neutrino decoupling region, both for mean energy neutrinos and for several neutrino energies. Finally, in Sect. 4 we summarize our results.

The optical depth
To evaluate the region where neutrinos decouple from matter, it is necessary to compute the optical depth for neutrinos at any point inside the simulation domain. The optical depth quantifies the global degree of opacity of matter to radiation between two points (A and B) along a certain path. If γ denotes the path between A and B, the optical depth is defined as where κ(s) is the opacity (corresponding to the inverse of the local mean free path) and ds = g i j dx i dx j is an infinitesimal displacement along the chosen path, with g i j being the local spatial metric. The physical interpretation of τ emerges from its definition: τ counts the average number of interactions that radiation particles experiences along γ . In radiation transport problems involving astrophysical objects, radiation is often produced in opaque regions from which it diffuses out on the diffusion timescale towards optically thin regions located at the boundary of the physical system. In this context, a statistical description of the radiation and of its global behavior is more relevant than the behavior of single radiation particles. Even though radiation is often produced isotropically on the microscopic scales, interaction with matter can change the propagation direction on macroscopic scales. On the one hand, if radiation particles move towards a region of increasing mean free path, they will more probably retain their direction and eventually move away freely from the production site. On the other hand, particles moving towards a region of decreasing mean free path will more likely interact with matter, changing their original propagation direction. This implies that, independently from the emission properties, macroscopically radiation moves preferentially towards regions of larger mean free path. Eventually, radiation particles reach regions where the local mean free path becomes larger than the relevant domain size and they will stream out of the system, practically without further interactions. If we are interested in the global radiation behavior, we consider any point inside the domain as starting point ( A) while the final point (B) can be any point on the boundary of the physical domain. According to the statistical interpretation presented above, among all the possible paths connecting A to the boundary, the most likely ways for radiation to escape are the paths that minimize the optical depth. We thus define the optical depth of a point x as where x b is any point of the boundary from which radiation escapes freely. The surface where radiation decouples from matter is defined as the region where τ ∼ 1 and it is referred as neutrino surface. If its curvature is not very pronounced, Eddington approximation applies and a neutrino surface is often referred as the surface where τ = 2/3. A common approach to effectively evaluate Eq. (2) is to select a certain number of direction moving away from x, evaluate τ along those, and then take the minimum value obtained across all directions. In the past this has been done already for Newtonian simulation, using cylindrical coordinates and 3-7 directions (see e.g., [32,35]). In our scheme, instead, we chose 17 directions, corresponding to 5 on-axis directions (we disregard all z − directions, since all simulations considered in this work are performed with z-symmetry), 8 planar diagonals and 4 full diagonals A representation of the directions in a Cartesian grid is given by Fig. 1. Even limiting the number of paths to 17, computing the integrals over the entire simulation domain is usually very expensive in large simulations covering thousands of km in each direction. However, the contribution to the optical depth given by very low-density matter is negligible. For this reason, when performing the integration, we reduce ourselves to a smaller region of the simulated volume centered around the remnant, where we expect to find the neutrino surfaces. In order to select this region, we define the characteristic length of the system l C as with r 11 being the average radius of a density isosurface at ρ th = 2 × 10 11 g cm −3 . In the simulations presented in this work r 11 ∼ 30 km. To test our approximation, we have computed τ with l C = 5, 6, 8 r 11 for a few selected cases, without finding significant differences in the thermodynamics properties of the decoupling region.

Neutrino opacities and BNS merger simulations
Neutrino opacities are provided by weak interaction processes and depends on the local matter properties, on the neutrino flavor and energy, and possibly on the neutrino radiation fields itself. Under the assumption of Nuclear Statistical Equilibrium (NSE), the thermodynamics state of matter is fully identified by its rest mass density, temperature, and electron fraction. In this work we distinguish between three independent neutrino species, namely electron neutrinos (ν e ), electron antineutrinos (ν e ) and a collective species for μ and τ (anti)neutrinos (ν x ). The most relevant neutrino-matter reactions considered in this work are listed in Table 1. We consider the corresponding absorptivity κ ab or scattering opacity κ sc , as they are defined in the Boltzmann transport equation.
On the one hand all processes equally contribute to the neutrino opacity relevant for the diffusion process. We thus define the diffusion opacity κ diff as Table 1 Weak reactions providing the neutrino opacities used in this work and references for their implementation. ν ∈ {ν e ,ν e , ν x } denotes a neutrino species with ν x referring to any heavy-lepton neutrino species. N ∈ {n, p} denotes a nucleon, A a generic nucleus (including α particles), e ± electrons and positrons

Reaction
References [70,74] ν +ν → e + + e − [41,55] N + N + ν +ν → N + N [75] where the indices r and s run over all the considered absorption and scattering reactions, respectively. On the other hand only a subset of reactions are effective in keeping the neutrino field in thermal and weak equilibrium with the plasma.
We define an equilibrium opacity as the geometrical mean between the diffusion opacity and the opacity only due to absorption processes (see e.g. [67] or [68] for analogous expressions): These two kinds of opacities provide two different optical depths, τ diff and τ eq , and two different kinds of neutrino surfaces. The diffusion surfaces mark the transition from the semi-transparent to the optically thin (free-streaming) regime. They can be identified as the last interaction surfaces, meaning that outside them neutrinos are likely not to interact anymore with matter, regardless of the nature of the reaction. On the other hand, the equilibrium surfaces identify the conditions at which neutrino radiation decouples from the background medium, in the sense that freeze-out from thermal and weak equilibrium occurs. This equilibrium is guaranteed by inelastic reactions that produce and absorb neutrinos. Since the latter form a subset of all reactions, κ diff ≥ κ eq and the equilibrium surfaces lie always inside the diffusion ones. In particular, if quasi-elastic scattering is efficient in providing opacity beyond the equilibrium surfaces, a scattering atmosphere, where neutrinos diffuse far from local equilibrium, can form. To compute energy-dependent neutrino opacities for neutrino absorption on nucleons and scattering off nucleons and nuclei we use the publicly available library NuLib [69]. In particular, we rely on expressions for the transport opacities and we include weak magnetism corrections [70], ionion correlations, form factor correction [71,72], and electron polarization correction [73] according to [74]. Differently from the standard code version, we do not include the effect of stimulated absorption in our calculations because it provides unphysically highν e opacities in cold, low-density matter, well below ρ th . For the inverse nucleon-nucleon bremsstrahlung and the neutrino-antineutrino pair annihilation we compute the reaction kernels according to [75] and [41,55], respectively. These kernels depend on the species and energies of both incoming neutrinos. Computations of the corresponding opacities would require the knowledge of the detailed neutrino distribution functions, since neutrinos are not only the colliding, but also the target particles. Since this information is in general not available, the absorption opacity for each neutrino species and energy is computed assuming the target neutrinos to be in equilibrium with matter (for example, in the case of ν e +ν e → e + +e − , we fix the ν e energy and we integrate over an equilibrium Fermi-Dirac distribution function forν e ). While this hypothesis is not correct in optically thin conditions, it is well verified in neutrino trapped regions. Since in this work we explore the intermediate, semi-transparent regime, the validity of this approach is a priori uncertain. In general, due to the pair nature of the reaction, the corresponding opacity drops quickly outside the neutrino surfaces even assuming equilibrium conditions for the target particles. Thus, even if not fully correct, their contributions to the absorption opacity in optically thin conditions in Eq. (2) is small and we do not expect a significant change in the calculations of the optical depths and in the determination of the neutrino surfaces. Moreover, these processes are relevant in the case of the ν x . Due to the presence of an extended scattering atmosphere, we expect heavy-lepton neutrinos to form a trapped gas also when pair processes become inefficient in keeping the radiation in equilibrium with matter. Thus, even if not accurate, we consider our approximation reasonable in the determination of the location of the neutrino surfaces for all neutrino species.
The thermodynamics conditions of matter used as input for the neutrino opacity calculations are taken from the output of (3+1)D numerical relativity simulations of BNS mergers. We consider equal mass binaries with M NS = 1.364 M (corresponding to a chirp mass of 1.188 M compatible with GW170817) and we make use of two different finite temperature, composition-dependent EOSs: DD2 [76,77] and SLy4 [78]. These tables assume NSE to determine the nuclear composition. We produce two distinct opacity tables for the two different simulations, using the appropriate EOS table as input. We solve the general relativistic initial data problem to produce irrotational BNS configurations in quasicircular orbit using the Lorene pseudo-spectral code [79]. The initial separation between the NS is set to 45 km, corresponding to ∼ 2−3 orbits before merger. Low temperature (T 0.1 MeV), ν-less weak equilibrium slices of the EOS are employed to construct the initial data. To correct for the presence of photons at low density we subtract their pressure contribution from the cold slices.
The (3+1)D evolutions are performed with the WhiskyTHC code [80][81][82][83], complemented by a leakage scheme to account for compositional and energy changes in the matter due to weak reactions involving ν e ,ν e , and ν x . Free-streaming neutrinos are emitted at an average energy and then evolved according to the M0 scheme introduced in [8,50]. All the technical details for these simulations are given in [30,84] which we refer to also for the employed grid setup. A detailed discussion of the thermodynamics properties of the matter during merger can be found in [30]. The DD2 simulation is the same discussed in [30], the SLy4 simulation is presented here for the first time. The domain covered by our simulation is a cube of size length 3024 km (assuming mirror symmetry along z) whose center is at the center of mass of the binary. It is resolved by a Cartesian grid hierarchy composed of 7 2:1 nested refinement levels. The finest refinement level, covering both NSs during the inspiral and the central remnant after merger, has a resolution of h 185 m. Since our optical depth calculation scheme works on uniform grids, we interpolate the simulation results onto one uniform grid of resolution h ∼ 0.74 km representing our standard resolution. In Appendix A we study the sensitivity of our results by running both higher and lower resolution post-processing analysis, and we find no significant differences. We repeat our analysis also using a different optical depth computation scheme, and we find again no significant difference, as shown in Appendix B. Due to the approximate nature of our neutrino transport scheme, the results we present in this paper carry a certain degree of uncertainty. Comparison with simulations performed with different neutrino schemes (e.g. [52], however, did not reveal significant differences in the density and temperature profiles. On the other hand, the electron fraction in the remnant is more sensitive to the neutrino transport, especially in the density interval where high-energy neutrino decoupling occurs. The two binary models have been selected in order to span the different outcomes of a BNS: the DD2 EOS model results in a massive NS remnant, which lives longer than the extent of the simulation, while the SLy4 EOS remnant collapses to BH around ∼ 12 ms after merger. We chose three time snapshots at which to evaluate the optical depth that were reached by both simulations. The timesteps correspond to 1, 10, and 20 ms after merger time (identified as the time corresponding to the maximum in the strain of the = m = 2 mode of the gravitational wave signal). Three dimensional profiles of the density, temperature, electron fraction, and spatial metric tensor were extracted for all the three snapshots and from both simulations.

Analysis strategy
Neutrino opacities have a significant dependence on the energy of the incoming neutrino. A crucial ingredient is thus the determination of the neutrino energy at which the opacity entering Eq. (2) should be evaluated.
Numerical simulations of BNS mergers including neutrino radiation provide values for the average energies of the neutrinos escaping to infinity. Despite the large variety of approximated schemes employed in these analysis, reported values usually agree within 10-20% [48,49,66].
As a first approach, we fix the neutrino energy to the following set of values compatible with those results: E ν e ≈ 9.34 MeV, Eν e ≈ 15.16 MeV, and E ν x ≈ 23.98 MeV, and we compute spectral optical depths and neutrino sur-faces for them. This approach assumes that the neutrino spectrum at infinity is mainly determined by the spectrum emerging from the neutrino surfaces. We determine the thermodynamics conditions of matter (density, temperature and electron fraction) at the neutrino surfaces and we use them to characterize the typical conditions at which the largest fraction of neutrinos decouples from matter. In our discretized domain, we identify the neutrino surface as the region where 0.5 ≤ τ ≤ 0.85. For each extracted quantity q we compute the volumetric mean q mean and to give an estimate of its distribution around the mean value we compute the corresponding standard deviation as: where the index i runs over all N cells belonging to the neutrino surface. Due to the actual distribution of conditions, we choose q to be log 10 ρ [g cm −3 ] , log 10 (T [MeV]) and Y e . As a second approach, we compute the optical depths and the corresponding neutrino surfaces for a large set of neutrino energies between 3 and 88.67 MeV, and analyze the decoupling conditions as a function of the neutrino energy. This energy range covers the most relevant part of the neutrino spectra emerging from a NS merger. Within this analysis we will investigate how different the decoupling thermodynamics conditions are for neutrinos of different energies, regardless of their importance in the emerging spectrum.
Finally, we compute the optical depths using energyintegrated opacities. This case is relevant since the few compact binary merger simulations with the M1 transport schemes for neutrinos in the literature are performed assuming a gray approximation [48,66]. The latter consists in prescribing a certain functional form for the neutrino distribution function, thus removing the 3N degrees of freedom represented by N energy groups for each species. When neutrinos are coupled to matter, weak and thermal equilibrium drives their distributions to an isotropic Fermi-Dirac distribution: where E ν is the neutrino energy, μ ν the neutrino chemical potential, k b is the Boltzmann constant, and T is the temperature of the fluid at absorption (or scattering) point. We evaluate the neutrino chemical potentials at equilibrium as where ν n , ν p and ν e are the relativistic neutron, proton and electron chemical potentials, respectively. We introduce a spectrum-averaged opacity defined as:κ By calculating this average, we compute the typical opacity locally experienced by neutrinos at equilibrium and diffusing from optically thick towards optically thin regions. Since this approach assumes that neutrinos are in equilibrium condition with matter, it is more appropriate for neutrinos in optically thick conditions and more in line with the equilibrium optical depth, τ eq . The average used in this work, Eq. (9), is slightly different from the average computed in some gray neutrino transport schemes, e.g. [66], where the factor E 2 is replaced by E 3 . In the latter approach, more emphasis is put on the energy transport (I ν ∝ E 3 ) rather than on the particle transport. We do not expect a qualitative difference between the two approaches, even if the energy transport opacities are usually larger than the particle transport ones.

Opacity variability
We start our analysis by exploring the opacities for the different neutrino species. In Fig. 2, we present the dependency of the absorption and scattering opacities (k abs and k scat , represented by solid and dashed lines, respectively) on rest mass density, temperature and electron fraction of the medium for the three neutrino species. We consider thermodynamics conditions close to the regions where we expect mean energy neutrinos to dynamically decouple from matter (ρ th = 2 × 10 11 g cm −3 , T th = 4 MeV and Y e,th = 0.1; see for example [35]). In particular we vary one of the three variables and keep the other two fixed. For concreteness, we fix the energy of the incoming neutrinos to the mean energies E ν reported in Sect. 2.3 and we use the DD2 EOS to compute nuclear properties, abundances and chemical potentials for matter in NSE.
Due to the dominant role of quasi-elastic scattering off free nucleons, scattering opacities present a very similar behavior among the three neutrino species. Since in very good approximation κ scat ∝ n b E 2 ν (with n b being the free baryon particle density), k scat is approximately proportional to ρ, while the larger (smaller) scattering opacity observed for ν x (ν e ) results from the larger (smaller) neutrino mean energy. Due to the quasi-elastic nature of the dominant scattering process and to the isotropic spatial distribution of all interacting particles, k scat is rather insensitive to variations in temperature and electron fraction. Weak magnetism, recoil and nucleon form factors introduce differences in the cross-sections possibly sensitive to the neutron-to-proton content, but on a scale much smaller than the variation induced by the change in nucleon number density and in the neutrino energy. The increase of The thick lines represent the full opacities, while the thin ones are the single contributes to the absorption opacity, including inverse nucleon-nucleon bremsstrahlung (NNB), ν-ν pair annihilation, and absorption on free nucleons k scat at low temperatures (T 1 MeV) shown in the central panel is related with the appearance of bound nuclei, instead of free baryons, that provide a more effective coherent scattering; see e.g. [41,74].
In the case of k abs , the dependence on the thermodynamics conditions is more complex, due to the presence of several processes characterized by different target particles. In Fig. 2, we present also the different contributions to the absorption opacity for the three neutrino species, including the inverse nucleon-nucleon bremsstrahlung, the absorption on free nucleons, and the annihilation of neutrino-antineutrino pairs. As visible in the left panel, for densities lower than ∼ 10 12 g cm −3 the neutron richness and the larger Q-value favor ν e captures on neutrons with respect toν e captures on protons. 1 For larger densities, Pauli blocking for the more degenerate final state particles drastically reduces the opacity due to charged-current absorption processes. For densities in excess of 10 13 g cm −3 , assuming again temperatures around a few MeV, the inverse nucleon-nucleon bremsstrahlung becomes the dominant absorption process. Depending on the neutrino chemical potentials at equilibrium, the abundance of the target neutrinos can significantly change. For example, for the explored conditionsν e are suppressed by degeneracy at high densities and this translates in a much smaller inverse bremsstrahlung opacity for ν e compared to ν x andν e . Variations in temperature have also a significant effect on k abs , as shown in the central panel of Fig. 2. The final state blocking effect that characterizes charged-current absorption reactions is reduced once the temperature increases between one and a few MeV, while for temperatures in excess of ∼ 10 MeV pair processes of thermal origin (like ν-ν annihilation) dominate. In particular, the corresponding inverse mean free path is roughly proportional to the Y ν Yν σ ν+ν , which for thermally produced neutrinos translate in a strong dependence on the matter temperature, ∝ T 8 .
Changes in the electron fraction also introduce significant dependencies on the k abs contributions. They are ultimately due to changes in the number densities of the target particles in the case of charged-current absorption processes, and to variations in the electron chemical potential for neutrino pair annihilation.

Conditions at the decoupling surfaces for mean neutrino energies
We move to the analysis of the neutrino surfaces and the corresponding thermodynamics conditions obtained for the average neutrino energies E ν .

Diffusion surfaces
In Figs. 3 and 4 we present the contours of the diffusion surfaces (τ diff = 2/3) at three different snapshots for the DD2 and SLy4 models, respectively. For each snapshot we show the rest mass density and temperature (respectively on the left and right half of each panel) in the meridional (top row) and equatorial (bottom row) planes. While around 10 ms after merger the system has already reached an almost axisymmetric configuration, in the first snapshot (∼ 1 ms after merger) the density and temperature profiles reveal a still highly dynamical evolution. The latter is characterized by the presence of spiral arms that rotate and shock matter in the mantle forming around the merging cores.
Due to the quadratic dependence on the neutrino energy of the elastic scattering off nucleons, the neutrino surfaces of heavy-flavor neutrinos are the most extended. At the same time, despite the lowest ν e mean energy, the more abundant free neutrons provide larger opacities to ν e than toν e , i.e.
Thus, the averageν e start to stream freely from (marginally) smaller radii. The last snapshots in Figs. 3 and 4 show two different remnants: while the DD2 massive NS still survives, in the SLy4 simulation a BH has formed. In Fig. 4 we highlight the apparent horizon (AH) with a white contour and black fill. As expected, when a BH has formed the neutrino surface shape changes. In the case of ν e andν e an outer and inner contour are both present, resulting in a toroidal-shaped neutrino surfaces. The heavy-flavor neutrinos instead only show an outer contour, while the inner contour was swallowed inside the AH. 2 In the DD2 case, where the massive NS is still surviving, the degree of axisymmetry and the disk compactness increase moving from 10 to 20 ms. In terms of neutrino surfaces, we notice that, for more compact disks, corresponding to late time configurations or to the softer SLy4 EOS, the distance between the outer electron neutrino and antineutrino surfaces decrease, due to the steeper density gradients.
The hierarchy observed in the neutrino surface locations is also visible in the mean densities and temperatures recorded at the diffusion surfaces and summarized in Table 2. The mean values obtained for log 10 (ρ) for each model and at any snapshots are well representative, since their spreads are rather small compared to the mean values (a few percent) and to the large range over which the matter density varies. Moreover, for each neutrino species they are rather insensitive to the model and to the time evolution (at least, within the uncertainty represented by their spreads). More variability is instead observed in the temperature distributions: the mean temperatures are significantly larger for the SLy4 model (5)(6) MeV) than for the DD2 model (3-4 MeV), as expected due to the larger temperatures that characterize soft EOS remnants. Moreover, in both cases the spreads around the means are of the order of 1/3 of the mean value just after merger and reduce down to ∼ 1 MeV at later time (corresponding to the 10 and 20 ms snapshots). This is a consequence of the evolution of matter distribution in the density-temperature plane; see [30] and the histograms below. In the early post-merger the gradient of the temperature around the diffusion surfaces can be quite steep, due to the hot tidal tails extending from the remnant into the forming accretion disk. This feature is present in both EOS models, even though in the DD2 model the tidal tails are more defined and develop stronger gradients with respect to the SLy4 model. This leads to a higher variability in the estimation of the mean temperature in the transient region at early stages of the remnant evolution. At ∼ 10 ms after merger, the outflow of matter from the remnant injects enough hot and dense matter into the disk such that the transient region for the diffusion surfaces is pushed further out and the tidal tails become less prominent. For times larger than 10 ms the remnants are more axisymmetric, and the hot tidal tails transfer heat to the material in the transient region more uniformly, resulting in almost stationary mean temperatures and lower spreads. This is especially clear in the SLy4 model, where the temperature spread in the decoupling volume is halved after 9 ms. A similar trend is also observed for Y e at the diffusion neutrino surfaces: remnants obtained by softer EOSs are characterized by larger Y e at the decoupling surfaces, while the more axisymmetric structure obtained at late times ( 10 ms after merger) reduces significantly the associated spreads.
To reveal possible correlations between the decoupling rest mass density and temperature, in Fig. 5 we present histograms describing how the volume around the neutrinosphere is distributed in terms of the matter density and temperature. We consider each of the three snapshots of the DD2 model. More specifically, we show the volumetric amount of matter at specific (ρ, T ) conditions normalized to the total volume of the neutrinosphere region V ns . These histograms display two important features. First, the initial large spreads in density and in temperature are initially weakly correlated, but the latter (for a fixed density) diminishes noticeably when moving further away from merger time. Second, at late times there is a clear correlation between the density and temperature conditions. In particular, the bulk of mass shifts towards lower densities and temperatures in the first 10 ms, while maintaining a noticeable spread in temperature   in the high density regions. This feature persists for the DD2 model even at 20 ms, while in the SLy4 case (not shown in the figure) the average temperature (at any given density) steadily increases with density, while the temperature spread becomes narrower. This difference is due to the presence of a massive remnant in the DD2 model, whose outflow in the inner disk region keeps a higher variability in density and temperature profiles in the high density part of the decoupling region. The correlation between the matter density and the temperature reflects the dominant dependency on ρ both for k abs and k sc (see Sect. 3.1), and it can be traced back to the correlation between density and temperature inside the remnant, which becomes narrower for increasing time after the merger [30].

Equilibrium surfaces
The surfaces where the freeze-out from weak and thermal equilibrium occurs for mean energy neutrinos (i.e. i.e. the region of the disk where neutrinos still diffuse due to scattering processes up to their diffusion surfaces, but out of local thermodynamics and weak equilibrium. Regardless of the EOS, we find an order inversion between the neutrino species with respect to the scattering surfaces: heavy-flavor neutrinos now decouple first, while the electron neutrinos remain longer in equilibrium. This is due to the different behavior of absorption processes involving different neutrino species. In the case of the ν e , the dominant equilibration process is the absorption on free neutrons. Due to the high neutron abundance inside the disk, this process effectively provides neutrino opacity almost up to the point where neutrino  Table 3. On the contrary, the relative paucity of free protons available to absorbν e locates their equilibrium decoupling surfaces at slightly larger densities and temperatures than the ones of the ν e . A qualitatively different behavior characterizes the ν x . In this case, thermal and weak equilibrium is guaranteed by pair processes, e.g. ν +ν → e + + e − . Due to their strong temperature dependence, the corresponding equilibrium surfaces shrink deeply inside the remnant, at densities ∼ 10 13 g cm −3 where temperatures are usually in excess of 8 MeV. We also notice that after BH formation (i.e. in the SLy4 model at 20 ms) there is no region where τ eq,ν x > 2/3, which means that even for E ν x 20 MeV the combination of scattering and pair processes that enters Eq. (5) requires matter with densities around ∼ 10 13 g cm −3 to reach high enough opacities for heavy-flavor neutrinos to be significantly absorbed. Since after BH formation the density in the disk drops below a few times 10 12 g cm −3 there is no region in the disk where equilibration can occur for heavy neutrinos with average energy.
We have also inspected the volumetric distribution of matter across the equilibrium neutrino surfaces in the (ρ, T ) plane and for the DD2 model. In the case of the ν e the distributions are practically identical to the one extracted for the scattering surfaces. For theν e , the correlation between ρ and T is even tighter than for the scattering surfaces, even if the surfaces are located at slightly larger densities. Qualitative difference are visible in the case of the ν x , where the more relevant role of the matter temperature in the decoupling conditions emerges from the distributions, in particular at late times.
Finally, the values of the temperature at the equilibrium decoupling surfaces allow us to test our initial assumptions about the values of the mean neutrino energies. In fact, neutrinos decoupling first are emitted with higher energies, while neutrinos decoupling at lower densities and temperatures will have a lower average energy. More precisely, the expected relation between the mean energies and the decoupling temperature, ν ∼ F 2 (0) T (where F 2 (0) ≈ 3.15 is the Fermi integral of order 2 evaluated for 0 degeneracy parameter), is approximately verified. The agreement is better for the DD2 model, while the softer SLy4 EOS produces a hotter remnant and thus possibly harder neutrino spectra emerging from the equilibrium decoupling surfaces.

Neutrino energy dependence of the decoupling conditions
We now turn to the study of dependency of the neutrino surfaces (and therefore of their thermodynamics properties) on the neutrino energy. To do that we have extracted the diffusion and the equilibrium optical depths at 20 ms after merger for eight different energy bins between 3 and 80.33 MeV. We have analyzed both models and for each of them we have obtained different sets of thermodynamics conditions for each neutrino species and optical depth kind.

Diffusion surfaces
In Figs. 8 and 9 we show mean densities, temperatures and electron fractions (with their respective standard deviation) around the diffusion surface as a function of the neutrino energy, for the DD2 and SLy4 model, respectively. Due to the ∼ n target E 2 ν dependence of the opacity for nucleon scattering and absorption on free baryons, lower energy neutrinos decouple at higher rest mass densities. Since inside the remnant matter the temperature increases for increasing density (at least, up to nuclear saturation density, [30]), the decoupling temperature also decreases monotonically with E ν . Matter deep inside the remnant is closer to weak equilibrium Y e , while the electron fraction inside the disk has been increased by neutrino emission and absorption processes up   Fig. 9 Same as Fig. 8 around the diffusion neutrino surfaces for the SLy4 model. In this case, the colored bands represent the neutrino energies for which no optically thick region is found to 0.30-0.35 in the outer regions of the disk. Thus, Y e at the decoupling surface for high-energy neutrinos is significantly larger than the one of low-energy neutrinos. The two orders of magnitude explored in E ν translate into four orders of magnitude in decoupling density. Deviations from this leading dependence are visible in the low-energy tail, due to final state Pauli blocking and to second order effects in the interaction cross-section, e.g. finite electron mass correction. Moreover, for a fixed neutrino energy, the more significant contribution to the total opacity provided by absorption processes, and in particular by the more abundant neutrons, results in a larger inverse mean free path for the ν e than for theν e and the ν x and, consequently, larger decoupling densities and temperatures.
While in the DD2 model (comprising a massive NS in the center) it is always possible to find a scattering decoupling region for all investigated neutrino energies, this is not always the case in the SLy4 model. After the BH has formed, the high density part of the remnant falls inside the AH. In Fig. 9 we do not find a diffusion surface for eitherν e or ν x with energies of 5 MeV or lower. For ν e we are able to find a small region that satisfies the condition 0.5 < τ diff < 0.85, but ∼ 90% of it comes from matter at τ < 2/3 located in a rather small volume.

Equilibrium surfaces
In Figs. 10 and 11 we show the result obtained for calculations at different energies of the equilibrium surfaces. Many features observed for the scattering surfaces can be observed also for the equilibrium ones, including the significant dependence on the matter rest mass density. Since k eq ≤ k diff , equilibrium decoupling happens always at larger densities and temperatures than the diffusion decoupling, for every neutrino energy. Differences in the decoupling densities are of the order of a few, while they are of 10-15% in the decoupling temperatures. Significant deviations from the ν e andν e behavior can be observed in the case of ν x equilibrium surfaces, for which the relevant absorption processes (neutrino pair annihilation and inverse nucleonnucleon bremsstrahlung) move the decoupling at systematically larger temperatures and densities. In the massive NS case, the decoupling conditions of the three neutrino species become similar for low-energy neutrinos. This is an indication that at those conditions pair absorption processes are the most relevant reactions also for ν e andν e . However, due to the low neutrino degeneracy, to the more abundant neutrons and to the larger abundance of electron antineutrinos for ρ 10 13 g cm −3 , T 10 MeV and Y e ∼ 0.07, electron neutrinos decouple again at slightly lower rest mass densities and temperature thanν e and ν x . In the BH remnant case, the lack of equilibrium decoupling surfaces for low-energy neutrinos is more severe. This holds especially true for ν x , whose equilibrium reactions require very high temperatures and rest mass densities to be able to keep low-energy neutrinos in thermal equilibrium. This is visible in the more extended shaded areas in Fig. 11.

Energy-integrated approach
Finally, to test the robustness and the coherence of the different approaches used in the calculations of the neutrino opacities, we perform the computation of the energy-integrated optical depthτ for both opacity prescriptions, Eqs. (4) and (5), and we compare them with the energy-dependent calculations evaluated for the mean neutrino energies at infinity, τ ( E ν ).
The energy-integrated approach assumes the neutrino radiation to be in equilibrium with matter. Thus we expect it to be more coherent with the τ eq case. Indeed, in Fig. 12 we show τ eq ( E ν ) alongsideτ eq . For concreteness we consider the DD2 model at each of the three snapshots. The location of the neutrino surfaces is remarkably close among the two approaches, in particular after the initial transient phase. We additionally note that the actual values of τ inside the neutrino surfaces can differ significantly among the two cases (and in particular for the ν e case). This is a consequence of the constant neutrino energy used in the energy-dependent approach, while in the energy-integrated one the relation between the neutrino energy and the local matter temperature and chemical potentials is more correctly taken into account. When considering the diffusion optical depth, the hypothesis of equilibrium between matter and radiation is not always verifies. Consequently in Fig. 13, where we show the energy-dependent and the energy-integrated results for the diffusion optical depth, we observe larger discrepancies between the two approaches. In particular, for the ν e the relevance of the absorption on free neutrons guarantees that the diffusion and equilibrium neutrino surfaces stay close, also in the energy-integrated case. However, in theν e case and, moreover, in the ν x case the assumption of matter-radiation equilibrium even forτ diff decreases the mean energies and lowers the opacity. Thus, the neutrino surfaces are apparently located deeper inside in comparison with the expected locations.

Summary
In this paper we have studied the thermodynamical properties of matter in the neutrino decoupling region of binary neutron star mergers remnants. We have considered two binary neutron star models, characterized by the same total mass (2.728 M ) and by two different nuclear EOSs (DD2 and SLy4), both compatible with present nuclear and astrophysical constraints. The stiffer DD2 EOS supports a long-lived remnant, while the softer SLy4 EOS produces a black hole roughly ∼ 10 ms after merger. In order to span the dependence on the merger dynamics and on the nature of the remnant we have considered three snapshots at times 1, 10, and 20 ms after merger. We have investigated two kinds of decoupling surfaces: the surfaces from where neutrinos stream freely (diffusion surfaces) and the ones from where thermal and weak equilibrium freezes out (equilibrium surfaces). To account for the dependence of on the neutrino energy, we have considered both neutrinos mean energies (as previously obtained by numerical simulations including neutrino transport) and a large range of binned neutrino energies. Our major findings are that the rest mass density and the neutrino energies are the most relevant quantities in determining the location of the decoupling surfaces and, consequently, the relevant thermodynamical conditions at decoupling for ν e andν e . This is due to the dominant role of quasielastic scattering and charged-current absorption reactions on free nucleons. For heavy-flavor neutrinos, weak and thermal equilibrium is guaranteed by pair processes, like inverse nucleon-nucleon bremsstrahlung and ν-ν annihilation. Due to the strong dependence of the target and incident particle densities on matter temperature, the latter becomes a relevant thermodynamics quantity as well.
For each model, the diffusion decoupling surfaces of mean energy neutrinos are characterized by similar conditions among the different neutrino species. This is due to a compensation effect between the larger mean energy associated to heavy-flavor neutrinos (k scat ∝ E 2 ν ) and the larger abundance of absorbing neutrons relevant for ν e . Typical rest mass densities at the corresponding diffusion surfaces are ∼ 10 11 g cm −3 , with slightly smaller values obtained for ν x . Temperatures can significantly differ based on the properties of the underlying nuclear EOS. Softer EOSs, like SLy4, produce hotter remnant than stiffer ones, like DD2. In the former case, the transition from diffusion to free-streaming regimes Qualitative differences are observed when moving from diffusion to equilibrium surfaces for mean energy neutrinos. On the one hand, due to the neutron richness of the plasma, thermodynamics conditions at the equilibrium and diffusion surfaces are very close for the ν e . On the other hand, the relative paucity of absorbing protons moves the equilibrium surfaces of theν e deeper inside the remnant, at rest mass densities of a few times 10 11 g cm −3 . Consequently, theν e equilibrium decoupling temperature increases to 5 MeV in the DD2 case and to 6 MeV in the SLy4 case. For heavyflavor neutrinos, matter temperature becomes the most relevant quantity and the equilibrium decoupling happens around 10 MeV. Inside the remnant profile this corresponds to a rest mass density in excess of 10 12 g cm −3 . Our results for the equilibrium decoupling temperatures broadly satisfy the expected blackbody emission relation, i.e. E ν ∼ 3.15 T , with a closer agreement in the case of the stiff EOS model DD2. Deviations from it largely depend on the neutrino emis-sion coming from the volume outside the neutrino surfaces (which is in fact more relevant for hotter remnants). We also notice that in the case of softer EOS differences between ν e andν e equilibrium decoupling temperatures reduce, as also observed in [38].
We have generalized our analysis by considering neutrinos with a broad range of energies, between 3 and 88 MeV, covering the most relevant part of the emitted spectrum, but irrespective of their relevance. Due to the strong dependence on E ν of all cross-sections, such a large energy range translates in broad ranges of decoupling radii and thermodynamics conditions. In particular, the rest mass density at the diffusion decoupling surface increases between a few times 10 8 g cm −3 for hard neutrinos to 10 13 g cm −3 for soft ones, with a systematic increase going from ν e toν e , and finally to ν x . The corresponding decoupling temperatures cover one order of magnitude, moving from ∼ 2 to ∼ 10 MeV. Similar trends are observed also for electron (anti)neutrinos in the case of the equilibrium surfaces, just for densities and temperatures larger by a factor of a few and 10−15%, respectively. A qual-itatively different behavior is observed for the equilibrium decoupling of heavy-flavor neutrinos, for which temperatures in the range 5-12 MeV are observed, usually increasing the equilibrium decoupling density by one order of magnitude. Important differences can be observed by comparing the DD2 and the SLy4 models. In addition to the larger decoupling temperatures and electron fractions that characterize the SLy4 model (which are a consequence of the hotter and less neutron-rich remnant emerging from softer EOS models), the presence of a massive NS or of a BH in the center largely affects the possibility for the low-energy part of the spectrum to thermalize. Indeed, after BH formation, lowenergy neutrinos do not have decoupling surfaces. This is a consequence of the reduction of the maximum rest mass density and temperature inside the remnant once an AH has formed. The effect is more severe forν e (E ν 10 MeV) and even more for ν x (E ν 40 MeV) than for ν e (E ν 5 MeV).
As a final check, we have computed the optical depths (and the corresponding decoupling surfaces) using an energyintegrated approach that assumes everywhere equilibrium between the neutrino field and the matter. We have verified that the equilibrium neutrino surfaces obtained with the energy-dependent approach (evaluated for the mean neutrino energies) and with the energy-integrated one lie very close, while relevant differences are observed in the case of the diffusion surfaces.
All the results we have reported are robust with respect to the optical depth calculation algorithm and to the grid resolution. However, all our simulations employ an approximate neutrino transport. While the present leakage+M0 treatment is expected to catch the dominant cooling features of neutrino emission, the lack of an explicit modeling of trapped neutrinos and uncertainties in the challenging radiation transport can affect more significantly the electron fraction inside the remnant. Thus, we expect the values of the density and temperature at the decoupling surfaces to be more robust than the ones of the electron fraction.
Our study represents one of the first systematic investigations of neutrino opacities in binary neutron star mergers and of the thermodynamics conditions occurring where neutrino decouple from matter. Our results have revealed how the large variety of thermodynamics conditions inside the remnant, combined with its disk geometry, translates in broad ranges of decoupling conditions in all the relevant thermodynamics quantities. These can influence the properties of the emitted spectrum and ultimately of many multimessenger observables through the evolution of the remnant and of the ejecta. Matter conditions predicted for BNS mergers share many similarities with the ones experienced by matter in the gravitational collapse of the iron core of massive stars. Indeed, also in CCSNe neutrinos of all flavors are copiously produced and become trapped at large enough densities on timescales much larger than the dynamical (and even the explosion) timescale. Due to the high sensitivity of CCSN explosions to weak reactions, neutrino-matter interactions in CCSNe have been deeply investigated in many studies, e.g. [56][57][58][59][60][61]. In particular, the densities and temperatures where neutrinos stream freely or decouple from matter are similar to the ones we found in BNS mergers. During the collapse before core bounce, Y e 0.3 and the maximum core temperature (and thus the neutrino mean energies) stay low, a few MeV. Under these conditions, decoupling occurs at higher densities (e.g., several times 10 11 g cm −3 for mean energy ν e ), lower temperatures and higher Y e . After core bounce, during the accretion phase, matter and neutrino temperatures increase at the surface of the forming PNS. As a consequence, the decoupling densities and temperatures of CCSNe become more similar to the ones expected in BNS mergers. In the case of very massive zero-age main sequence stars ( 25M ), even the maximum temperatures are comparable to the one observed in BNS mergers. During the subsequent PNS cooling phase, decoupling temperatures progressively decrease and the neutrino surfaces sink again towards larger densities. However, the most noteworthy difference between neutrino surfaces in CCSNe and BNS mergers resides in the surface geometry: while in CCSNe they are all spheroidal, in BNS merger remnants containing a massive NS they are close to spherical geometry for low-energy neutrinos and increase to a disk-like shape inside the disk for larger neutrino energies.
In the light of our results, the accurate modeling of neutrinos and of their emitted spectrum in BNS mergers requires the consistent inclusion of neutrino-matter interactions for a broad range of nuclear conditions, characterized by more than four orders of magnitude in rest mass densities, one order of magnitude in temperatures, and a factor of a few in electron fraction. This would require the usage in merger simulations of more consistent and accurate weak reaction rates, e.g. [85][86][87][88]. This represents a great challenge in nuclear and neutrino physics for the years to come. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A scheme convergence
In order to check the robustness of our scheme, we extracted the optical depth over grids with different resolutions (interpolated from the same simulation results). We call standard resolution (SR) the one used throughout the paper (h ∼ 0.74 km); we then picked as our best resolution (HR) a grid with h ∼ 0.49 km and as our worst resolution (LR) a grid with twice the spacing of our SR (h ∼ 1.48 km). In order to make a more informed convergence study, we also considered a medium resolution (MR) grid, with h ∼ 1.11 km. Table 4 In this table we present the mean values of density ρ, temperature T and electron fraction Y e in the regions with opacity 0.5 < τ < 0.85 using Eq. (4) for four different resolutions. The considered model is the DD2 at 20 ms after merger. Densities are given as log 10 Table 5 In this table we present the mean values of density ρ, temperature T and electron fraction Y e in the regions with opacity 0.5 ≤ τ s ≤ 0.85 using Eq. (4) for both models at different stages of the post-merger evolution. In the Sly4 case, at t ∼ 20 ms, the remnant already collapsed to BH. Densities are given as log 10 (ρ [g cm −3 ]) and temperatures as T [MeV]. For each quantity we also give its variance in the distribution DD2 SLy4 t ∼ 10 ms t ∼ 20 ms t ∼ 10 ms t ∼ 20 ms log 10 (ρ νe ) 11.17 ± 0.27 11.21 ± 0.26 11.04 ± 0. 25  Y νx e 0.14 ± 0.03 0.14 ± 0.01 0.24 ± 0.02 0.18 ± 0.02 In Table 4 we give the mean values for the main hydrodynamic quantities in the transition region for the three resolutions considered. We can see how the increase in resolution translates into an increase in the mean values of rest mass density and temperature in the transient region. We expect this to happen due to how optical depth is evaluated. The higher the number of cells in the region where rest mass density varies quickly means that a cell that previously occupied a volume V with density ρ now is split into eight cells, increasing the possibility that at least one of them has a density ρ 1 < ρ that is just enough to push the value of τ for that cell below the threshold. This only works towards decreasing τ , since we always pick the minimum optical depth between all directions considered, and therefore is sufficient that this happens in one direction to decrease τ , while for this mechanism to increase τ the new grid should have all the neighboring cells to end up with higher ρ. This is unlikely given the distribution of density in the system.

B Scheme comparison
To test the robustness of our results with respect to the optical depth computations, we repeat the calculation of the decoupling thermodynamics conditions for mean energy neutrinos using MODA [65]. For concreteness, we perform the analysis for the last two simulation snapshots (10 and 20 ms after merger) and we focus on the diffusion optical depth. In Table 5 we summarize the thermodynamic decoupling conditions. All the results are compatible with our analysis within uncertainties.