Neutrino emission from binary neutron star mergers: characterising light curves and mean energies

Abstract

Neutrinos are copiously emitted by neutron star mergers, due to the high temperatures reached by dense matter during the merger and its aftermath. Neutrinos influence the merger dynamics and shape the properties of the ejecta, including the resulting r-process nucleosynthesis and kilonova emission. In this work, we analyse neutrino emission from a large sample of binary neutron star merger simulations in Numerical Relativity, covering a broad range of initial masses, nuclear equation of state and viscosity treatments. We extract neutrino luminosities and mean energies, and compute quantities of interest such as the peak values, peak broadnesses, time averages and decrease time scales. We provide a systematic description of such quantities, including their dependence on the initial parameters of the system. We find that for equal-mass systems the total neutrino luminosity (several $10^{53}{\hbox {erg}~{\hbox {s}}^{-1}}$) decreases as the reduced tidal deformability increases, as a consequence of the less violent merger dynamics. Similarly, tidal disruption in asymmetric mergers leads to systematically smaller luminosities. Peak luminosities can be twice as large as the average ones. Electron antineutrino luminosities dominate (initially by a factor of 2-3) over electron neutrino ones, while electron neutrinos and heavy flavour neutrinos have similar luminosities. Mean energies are nearly constant in time and independent on the binary parameters. Their values reflect the different decoupling temperature inside the merger remnant. Despite present uncertainties in neutrino modelling, our results provide a broad and physically grounded characterisation of neutrino emission, and they can serve as a reference point to develop more sophisticated neutrino transport schemes.

1 Introduction

BNS mergers represent one of the main research topics in modern astrophysics. Due to the wide range of densities and temperatures required to study the dynamics of these events [1], their understanding connects several branches of physics spanning from nuclear physics to relativistic hydrodynamics in strong-field conditions. They can be considered natural laboratories to investigate the behaviour of matter at extreme densities, which cannot be produced in Earth-based facilities (see e.g. [2,3,4] for recent reviews).

BNS mergers are prominent sources of gravitational waves (GWs) [5, 6], and a primary target for ground-based GW detector facilities such as LIGO [7], VIRGO[8] and KAGRA [9]. Furthermore, they have long been considered one of the most likely progenitors of high-energy astronomical signals such as short gamma ray bursts (sGRBs) [10,11,12,13] and kilonovae [14, 15], see also [16, 17] for recent reviews. Kilonovae (sometimes also referred to as macronovae) are powered by the decay of radioactive heavy elements that are synthesised in the ejecta of binary neutron star (BNS) mergers [see e.g. [18, 19], and references therein]. This aspect links these systems to open issues regarding the evolution of the chemical composition of the Galaxy and of the Cosmos. BNS mergers have indeed emerged as sites (perhaps the main ones) of production of heavy elements in the Universe [20,21,22,23,24].

All these hypotheses recently received a direct confirmation by the first multimessenger detection of a BNS merger. This event (hereafter referred to as GW170817) was observed as a GW signal [25], followed by a sGRB (GRB170817A) and, finally, by a kilonova lasting from a few hours to several days after the merger [21, 24, 26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47]. This detection opened the era of multimessenger astronomy from compact binary mergers. A second detection of a GW signal from a BNS merger, GW190425 [48], was observed a couple of years later, but without the firm identification of associated electromagnetic counterparts.

BNS mergers produce copious amounts of neutrinos, starting from the latest moments of the inspiral until the merger remnant collapses or cools down. This emission is a key element in the dynamics of the system. On one hand, neutrinos are thought to play a significant role in the jet-launching mechanism that powers sGRBs [e.g. [10, 49,50,51]]. Neutrino absorption and energy deposition in the funnel above the poles of the merger remnant could contribute to clean this region, reducing its baryon density and allowing the launch of a relativistic jet [52, 53]. It has also been suggested that neutrino/antineutrino pair annihilation could deposit an amount of energy compatible with the one necessary to explain sGRBs [see e.g. [10, 54,55,56,57,58,59,60]]. Neutrino absorption is also likely one of the mechanisms for matter ejection from BNS mergers, in association to the production of neutrino-driven winds on time scales of $\sim 100 {\hbox { ms}}$ after the merger [e.g. [50, 57, 60, 61]] . Even more importantly, neutrino-matter interactions affect the composition of the ejecta, by driving the evolution of the relative abundance of neutrons and protons, starting from the decompression of beta-equilibrated, cold neutron star (NS) matter. The neutron richness in the ejecta directly impacts the outcome of the r-process nucleosynthesis and of the resulting kilonova signal [62,63,64]. It was shown that the neutrino transport used in the simulations influences essential ejecta properties like the radial speed, the electron fraction and the entropy [65,66,67]. To reliably model these phenomena it is therefore of the utmost importance to characterise the properties of neutrino emission in BNS mergers.

BNS mergers are intrinsically multi-dimensional events. Moreover, their thermodynamic conditions are such that the neutrino optical depth decreases by several orders of magnitude from the optically thick central remnant to the optically thin accretion disc [68]. The quantitative modelling of neutrino production and diffusion in BNS mergers is, thus, a non-trivial task that has only been made possible by the advent of sophisticated numerical simulations in three spatial dimensions. The employed transport methods range from light bulb models in Newtonian spacetime, to moment schemes, and even to Monte Carlo schemes in full general relativity (GR) [e.g. [54, 69,70,71,72,73,74,75,76,77,78,79]]. Our understanding of neutrino physics and transport in BNS mergers largely benefits from Core Collapse-Supernova (CCSN) modelling [80, 81], and references therein]. However, compared to the wealth of literature regarding neutrinos in CCSNe, only few studies in the past have examined neutrino luminosities and mean energies in BNS mergers [54, 65, 69, 79, 82,83,84,85,86,87]. From these seminal studies, a few robust features emerged. Due to the initial neutron richness, electron antineutrinos dominate over the other flavours. Moreover, heavy flavour neutrinos are more energetic, since they decouple deeper inside the remnant. Additionally, more compact BNSs produce more violent mergers, resulting in larger neutrino luminosities. Despite the general consensus about these features, quantitative differences have emerged, such that both the absolute and the relative importance of the different neutrino species, as well as their temporal evolution during the transition between the merger and the remnant cooling phase still remain largely unexplored. One of the main reasons behind these limitations is that neutrino luminosities are only studied for a few milliseconds, while neutrino cooling is relevant during the entire cooling phase, lasting up to tens of seconds.

Table 1 Weak reactions accounted for in the neutrino transport schemes. The following notation is used: $\text {N}\in \{\text {n}, \text {p}\}$ denotes a free nucleon, $\text {A}$ a nucleus, $\nu \in \{\nu _{\text {e}},{\bar{\nu }}_{\text {e}},\nu _{\text {x}}\}$ a neutrino. The “Reference” column accounts for the corresponding rate implementation

In this work, we consider BNS simulations spanning a wide range in total mass, mass ratio, and dense matter equation of state (EOS). Moreover, we consider some of the longest BNS merger simulations in 3+1 numerical relativity (NR). We also consider the effects of the inclusion of physical viscosity of magnetic origin in our simulations. Based on this ample trove of data, we endeavour to find patterns, trends and commonalities in the temporal evolution of the neutrino luminosities and mean energies. We strive to identify in neutrino data universal relations, i.e., relations between parameters describing neutrino emission and quantities characterising BNS models that are EOS independent. Similar relations have been found in the context of NS structure and GW emission [88,89,90,91,92,93,94]. The broad scope of our data sample, which allows us to avoid as much as possible being biased towards a too specialised subset of BNS merger configurations, represents a major innovation of this work.

All the simulations considered in this work, in addition to being homogeneous with respect to the general numerical setup, share the same neutrino physics input and transport scheme. In particular, the minimal set of necessary neutrino reactions has been included (see the main text and Table 1 for details). Moreover, neutrino transport is taken in account using the combination of a leakage scheme and a so-called M0 scheme. These schemes attempt to strike a balance between computational cost and physical realism. In our setup, neutrinos are assumed to be massless and we neglect neutrino oscillations.

The paper is organised as follows: in Sect. 2 we summarise the numerical methods employed to perform the simulations, which we base our analysis on; Sect. 3 describes our simulation sample, the overall properties of neutrino emission, and the analysis strategy that we follow; Sect. 4 contains the main results of this work, in the form of a detailed analysis of the properties of neutrino emission in BNS mergers and their likely explanation in terms of the system dynamics. We discuss our results in the context of multimessenger astrophysics in Sect. 5. We finally summarise our findings and discuss their implications in Sect. 6. Several appendices provide additional details on our analysis, including information about each simulation in our sample.

2 Methods

2.1 Numerical setup

We base our analysis on results collected from a large sample of BNS mergers simulations in NR. All simulations share the same numerical setup and evolution scheme. In the following, we summarise them and we briefly introduce the codes used to produce our data. More details can be found in Ref. [95].

The BNS initial data are evolved with the infrastructure provided by the Einstein Toolkit [96,97,98]. The hyperbolic sector of Einstein’s field equations is evolved with the Z4c formalism [99], implemented in the CTGamma solver [100, 101]. Moreover, general relativistic hydrodynamics is handled by the WhiskyTHC code [95, 102,103,104,105]. The code solves Euler’s equations for the balance of energy and momentum:

$$\begin{aligned} \nabla _\nu T^{\mu \nu } = Q u^{\mu }\,, \end{aligned}$$

(1)

where $T^{\mu \nu }$ is the stress-energy tensor and Q is the net energy deposition rate due to the absorption and emission of neutrinos and antineutrinos (see Sect. 2.3.1). WhiskyTHC evolves neutron and proton number densities separately as:

$$\begin{aligned} \nabla _\mu \left( n_{\text {p,n}} u^\mu \right) = R_{\text {p,n}}, \end{aligned}$$

(2)

where $n_{\text {p,n}}$ are the proton and neutron number densities, respectively, $u^\mu$ is the fluid four-velocity and $R_{\text {p,n}}$ is the net lepton number exchange rate due to the absorption and emission of electron flavour neutrinos and antineutrinos. Due to charge neutrality the electron fraction is directly related to the proton number density, i.e. $Y_{\text {e}} \equiv n_{\text {e}}/(n_p+n_n) = n_{\text {p}}/(n_p+n_n)$. Neutrino emission and cooling are handled with a leakage scheme, while neutrino absorption and heating in optically thin conditions are treated with the so-called M0 scheme (see Sect. 2.3.1). Equations (1) and (2) are closed by a finite-temperature, composition dependent, nuclear EOS (see Sect. 2.3.2). The code also implements the general-relativistic large-eddy simulation (GRLES) method to account for turbulent viscosity of magnetic origin (see Sect. 2.3.3).

The computational domain of the simulations is a cube of side $\sim 3024{\hbox { km}}$ centred on the binary’s centre of mass. The code uses a box-in-box Berger-Oliger adaptive mesh refinement (AMR) scheme with refluxing [106, 107] provided by the Carpet module of the Einstein Toolkit, and composed of seven refinement levels. The finest refinement level covers both NSs during the inspiral and the remnant after the merger, and it has a resolution of $h\approx 246\,{\textrm{m}}$ (for grid setup named here low resolution; LR), $h\approx 185\,{\textrm{m}}$ (standard resolution; SR) or $h\approx 123\,{\textrm{m}}$ (high-resolution; HR) [see also [108]]

2.2 Relevant simulation parameters

Each BNS is characterised by the gravitational masses of the two NSs at infinity, $M_{A,B}$¹. The total gravitational mass and mass ratio are defined as $M_{\text {tot}}=M_A+M_B$ and $q=M_A/M_B$, respectively. A further characterisation system is provided by the dimensionless reduced tidal deformability ${\tilde{\varLambda }}$, since it also depends on the stars’ EOS. It is a weighted average of the dimensionless tidal deformabilities $\varLambda _i$, $i\in {A,B}$, of the two NSs, defined as [109]:

