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

Neutrino-matter interactions play a key role in binary neutron star mergers. Thermodynamics conditions at the surfaces where neutrinos decouple from matter influence neutrino spectra, ultimately affecting the evolution of the remnant and the properties of the ejecta. In this work, we post-process results of general relativistic merger simulations employing microphysical equations of state and approximate neutrino transport to investigate the thermodynamics conditions at which weak and thermal equilibrium freezes out (equilibrium surfaces), as well as conditions at which the transition between diffusion and free-streaming regime occurs (diffusion surfaces). We find that the rest mass density and the neutrino energy are the most relevant quantities in determining the location of the decoupling surfaces. For mean energy neutrinos ($\langle E_{{\nu}_e} \rangle \approx 9~{\rm MeV}$, $\langle E_{{\bar{\nu}}_e} \rangle \approx 15~{\rm MeV}$, $\langle E_{{\nu}_{\mu,\tau}} \rangle \approx 25~{\rm MeV}$), diffusion surfaces are located around $10^{11}{\rm g~cm^{-3}}$ for all neutrino species, while equilibrium surfaces for heavy flavor neutrinos are significantly deeper (several $10^{12}{\rm g~cm^{-3}}$) than the ones of $\bar{\nu}_e$ and $\nu_e$ ($\gtrsim 10^{11}{\rm g~cm^{-3}}$). The resulting decoupling temperatures are in good agreement with the average neutrino energies ($\langle E_{\nu} \rangle \sim 3.15~T$), with the softer equation of state characterized by systematically larger decoupling temperatures ($\Delta T \lesssim 1~{\rm MeV}$). Neutrinos streaming at infinity with different energies come from very different regions of the remnant. The presence of a massive NS or of a BH in the remnant influences the neutrino thermalization process.


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 immedi-Send offprint requests to: e-mail address of corresponding author ately synthetize 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 gammaray 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 kilonova (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. Modelling 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 one hundred 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 lenghtscale over which termodynamics 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 neutron-to-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. [46,47,37,48,49,50,51,52,53,8,54]. In most of the cases, grey (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 and proto-neutron stars 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 these works, the location of grey 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 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,66,49,48].
In this work, we present a first quantitative study of the the thermodynamics conditions of matter at the neutrinosphere in BNS remnants. Neutrinos opacities and optical depths are calculated by postprocessing 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 Sec. 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 Sec. 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 Sec. 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 ij dx i dx j is an infinitesimal displacement along the chosen path, with g ij 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 behaviour is more relevant than the behaviour 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 to 7 directions (see e.g., [32] [35]). In our scheme, instead, we chose 17 directions, corresponding to 5 on-axis directions x + , x − , y + , y − , z + , (we disregard all z − directions, since all simulations considered in this work are performed with z−symmetry), 8 planar diagonals x + y + , x + y − , x − y + , x − y − , x + z + , x − z + , y + z + , y − z + , 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 directions. 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 l C = 5r 11 , 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 ∼ 30km. 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 anti-neutrinos (ν 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 where the indexes 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, while the equilibrium surfaces identify the conditions at which neutrino radiation decouples from the background medium. 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], ion-ion 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 dense matter, well below ρ th . For the inverse nucleon-nucleon bremsstrahlung and the neutrinoantineutrino 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 ν x 's. 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 quasi-circular 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 [50,8]. All the technical details for these simulations are given in [84,30] 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 3,024 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.74km representing our standard resolution. In Appendix A we study the sensitivity of our results by running both higher and lower resolution postprocessing 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 visible in Appendix B.
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 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
Ref.
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% [66,49,48] 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 surfaces 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 energy-integrated 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 [66,48]. 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 spectrumaveraged 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 approch 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 Sec. 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 k scat at low temperatures (T 1 MeV) visible 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 nucleonnucleon 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 nucleonnucleon 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 ok k abs , as visible in the central panel of Fig. 2. The fi-nal 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 dependences 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 Fig. 3 and Fig. 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. (1 − Y e ) E νe 2 Y e Eν e 2 . Thus, the averageν e 's start to stream freely from (marginally) smaller radii. The last snapshots in Fig. 3 and Fig. 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 contourare 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 the histograms below. In the early post-merger the gradient of 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 den-  sities and temperatures in the first 10 ms, while mantaining 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 Sec. 3.1), and trace back to the correlation between density and temperature inside the remnant, which becomes narrower for increasing time after merger [30].

Equilibrium surfaces
The surfaces where the freeze-out from weak and thermal equilibrium occurs for mean energy neutrinos (i.e. τ eq ( E ν ) = 2/3) are shown in Fig. 6 and Fig. 7. The difference in the location between the scattering and equilibrium surfaces defines the extension of the diffusion atmosphere, 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 dif-ferent neutrino species. In the case of ν e 's, 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 free streaming occurs. As a consequence, the thermodynamics conditions at the ν e equilibrium decoupling surfaces are close to the ones obtained at the scattering decoupling surfaces, as reported in Table 3. On the contrary, the relative paucity of free protons available to absorbeν e locates their equilibrium decoupling surfaces at slightly larger densities and temperatures than the ones of ν e 's. A qualitatively different behavior characterizes ν x 's. 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 ν e 's the distributions are practically identical to the one extracted for the scattering surfaces. Forν e 's, 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 ν x 's, 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, 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 8 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 Figures 8-9 we show mean densities, temperatures and electron fractions (with their respective standard deviation) around the diffusion surface as 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 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 cold weak equilibrium Y e , while the electron fraction inside the disk has been increased by neutrino emission and absorption processes up 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 in 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 signif-icant 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 ν e 's than forν e 's and ν x 's 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 Fig. 10-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 denisties 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 nucleon-nucleon bremsstrahlung) move the decoupling at systematically larger temperatures and densities. In the massive NS case, the decouplig 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, Eq. (4) and Eq. (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 ν e 's 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, even more, in the ν x case the assumption of matterradiation equilibrium even forτ diff decreases the mean energies     Fig. 11. Same as Fig. 10 around the equilibrium neutrino surfaces for the SLy4 model. In this case, the colored bands represent the neutrino energies for which no optically thick region is found.

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 consider 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 quasi-elastic scattering and charged current absorption reactions on free nucleons. For heavy-flavor neutrinos, weak and thermal equilibrium is guarateed by pair processes, like inverse nucleon-nucleon bremsstrahlung and ν-ν annihilation. Due to the strong dependence of the target and incident particle densities on mat- ter 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 for mean energy neutrinos happens between 4 and 5 MeV, while in the latter between 3 and 4 MeV.
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 ν e 's. On the other hand, the relative paucity of absorbing protons moves the equilibrium surfaces ofν e 's 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 heavy flavor 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 equiliubrim decoupling temperatures broadly satisfy the expected blackbody emission relation, i.e. E ν ∼ 3.15T , with a closer agreement in the case of the stiff EOS model DD2. Deviations from it largely depend on the neutrino emission 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 a softer EOS differences between ν e andν e equilibrium decoupling temperatures reduces, 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 trend 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 qualitatively 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, low energy 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 energy integrated 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.
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.
In light of our results, the accurate modelling 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 order of magnitudes 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.

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.
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 8 cells, increasing the possiblity 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 neighbouring 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 Table 4. In this table we present the mean values of density ρ, temperature T and electron fraction Ye 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 (ρ [g/cm 3 ]) and temperatures are expressed in MeV. For each quantity we also give its variance in the distribution. 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.