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

TL;DR

This study investigates the thermodynamic conditions at neutrino decoupling surfaces in neutron star mergers, revealing how density, neutrino energy, and remnant composition influence neutrino spectra and thermalization.

Contribution

It provides detailed analysis of decoupling and diffusion surfaces using post-processed merger simulation data with microphysical equations of state.

Findings

Diffusion surfaces are around 10^{11} g/cm^3 for all neutrinos.

Heavy flavor neutrinos decouple at deeper layers (~10^{12} g/cm^3).

Decoupling temperatures correlate with average neutrino energies and EOS softness.

Abstract

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.

Figures

Tables

Table 1. Table 1: Weak reactions providing the neutrino opacities used in this work and references for their implementation. ν ∈ { ν e , ν ¯ e , ν x } 𝜈 subscript 𝜈 𝑒 subscript ¯ 𝜈 𝑒 subscript 𝜈 𝑥 \nu\in\{\nu_{e},\bar{\nu}_{e},\nu_{x}\} denotes a neutrino species with ν x subscript 𝜈 𝑥 \nu_{x} referring to any heavy-lepton neutrino species. N ∈ { n , p } 𝑁 𝑛 𝑝 N\in\{n,p\} denotes a nucleon, A 𝐴 A a generic nucleus (including α 𝛼 \alpha particles), e ± superscript 𝑒 plus-or-minus e^{\pm} electrons and positrons.

Reaction	Ref.
${\nu }_e+n\to p+e^{-}$	Burrows:2004vq ; Horowitz:2001xf
${\overline{\nu }}_e+p\to n+e^{+}$	Burrows:2004vq ; Horowitz:2001xf
$\nu +\overline{\nu }\to e^{+}+e^{-}$	Bruenn:1985en ; Mezzacappa:1993gn
$N+N+\nu +\overline{\nu }\to N+N$	Hannestad:1997gc
$\nu +N\to \nu +N$	Burrows:2004vq
$\nu +A\to \nu +A$	Burrows:2004vq

Table 2. Table 2: In this table we present the mean values of density ρ 𝜌 \rho , temperature T 𝑇 T and electron fraction Y e subscript 𝑌 𝑒 Y_{e} in the regions with opacity 0.5 ≤ τ s ≤ 0.85 0.5 subscript 𝜏 𝑠 0.85 0.5\leq\tau_{s}\leq 0.85 using Eq. ( 4 ) for both models at different stages of the postmerger evolution. In the Sly4 case, at t ∼ 20 similar-to 𝑡 20 t\sim 20 ms, the remnant already collapsed to BH. Densities are given as log 10 ⁡ ( ρ ​ [ gcm − 3 ] ) subscript 10 𝜌 delimited-[] superscript gcm 3 \log_{10}(\rho\,[\mathrm{gcm}^{-3}]) and temperatures as log 10 ⁡ ( T ​ [ MeV ] ) subscript 10 𝑇 delimited-[] MeV \log_{10}(T\,[\mathrm{MeV}]) . For each quantity we also give its variance in the distribution.

		DD2			SLy4
	$t\sim 1$ ms	$t\sim 10$ ms	$t\sim 20$ ms	$t\sim 1$ ms	$t\sim 10$ ms	$t\sim 20$ ms
${\log }_{10}\left({\rho }^{{\nu }_e}\right)$	$10.89\pm 0.33$	$10.94\pm 0.23$	$10.94\pm 0.23$	$11.36\pm 0.22$	$10.97\pm 0.22$	$11.33\pm 0.24$
${\log }_{10}\left(T^{{\nu }_e}\right)$	$0.49\pm 0.14$	$0.58\pm 0.09$	$0.58\pm 0.07$	$0.81\pm 0.20$	$0.75\pm 0.06$	$0.75\pm 0.07$
$Y_e^{{\nu }_e}$	$0.14\pm 0.05$	$0.14\pm 0.02$	$0.14\pm 0.01$	$0.20\pm 0.08$	$0.25\pm 0.02$	$0.18\pm 0.01$
${\log }_{10}\left({\rho }^{{\overline{\nu }}_e}\right)$	$11.07\pm 0.28$	$11.07\pm 0.17$	$11.07\pm 0.19$	$11.35\pm 0.20$	$10.97\pm 0.18$	$11.30\pm 0.22$
${\log }_{10}\left(T^{{\overline{\nu }}_e}\right)$	$0.52\pm 0.16$	$0.62\pm 0.07$	$0.61\pm 0.06$	$0.82\pm 0.17$	$0.75\pm 0.06$	$0.74\pm 0.06$
$Y_e^{{\overline{\nu }}_e}$	$0.13\pm 0.05$	$0.15\pm 0.02$	$0.14\pm 0.01$	$0.21\pm 0.08$	$0.25\pm 0.02$	$0.18\pm 0.01$
${\log }_{10}\left({\rho }^{{\nu }_x}\right)$	$10.76\pm 0.25$	$10.80\pm 0.18$	$10.80\pm 0.19$	$11.20\pm 0.19$	$10.82\pm 0.18$	$11.11\pm 0.21$
${\log }_{10}\left(T^{{\nu }_x}\right)$	$0.48\pm 0.14$	$0.55\pm 0.09$	$0.55\pm 0.06$	$0.82\pm 0.17$	$0.71\pm 0.05$	$0.69\pm 0.05$
$Y_e^{{\nu }_x}$	$0.14\pm 0.05$	$0.14\pm 0.02$	$0.13\pm 0.01$	$0.24\pm 0.09$	$0.24\pm 0.02$	$0.18\pm 0.02$