$$\begin{aligned} {\tilde{\varLambda }}= \frac{16}{3}\frac{\left( M_A+12M_B\right) M_A^4\varLambda _A}{M_\text {tot}^5} + (A\leftrightarrow B)\,. \end{aligned}$$

(3)

In Eq. (3) the notation $(A\leftrightarrow B)$ indicates a second term identical to the first except that the indices A and B are exchanged. The dimensionless tidal deformabilities in turn are related to the quadrupolar Love number, $k_2$, describing the static quadrupolar deformation of a star in the gravitoelectric field of the companion [110], by:

$$\begin{aligned} \varLambda _i=\frac{2}{3} k_2 C_i^{-5}\,, \end{aligned}$$

(4)

where $C_i=GM_i/c^2R_i$ is the NS compactness and $R_i$ is the areal radius prior to deformation.

The initial data for all the selected simulations are constructed by solving for irrotational stars of varying masses and different EOSs, using the spectral elliptic solver LORENE [111]. The binaries are set to quasi-circular orbits at an initial separation which, in most cases, is 45 km. This orbital separation corresponds to an inspiral phase of 2–3 orbits before merger. Note that our results do not depend sensitively on the initial separation or the number of orbits before merger, since neutrino emission is linked to the dynamics of the system in the post-merger phase. The EOS used in solving for the initial data are the minimum temperature slice of the EOS table used for the evolution composition fixed assuming neutrinoless beta-equilibrium.

In the following, we use the term model to describe a BNS system with a given combination of initial masses and EOS. For each model, we can have multiple realisations of it, i.e. simulations, which differ from one another by having been run at different resolution, or by including or not a model of the magnetic viscosity of turbulent origin.

2.3 Input physics

2.3.1 Neutrino transport

Since the focus of the present work are the properties of neutrino emission, we provide here a brief, yet fairly detailed, description of the methods of neutrino transport implemented in the simulations that we use. These methods (a leakage scheme and the so-called M0 scheme) are described in detail in Refs. [71, 105] and references therein. They are both “grey” schemes, i.e. schemes in which the dependence of various quantities on the energy of the neutrinos is not explicitly taken into account: instead, energy-averaged quantities are considered. They account for three distinct neutrino species: electron neutrinos, $\nu _{\text {e}}$; electron antineutrinos, ${\bar{\nu }}_{\text {e}}$; and a collective species for heavy neutrinos, $\nu _{\text {x}}$. The last one models muonic and tauonic neutrinos and antineutrinos as a single species of statistical weight 4.

Neutrino emission. The emission of neutrinos from the fluid and the subsequent loss of energy is described by a neutrino leakage scheme (NLS). It is based on the method outlined in Ref. [112], where the local thermodynamical equilibrium chemical potential is used everywhere for all neutrino species while computing opacities as in Ref. [55]. Table 1 lists the reactions taken into account by this scheme to compute the neutrino production free rates, $R_\nu ^{\text {free}}$, $\nu \in \{\nu _{\text {e}},{\bar{\nu }}_{\text {e}},\nu _{\text {x}}\}$, the free energy release, $Q_\nu ^{\text {free}}$, and the neutrino absorption, $\kappa _{\nu ,\text {a}}$, and scattering, $\kappa _{\nu ,\text {s}}$, opacities. These reactions include charged current absorption reactions on free nucleons, namely electron neutrino and antineutrino absorption on free neutrons and protons, respectively; and their inverse reactions. The direct ones are the main responsible for the absorption of $\nu _{\text {e}}$ and ${\bar{\nu }}_{\text {e}}$ both in optically thick and thin conditions, and they provide a relevant contribution to neutrino opacity. The inverse ones are the main processes responsible for the production of electron neutrinos and antineutrinos in hot and dense matter. Additionally, we consider the production of neutrino pairs of all flavours through electron-positron annihilation, nucleon-nucleon bremsstrahlung and plasmon decay. The first one is expected to be the most relevant source of $\nu _{\text {x}}$’s in mildly and non-degenerate matter conditions, while the second one at very high density [113, 114]. We neglect their explicit contribution to the absorption opacity, since we expect it to be subdominant due to the pair nature of the inverse reactions, while their thermalisation effect is implicitly taken into account inside a NLS. Neutrino scattering off free nucleons is included as a major source of scattering opacity for neutrinos of all flavours and it is treated in the elastic approximation. In the case of $\nu _{\text {e}}$’s and ${\bar{\nu }}_{\text {e}}$’s, this opacity contribution is comparable to the one of absorption reactions, while in the case of $\nu _{\text {x}}$ this is the dominant one [see e.g. [68]]. Coherent scattering off nuclei is also included, even if the paucity of nuclei makes its impact negligible in the context of BNS mergers. It is important to recall that, at leading order, both the absorption and the scattering opacity off free nucleons depends quadratically on the energy of the incoming neutrinos. This quadratic dependence is taken into account when computing absorption opacities for the M0 scheme.

The scheme distinguishes number density weighted opacities, $\kappa ^0_{\nu ,\text {a}}$ and $\kappa ^0_{\nu ,\text {s}}$, that determine the rate at which neutrinos diffuse out of the material, from energy density weighted opacities, $\kappa ^1_{\nu ,\text {a}}$ and $\kappa ^1_{\nu ,\text {s}}$, that determine the rate at which energy is released due to the loss of neutrinos. The neutrino optical depth $\tau _\nu$ is evolved in time following the scheme presented in [115], which allows the optical depth profile to adapt to the complex geometry of the system. In particular, the optical depth evolves as:

$$\begin{aligned} \tau _\nu ^{n+1}=\max ( (k_{\nu ,\text {s}} + k_{\nu ,\text {a}})\text {d}l + \tau _\nu ^n)\,, \end{aligned}$$

(5)

where $\text {d}l$ is a local displacement of one grid point and the maximum is taken over all spatial directions.

The optical depth is used to define the effective emission rates:

$$\begin{aligned} R_\nu ^\text {eff} =\frac{R_\nu ^{\text {free}}}{1+t^0_\text {diff}(t^0_\text {loss})^{-1}}\,, \end{aligned}$$

(6)

where $t_\text {diff}$ is the effective diffusion time

$$\begin{aligned} t_\text {diff}^0 = {\mathcal {D}} \frac{(\tau ^0_\nu )^2}{\kappa ^0_{\nu ,\text {a}}+\kappa ^0_{\nu ,\text {s}}}\,, \end{aligned}$$

(7)

and $t_\text {loss}$ is the neutrino emission time scale

$$\begin{aligned} t_\text {loss}^0 = \frac{R_\nu ^{\text {free}}}{n_\nu }\,. \end{aligned}$$

(8)

In Eq. (7), ${\mathcal {D}}$ is a (dimensionless) tuning parameter set to 6², and $n_\nu$ in Eq. (8) is the neutrino number density computed assuming thermal and weak equilibrium. The effective energy emission rates $Q_\nu ^\text {eff}$ are computed with the same procedure as $R_\nu ^\text {eff}$, but using the appropriate opacities and optical depths. This method of computing effective rates provides a smooth interpolation between an estimate of the diffusion rate in optically thick condition and the local production rate in optically thin conditions, based on the optical depth.

Neutrino transport and absorption in optically thin conditions. Neutrino transport and absorption in optically thin conditions is accounted for by the moment scheme introduced in [105], called M0 scheme. Neutrinos are split into two components: a free-streaming one, $n_\nu ^\text {fs}$, and a trapped one, $n_\nu ^\text {trap}$, which is treated with the NLS previously described. The M0 scheme evolves the zeroth moment of the distribution function of free streaming neutrinos, and allows to compute their number densities and average energies on a polar grid.

This scheme assumes that neutrinos propagate radially at the speed of light along four-vectors:

$$\begin{aligned} k^\alpha = u^\alpha + r^\alpha \,, \end{aligned}$$

(9)

where $r^\alpha$ represents the spatial direction of propagation orthogonal to the fluid four-velocity $u^\alpha$. This assumption implies that the neutrino number current $J^\alpha$ equals $n_\nu ^\text {fs}k^\alpha$. Under these assumptions it is possible to show that the free-streaming neutrino number density, $n_\nu ^\textrm{fs}$, satisfies:

$$\begin{aligned} \nabla _\alpha (n_\nu ^\text {fs}k^\alpha ) = R^\text {eff}_\nu - \kappa _{\nu ,\text {a}} n_\nu ^\text {fs} \,, \end{aligned}$$

(10)

where $\kappa _{\nu ,\text {a}}$ is the absorption opacity. This finally results in an evolution equation for the neutrino number density, namely:

$$\begin{aligned}&\partial _t(\sqrt{-g}n^\text {fs}_\nu k^t) + \partial _r(\sqrt{-g}n^\text {fs}_\nu k^r)\nonumber \\&\quad =\sqrt{-g}(R^\text {eff}_\nu -\kappa ^\text {eff}_\nu n^\text {fs}_\nu )\,, \end{aligned}$$

(11)

where g is the four-metric determinant in spherical coordinates. This equation is solved on a series of independent radial rays using a first order, fully-implicit, finite volume method.

Free-streaming neutrino mean energies are estimated under the additional assumption of a stationary spacetime. Accordingly, $t^\alpha :=(\partial _t)^\alpha$ is assumed to be a Killing vector so that $p^\alpha _\nu (\partial _t)_\alpha$, with $p^\alpha$ being the neutrino four-momentum, is conserved. Therefore the quantity $\varepsilon _\nu = -p^\alpha t_\alpha$ represents the energy of neutrinos as seen by the “coordinate observer” (a non-physical observer with four-velocity $t^\alpha$), and can be rewritten as $\varepsilon _\nu = E_\nu \chi$, with $E_\nu$ the neutrino energy as measured by an observer comoving with the fluid and $\chi =-k_\alpha t^\alpha$. Within this approximation, the evolution equation for the average neutrino energy is written as:

$$\begin{aligned} n_\nu ^\text {fs}k^t\partial _t \varepsilon _\nu + n_\nu ^\text {fs}\partial _r k^r \varepsilon _\nu =\big (\chi Q_\nu ^\text {eff}-\varepsilon _\nu R^\text {eff}_\nu \big )\,, \end{aligned}$$

(12)

where $Q_\nu ^\text {eff}$ and $R_\nu ^\text {eff}$ are the effective neutrino energy and number emission rates taken from the NLS. This equation is solved using a fully-implicit upwind 1st order finite-difference method.

The coupling with hydrodynamics is handled by interpolating quantities from/to the standard Cartesian AMR grid at every timestep, by means of trilinear interpolation. In the setup of our sample of simulations, the M0 grid consists of 2048 rays uniformly spaced in latitude and longitude with a radial resolution $\varDelta r\approx 244\,{\textrm{m}}$.

The neutrino number and energy rates computed by the combined leakage and M0 schemes appear as sources in the Euler equations for the NS matter, see Sect. 2.1. The coupling in this case is handled, at every timestep, by first advancing the hydrodynamic quantities in time disregarding neutrino contributions; neutrino sources are then added to the Euler equations with a semi-implicit first-order method, in an operator split approach.

2.3.2 Equations of state

In our simulation sample, we consider six finite temperature, composition dependent EOSs, namely: LS220 [122], SLy4 [123, 124], DD2 [125, 126], SFHo [127], BHB$\varLambda \phi$ [128], and BLh [129]. They are widely used in the literature on BNS mergers and are broadly consistent with current constraints, including astrophysical constraints derived from GW observations [25, 48, 130,131,132]. The above EOSs satisfy properties of symmetric nuclear matter at saturation density. They also provide values for the symmetry energy and its slope in agreement with recent experimental estimates [133, 134], with the possible exception of PREX II results [135] that reported a quite large value of the slope of the symmetry energy at saturation density. The matter modelled by these EOSs is composed of neutron, protons, electrons, positrons and photons. One of them, namely BHB$\varLambda \phi$, also includes $\varLambda$-hyperons. In all our EOSs we do not take into account the presence of muons. They would lead to a slight softening of the EOS and their correct inclusion in the EOS may be important to describe the emission spectrum of neutrinos in a more accurate way. This task is left for future work.

The LS220 and SLy4 EOSs are based on a non-relativistic liquid drop model with a Skyrme-like interaction. This model includes surfaces effects and considers in the low density region an ideal classical gas formed by $\alpha$ particles and heavy nuclei. The latter are treated using the single nucleus approximation (SNA). The SLy4 EOS employed in this work is constructed on the original Skyrme parametrisation proposed in Ref. [123] for cold nuclear matter. It is extended to finite temperature [124], employing an improved version of the LS220 framework that includes non-local isospin asymmetric terms, a better treatment of nuclear surface properties, and a more consistent treatment of heavy nuclei sizes.

The DD2, SFHo, and BHB$\varLambda \phi$ EOSs are based on relativistic mean field (RMF) models. Besides single nucleons, their composition includes light nuclei (such as deuterium, tritium, helium) as well as heavy nuclei in nuclear statistical equilibrium (NSE). The Lagrangian that models the mean-field nuclear interaction is parametrised differently for the three EOSs. While DD2 and BHB$\varLambda \phi$ use density dependent coupling constants, the SFHo parametrisation employs constant couplings adjusted to reproduce NS radius measurements from low-mass X-ray binaries. In all three cases, the resulting RMF equations are solved in Hartree’s approximation.

The BLh EOS is a microscopical, finite temperature EOS obtained as an extension to the zero-temperature BL EOS [136]. At densities larger than 0.05 fm$^{-3}$, the latter was derived in the framework of the non-relativistic many-body Brueckner-Hartree-Fock (BHF) approach. The nucleon-nucleon interactions are described through a potential derived in the context of chiral effective theory [137]. They include two-body interactions [138] calculated up to next-to-next-to-next-to-leading order (N3LO), and an effective treatment of three-body interaction up to next-to-next-to-leading order (N2LO) [139]. Both interactions include contributions from $\varDelta$-excitation in the intermediate states of the nucleon-nucleon and three-nucleon interactions. Finite temperature and arbitrary nuclear composition effects are calculated using the finite temperature extension of the Brueckner–Bethe–Goldstone quantum many-body theory in the BHF approximation. At low densities the BLh EOS is smoothly connected to the SFHo EOS.

The EOSs employed in this work have been chosen in order to cover a broad range of stiffness. The stiffest EOS is the DD2 EOS, while the softest is the SLy4 EOS. These two EOSs support cold, non-rotating NSs maximum masses of $2.42{\hbox {M}_\odot }$ and $2.06{\hbox {M}_\odot }$, respectively. Operating on a broad stiffness range is important on one hand to avoid as much as possible any bias in our analysis, and on the other to allow us to look for universal relations in our data.

2.3.3 Viscosity

Slightly more than one third of the models analysed in this work employs the GRLES method of Ref. [140] to investigate the impact of turbulent viscosity on the merger dynamics [see also [141], for an alternative version of this formalism].

In essence, the GRLES method consists in taking into account that, due to finite resolution, any simulation deals only with a coarse-grained version of the hydrodynamics equations. Formally, this means introducing a linear filtering operator on the hydrodynamics variables that removes features at small scales (in our case this is simply the cell-averaging of the finite-volume discretization of the equations).

In turn, this implies that applying the filtering to the hydrodynamics equations requires the introduction of closure terms. In the resulting equations, the turbulent viscosity, $\nu _\text {T}$, is expressed in terms of the sound speed, $c_\text {s}$, and a free parameter, $\ell _\text {mix}$, that sets the characteristic length at which the turbulence operates, as $\nu _\text {T}=\ell _\text {mix} c_\text {s}$. In the simulations that we consider, $\ell _\text {mix}$ is estimated as a function of the rest mass density by fitting the results of very high resolution magnetohydrodynamics BNS merger simulations [142, 143].

3 Overview of simulations and analysis methods

3.1 Simulations sample

For our analysis, we consider a subset of the simulations presented in Ref. [1], whose setup is generic and not targeted to model a specific BNS configuration. In addition, we consider a subset of the simulations targeted to GW170817 and extensively discussed in Refs. [108, 144, 145], and data extracted from more recent simulations targeted to GW190425 [146]. Finally, we include also eight simulations which have not been published in earlier works but are presented for the first time in this article. In summary, we work on a sample of 66 simulations of 51 models of BNS mergers. The range of total gravitational mass $M_{\text {tot}}$ spanned by these models is $\left[ 2.600,\,3.438\right] M_{\odot }$ and the range of mass ratio q is $\left[ 1.0,\,1.82\right]$.

The reduced dimensionless tidal deformability of our set of models spans the wide range ${\tilde{\varLambda }}\in \left[ 90,\,1108\right]$. By comparison, data from the only two detected GW signals compatible with BNS mergers, namely GW170817 and GW190425, suggest that for those systems ${\tilde{\varLambda }}<700$ at the 90% confidence level [147]. However, we remind the reader that ${\tilde{\varLambda }}$ depends on the masses and mass ratio of the stars, so future events could also have larger ${\tilde{\varLambda }}$.

Regarding resolution, the sample consists of 15 LR simulations, 49 SR simulations, and 2 HR simulations. Where possible, we have decided to work with SR simulations because these tend to offer a better balance between accuracy and time extent of the post-merger data. Finally, 25 simulations out of the 66 employ the GRLES method described in Sect. 2.3.3 to account for viscous effects.

For each simulation we consider the neutrino energy luminosities and mean energies as extracted at the edge of the M0 computational domain and integrated over the outermost coordinate sphere. The luminosities and the mean energies are given in retarded time with respect to the time of merger (computed as the instant where the amplitude of the strain of the GW is maximum). The main properties of our sample of BNS simulations are summarised in Tables 2 and 3 of Appendix C.

3.2 Neutrino emission: a qualitative overview

Fig. 1

Evolution of the maximum density and temperature (top panels), normalised GW strain and GW luminosity (middle panels), and neutrino luminosity for the three neutrino species (bottom panels) for four models representative of the considered simulation categories. a PC simulation with LS220 EOS, NS masses of $1.772M_{\odot }$ and $1.065M_{\odot }$, and no viscosity; b VSL simulation of an equal mass binary ($1.364M_{\odot }$) with SFHo EOS and no viscosity; c DC simulation of an equal mass binary ($1.364M_{\odot }$) with SLy4 EOS and no viscosity; d LL simulation with DD2 EOS, NS masses of $1.509M_{\odot }$ and $1.235M_{\odot }$, with viscosity. The dashed vertical line in panel d indicates the change from linear to logarithmic scale in the time axis

We first present an overview of the observed properties of neutrino emission that are common to large subsets of models and simulations in our sample. According to the remnant fate, we distinguish our simulations into four categories: prompt collapse (PC), very short lived (VSL), delayed collapse (DC) and long lived (LL). We define PC simulations as the ones for which, at the time of merger, the minimum of the lapse function over the computational domain decreases monotonically. These conditions provide a proxy for detecting the collapse of the central object to a black hole (BH), i.e. for all these simulations the remnant is too heavy to sustain the formation of an massive neutron star (MNS). In VSL simulations, the merger remnant does not collapse promptly, but within 5 ms from the merger. DC simulations are those for which the collapse happens between 5 ms and the end of the simulation. Finally, in LL simulations no BH are observed until the simulation end.

Before discussing the main qualitative features, it is useful to summarise the origin of the neutrino emission. Neutrinos are emitted mostly from three sites: 1) from matter expanding from the contact interface between the two stars at merger and soon after it; 2) from the merger remnant, before collapse; 3) from the innermost, hot part of the post-merger accretion disc. The relative importance of these three sites varies depending on the dynamics of the system. For example, in PC simulations the remnant collapses immediately and the accretion disc has very low mass, making the contact interface during merger the only significant source of neutrinos. In LL models with high-mass ratio, the disc can be rather massive and contribute more to the neutrino emission than in near equal-mass models. The definitions of the above simulation classes therefore are also motivated by the mechanics of neutrino emission, since every simulation within one of these groups has similar properties and behaviour regarding neutrino luminosities and mean energies.

With reference to Fig. 1, we observe that in all cases the neutrino luminosity increases just before the merger. During the inspiral, tidal interaction heats up the two NSs, however this effect is expected to be small, $T \lesssim 1{\hbox {MeV}}$ [see e.g. [148]], and not accompanied by an intense neutrino emission. However, a non-negligible luminosity is observed in our simulations also during the inspiral. This is due to a spurious numeric increase in temperature ($T \lesssim 10{\hbox {MeV}}$) at the NS surfaces resulting from the fast NS motion inside a much more dilute atmosphere. Note that this has a negligible effect on the ejecta composition, since the emitting matter represents a small fraction both of the total mass and of the ejecta. A significant increase is observed around merger, due to the direct contact between the NS surfaces. This process continues during the merger and its immediate aftermath, causing the neutrino luminosity to peak at this time to typical values around $10^{53}\,{\hbox {erg}~{\hbox {s}}^{-1}}=100\,{\hbox {Bethe}}$. This is primarily due to the rapid increase in matter temperature (up to several tens of MeV) due to the NS collision and core fusion, two processes in which kinetic bulk energy is efficiently converted into thermal energy available to be radiated in neutrinos.

PC simulations present a single, relatively low peak generally between the merger and 1 ms after it. This is due to the main source of neutrinos, the merger remnant, being cut off by its collapse. In VSL simulations this peak is also present, but typically a few times higher than for PC ones. By contrast, simulations have, typically, between 3 to 4 well defined luminosity peaks in the first 10-15 ms after the merger for each neutrino flavour. We notice that what we consider as the “first peak” is always the highest one, and we disregard smaller, secondary peaks in the luminosity before it. While these secondary peaks are likely physical in origin, they cannot be modelled robustly in our simulations. In particular, their number, position and width vary with resolution (see Appendix B). While this statement can apply to the highest peak as well, it can still be unambiguously defined. Therefore we focus on its analysis. These luminosity peaks are likely related to the oscillations of the MNS in the early post merger. In this phase, the contractions and expansions of the merger remnant as it evolves towards a more stable configuration drive shock waves outwards through the remnant itself and the surrounding matter, raising its temperature via shock heating and therefore enhancing neutrino emission. Additionally, matter compressed at the NS collision interface and between the two merging cores is heated up and expelled from the centre of the remnants, expanding and decreasing its density inside the forming accretion disc. It is however non-trivial to link neutrino luminosity peaks to, e.g., features in the density evolution of the MNS or in the GW signal, as Fig. 1 illustrates. This is due to the fact that neutrinos can escape the system only when produced or transported outside the neutrinosphere, which is located $\sim 20\,{\hbox {km}}$ from the remnant and is itself evolving and growing in radius [68], making it very difficult to look for time coincidences.

Most VSL remnants approach a BH-torus configuration shortly after merger. We observe that after this point the neutrino luminosity decreases very rapidly, even if it does not drop to zero, as the inner, hot parts of the remaining disc are still neutrino sources. A similar behaviour can be seen in the DC case, but the drop in luminosity is not as steep as in VSL simulations and the post collapse luminosity is $\lesssim 50\%$ of the one before merger. This is due to the fact that the accretion disc mass is usually larger if the system is less massive (i.e. less prone to a fast collapse) or asymmetric (i.e. more prone to a tidal deformation of the secondary). Indeed, since the disc formation process lasts for several milliseconds after merger [95, 145], a faster collapse of the central MNS prevents the formation of a massive disc that can sustain a significant luminosity also after the MNS collapse. We recall, in this respect, that the collapse of the MNS drags inside the apparent horizon roughly half of the disc mass, corresponding to the innermost, hotter portion of the disc.