Table 3. Table 3: Same as Table  2 but for the optical depth τ e ​ q subscript 𝜏 𝑒 𝑞 \tau_{eq} , extracted with Eq. ( 5 ).

		DD2			SLy4
	$t\sim 1$ ms	$t\sim 10$ ms	$t\sim 20$ ms	$t\sim 1$ ms	$t\sim 10$ ms	$t\sim 20$ ms
${\log }_{10}\left({\rho }^{{\nu }_e}\right)$	$11.02\pm 0.39$	$11.02\pm 0.25$	$11.04\pm 0.27$	$11.45\pm 0.23$	$11.04\pm 0.23$	$11.44\pm 0.25$
${\log }_{10}\left(T^{{\nu }_e}\right)$	$0.48\pm 0.17$	$0.60\pm 0.10$	$0.60\pm 0.08$	$0.79\pm 0.22$	$0.76\pm 0.07$	$0.78\pm 0.07$
$Y_e^{{\nu }_e}$	$0.13\pm 0.05$	$0.14\pm 0.02$	$0.14\pm 0.01$	$0.18\pm 0.09$	$0.25\pm 0.02$	$0.18\pm 0.01$
${\log }_{10}\left({\rho }^{{\overline{\nu }}_e}\right)$	$11.79\pm 0.31$	$11.30\pm 0.20$	$11.34\pm 0.22$	$11.63\pm 0.26$	$11.14\pm 0.21$	$11.53\pm 0.21$
${\log }_{10}\left(T^{{\overline{\nu }}_e}\right)$	$0.61\pm 0.17$	$0.70\pm 0.08$	$0.69\pm 0.07$	$0.79\pm 0.24$	$0.79\pm 0.06$	$0.81\pm 0.06$
$Y_e^{{\overline{\nu }}_e}$	$0.11\pm 0.05$	$0.16\pm 0.02$	$0.16\pm 0.01$	$0.16\pm 0.08$	$0.25\pm 0.02$	$0.18\pm 0.01$
${\log }_{10}\left({\rho }^{{\nu }_x}\right)$	$12.97\pm 0.25$	$12.40\pm 0.38$	$12.74\pm 0.22$	$12.54\pm 0.27$	$12.45\pm 0.19$	$-$
${\log }_{10}\left(T^{{\nu }_x}\right)$	$0.96\pm 0.15$	$0.95\pm 0.06$	$0.95\pm 0.05$	$1.01\pm 0.22$	$0.98\pm 0.10$	$-$
$Y_e^{{\nu }_x}$	$0.08\pm 0.03$	$0.10\pm 0.02$	$0.08\pm 0.01$	$0.12\pm 0.07$	$0.11\pm 0.03$	$-$

Table 4. Table 4: In this table we present the mean values of density ρ 𝜌 \rho , temperature T 𝑇 T and electron fraction Y e subscript 𝑌 𝑒 Y_{e} in the regions with opacity 0.5 < τ < 0.85 0.5 𝜏 0.85 0.5<\tau<0.85 using Eq. ( 4 ) for four different resolutions. The considered model is the DD2 at 20 ​ ms 20 ms 20\,\mathrm{ms} after merger. Densities are given as log 10 ⁡ ( ρ ​ [ g / cm 3 ] ) subscript 10 𝜌 delimited-[] g superscript cm 3 \log_{10}(\rho\,[\mathrm{g/cm}^{3}]) and temperatures are expressed in MeV. For each quantity we also give its variance in the distribution.

	LR $\left(1.48\mathrm{km}\right)$	MR $\left(1.11\mathrm{km}\right)$	SR $\left(0.74\mathrm{km}\right)$	HR $\left(0.49\mathrm{km}\right)$
${\log }_{10}\left({\rho }^{{\nu }_e}\right)$	$11.13\pm 0.27$	$11.16\pm 0.28$	$11.18\pm 0.29$	$11.24\pm 0.33$
${\log }_{10}\left(T^{{\nu }_e}\right)$	$0.62\pm 0.07$	$0.63\pm 0.07$	$0.63\pm 0.07$	$0.64\pm 0.08$
$Y_e^{{\nu }_e}$	$0.13\pm 0.01$	$0.13\pm 0.01$	$0.13\pm 0.01$	$0.13\pm 0.01$
${\log }_{10}\left({\rho }^{{\overline{\nu }}_e}\right)$	$11.16\pm 0.18$	$11.17\pm 0.18$	$11.21\pm 0.21$	$11.22\pm 0.21$
${\log }_{10}\left(T^{{\overline{\nu }}_e}\right)$	$0.63\pm 0.05$	$0.63\pm 0.06$	$0.64\pm 0.05$	$0.64\pm 0.06$
$Y_e^{{\overline{\nu }}_e}$	$0.14\pm 0.01$	$0.14\pm 0.01$	$0.15\pm 0.01$	$0.14\pm 0.01$
${\log }_{10}\left({\rho }^{{\nu }_x}\right)$	$10.86\pm 0.17$	$10.90\pm 0.19$	$10.92\pm 0.19$	$10.94\pm 0.20$
${\log }_{10}\left(T^{{\nu }_x}\right)$	$0.56\pm 0.05$	$0.57\pm 0.05$	$0.57\pm 0.05$	$0.58\pm 0.05$
$Y_e^{{\nu }_x}$	$0.13\pm 0.01$	$0.13\pm 0.01$	$0.13\pm 0.01$	$0.13\pm 0.01$

Table 5. Table 5: In this table we present the mean values of density ρ 𝜌 \rho , temperature T 𝑇 T and electron fraction Y e subscript 𝑌 𝑒 Y_{e} in the regions with opacity 0.5 ≤ τ s ≤ 0.85 0.5 subscript 𝜏 𝑠 0.85 0.5\leq\tau_{s}\leq 0.85 using Eq. ( 4 ) for both models at different stages of the postmerger evolution. In the Sly4 case, at t ∼ 20 similar-to 𝑡 20 t\sim 20 ms, the remnant already collapsed to BH. Densities are given as log 10 ⁡ ( ρ ​ [ g ​ cm − 3 ] ) subscript 10 𝜌 delimited-[] g superscript cm 3 \log_{10}(\rho\,[\mathrm{g~{}cm^{-3}}]) and temperatures as T ​ [ MeV ] 𝑇 delimited-[] MeV T\,[\mathrm{MeV}] . For each quantity we also give its variance in the distribution.

	DD2		SLy4
	$t\sim 10$ ms	$t\sim 20$ ms	$t\sim 10$ ms	$t\sim 20$ ms
${\log }_{10}\left({\rho }^{{\nu }_e}\right)$	$11.17\pm 0.27$	$11.21\pm 0.26$	$11.04\pm 0.25$	$11.40\pm 0.25$
${\log }_{10}\left(T^{{\nu }_e}\right)$	$0.67\pm 0.21$	$0.66\pm 0.18$	$0.77\pm 0.16$	$0.77\pm 0.15$
$Y_e^{{\nu }_e}$	$0.15\pm 0.02$	$0.15\pm 0.01$	$0.24\pm 0.03$	$0.18\pm 0.01$
${\log }_{10}\left({\rho }^{{\overline{\nu }}_e}\right)$	$10.91\pm 0.24$	$10.93\pm 0.21$	$11.01\pm 0.21$	$11.35\pm 0.22$
${\log }_{10}\left(T^{{\overline{\nu }}_e}\right)$	$0.59\pm 0.17$	$0.58\pm 0.14$	$0.76\pm 0.24$	$0.76\pm 0.13$
$Y_e^{{\overline{\nu }}_e}$	$0.14\pm 0.02$	$0.14\pm 0.02$	$0.25\pm 0.02$	$0.18\pm 0.01$
${\log }_{10}\left({\rho }^{{\nu }_x}\right)$	$10.46\pm 0.32$	$10.46\pm 0.32$	$10.86\pm 0.20$	$11.15\pm 0.25$
${\log }_{10}\left(T^{{\nu }_x}\right)$	$0.50\pm 0.19$	$0.50\pm 0.12$	$0.72\pm 0.13$	$0.70\pm 0.13$
$Y_e^{{\nu }_x}$	$0.14\pm 0.03$	$0.14\pm 0.01$	$0.24\pm 0.02$	$0.18\pm 0.02$

Equations

$${\tau }_{\gamma :A\to B}={\int }_{\gamma :A\to B}\kappa (s)\mathrm{d}s$$

$$\tau (\mathbf{x})=\underset{\{\gamma |\gamma :\mathbf{x}\to {\mathbf{x}}_{\mathrm{b}}\}}{\min }{\int }_{\gamma }\kappa (s)\mathrm{d}s,$$

$$x^{+},x^{-},y^{+},y^{-},z^{+},$$

$$x^{+}{y}^{+},x^{+}{y}^{-},x^{-}{y}^{+},x^{-}{y}^{-},x^{+}{z}^{+},x^{-}{z}^{+},y^{+}{z}^{+},y^{-}{z}^{+},$$

$$x^{+}{y}^{+}{z}^{+},x^{+}{y}^{-}{z}^{+},x^{-}{y}^{+}{z}^{+},x^{-}{y}^{-}{z}^{+}.$$