In LL simulations, after the first oscillatory phase, the neutrino luminosity decreases exponentially in time at a much smaller and steady rate, remaining comparable to the luminosity observed in the first milliseconds after merger on time scales even of hundreds of milliseconds, i.e. comparable to the MNS lifetime. In a MNS+disc configuration, both the central object and the disc significantly contribute to the neutrino emission. The cooling of the central object and the release of gravitational energy inside the accretion disc are both active mechanisms in sustaining the neutrino emission over the longer cooling and accretion time scales. In particular, the optical depth for the most relevant neutrino energies inside the disc is of the order of a few, while it is two to three orders of magnitude larger inside the central MNS. As a consequence, the cooling time scale of the disc is a few ms and its luminosity is sustained until accretion takes place, while the cooling time scale of the MNS is of several seconds and the corresponding luminosity lasts until the central object is hot enough [see e.g. [59]]. In Ref. [149], it was estimated that a LL remnant should liberate $\sim 0.08 M_{\odot } c^2$ in its cooling phase. This corresponds to $\sim 1.4 \times 10^{53} \textrm{erg}$. This is compatible with a total neutrino luminosity of the order of $10^{53}$–$10^{52}\mathrm{erg~s^{-1}}$, lasting for a few seconds.

Regarding the relative abundance of neutrino species, during and after the merger positron captures on free neutrons are favoured since matter is initially extremely neutron rich ($Y_{\text {e}}\sim 0.1$) and hot $(T\sim$ 10–50 MeV). Therefore the electron antineutrino luminosity is dominant in every model. For electron neutrinos the most relevant production reaction is the capture of electrons on free protons. Due to the relative paucity of protons, $\nu _{\text {e}}$ are emitted in a subdominant fashion with respect to ${\bar{\nu }}_{\text {e}}$. Moreover $\nu _{\text {e}}$ are also more easily absorbed in typically thin conditions in their way out from the remnant. Around the time of merger heavy flavour neutrinos are emitted with a luminosity comparable to that of electron neutrinos. These heavy flavour neutrinos are produced by very hot matter ($T \sim$ tens of MeV) initially expelled from the bouncing remnant and rapidly expanding in optically thin conditions. Electron-positron annihilation and plasmon decay (which are the most relevant reactions producing $\nu _{\text {x}}$’s) have an extreme dependence on temperature (with production rates $Q_{\nu _{\text {x}}}\propto T^9$, [see e.g. [71]]). Once the remnant has settled on a quasi-stationary configuration, $\nu _{\text {x}}$’s emission mostly reduces to the thermal diffusion from the optically thick central remnant. On the other hand electron (anti-)neutrinos are mostly produced via electron/positron captures on nucleons, reactions with a milder dependence on the temperature and happening also inside the accretion disc. Therefore, as the system stabilises and cools, the heavy neutrino production is significantly reduced with respect to the other neutrino flavours. In the case of LL simulations, we also note that with time the difference in luminosity between $\nu _{\text {e}}$ and ${\bar{\nu }}_{\text {e}}$ tends to decrease, such that for all long lasting simulations we observe that $L_{{\bar{\nu }}_{\text {e}}} \sim L_{\nu _{\text {e}}}$. This is due to the matter being leptonised, reducing the dominance of the ${\bar{\nu }}_{\text {e}}$’s emission mechanisms.

The neutrino mean energies present a different pattern with respect to the neutrino luminosities. In the first few milliseconds after merger, we observe that they oscillate wildly and rapidly. However, this might be an artefact due to the approximate character of the neutrino transport schemes we rely on. We therefore do not attempt to characterise this phase any further. After this oscillatory phase the neutrino mean energies show a much more stable behaviour, in fact they are nearly constant until the end of simulation or the collapse to BH of the merger remnant. Clearly this second phase is only present in DC and LL simulations. This behaviour is related to the thermodynamic conditions of matter around the surfaces of neutrino decoupling. Neutrinos leave the system if emitted outside the neutrinosphere, and their energy distribution is strongly influenced by the temperature of the emitting medium at the density where thermal and weak decoupling between neutrinos and matter occurs. In the aftermath of BNS mergers, the neutrinospheres for each flavour and neutrino energy are mostly determined by the density profile inside the disc [68], and the latter changes very slowly, only over the accretion time scale. This in turn implies that the neutrinos are emitted by matter whose thermodynamic conditions do not significantly vary within the analysed time.

3.3 Analysis strategy

Based on the general features summarised in Sect. 3.2, we focus our analysis on neutrino luminosities $L_\nu$ and mean energies $E_\nu$ for all three flavours, i.e. for $\nu \in \{\nu _{\text {e}},{\bar{\nu }}_{\text {e}},\nu _{\text {x}}\}$.

For all simulations we consider the peak luminosity $L_{\text {peak},\nu }$, which is simply the highest peak for a given simulation. We also examine the full width at half maximum (FWHM) $\varGamma$ of the peak by fitting the neutrino luminosity in a window of width 1 ms centred on the peak luminosity time $t_{\text {peak}}$. As a fitting function, we employ a Gaussian function:

$$\begin{aligned} L=L_\text {peak}\exp {\left( -\frac{\left( t-t_\text {peak} \right) ^2}{2\sigma ^2}\right) }\,, \end{aligned}$$

(13)

where the amplitude and peak centre position are fixed as the peak luminosity and time, respectively, while the peak width $\sigma$ is the fitting parameter. Finally the FWHM is related to $\sigma$ as:

$$\begin{aligned} \varGamma = 2\sqrt{2\ln {2}}\sigma \,. \end{aligned}$$

(14)

For DC and LL simulations, we also analyse the values of the time-averaged luminosity $\langle L_\nu \rangle$ and the time-averaged neutrino mean energy $\langle E_\nu \rangle$. Explicitly, the time average of a quantity $X_\nu$ is computed as:

$$\begin{aligned} \langle X_\nu \rangle = \frac{1}{t_{\text {stop}}-t_{\text {merg}}} \int ^{t_{\text {stop}}}_{t_\text {merg}}X_\nu (t)dt\,, \end{aligned}$$

(15)

where $t_\text {merg}$ is the time of merger and $t_\text {stop}$ is a suitable final time. To these time-averages we associate their standard deviations, computed as:

$$\begin{aligned} \sigma _{X_\nu } = \sqrt{\langle X_\nu ^2 \rangle -\langle X_\nu \rangle ^2}\,, \end{aligned}$$

(16)

where $\langle X_\nu ^2 \rangle$ is the average of $X^2$. For the luminosity, the time average is computed using a window starting at the time of merger and extending either to 10 ms after merger or until BH formation. This window has been chosen to be long enough so that computing the average is meaningful, but not so long that in LL simulations the final value is influenced by the late time decrease. For the time-averaged neutrino mean energies we select a different time window, extending from the point at which the neutrino mean energies begin to stabilise (typically $2-5\,{\hbox {ms}}$ after the merger), to either the end of the simulation or BH formation. In this case too the window has been chosen to be long enough to get a meaningful average energy. Differently from the case of the time-averaged luminosity however, the final computed value is not sensitive to the end point of the window, because as noted above mean energies are essentially constant until collapse or the end of the simulation.

Fig. 2

Peak luminosity $L_\text {peak}$ plotted against the reduced dimensionless tidal deformability ${\tilde{\varLambda }}$ for electron neutrinos (panel a), electron antineutrinos (b) and heavy lepton neutrinos (c). Colour indicates the BNS mass ratio. Note the different abscissa scales in the three panels

4 Results

4.1 Luminosity peak and peak broadness

We start by exploring the peak luminosities for the different neutrino species. Fig. 2 displays their dependency on the tidal deformability and mass ratio for our BNS models. The peak luminosities approximately span the range $1\cdot 10^{52}-5.5\cdot 10^{53}\,{\hbox {erg}~\hbox {s}^{-1}}$ for electron antineutrinos, while the other two flavours do not go beyond $\sim 1.7 \cdot 10^{53}\,{\hbox {erg}~\hbox {s}^{-1}}$ even in the most extreme cases. The extreme neutron richness and high temperatures of the MNS matter enhance the production of electron antineutrinos, hence the differences in the peak strengths. Within the observed ranges we notice that the peak luminosity values follow very similar trends in different neutrino species. There is a roughly constant factor of $\sim 3$ between $\nu _{\text {e}}$ and $\nu _{\text {x}}$ neutrinos with respect to ${\bar{\nu }}_{\text {e}}$ ones. This similarity can be understood by noting that the qualitative behaviour of neutrino emission in this phase is influenced more by the bulk dynamics of matter than the specifics of neutrino interactions.

PC simulations have very low peak luminosities, up to six times lower than other models. For symmetric systems, this is due to two related phenomena. The merger remnant collapses right after merger and a massive disc cannot form since most of the matter is caught in the collapse. While equal-mass PC simulations cluster at low values of ${\tilde{\varLambda }}$ (bottom lower part of Fig. 2), high-q models with higher ${\tilde{\varLambda }}$ can also result in a prompt collapse. With respect to q (and thus to ${\tilde{\varLambda }}$) we observe a slightly upward trend, which can be understood by noting that the lighter object is more easily tidally disrupted as q increases, allowing for a more massive disc that contributes to neutrino emission.

The remaining simulation categories show a different and much stronger dependence on ${\tilde{\varLambda }}$. Equal- or nearly equal-mass DC and LL models generally have higher peak luminosity than their asymmetric counterparts and the luminosity peak values present a downward trend with respect to ${\tilde{\varLambda }}$. Systems characterised by a higher tidal deformability contain less compact stars, which collide in a less violent fashion. Under these conditions, shock heating is less prominent and neutrino emission is correspondingly smaller. We observe the largest peak luminosities for $380 \lesssim {\tilde{\varLambda }}\lesssim 420$; note however that the limits of this interval depend on the sample of EOSs and masses that we consider, and might change by considering a wider sample. Systems with higher mass ratio fall on a second branch, because the tidal disruption of one of the two stars leads to less violent coalescences. Because an increased tidal disruption also tends to increase the disc mass and its contribution to the neutrino emission, this trend is not strictly monotonic with respect to the mass ratio. This behaviour can be contrasted to the analogous one of the time-averaged neutrino luminosity, where both branches (the equal- and unequal-mass ones) show a much more well delineated trend with respect to ${\tilde{\varLambda }}$ (see Sect. 4.2).

Note finally that the VSL simulations provide a sort of transition between the $q \approx 1$ maximum of the DC+LL sample and the PC $q \approx 1$ branch.

We then consider the broadness of the first peak of the neutrino luminosity, computed as detailed in Sect. 3.3. To measure the goodness of the fit, we consider the relative residuals between the data and the fit at fixed times. We observe that, for every flavour and every simulation, they do not exceed 5% at any point in the fit interval. Furthermore the coefficient of determination of the fit $R^2$ is $\sim 0.99$ in all cases. To further test the goodness of the fit we also compare the fitted $\varGamma$ values with the ones calculated directly from the light curves as the FWHM (when the peak shapes allow this calculation, i.e. for PC and VSL simulations). The differences in the results of these two procedures does not exceed the $20\%$ level in most cases, barring two PC outliers, with very low and broad peaks and a relative difference of $30\%$. Fig. 3 presents the dependence of $L_\text {peak}$ to $\varGamma$. Here too we observe two trends, one for the PC simulations and one for the other three categories. In the VSL, DC and LL simulations, as the peak luminosity decreases with increasing ${\tilde{\varLambda }}$, the peak broadness increases instead. The peaks of the neutrino luminosity arise as the results of shock waves generated by the oscillations of the merger remnant. The time scale of these oscillations is $\sim 1\,{\hbox {ms}}$, which is indeed the typical value of $\varGamma$. The time scale of the oscillations is related to the free-fall time scale of the remnant, which scales as $t_{\text {ff}}\propto \langle {\rho }\rangle ^{-1/2}$, where $\langle {\rho }\rangle$ is the mean density of the MNS. Since stars with higher deformability have generally lower $\langle \rho \rangle$, their oscillations time scale is longer, and $\varGamma$ is broader. This observation can be recast in a way which is physically more meaningful. We note that for the VSL+DC+LL branch, $L_{\text {peak}}$ and $\varGamma$ are loosely inversely proportional to each other, ans their product roughly constant (see Fig. 3). The time integral of the Gaussian we employed as fitting function,

$$\begin{aligned} E_{\text {peak}}=\frac{1}{2}\sqrt{\frac{\pi }{\ln 2}} ~L_{\text {peak}}\varGamma \,, \end{aligned}$$

(17)

represents an estimate of the energy released by the first neutrino peak. We plot this quantity in Fig. 4. Clearly $E_{\text {peak}}$ is broadly constant, with typical values of $0.75\cdot 10^{50}\,{\hbox {erg}}$ for $\nu _{\text {e}}$’s and $\nu _{\text {x}}$’s, and of $2.25\cdot 10^{50}\,{\hbox {erg}}$ for ${\bar{\nu }}_{\text {e}}$’s with a maximum deviation of a factor of 2. This allows us to provide a very concise characterisation of the first neutrino luminosity peak: as long as the remnant does not collapse promptly after merger, the first luminosity peak releases a roughly constant amount of energy of $\approx 6 \times 10^{50}{\hbox {erg}}$.

Finally, it is also clear that the argument outlined above does not apply to PC simulations, which due to immediate collapse have not only very low $L_\text {peak}$, but also very low $E_{\text {peak}}$. Furthermore no time scale argument can apply to a collapsed remnant since it does not emit neutrinos. Indeed while the typical values of $\varGamma$ are the same for PC simulations too, they do not follow any particular trend with respect to either ${\tilde{\varLambda }}$ or q.

Fig. 3

Peak luminosity, $L_\text {peak}$, as a function of the FWHM of the peak, $\varGamma$, for electron neutrinos (panel a), electron antineutrinos (b) and heavy lepton neutrinos (c). Colour indicates the BNS mass ratio. Note the different abscissa scales in the three panels

Fig. 4

Energy emitted by the peak $E_{\text {peak}}$ as a function of the reduced dimensionless tidal parameter ${\tilde{\varLambda }}$ for electron neutrinos (panel a), electron antineutrinos (b) and heavy lepton neutrinos (c). Colour indicates the BNS mass ratio

4.2 Time-averaged luminosities

Next we examine the average neutrino luminosity for different neutrino species, showed in Fig. 5. We limit the analysis to the DC and LL classes, since for the other two the average luminosity is not well defined. The values span the range $0.6\cdot 10^{53}-1.4\cdot 10^{53}\,{\hbox {erg}~\hbox {s}^{-1}}$ for electron antineutrinos. The other two flavours mostly vary in the range $0.2\cdot 10^{53}-0.6\cdot 10^{53}\,{\hbox {erg}~\hbox {s}^{-1}}$. Similarly to the peak luminosities, different neutrino species follow very similar trends, differing in this case by a roughly constant factor between 2 and 2.5 . The physical explanation of this trend outlined in the previous section applies here too. The reduction of the scaling factor with respect to $L_{\textrm{peak}}$ is due to the unbalanced ${\bar{\nu }}_{\text {e}}$ emission, which leptonises the remnant, partially suppressing its own emission mechanism.

Fig. 5

Average luminosity $\langle L_\nu \rangle$ as a function of the reduced dimensionless tidal deformability ${\tilde{\varLambda }}$ for electron neutrinos (panel a), electron antineutrinos (b) and heavy lepton neutrinos (c) for DC and LL simulations. Colour indicates the BNS mass ratio while grey bars the standard deviation of the values

To corroborate these observations we consider the dependence of the average luminosities for a flavour on the other two, as shown in Fig. 6. Clearly there is a linear correlation between the average luminosities of any two flavours, with a proportionality factor of $\sim 2.5$ between electron neutrinos and antineutrinos, and a slightly smaller factor between electron antineutrinos and heavy neutrinos (however we refrain from fitting a straight line trough our data points, judging their quality too poor to warrant such an analysis).

Fig. 6

Average luminosity $\langle L_\nu \rangle$ of two neutrino flavours plotted against each other for DC and LL simulations. Panel a: $\langle L_{{\bar{\nu }}_{\text {e}}}\rangle$ vs. $\langle L_{\nu _{\text {e}}}\rangle$; panel b: $\langle L_{{\bar{\nu }}_{\text {e}}}\rangle$ vs. $\langle L_{\nu _{\text {e}}}\rangle$. Colour indicates the BNS mass ratio and the dashed lines have slopes of 2.5, for comparison

In Fig. 5 we also see that equal-mass models values decrease with increasing tidal deformability, and in this case too the explanation outlined in Sect. 4.1 holds true. Also in this case varying the mass ratio creates a second branch, with generally smaller average luminosities than equal-mass binaries. It is however much more prominent in the case of average luminosities and our data suggests it is monotonically increasing with respect to ${\tilde{\varLambda }}$, at least for ${\tilde{\varLambda }}\lesssim 700$.

The explanation of the differences between the peak and the average luminosities is in the act of taking a time average. Peak luminosities are associated to a transient and quite violent phase, whose properties cannot be satisfactorily described with a single parameter such as ${\tilde{\varLambda }}$ or q. Therefore it is to be expected for the peak luminosities to show a larger variability. On the other hand taking an average value can help to better isolate a trend present in the data, as shown in Fig. 5. We find further support for this line of reasoning by looking at the grey bars in Fig. 5, representing the time variability of the data around the average values. The bars are quite wide, spanning a range that in some cases is as wide as the value of the average value to which they are associated: stated differently, the neutrino luminosities oscillate rather widely as a function of time. Note that the origin of this variability is physical, being linked to e.g. the oscillations of the central object. Moreover their width also shows a trend with ${\tilde{\varLambda }}$: BNS mergers characterised by smaller ${\tilde{\varLambda }}$ and $q \approx 1$ present a more significant variability between the peaks and the valleys in the luminosity behaviour, reflecting the more violent dynamics of the merger.

Different resolutions and/or the inclusion of physical viscosity in the simulations do not seem to have a significant impact on the major results concerning the peak and average luminosities. A more detailed discussion about these points is documented in Appendices A and B.

4.3 Long term behaviour of the luminosity

In order to better characterise the time evolution of the neutrino luminosity over longer time scales, we focus on LL remnants and only select simulations that extend further than 20 ms after merger, for a total of 10 simulations that last between 40 and 110 ms. We choose the simple model

$$\begin{aligned} L_{\nu }(t)=L_{0,\nu }\exp \left( \frac{\alpha _{\nu }}{t}-\frac{t}{\beta _{\nu }} \right) \text { for }\nu \in {(\nu _{\text {e}},{\bar{\nu }}_{\text {e}},\nu _{\text {x}})}\,, \end{aligned}$$

(18)

and fit $L_{0,\nu }$, $\alpha _{\nu }\ge 0$ and $\beta _{\nu }>0$ to the neutrino luminosity curves, starting from the time at which all neutrino flavours monotonically decrease until the end of the simulation. Equation (18) is an exponential decay, augmented by a term which allows for deviations from a purely decaying exponential at early times.

The relative residuals between the data and the fit vary by up to 15% for heavy lepton neutrinos, and up to 10% for other flavours. We also observe that the largest residuals are observed at early times, when the luminosity is still characterised by residual oscillations. We also compute the coefficient of determination $R^2$, which equals $\sim 0.93$ for heavy lepton neutrinos and $\sim 0.99$ for the other two flavours. We conclude that Eq. (18) is a good description of the long term evolution of neutrino luminosities.

Fig. 7

Neutrino luminosity decay time scale $\beta _\nu$ as a function of the reduced dimensionless tidal deformability ${\tilde{\varLambda }}$ for electron neutrinos (panel a), electron antineutrinos (b) and heavy lepton neutrinos (c) for the longest LL simulations. Colour indicates the BNS mass ratio and the dashed line in panel c indicates the upper limit of the other two panels for ease of comparison

We focus on the coefficients $\beta _\nu$, i.e. the time scale over which the luminosity drops, shown in Fig. 7. This quantity does not seem to correlate with either ${\tilde{\varLambda }}$ or q, but a few interesting observations are possible. Typical values of $\beta _{\nu }$ for electron neutrinos and antineutrinos are of the order of 100 ms. Barring a few outlying points, the corresponding value for heavy neutrinos is between 100 and 400 ms. These are rather long time scales, compared to the dynamical time scales associated with the MNS ($\sim 1\,{\hbox {ms}}$). Clearly the decline of neutrino emission reactions is a steady and relatively slow process, associated with the cooling of matter in the remnant, and indeed a time scale of several hundreds of milliseconds is more in line with both the cooling time scale of the MNS and with the accretion time scale of the disc [see e.g. [59, 150], and Appendix A].

The difference in the decrease rate between $\nu _{\text {e}}$/${\bar{\nu }}_{\text {e}}$’s and $\nu _{\text {x}}$’s is related to two causes. First, the neutrino origin: $\nu _{\text {e}}$’s and ${\bar{\nu }}_{\text {e}}$’s are both emitted by the accretion disc and the central MNS, while $\nu _{\text {x}}$’s mostly by the latter. Second, the different mean energies at the decoupling surfaces: $\nu _{\text {x}}$’s decouple deeper inside the remnant and their spectrum is significantly harder (see next section). These hotter neutrinos still diffuse between the equilibrium decoupling surface and the last scattering surface, due to the opacity provided by quasi-elastic scattering off free baryons. Since the cross section for this process depends quadratically on the neutrino energy, the opacity for $\nu _{\text {x}}$’s (and consequently also its cooling time scale) is significantly larger and the cooling of the deepest layers proceeds at a slower pace.

Extrapolating Eq. (18) to late times, the total emitted energy would be a few times $10^{52}\,{\hbox {erg}}$, i.e. almost one order of magnitude smaller than expected. We speculate that the exponential decrease we observe for $\nu _{\text {e}}$’s and ${\bar{\nu }}_{\text {e}}$’s is mostly due to the evolution of the accretion luminosity. However, once a significant portion of the disc has been consumed, the luminosity coming from the cooling of the central object will take over and it will likely decrease with a different time scale, which our fit over a limited time window cannot account for.

Fig. 8

Neutrino average mean energy $\langle E_\nu \rangle$ plotted against the reduced tidal deformability ${\tilde{\varLambda }}$ for electron neutrinos (panel a), electron antineutrinos (b) and heavy lepton neutrinos (c) for DC and LL simulations. Grey bars indicate the standard deviation in time from the average value

4.4 Time-averaged mean energies

The neutrino mean energies of DC and LL simulations plotted in Fig. 8 present a radically different behaviour compared to the luminosities. The typical energy values are $\sim 10\,{\hbox {MeV}}$, $\sim 14\,{\hbox {MeV}}$ and $\sim 23\,{\hbox {MeV}}$ for electron neutrinos, electron antineutrinos and heavy neutrinos, respectively (note that these are the same values reported in Ref. [68] and references therein). This hierarchy can be explained in relation to the properties of the neutrino decoupling regions. Of relevance here are the equilibrium surfaces, where neutrinos decouple from the fluid but are not yet free-streaming. It has been shown (for long-lived remnants) that these surfaces lie at increasing radii further away from the remnant for heavy neutrinos, electron antineutrinos and electron neutrinos, in this order. As temperature also decreases further away from the remnant, this explains the energy hierarchy between neutrino flavours. Furthermore the grey bars, representing the time variability of the mean energies (cf. Sect. 3.3), are extremely small, not being even visible in the leftmost two panels. Stated differently the neutrino mean energies are constant in the early post merger phase. This can be explained by noting that the thermodynamic conditions of matter at the surface of neutrino decoupling are constant in time in the later part of the evolution, since the location of the neutrinosphere itself does not evolve significantly at this point [68].

A similar observation also explains why the neutrino mean do not depend on the masses or EOS. The thermodynamics condition at the neutrinosphere are not only constant in time, but being located at rather large radii ($\sim 20\,{\hbox {km}}$), they are also rather insensitive to variations in the bulk dynamics of the system. We speculate that this could result from two reasons. On one hand, since the location of the equilibrium decoupling surface depends at leading order on the matter density, it is likely that the properties of the accretion disc (and in particular of the density-temperature profile) are rather independent from the specific binary system, especially once the disc has reached a quasi-stationary state and a high degree of axisymmetry. On the other hand, matter temperature also influences the neutrino opacity, mostly through the energy of the diffusion of thermal neutrinos. If a disc is hotter, the larger temperatures increase the opacity inside the disc, moving the decoupling surfaces at larger radii and, thus, lower temperatures. Clearly, these two effects tends to compensate each other, providing similar decoupling temperatures in all cases.