$$l_C=5r_{11},$$

$${\kappa }_{\mathrm{diff}}=\sum _r{\kappa }_{\mathrm{ab},r}+\sum _s{\kappa }_{\mathrm{sc},s},$$

$${\kappa }_{\mathrm{eq}}=\sqrt{\left(\sum _r{\kappa }_{\mathrm{ab},r}\right){\kappa }_{\mathrm{diff}}}.$$

$${\sigma }_q=\frac{1}{N}\sum _{i=1}^N\sqrt{{\left(q_i-q_{\mathrm{mean}}\right)}^2},$$

$$f_{\nu }(E_{\nu })={\left(e^{(E_{\nu }-{\mu }_{\nu })/k_{b}T}+1\right)}^{-1},$$

$${\mu }_{{\nu }_e}$$

$${\mu }_{{\bar{\nu }}_e}$$

$${\mu }_{{\nu }_x}$$

$${\tilde{\kappa }}_{\nu }=\frac{{\int }_0^{\infty }{f}_{\nu }(E){\kappa }_{\nu }(E)E^{2}dE}{{\int }_0^{\infty }{f}_{\nu }(E)E^{2}dE}.$$

Full text

††offprints: andrea.endrizzi@uni-jena.de

1 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 Rosswog:2015nja ; Baiotti:2016qnr ; Barack:2018yly 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.Hotokezaka:2011dh ; Bauswein:2013jpa ; Hotokezaka:2013iia ; Zappa:2017xba ; Radice:2018pdn ; Dietrich:2018phi ; Koppel:2019pys . Due to its high neutron content, the ejecta expelled in the interstellar medium can immediately synthetize heavy elements via the so-called rapid neutron capture process ($r$-process) nucleosynthesis Lattimer:1974a ; Symbalisty:1982a ; Thielemann:2017acv . 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 Li:1998bw ; Kulkarni:2005jw ; Metzger:2016pju . 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 Paczynski:1986px ; Eichler:1989ve .

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 TheLIGOScientific:2017qsa . 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, Mooley:2018dlz ; Ghirlanda.etal:2019 ) and of a kilonova (AT2017gfo) produced by the merger Monitor:2017mdv ; GBM:2017lvd . 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. Kasen:2017sxr ; Villar:2017wcc ; Perego:2017wtu ; Wollaeger:2017ahm ; Kawaguchi:2018ptg ; Miller:2019dpt . 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 $\sim 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, Perego:2019adq ) 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. Eichler:1989ve ; Ruffert:1996by ; Rosswog:2003rv ; Perego:2017fho . Deep inside the forming hot remnant, for densities in excess of $10^{12}{\rm 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 ($\sim 10^{11}{\rm 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. Dessart:2008zd ; Perego:2014fma ; Just:2014fka ; Foucart:2015gaa ; Sekiguchi:2016bjd ; Martin:2017dhc .

Neutrino cross sections depend significantly on the incoming neutrino energy: for energies in the range 1-100 MeV the dependence is approximately quadratic Dicus:1972yr ; Bruenn:1985en ; Formaggio:2013kya . 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. * Martin:2015hxa ; Metzger:2014ila ; Metzger:2017wot .

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. Galeazzi:2013mia ; Palenzuela:2015dqa ; Foucart:2015gaa ; Foucart:2016rxm ; Perego:2015agy ; Radice:2016dwd ; Bovard:2017mvn ; Perego:2017xth ; Foucart:2018gis ; Radice:2018pdn ; Ardevol-Pulpillo:2018btx . 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 (CCSNe) and proto-neutron stars (PNSs) many sensitivities studies on the relevant neutrino physics have been performed, e.g. Mezzacappa:1993gn ; Roberts.etal:2012 ; Lentz:2012xc ; Abdikamalov:2014oba ; Melson:2015 ; Fischer:2016 ; OConnor:2018sti ; Pan:2018vkx ; Cabezon:2018lpr ; Just.etal:2018 , 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 Ruffert:1996by ; Rosswog:2003rv . 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 Perego:2014fma , 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. Perego:2014qda , 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 Galeazzi:2013mia ; Foucart:2015vpa ; Perego:2015agy ; Foucart:2016rxm .

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.

2 Method

2.1 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 $\gamma$ denotes the path between $A$ and $B$, the optical depth is defined as:

[TABLE]

where $\kappa(s)$ is the opacity (corresponding to the inverse of the local mean free path) and ${\rm d}s=\sqrt{g_{ij}{\rm d}x^{i}{\rm d}x^{j}}$ is an infinitesimal displacement along the chosen path, with $g_{ij}$ being the local spatial metric. The physical interpretation of $\tau$ emerges from its definition: $\tau$ counts the average number of interactions that radiation particles experiences along $\gamma$. 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 $\mathbf{x}$ as

[TABLE]

where $\mathbf{x}_{\rm 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 $\tau\sim 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 $\tau=2/3$.

A common approach to effectively evaluate Eq. (2) is to select a certain number of direction moving away from $\mathbf{x}$, evaluate $\tau$ 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., Rosswog:2003rv Perego:2014fma ). In our scheme, instead, we chose $17$ directions, corresponding to $5$ on-axis directions

[TABLE]

(we disregard all $z^{-}$ directions, since all simulations considered in this work are performed with $z-$symmetry), $8$ planar diagonals

[TABLE]

and $4$ full diagonals

[TABLE]

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

[TABLE]

with $r_{11}$ being the average radius of a density isosurface at $\rho_{\rm th}=2\times 10^{11}\,\mathrm{g~{}cm^{-3}}$. In the simulations presented in this work $r_{11}\sim 30\mathrm{km}$. To test our approximation, we have computed $\tau$ 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.

2.2 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 ($\nu_{e}$), electron anti-neutrinos ($\bar{\nu}_{e}$) and a collective species for $\mu$ and $\tau$ (anti)neutrinos ($\nu_{x}$). The most relevant neutrino-matter reactions considered in this work are listed in Table 1. We consider the corresponding absorptivity $\kappa_{{\rm ab}}$ or scattering opacity $\kappa_{{\rm 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 $\kappa_{\rm diff}$ as

[TABLE]

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. Shapiro:1983du or Raffelt:2001kv for analogous expressions):

[TABLE]

These two kinds of opacities provide two different optical depths, $\tau_{\rm diff}$ and $\tau_{\rm 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, $\kappa_{\rm diff}\geq\kappa_{\rm 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 OConnor:Thesis:2012 . In particular, we rely on expressions for the transport opacities and we include weak magnetism corrections Horowitz:2001xf , ion-ion correlations, form factor correction Tubbs:1975jx ; Burrows:1981zz , and electron polarization correction Leinson:1997ns according to Burrows:2004vq . Differently from the standard code version, we do not include the effect of stimulated absorption in our calculations because it provides unphysically high $\bar{\nu}_{e}$ opacities in cold, low dense matter, well below $\rho_{\rm th}$. For the inverse nucleon-nucleon bremsstrahlung and the neutrino-antineutrino pair annihilation we compute the reaction kernels according to Hannestad:1997gc and Bruenn:1985en ; Mezzacappa:1993gn , 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 $\nu_{e}+\bar{\nu}_{e}\rightarrow e^{+}+e^{-}$, we fix the $\nu_{e}$ energy and we integrate over an equilibrium Fermi-Dirac distribution function for $\bar{\nu}_{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 $\nu_{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_{\rm NS}=1.364~{}M_{\odot}$ (corresponding to a chirp mass of $1.188~{}M_{\odot}$ compatible with GW170817) and we make use of two different finite temperature, composition dependent EOSs: DD2 Typel:2009sy ; Hempel:2009mc and SLy4 daSilvaSchneider:2017jpg . 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 Gourgoulhon:2000nn . The initial separation between the NS is set to $45\,{\rm km}$, corresponding to $\sim 2{-}3$ orbits before merger. Low temperature ($T\lesssim 0.1$ MeV), $\nu$-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 Radice:2012cu ; Radice:2013hxh ; Radice:2013xpa ; Radice:2015nva , complemented by a leakage scheme to account for compositional and energy changes in the matter due to weak reactions involving $\nu_{e}$, $\bar{\nu}_{e}$, and $\nu_{x}$. Free-streaming neutrinos are emitted at an average energy and then evolved according to the M0 scheme introduced in Radice:2016dwd ; Radice:2018pdn . All the technical details for these simulations are given in Radice:2018xqa ; Perego:2019adq 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 Perego:2019adq . The DD2 simulation is the same discussed in Perego:2019adq , 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\simeq 185\,{\rm m}$. Since our optical depth calculation scheme works on uniform grids, we interpolate the simulation results onto one uniform grid of resolution $h^{\prime}\sim 0.74\mathrm{km}$ 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. 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. Perego:2017xth , 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 $\sim 12\,\mathrm{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\,\mathrm{ms}$ after merger time (identified as the time corresponding to the maximum in the strain of the $\ell=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.

2.3 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% Foucart:2015vpa ; Perego:2015agy ; Foucart:2016rxm

As a first approach, we fix the neutrino energy to the following set of values compatible with those results: $\langle E_{\nu_{e}}\rangle\approx 9.34~{}{\rm MeV}$, $\langle E_{\bar{\nu}_{e}}\rangle\approx 15.16~{}{\rm MeV}$, and $\langle E_{\nu_{x}}\rangle\approx 23.98~{}{\rm 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\leq\tau\leq 0.85$. For each extracted quantity $q$ we compute the volumetric mean $q_{\rm mean}$ and to give an estimate of its distribution around the mean value we compute the corresponding standard deviation as:

[TABLE]

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}\left({\rho}~{}[{\rm g~{}cm^{-3}}]\right)$, $\log_{10}{\left(T~{}[{\rm MeV}]\right)}$ 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  Foucart:2015vpa ; Foucart:2016rxm . 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:

[TABLE]

where $E_{\nu}$ is the neutrino energy, $\mu_{\nu}$ 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

[TABLE]

where $\nu_{n}$, $\nu_{p}$ and $\nu_{e}$ are the relativistic neutron, proton and electron chemical potentials, respectively. We introduce a spectrum-averaged opacity defined as:

[TABLE]

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, $\tau_{\rm eq}$. The average used in this work, Eq. (9), is slightly different from the average computed in some gray neutrino transport schemes, e.g. Foucart:2015vpa , where the factor $E^{2}$ is replaced by $E^{3}$. In the latter approach, more emphasis is put on the energy transport ($I_{\nu}\propto 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.

3 Results

3.1 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_{\rm abs}$ and $k_{\rm 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 ($\rho_{\rm th}=2\times 10^{11}~{}{\rm g~{}cm^{-3}}$, $T_{\rm th}=4~{}{\rm MeV}$ and $Y_{e,{\rm th}}=0.1$, see for example Perego:2014fma ). 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 $\langle E_{\nu}\rangle$ 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 $\kappa_{\rm scat}\propto n_{\rm b}E_{\nu}^{2}$ (with $n_{\rm b}$ being the free baryon particle density), $k_{\rm scat}$ is approximately proportional to $\rho$, while the larger (smaller) scattering opacity observed for $\nu_{x}$ ($\nu_{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_{\rm 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_{\rm scat}$ at low temperatures ($T\lesssim 1~{}{\rm 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. Bruenn:1985en ; Burrows:2004vq .

In the case of $k_{\rm 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 $\sim 10^{12}{\rm g~{}cm^{-3}}$ the neutron richness and the larger $Q$-value favor $\nu_{e}$ captures on neutrons with respect to $\bar{\nu}_{e}$ captures on protons111For $\rho\lesssim 10^{10}{\rm g~{}cm^{-3}}$, the dominant role of neutrino pair annihilation is a spurious effect due to the assumption of having a neutrino gas even at low densities. We notice that, even for a moderately high temperature, the values for the opacities are relatively low and thus they do not impact on the determination of the decoupling conditions.. 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}{\rm 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 $\bar{\nu}_{e}$ are suppressed by degeneracy at high densities and this translates in a much smaller inverse bremsstrahlung opacity for $\nu_{e}$ compared to $\nu_{x}$ and $\bar{\nu}_{e}$. Variations in temperature have also a significant effect ok $k_{\rm abs}$, as visible 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 $\sim 10$ MeV pair processes of thermal origin (like $\nu$-$\bar{\nu}$ annihilation) dominate. In particular, the corresponding inverse mean free path is roughly proportional to the $Y_{\nu}Y_{\bar{\nu}\sigma_{\nu+\bar{\nu}}}$, which for thermally produced neutrinos translate in a strong dependence on the matter temperature, $\propto T^{8}$.

Changes in the electron fraction also introduce significant dependencies on the $k_{\rm 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.

3.2 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 $\langle E_{\nu}\rangle$.

3.2.1 Diffusion surfaces

In Fig. 3 and Fig. 4 we present the contours of the diffusion surfaces ($\tau_{\rm 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~{}{\rm ms}$ after merger the system has already reached an almost axisymmetric configuration, in the first snapshot ($\sim 1~{}{\rm 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 $\nu_{e}$ mean energy, the more abundant free neutrons provide larger opacities to $\nu_{e}$ than to $\bar{\nu}_{e}$, i.e. $(1-Y_{e})\langle E_{\nu_{e}}\rangle^{2}\gtrsim Y_{e}\langle E_{\bar{\nu}_{e}}\rangle^{2}$. Thus, the average $\bar{\nu}_{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 $\nu_{e}$ and $\bar{\nu}_{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 AH222We note that the region within the AH is not considered in our evaluation of hydrodynamics properties in the transient region.. 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\,\mathrm{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}{\left(\rho\right)}$ 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 $\sim 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 Perego:2019adq and 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 $\sim 10\,\mathrm{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\,\mathrm{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 ($\gtrsim 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 $(\rho,T)$ conditions normalized to the total volume of the neutrinosphere region $V_{\rm 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\,\mathrm{ms}$, while maintaining a noticeable spread in temperature in the high density regions. This feature persists for the DD2 model even at $20\,\mathrm{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 $\rho$ both for $k_{\rm abs}$ and $k_{\rm 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 Perego:2019adq .

3.2.2 Equilibrium surfaces

The surfaces where the freeze-out from weak and thermal equilibrium occurs for mean energy neutrinos (i.e. $\tau_{\rm eq}(\langle E_{\nu}\rangle)=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 different neutrino species. In the case of $\nu_{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 $\nu_{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 absorb $\bar{\nu}_{e}$ locates their equilibrium decoupling surfaces at slightly larger densities and temperatures than the ones of $\nu_{e}$’s. A qualitatively different behavior characterizes $\nu_{x}$’s. In this case, thermal and weak equilibrium is guaranteed by pair processes, e.g. $\nu+\bar{\nu}\rightarrow e^{+}+e^{-}$. Due to their strong temperature dependence, the corresponding equilibrium surfaces shrink deeply inside the remnant, at densities $\sim 10^{13}{\rm g~{}cm^{-3}}$ where temperatures are usually in excess of $8~{}{\rm MeV}$. We also notice that after BH formation (i.e. in the SLy4 model at $20~{}{\rm ms}$) there is no region where $\tau_{\rm eq,\nu_{x}}>2/3$, which means that even for $E_{\nu_{x}}\geqslant 20\,\mathrm{MeV}$ the combination of scattering and pair processes that enters Eq. (5) requires matter with densities around $\sim 10^{13}~{}\mathrm{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}{\rm 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 $(\rho,T)$ plane and for the DD2 model. In the case of $\nu_{e}$’s the distributions are practically identical to the one extracted for the scattering surfaces. For $\bar{\nu}_{e}$’s, the correlation between $\rho$ 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 $\nu_{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, $\langle\epsilon_{\nu}\rangle\sim F_{2}(0)\,T$ (where $F_{2}(0)\approx 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.

3.3 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.

3.3.1 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 $\sim n_{\rm target}E_{\nu}^{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, Perego:2019adq ), the decoupling temperature also decreases monotonically with $E_{\nu}$. 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 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_{\nu}$ 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 $\nu_{e}$’s than for $\bar{\nu}_{e}$’s and $\nu_{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 $\bar{\nu}_{e}$ or $\nu_{x}$ with energies of $5\,\mathrm{MeV}$ or lower. For $\nu_{e}$ we are able to find a small region that satisfies the condition $0.5<\tau_{\rm diff}<0.85$, but $\sim 90\%$ of it comes from matter at $\tau<2/3$ located in a rather small volume.

3.3.2 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_{\rm eq}\leq k_{\rm 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 $\nu_{e}$ and $\bar{\nu}_{e}$ behavior can be observed in the case of $\nu_{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 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 $\nu_{e}$ and $\bar{\nu}_{e}$. However, due to the low neutrino degeneracy, to the more abundant neutrons and to the larger abundance of electron antineutrinos for $\rho\gtrsim 10^{13}{\rm g~{}cm^{-3}}$, $T\gtrsim 10~{}{\rm MeV}$ and $Y_{e}\sim 0.07$, electron neutrinos decouple again at slightly lower rest mass densities and temperature than $\bar{\nu}_{e}$ and $\nu_{x}$. In the BH remnant case, the lack of equilibrium decoupling surfaces for low energy neutrinos is more severe. This holds especially true for $\nu_{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.

3.4 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 $\bar{\tau}$ 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, $\tau(\langle E_{\nu}\rangle)$.

The energy integrated approach assumes the neutrino radiation to be in equilibrium with matter. Thus we expect it to be more coherent with the $\tau_{\rm eq}$ case. Indeed, in Fig. 12 we show $\tau_{\rm eq}(\langle E_{\nu}\rangle)$ alongside $\bar{\tau}_{\rm 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 $\tau$ inside the neutrino surfaces can differ significantly among the two cases (and in particular for the $\nu_{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 $\nu_{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 $\bar{\nu}_{e}$ case and, even more, in the $\nu_{x}$ case the assumption of matter-radiation equilibrium even for $\bar{\tau}_{\rm diff}$ decreases the mean energies and lowers the opacity. Thus, the neutrino surfaces are apparently located deeper inside in comparison with the expected locations.

4 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_{\odot}$) 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 $\sim 10~{}\mathrm{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\,\mathrm{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 $\nu_{e}$ and $\bar{\nu}_{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 guaranteed by pair processes, like inverse nucleon-nucleon bremsstrahlung and $\nu$-$\bar{\nu}$ 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_{\rm scat}\propto E_{\nu}^{2}$) and the larger abundance of absorbing neutrons relevant for $\nu_{e}$. Typical rest mass densities at the corresponding diffusion surfaces are $\sim 10^{11}{\rm g~{}cm^{-3}}$, with slightly smaller values obtained for $\nu_{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\,\mathrm{MeV}$, while in the latter between $3$ and $4\,\mathrm{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 $\nu_{e}$’s. On the other hand, the relative paucity of absorbing protons moves the equilibrium surfaces of $\bar{\nu}_{e}$’s deeper inside the remnant, at rest mass densities of a few times $10^{11}{\rm g~{}cm^{-3}}$. Consequently, the $\bar{\nu}_{e}$ equilibrium decoupling temperature increases to $5\,\mathrm{MeV}$ in the DD2 case and to $6\,\mathrm{MeV}$ in the SLy4 case. For heavy flavor neutrinos, matter temperature becomes the most relevant quantity and the equilibrium decoupling happens around $10~{}{\rm MeV}$. Inside the remnant profile this corresponds to a rest mass density in excess of $10^{12}{\rm g~{}cm^{-3}}$. Our results for the equilibrium decoupling temperatures broadly satisfy the expected blackbody emission relation, i.e. $\langle E_{\nu}\rangle\sim 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 $\nu_{e}$ and $\bar{\nu}_{e}$ equilibrium decoupling temperatures reduces, as also observed in Sekiguchi:2016bjd .

We have generalized our analysis by considering neutrinos with a broad range of energies, between $3$ and $88\,\mathrm{MeV}$, covering the most relevant part of the emitted spectrum, but irrespective of their relevance. Due to the strong dependence on $E_{\nu}$ 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}{\rm g~{}cm^{-3}}$ for hard neutrinos to $10^{13}{\rm g~{}cm^{-3}}$ for soft ones, with a systematic increase going from $\nu_{e}$ to $\bar{\nu}_{e}$, and finally to $\nu_{x}$. The corresponding decoupling temperatures cover one order of magnitude, moving from $\sim 2$ to $\sim 10\,\mathrm{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 qualitatively different behavior is observed for the equilibrium decoupling of heavy flavor neutrinos, for which temperatures in the range $5$-$12\,\mathrm{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 $\bar{\nu}_{e}$ ($E_{\nu}\lesssim 10\,\mathrm{MeV}$) and even more for ${\nu}_{x}$ ($E_{\nu}\lesssim 40\,\mathrm{MeV}$) than for ${\nu}_{e}$ ($E_{\nu}\lesssim 5\,\mathrm{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.

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 emply 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 modelling 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 time scales much larger than the dynamical (and even the explosion) time scale. Due to the high sensitivity of CCSN explosions on weak reactions, neutrino-matter interactions in CCSNe have been deeply investigated in many studies, e.g. Roberts.etal:2012 ; Lentz:2012xc ; Abdikamalov:2014oba ; Melson:2015 ; Fischer:2016 ; OConnor:2018sti . 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}\gtrsim 0.3$ and the maximum core temperature (and thus the neutrino mean energies) stay low, $\lesssim$ a few MeV. Under these conditions, decoupling occurs at higher densities (e.g., several times $10^{11}{\rm g~{}cm^{-3}}$ for mean energy $\nu_{e}$’s), 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 ($\gtrsim 25M_{\odot}$), 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 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 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. Bacca:2011qd ; Rrapaj:2014yba ; Roberts:2016mwj ; Guo:2019cvs . This represents a great challenge in nuclear and neutrino physics for the years to come.

Acknowledgments

The authors thank Bernd Brügmann and Domenico Logoteta for useful discussions. The authors thank the organizers and participants of the INT Program INT-18-72R “First Multi-Messenger Observation of a Neutron Star Merger and its Implications for Nuclear Physics INT workshop” held at (Seattle, March 2018), of the ExtreMe Matter Institute’s rapid task force meeting at GSI/FAIR (Darmstadt, June 2018),of the GWEOS workshop (Pisa, February 2019), for stimulating discussions. SB acknowledges support by the EU H2020 under ERC Starting Grant, no. BinGraSp-714626. DR acknowledges support from a Frank and Peggy Taplin Membership at the Institute for Advanced Study and the Max-Planck/Princeton Center (MPPC) for Plasma Physics (NSF PHY-1804048). Computations were performed on the supercomputer SuperMUC at the LRZ Munich (Gauss project pn56zo), on supercomputer Marconi at CINECA (ISCRA-B project number HP10B2PL6K and HP10BMHFQQ); on the supercomputers Bridges, Comet, and Stampede (NSF XSEDE allocation TG-PHY160025); on NSF/NCSA Blue Waters (NSF AWD-1811236). AE, AP, and BG acknowledge computational support also by the INFN initiative TEONGRAV. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

Appendix 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^{\prime}\sim 0.74\,\mathrm{km}$); we then picked as our best resolution (HR) a grid with $h^{\prime}\sim 0.49\,\mathrm{km}$ and as our worst resolution (LR) a grid with twice the spacing of our SR ($h^{\prime}\sim 1.48\,\mathrm{km}$). In order to make a more informed convergence study, we also considered a medium resolution (MR) grid, with $h^{\prime}\sim 1.11\,\mathrm{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 $\rho$ now is split into $8$ cells, increasing the possibility that at least one of them has a density $\rho_{1}<\rho$ that is just enough to push the value of $\tau$ for that cell below the threshold. This only works towards decreasing $\tau$, since we always pick the minimum optical depth between all directions considered, and therefore is sufficient that this happens in one direction to decrease $\tau$, while for this mechanism to increase $\tau$ the new grid should have all the neighbouring cells to end up with higher $\rho$. This is unlikely given the distribution of density in the system.

Appendix 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 Perego:2014qda . 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.