5 Discussion

5.1 Comparison with GW luminosities

Neutrinos provide the most relevant radiation loss from merging BNSs on the cooling time scale of the remnant, but the inspiral and the early post-merger ($\lesssim 20\,{\hbox {ms}}$) are GW-dominated [151]. In Fig. 9 we present a comparison between the GW and total neutrino peak luminosities, the former being computed as the first peak that the GW luminosity reaches during the merger. We recognise three different regimes. For VSL, DC and LL near-equal mass ($q \lesssim 1.25$) systems, there is a correlation between the luminosity in GWs and $\nu$’s. This is due to the fact that neutrino radiation is emitted by the same matter that produces also the GW emission. Since the binary properties that boost the GW emission [see, e.g. [152]] are the same ones that increase the remnant temperature, the two luminosities increase together. If the mass ratio becomes significantly higher than 1, ${\tilde{\varLambda }}$ decreases and both $L_{\textrm{peak,GW}}$ and $L_{\textrm{peak},\nu }$ decrease, but the reduction in GWs is less significant. This is due to the fact that the strong-field behaviour for $L_{\textrm{GW,peak}}$ is not precisely captured by ${\tilde{\varLambda }}$, but by the so-called $\kappa _2^L$ parameter [152]. In particular $\kappa _2^L$ is the perturbative parameter that enters the 5th order post-Newtonian term related to the tidal effects in the binary evolution. Finally, in the PC cases the two luminosities follow opposite trends: GW emission is the brightest for PCs resulting from symmetric BNSs merger, but these are the systems for which $L_{\mathrm{peak,\nu }}$ is the lowest (see Sect. 4.1). This effect is partially mitigated by the tidal disruption of the secondary happening in the high-q cases.

Fig. 9

Neutrino peak luminosity $L_{\nu \text {,peak}}$ as a function GW peak luminosity $L_{\text {GW,peak}}$. Colour indicates the BNS mass ratio

5.2 Influence on the electron fraction and kilonova colour

Neutrino interactions change the electron fraction, $Y_{\text {e}}$, of matter through charged current reactions, including electron, positron, $\nu _{\text {e}}$ and ${\bar{\nu }}_{\text {e}}$ captures on free neutrons and protons. All these reactions are relevant inside the neutrino surfaces to change $Y_{\text {e}}$ from cold, neutrino-less, $\beta$-equilibrium conditions ($Y_{\text {e}}\sim 0.05$ for the relevant densities) to finite temperature, neutrino trapped equilibrium conditions. Additionally, neutrino emission and irradiation can further change $Y_{\text {e}}$ also outside the neutrino surface in out-of-equilibrium conditions. Simulations including neutrino transport can follow in detail the evolution of the ejecta properties. The contribution of the different processes and their outcome can also be analysed in post-processing [see e.g. [153]]. Here we want to focus on a simpler question: how does the variation in the neutrino luminosity observed in our simulation sample possibly translate in a variation of $Y_{\text {e}}$ for a representative expanding fluid elements?

To answer this question, we assume a simplified model for the evolution of the $\nu _i$ (with $i=\nu _{\text {e}},{\bar{\nu }}_{\text {e}}$) luminosities:

$$\begin{aligned} L_{\nu _i}(t) = {\left\{ \begin{array}{ll} \langle L_{\nu _i} \rangle &{} 0< t < t_{\nu _i} \, , \\ \langle L_{\nu _i} \rangle \exp \left( -\frac{t-t_{\nu _i}}{\beta _{\nu }} \right) &{} t \ge t_{\nu _i} \,, \end{array}\right. } \end{aligned}$$

(19)

where $\langle L_{\nu _i} \rangle$ are the average luminosity presented in Sect. 4.2 and t the time after the escape of the fluid element from the neutrino surface. We set $\beta _{{\bar{\nu }}_{\text {e}}}=\beta _{\nu _{\text {e}}}=100\,{\hbox {ms}}$, based on Fig. 7. For ${\bar{\nu }}_{\text {e}}$’s we assume $t_{{\bar{\nu }}_{\text {e}}}=10\,{\hbox {ms}}$ while $t_{\nu _{\text {e}}}$ is fixed by the condition

$$\begin{aligned} L_{\nu _{\text {e}}}(t_{\nu _{\text {e}}}) = L_{{\bar{\nu }}_{\text {e}}}(t_{\nu _{\text {e}}})\,, \end{aligned}$$

(20)

meaning that the two luminosities are the same on the time scale set by $\beta _{\nu }$, as visible in the long term evolution of our LL models.

We further consider constant mean energies, equal to the average ones extracted from the simulations and presented in Sect. 4.4. We compute the evolution of $Y_{\text {e}}$ based on the equation

$$\begin{aligned} \frac{\text {d}Y_e}{\text {d}t} = \lambda _{\nu _{\text {e}}} (1-Y_{\text {e}}) - \lambda _{{\bar{\nu }}_{\text {e}}} Y_{\text {e}}\,, \end{aligned}$$

(21)

where $\lambda _{\nu _{\text {e}}}$ and $\lambda _{{\bar{\nu }}_{\text {e}}}$ are the $\nu _{\text {e}}$ and ${\bar{\nu }}_{\text {e}}$ capture rates, respectively. The expressions of $\lambda _{\nu _{\text {e}}}$ and $\lambda _{{\bar{\nu }}_{\text {e}}}$ are taken from equations (C.4)–(C.10) and (3) of [153]. To better focus on the role of luminosities, we neglect the impact of electron and positron captures outside the neutrino surfaces. This approximation is valid as long as the temperature in the ejecta expanding outside the neutrino surfaces is below a few MeV [154]. According to Ref. [68], typical temperatures at the relevant outermost $\nu _e$ surfaces are 3-5 MeV, depending on the EOS stiffness. According to Ref. [153], $e^{\pm }$ captures alone combine in such a way that they do not change significantly $Y_{\text {e}}$ (see their “capture” case in Figs. 6 or 7; however, see [82] for different conclusions, possibly due to hotter ejecta). From equation (3) the neutrino flux depends on the radial distance and its evolution. We consider $R(t) = v t + R_0$ where $R_0$ is the typical radial distance of the neutrino surface and v the ejecta speed. We further know that neutrino emission is not isotropic, due to the shadow effect provided by dense matter in the disc along the equatorial plane. We then consider two possible directions identified by the polar angle $\theta$, namely $\theta =0$ (polar direction) and $\theta =\pi /2$ (orbital plane), and the angular dependence implied by equation (3) in [153], assuming $\alpha = 2$, which corresponds to a polar flux three times larger than the equatorial one.

Fig. 10

Estimated electron fraction $Y_{\text {e}}$ of an ejected fluid element as a function of the reduced tidal deformability ${\tilde{\varLambda }}$ for dynamical ejecta along the polar axis (top left panel, $Y_{\text {e}}(t=0)=0.2$) and orbital plane (top right, $Y_{\text {e}}(t=0)=0.05$), and for a disc wind along the polar axis (bottom left, $Y_{\text {e}}(t=0)=0.2$) and orbital plane (bottom right, $Y_{\text {e}}(t=0)=0.2$) for LL simulations. Colours indicate the BNS mass ratio while dashed lines indicate the threshold below which the mass fraction of lanthanides and actinides produced in the r-process nucleosynthesis increases above 10% for typical ejecta conditions ($Y_{\text {e}}=0.22$). The grey triangles represent, for each simulation, the corresponding equilibrium $Y_e$. See the text for more details

We consider two kinds of ejecta: the dynamical and the disc wind ejecta. The dynamical ejecta [see e.g. [55, 65, 82, 87, 95, 105, 116, 118, 119, 145, 155,156,157,158,159,160,161,162,163]] are the matter expelled within a few dynamical time scales after merger ($\lesssim 5\,{\hbox {ms}}$), with typical average speeds ranging between 0.1-0.3c, by tidal torques and shock waves propagating inside the remnant. We compute the speed of the ejecta as a function of ${\tilde{\varLambda }}$ and q, based on the fitting formula equation (6) presented in [67], using in particular results from the M0RefSet dataset. This fit is a second order polynomial that predicts the largest speeds either for $q \lesssim 1.2$ and ${\tilde{\varLambda }} \lesssim 400$ (corresponding to very violent PC mergers) or for ${\tilde{\varLambda }} > rsim 1000$, in which tidal ejection is very effective. Disc winds [see e.g. [52, 57, 59, 62,63,64, 144, 150, 164,165,166,167,168,169,170,171,172,173,174,175]] are possibly expelled on the disc evolution time scale ($\sim$ 10 ms–1 s) by a variety of mechanisms, including neutrino absorption itself, nuclear recombination following viscous spreading of the disc, spiral wave triggered by long-standing $m=1$ bar modes in the remnant, magnetic processes. In this case, the ejection speed is expected to be $\sim$ 0.05–0.1 c. In our calculation, we consider a representative value of 0.08c³. However, in order to account for the fact that disc winds are produced on rather long time scales, ${\mathcal {O}}(100\,{\hbox {ms}})$, we actually consider a lower velocity value of 0.008c for the first 100 ms of evolution, necessary for a fluid element to reach the typical nuclear recombination radius (250–300 km) inside the disc within 100 ms [see e.g. [59, 149, 150, 165, 168]]. Finally, we smoothly connect the two values. Therefore for the wind velocity, we impose: $v_{\text {wind}} = 0.008\cdot (1 - f) + 0.08\cdot f$, where $f = ( 1 + \exp {\left( -(t - 100 \,{\hbox {ms}})/2\,{\hbox {ms}} \right) })^{-1}$. For the ejecta expelled in the orbital plane we assume $R_0=20\,{\hbox {km}}$, corresponding to the typical radius of the most relevant $\nu _{\text {e}}$ and ${\bar{\nu }}_{\text {e}}$ neutrino surfaces inside the disc [68]. For the ejecta expelled along the polar axis, we consider $R_0=15\,{\hbox {km}}$, corresponding to the radius of the MNS.

This analysis critically relies on the initial $Y_{\text {e}}$. A distributions of $Y_{\text {e}}$ at the neutrino surface is expected. However, here we rely on representative values just to focus on the impact of the luminosity variation. For the dynamical ejecta of tidal origin moving across the equator, we take $Y_{\text {e}}(t=0) = 0.05$. This matter is indeed expected not to be significantly reprocessed by weak processes [e.g. [20, 108]]. Dynamical ejecta expanding close to the poles are more significantly re-processed by shock-heating and we take $Y_{\text {e}}(t=0) = 0.2$ [82]. In the case of the disc wind ejecta, since the ejection happens on the longer viscous time scale, $e^{\pm }$-captures have time to act and we consider again $Y_{\text {e}}(t=0) = 0.2$, irrespectively of the direction [176].

In Fig. 10 we present the final results of our $Y_{\text {e}}$ calculations (at 1 second after merger), for LL simulations only. In the top (bottom) panels, we collect results for the dynamical (wind) ejecta, while in the left (right) panels, along the polar axis (equatorial plane). The grey triangles mark the equilibrium $Y_{\text {e}}$, $Y_{\text {e,eq}}$, defined as the value of $Y_{\text {e}}$ obtained by assuming $\text {d}Y_e/\text {d}t=0$ in Eq. (21) and no evolution of the radius [see e.g. [154]]. In practice, it is the value of the electron fraction that the fluid element would reach if the neutrino absorption time scales were significantly smaller than the expansion time scales. We observe that $0.28 \lesssim Y_{\text {e,eq}} \lesssim 0.38$, without any clear trend with ${\tilde{\varLambda }}$. Due to the fast expansion, the final $Y_{\text {e}}$ computed by integrating Eq. (21) is smaller than $Y_{\text {e,eq}}$ for the dynamical ejecta, but comparable for disc winds, which expand more slowly. Additionally, the different flux intensities produce an appreciable difference between the polar and the equatorial directions, that it is more pronounced in the case of the dynamical ejecta. Finally, in all cases there is a weak trend both with respect to ${\tilde{\varLambda }}$ and q: the change in $Y_{\text {e}}$ is smaller for BNSs with higher tidal deformability and for more asymmetric binaries: this is consistent with the variation of the luminosities observed in Sect. 4. To asses possible systematics, we repeat our calculations using $Y_{\text {e}}(t=0)=0.15$ for all ejecta types and directions. Despite the fact that the differences in the final $Y_{\text {e}}$ decreases among the different cases, we found qualitatively similar results.

While appropriate to study general and robust trends, we stress that a detailed evolution requires to extract $Y_e$ from the simulations. We notice, however, that our results are in good qualitative agreement with simulations results. In particular, the polar irradiation is effective in increasing $Y_{\text {e}}$ in all possible configurations, due to the larger radiation flux. By comparing the calculated values of $Y_{\text {e}}$ in the different cases with the equilibrium ones, we can however conclude that the dependence on the final $Y_{\text {e}}$ on ${\tilde{\varLambda }}$ and q is rather weak.

Additionally, in Fig. 10 we highlight $Y_{\text {e}}= 0.22$, corresponding to the value of $Y_{\text {e}}$ above which the mass fraction of synthesised lanthanides and actinides drops below 10% [see e.g. [177]]. Ejecta with $Y_{\text {e}}$ above or around this value is more prone to power a blue kilonova, while for the ejecta whose electron fraction is below that value the production of lanthanides and actinides provides larger opacities to photons, resulting in a redder kilonova peaking at later times. Our results confirm previous findings: equatorial ejecta tend to produce red kilonovae in all configurations, while polar dynamical ejecta produced in equal mass and more compact mergers blue ones.

5.3 Comparison with previous results

The major outcomes of our work are in good qualitative agreement with previous works. For example, a comparison between the neutrino luminosities produced by BNS mergers with different masses and/or different EOSs in Numerical Relativity was carried out in Refs. [82, 83, 160]. The reported qualitative behaviours are similar to what we find, with peak luminosities of the order of several $10^{53}\,{\hbox {erg}~{\hbox {s}}^{-1}}$, dominant ${\bar{\nu }}_{\text {e}}$ emission, and an oscillatory phase lasting 10-15 ms post merger followed by a slow decay. As in our analysis, the softer SFHo EOS (resulting in smaller ${\tilde{\varLambda }}$’s) provides systematically larger luminosities. A relevant difference is the relative importance between $\nu _{\text {e}}$’s and $\nu _{\text {x}}$’s, whose luminosities are comparable in our simulations and in simulations from Ref. [83], while $\nu _{\text {e}}$ luminosities from Refs. [82, 160] are smaller but closer to the ${\bar{\nu }}_{\text {e}}$ ones. This difference is likely related to the different implementation details of the neutrino treatment.

Compared with Newtonian simulations, as for example the ones presented in Refs. [54, 69, 158], we see again a qualitative agreement, but some quantitative differences. In these cases, the lower neutrino luminosities were probably a consequence of the lower remnant temperature observed in the less violent merger dynamics that characterise Newtonian gravity simulations employing stiff EOSs. It is interesting to note that values of the luminosities intermediate between ours and the ones obtained in Newtonian simulations were obtained in Ref. [86], using a Smoothed Particle Hydrodynamics code with conformally flat spacetime approximation coupled with a leakage scheme. The duration of the oscillatory phase were in these models also shorter than ours, probably resulting from a possibly different post merger dynamics of the remnant. However, the hierarchy and the numerical values of the mean energies were very compatible with ours and rather independent on the BNS properties. Also the relative importance of $\nu _{\text {e}}$ and $\nu _{\text {x}}$ is closer to our results. Finally, we compare our results with the ones reported in [65] for a light BNS merger (1.2–1.2 $M_{\odot }$) employing the LS220 EOS. The rather low luminosities obtained in this case are in overall agreement with the fact that such a system is characterised by a relatively large value of ${\tilde{\varLambda }}$. $\nu _{\text {x}}$ luminosities are however more relevant here than in our results. This difference is partially explained by the larger $\langle E_{\nu _{\text {x}}} \rangle$ obtained in that analysis. A more important difference is represented by the different evolution of the luminosities with time. While also in this case one can see fast oscillations in all neutrino luminosities on the dynamical time scale, the presence of a strong peak in the very first post-merger phase is not present in these simulations. On the contrary, all luminosities tend to increase up to the end of the simulation. The relatively short duration of the simulation and the need of considering the neutrino time of flight make the comparison harder in this case.

5.4 Limitations of the present analysis

It is important to note the several limitations that affect the present analysis. First of all, since the data we work on has been generated by numerical simulations, the usual caveats that apply in this context apply in our case as well, namely the loss of accuracy due to finite resolution and the difficulty of obtaining proper convergence in the post-merger phase. In addition, a more serious limitation concerns the algorithms for neutrino transport that we rely on. They attempt to strike a balance between computational cost and physical realism, but in doing so neglect some of the finer details of neutrino dynamics. One such example is the assumption of purely radial propagation of neutrinos in the M0 scheme, which only approximately reflects the complex geometry of BNS systems. Moreover, a grey NLS is not a proper transport scheme, since it only approximates the diffusion regime through time scale arguments.

Detailed comparisons between different neutrino treatments in the context of CCSNs and BNS mergers recently addressed the problem of the accuracy of approximate neutrino transport schemes in astrophysical environments. The general outcome is that, while well gauged leakage schemes can still provide a qualitatively correct picture, the comparison with more sophisticated transport scheme revels possibly relevant differences at a quantitative level [e.g. [57, 78, 178,179,180,181,182]]. In the context of CCSNs (for which more detailed models are available and the geometry of the problem is simpler), the accuracy in the neutrino luminosity and mean energy provided by leakage schemes can be even of the order of 20–30$\%$, once directly compared with moment schemes or even with Boltzmann transport [e.g. [178, 181, 182]]. In the more complex and less studied case of BNS mergers, the discrepancy possibly increases up to a factor of a few. A recent direct comparison between the leakage+M0 scheme (used in this work) and a new M1 scheme [79], both implemented within the WhiskyTHC code and using the same progenitors and the same microphysics, revealed that the two schemes provide qualitatively similar features, but the former tends to overestimate the luminosity by a factor of $\sim$2. Additionally, the non-trivial angular dependence also introduces additional uncertainties [see, e.g., [78]]. Because of these reasons, we have decided to focus mostly on peak and integrated quantities, stressing in particular trends with respect to global binary properties and to the neutrino flavours.

Furthermore the neutrino treatment employed in this work uses what we think is the minimal set of neutrino-matter reactions necessary to account for in BNS merger scenarios, both in terms of reactions and reaction rate implementations. However a detailed analysis of the role and impact of these and other missing reactions is presently lacking. One of the main reason is that the large uncertainties that still plague neutrino transport in BNS merger simulations do not allow to robustly address this problem. In parallel to the improvement of transport schemes, it would be desirable also to improve the level of microphysics in the simulations, for example by extending the set of reactions and by implementing more detailed reactions rates and opacities, more consistent with nuclear matter properties [see e.g. [183,184,185]].

Finally, neutrinos are expected to undergo flavour conversions due to their small, but non-zero, masses. Neutrino oscillations will occur for the neutrinos emitted during and after a BNS merger. These oscillations will certainly happen in vacuum and due to matter interaction, in both cases relatively far from the merger remnant. However, collective and resonant neutrino oscillations could also happen closer to the neutrino surfaces and above the remnant [see e.g. [186,187,188,189,190]]. These effects are not included in our simulation setup, but since we are mainly interested in characterising the energy loss from the remnant this should not be a major limitation. The possible impact on the ejecta composition and on the neutrino-antineutrino annihilation is possibly more relevant and requires dedicated studies, that at the moment are usually done in a post-processing fashion [85, 86, 191].

6 Conclusions

In this paper, we have analysed the behaviour of the neutrino luminosities and mean energies produced in the merger of two NSs and during the first tens of ms after it. We have considered the outcome of 66 BNS merger simulations in Numerical Relativity, exploring 51 distinct models. The various models differ because of the masses of the binary constituents and the employed NS EOS. Each model can correspond to multiple simulations due to the inclusion of viscosity of physical origin and because of the numerical grid’s resolution. The simulation sample is homogeneous in terms of numerical setup and all simulations include neutrino emission through a grey neutrino leakage scheme [71] coupled to an M0 scheme for the propagation of neutrino radiation in optically thin conditions [105].

Despite the large variety of conditions, we found that the behaviour of the neutrino luminosities presents qualitative similarities, mostly depending on the fate of the remnant. Assuming that the central remnant does not collapse to a BH within the first 20 ms after merger (i.e. DC and LL cases in our classification), the early post merger phase is characterised by an intense neutrino emission (with total luminosities in excess of several times $10^{53} \,{\hbox {erg}~{\hbox {s}}^{-1}}$), showing a more pronounced first peak (usually occurring within the first 2-3 ms after merger), followed by ample oscillations whose period is comparable to the dynamical time scale of the merger remnant. After 10–15 ms, the oscillations subside and the luminosities enter an exponentially decreasing phase.

The bulk properties of the remnant, and in particular the matter temperature, determine the intensity of the emission. More symmetric and compact BNSs, resulting in more violent mergers and hotter remnants, have larger luminosities. The formation of shock waves produced by the bouncing central MNS and their propagation through the remnant up to the neutrino surfaces produce this characteristic peak structure. We additionally find that, unless the merger results in a PC, the neutrino luminosity correlates with the GW luminosities, since they both are enhanced by the same BNS properties.

Neutrinos come both from the cooling of the optically thick central MNS and from the innermost part of the accretion disc. The formation of a BH in the centre removes not only the MNS, but also a significant fraction of the disc. Then neutrino luminosities are significantly reduced when a BH forms. If the gravitational collapse happens within the first 5 ms (VSL simulations), only the first peak is present. If it happens promptly, i.e. without the formation of a MNS (PC simulations), only a weak and broad peak is observed. We then studied the dependence of both the peak and average luminosity (where the latter is computed over the oscillatory phase) on the reduced tidal deformability parameter ${\tilde{\varLambda }}$. We found that for equal or nearly equal BNS mergers that do not collapse too quickly to a BH (i.e. LL and DC cases), the luminosity significantly increases (up to a factor of 3) as ${\tilde{\varLambda }}$ decreases, down to ${\tilde{\varLambda }}\approx 380$. BNS mergers characterised by similar ${\tilde{\varLambda }}$, but with mass ratios significantly different from 1 produce a less intense neutrino emission, due to the less violent nature of the tidally dominated merger dynamics. PC cases populate the low-${\tilde{\varLambda }}$, low-$L_{\nu }$ portion of the result space, with an increasing trend in both quantities for more asymmetric binaries.

We additionally focused on the main luminosity peak. We found that, once the PC cases are excluded, in all cases the peak intensity anti-correlates with the peak width $\varGamma$: stronger peaks last less than weaker ones. Indeed, the energy emitted by this peak, $E_{\textrm{peak}} \sim L_{\nu ,\textrm{peak}} \varGamma$ does not show any clear trend for non-PC models.

All the trends described above apply to all modelled neutrino species, i.e. $\nu _{\text {e}}$, ${\bar{\nu }}_{\text {e}}$, and $\nu _{\text {x}}$, the latter being a collective species for heavy flavour (anti)neutrinos. Due to the neutron richness of the system and to the tendency of neutron rich matter to leptonise when decompressed and heated up, ${\bar{\nu }}_{\text {e}}$ emission dominates over $\nu _{\text {e}}$ and $\nu _{\text {x}}$, at least during the early post-merger phase. In particular, $L_{\text {peak},{\bar{\nu }}_{\text {e}}}\sim 3L_{\text {peak},\nu _{\text {e}}}$ and $L_{\text {peak},\nu _{\text {e}}}\sim L_{\text {peak},\nu _{\text {x}}}$. Similar relations hold for the luminosity averaged over the first 10 ms, even if $\langle L_{{\bar{\nu }}_{\text {e}}} \rangle / \langle L_{\nu _{\text {e}}} \rangle \sim$ 2–2.5. The reduction of the difference between the $\nu _{\text {e}}$ and ${\bar{\nu }}_{\text {e}}$ luminosities becomes more evident at later times, during the steadily declining phase, as a consequence of the early remnant leptonisation, driving its neutron-to-proton content towards a new equilibrium.

Finally, we investigated the value of the neutrino mean energies and their dependence on the BNS parameters. We found that $\langle E_{\nu } \rangle$ is the least sensitive quantity, for all neutrino flavours, with $\langle E_{\nu _{\text {e}}} \rangle \approx 10\,{\hbox {MeV}}$, $\langle E_{{\bar{\nu }}_{\text {e}}} \rangle \approx 14$–15 MeV and $\langle E_{\nu _{\text {x}}} \rangle \approx 20-25\,{\hbox {MeV}}$. This hierarchy can be easily explained in terms of the different location of the neutrino surfaces [68].

Building on a simplified, yet physically motivated model for the neutrino luminosities based on our results, we have studied the potential impact of neutrino irradiation on the electron fraction of the matter expelled from a BNS merger. Our results agree with previous findings: for example, $Y_e$ increases more significantly due to $\nu _{\text {e}}$ absorption in the polar direction and for slower disc winds. We further found that the range in luminosities has a weaker effect than the expansion time scale and the remnant geometry.

The major strengths of this work are the wide sample of models we employed and their relative homogeneity in terms of numerical setup. They qualify it as the first systematic study of the properties of the neutrino emission over a wide sample of BNS models available at present. As mentioned in Sect. 5.4, there are several areas in which our approach could be improved. Yet we believe that the results presented in this work are relevant and possibly very useful. This stems chiefly from two considerations. First of all, while obtaining more precise, accurate and realistic data is indeed desirable, it is important to start building a phenomenological and theoretical picture from the data as they are available at present. Secondly, while more realistic neutrino treatments and overall improvements in simulation machinery will undoubtedly provide quantitative corrections to the data we collected and presented here, we believe that our approach captures the fundamental aspects of neutrino emission in BNS mergers. Moreover, our characterisation of neutrino emission will likely work and find usefulness also as a reference point, to gauge the accuracy, performance and overall behaviour of the aforementioned advanced schemes.

Our analysis could also serve as input to study the detectability of neutrinos produced in a BNS merger [see e.g. [119]]. Due to their small cross sections, it will be impossible to detect thermal MeV-neutrinos produced by a merger at the typical distance of several tens of Mpc (or even more) we usually expect to observe them. However, in the very unlikely case of a Galactic BNS merger, Hyper-Kamiokande [192] will be able to detect several tens of thousands neutrinos, similar to the case of a CCSN or even larger due to the larger neutrino luminosities, especially for ${\bar{\nu }}_{\text {e}}$’s. A BNS merger occurring in the outskirt of our Galaxy (where it is more plausible to happen rather than inside the Galactic disc) will still result in a few thousands events. A handful of neutrinos could possibly be detected also if the merger happens in a nearby galaxy, up to a distance of a few times $10^3 \textrm{kpc}$. Our analysis could also be expanded towards the study of the spatial dependence of neutrino emission, as well as the its late post-merger properties. This information will be key to study, for example, the role of neutrino flavour conversions. However we leave these topics for future works.

Appendices

Appendix A: Influence of viscosity treatment

Part of our simulation sample employs an implementation of the GRLES method to effectively model the viscosity that results from the amplification of magnetic fields and the magneto-rotational instability (MRI) in the post-merger. In this appendix, we explore the impact that the inclusion or lack of viscosity has on the properties of neutrino emission. To this end Fig. 11 presents the same data that has been analysed in Secs. 4.1, but separating the simulations employing the GRLES technique from those that do not employ it.

Fig. 11

Same as the middle panel of Fig. 2, but distinguishing simulations which take into account the fluid viscosity by means of the GRLES method from those that do not

Fig. 12

Comparison of the long term behaviour of the luminosities for the three different neutrino flavours (top: $\nu _{\text {e}}$, middle: ${\bar{\nu }}_{\text {e}}$, bottom: $\nu _{\text {x}}$) for two equal mass simulations with $M_A=M_B=1.364 M_{\odot }$ employing the DD2 (orange lines) and the BLh (green lines) EOS. Solid (dashed) lines refer to the simulation with (without) GRLES viscosity

Rather interestingly, the figure highlights how the viscosity has essentially no impact on the neutrino emission, at least during the first 10 ms after merger A similar behaviour is observed also for the time average luminosities and mean energies. The explanation is related to the fact that many properties of neutrino emission are influenced primarily by the bulk dynamics of BNS systems and /or by the thermodynamic conditions of matter at the surface of neutrino decoupling. Small scale effects due to viscosity can happen on a time scale comparable to the disc dynamical time scale

$$\begin{aligned} t_{\textrm{dyn}}\sim & {} \frac{2 \pi }{\varOmega _{\textrm{K}}} \approx 0.011~\textrm{s} \left( \frac{M_{\textrm{rem}}}{2.7 M_{\odot }} \right) ^{-1/2} \left( \frac{R_{\textrm{disc}}}{100 \textrm{km}} \right) ^{3/2},\nonumber \\ \end{aligned}$$

(A.1)

(where $\varOmega _{\textrm{K}}$ is the Keplerian angular speed, $M_\textrm{rem}$ the central remnant mass, and $R_{\textrm{disc}}$ the typical disc extension) only on length scales comparable or smaller than the mixing length. Since the GRLES model was calibrated on MRI data for which the mixing length is $\ell _{\textrm{mix}} \lesssim 25\,{\hbox {m}}$, and only in a narrow density band (for both higher and lower densities, it decreases rapidly [142]), the inclusion of viscous effects have almost no impact on the bulk motion inside the remnant ($R_{\textrm{rem}} \sim 15\,{\hbox {km}} \gg \ell _{\textrm{mix}}$) and inside the disc ($R_{\textrm{disc}} \sim 100\,{\hbox {km}} \gg \ell _{\textrm{mix}}$) during the first ms after merger.

Fig. 13

Resolution comparison for the neutrino luminosity evolution of the three neutrino flavours ($\nu _{\text {e}}$, ${\bar{\nu }}_{\text {e}}$ and $\nu _{\text {x}}$ in panels a, b and c, respectively) and for the total neutrino luminosity (panel d). The model is a LL simulation of an equal mass system ($M_A=M_B=1.364\,{{\hbox {M}}_{\odot }}$) employing the BLh EOS and GRLES viscosity. Filled regions cover the range containing all three resolutions

On the other hand, viscosity induces matter accretion on the longer viscous time scale. The latter can be estimated as $t_{\textrm{vis}} \sim \nu _{\textrm{T}}/R_{\textrm{disc}}$ where $\nu _{\textrm{T}}$ is the viscosity coefficient. For a Keplerian, Shakura-Sunyaev disc [193] whose viscosity coefficient is parametrised in terms of a dimensionless $\alpha$ parameter,

$$\begin{aligned} t_{\textrm{vis}}\sim & {} \alpha ^{-1} \left( \frac{H}{R} \right) ^{-1} \varOmega _{\textrm{K}}^{-1} \approx 12.8~\textrm{s} \left( \frac{\alpha }{5\cdot 10^{-4}} \right) ^{-1} \nonumber \\{} & {} \times \left( \frac{H/R}{1/2} \right) ^{-1} \left( \frac{M_{\textrm{rem}}}{2.7 M_{\odot }} \right) ^{-1/2} \left( \frac{R_{\textrm{disc}}}{100\,{\hbox {km}}} \right) ^{3/2},\quad \end{aligned}$$

(A.2)

where H/R is the disc aspect ratio. Note that in this formula we used the estimate $\alpha = \ell _{\textrm{mix}} H$ obtained by considering a Shakura-Sunyaev disc with $\ell _\textrm{mix}=25\,{\hbox {m}}$. However BNS discs are not thin as in the Shakura-Sunyaev model, so this is only a qualitative estimate. More realistic values inside the disc are $\alpha \sim 0.01$, see [143], and the resulting accretion time scale are $\sim {\mathcal {O}}(1~\textrm{s})$. On such a time scale, simulations employing physical viscosity provide larger neutrino luminosities (especially for $\nu _{\text {e}}$’s and ${\bar{\nu }}_{\text {e}}$’s) due to the enhanced accretion rate, as visible in Fig. 12.

Appendix B: Resolution dependence

The simulations used in this study make use of a box-in-box AMR grid with three possible resolutions, namely: LR, SR and HR. Here we study the effect of different resolutions on the neutrino emission by considering one model for which all the three resolutions are available as representative. In each panel of Fig. 13 we present the luminosities obtained by the different simulations for each of the three neutrino species, alongside their sum (bottom right panel). The shaded area represents the maximum variability between resolutions.

Table 2 List of PC and VSL simulations. Columns from left to right provide: the mass of the heaviest star; the mass of the lightest star; the mass ratio; the reduced dimensionless tidal deformability; whether the simulation employs GRLES viscosity; the resolution; the time of collapse (and the time of the end of the simulation); the peal luminosity for the 3 neutrino flavors; the work in which the simulation was first presented

On one hand, within the first $\sim 10\,{\hbox {ms}}$ after merger, the three resolutions can differ up to a factor $\sim 2$ at corresponding times. This is due to the fact that the neutrino luminosity oscillates very rapidly and widely, as a consequence of the complex remnant’s dynamics. Clearly point-wise differences become comparable to the oscillation amplitude as soon as the remnant’s dynamics (characterised by bounces and sound/shock waves, and strongly dependent on the resolution) accumulates a difference comparable to the luminosity oscillation periods. Enlarging our view on the whole 0-10 ms interval, we recognise that different resolutions produce a very similar global behaviour, even if the HR simulation tends to have higher maximum peaks and lower secondary peaks with respect to the LR one. On the other hand, in the exponentially decreasing phase (after the $10{-}15\,{\hbox {ms}}$ mark), the three resolutions show a much closer behaviour.

Table 3 List of DC and LL simulations. Columns from left to right provide: the mass of the heaviest star; the mass of the lightest star; the mass ratio; the reduced dimensionless tidal deformability; whether the simulation employs GRLES viscosity; the resolution; the time of collapse (and the time of the end of the simulation); the peal luminosity for the 3 neutrino flavors; the work in which the simulation was first presented

Finally resolution mainly affects the luminosity peak, while the average luminosities and mean energies are less influenced. In order to quantify their variations, we average the values of $L_{\text {peak}}$, $\langle L \rangle$ and $\langle E\rangle$ over the three available resolutions and consider the maximum relative deviations from these values. While peak luminosities can vary by up to 20% from their average value, for average luminosities and average mean energies this figure is reduced to 15% and 8%, respectively.

This analysis suggests that while the precise values we quote in our results could of course be improved if we had access to higher-resolution simulations, the trends we identify are robust and valid.

Table 4 List of DC and LL simulations. Columns from left to right provide: the mass of the heaviest star; the mass of the lightest star; the reduced dimensionless tidal deformability; whether the simulation employs GRLES viscosity; the resolution; the average luminosity for the 3 neutrino flavors and respective standard deviations; the average mean energy for the 3 neutrino flavors and respective standard deviations

Appendix C: Data tables

In this section we collect detailed data pertaining to our simulation sample and our results. Table 2 lists PC and VSL simulations, providing details on the initial conditions, EOS and the value of peak luminosities. The same data is provided for DC and LL simulations in Table 3. In Table 4, we list the values of average luminosities and average mean energies, and their standard deviations for DC and LL simulations.