Calibration of the Advanced Spectral Leakage scheme for neutron star merger simulations, and extension to smoothed-particle hydrodynamics Free

ABSTRACT

1 INTRODUCTION

The combined detection of gravitational and electromagnetic waves from the neutron star merger event GW170817 (Abbott et al. 2017a, c) has solved several decades-old puzzles. For example, the detection of the short Gamma-Ray Burst (sGRB) GRB170817A ∼1.74 s after the merger (Goldstein et al. 2017; Savchenko et al. 2017) confirmed the long suspected association between sGRBs and compact binary mergers (Paczynski 1986; Eichler et al. 1989; Narayan, Paczynski & Piran 1992). Moreover, the detection of the macronova/kilonova transient AT2017gfo (Arcavi et al. 2017; Chornock et al. 2017; Coulter et al. 2017; Evans et al. 2017; Kasliwal et al. 2017; Kilpatrick et al. 2017; Lipunov et al. 2017; Pian et al. 2017; Smartt et al. 2017; Soares-Santos et al. 2017) confirmed binary neutron star mergers as a major, and possibly dominant, r-process nucleosynthesis site (Drout et al. 2017; Kasen et al. 2017; Rosswog et al. 2018), thus confirming earlier theoretical predictions on the subject (Lattimer & Schramm 1974; Eichler et al. 1989; Freiburghaus, Rosswog & Thielemann 1999; Rosswog et al. 1999). The event GW170817 also allowed to place important constraints on the Equation of State (EoS) of nuclear matter (Bauswein et al. 2017; Abbott et al. 2018; De et al. 2018; Most et al. 2018; Radice et al. 2018a; Coughlin et al. 2019; Jiang et al. 2019, 2020; Kiuchi et al. 2019; Radice & Dai 2019) and provided an independent measurement of the Hubble parameter (Abbott et al. 2017b; Dhawan et al. 2020).

The luminous blue component of the macronova transient AT2017gfo has also emphasized the role played by neutrinos in changing the electron fraction Y_e of the ejected matter. Weak interactions modify the electron fraction in a fair portion of the ejecta from initial values of Y_e ∼ 0.05 to values above Y_e = 0.25, where the resulting nucleosynthesis changes abruptly (Korobkin et al. 2012; Lippuner & Roberts 2015) and no more lanthanides and heavier nuclei are produced. This increase in Y_e also drastically reduces the optical opacities (Kasen, Badnell & Barnes 2013; Tanaka & Hotokezaka 2013) and leads to blue (rather than red) transients after about a day rather than a week. Fits to the blue component of the light curve from macronova models suggest ejecta of several 10⁻² M_⊙ with velocity ∼0.3c. There is a general agreement that secular mass ejection from the remnant disc is needed to achieve such large ejecta amounts (e.g. Radice et al. 2018b; Ciolfi & Kalinani 2020; Nedora et al. 2020), but it is still a matter of debate to which extent the different proposed secular ejection channels contribute to the blue macronova. The ejection channels include winds driven by magnetic fields (Ciolfi et al. 2017; Ciolfi & Kalinani 2020), viscosity (Fernández & Metzger 2013; Siegel & Metzger 2017; Fujibayashi et al. 2018, 2020; Fernández et al. 2019), neutrinos (Dessart et al. 2009; Perego et al. 2014a; Martin et al. 2015), and spiral-waves (Nedora et al. 2019, 2020). Our motivation here are neutrino-driven winds, for which no systematic study in dynamical simulations of binary neutron star mergers exists to date. For this reason, we have recently implemented an extension to the original Advanced Spectral Leakage (ASL) scheme (Perego, Cabezón & Käppeli 2016) with the purpose of modelling the neutrino absorption responsible for launching the winds in merger simulations (Gizzi et al. 2019).

In this paper, we extend the analysis of the ASL scheme by performing a careful calibration for the three parameters that enter the scheme. These parameters have an impact on both neutrino emission and absorption, and must therefore be gauged carefully. Two of these parameters were already introduced in Perego et al. (2016), but calibrated for core-collapse supernovae simulations. They define the neutrino emission and represent two physical effects. The first is a blocking parameter, accounting for both Pauli blocking effects due to the fermionic nature of neutrinos and for inward neutrino fluxes when calculating the neutrino emission above the decoupling region, while the second is a thermalization parameter, describing the number of weak interactions needed to thermalize neutrinos. We added a third parameter in Gizzi et al. (2019) specifically for modelling the neutrino absorption responsible for launching the neutrino-driven winds, and we gauged it by comparing the outcome of the simulations with the two-moment scheme (M1) approach implemented in flash (Fryxell et al. 2000; O’Connor 2015; O’Connor & Couch 2018). Here, we use both the M1 and the Monte Carlo neutrino transport code sedonu (Richers et al. 2015) as sources of comparison. We use neutrino quantities directly impacted by each parameter to perform our calibration and we extract them by post-processing a number of snapshots of binary neutron star remnants.

When using the ASL scheme, we present a new, mesh-free implementation for calculating spectral, multigroup optical depths for smoothed-particle hydrodynamics (SPH) (Monaghan 1992; Monaghan 2005; Rosswog 2009, 2015a, b, 2020) simulations. This makes our ASL fully mesh-free and ideal for dynamical SPH simulations of merging neutron stars with neutrino effects.

The paper is structured as follows. We begin with a description of the neutrino transport codes in Section 2. We then describe a new particle-based algorithm to calculate the optical depth in Section 3, where we also include a comparison with grid-based calculations. Section 4 describes the set of snapshots that we use, as well as the calibration strategy. The results of the calibration are discussed in Section 5. In Section 6, we summarize and draw our conclusions.

2 NEUTRINO TRANSPORT METHODS

The neutrino transport approaches we use for our analysis are: the Monte Carlo code sedonu, the two-moment scheme (M1) developed in the Eulerian hydrodynamics code flash, and the ASL scheme. All of them are spectral (i.e. neutrino energy is treated as an independent variable and neutrino–matter interactions are calculated for different neutrino energies), an important ingredient given the energy-squared dependence entering neutrino cross-sections (see e.g. Burrows & Thompson 2004 and references therein). Moreover, we choose the same set of weak interactions: production and absorption of electron neutrinos and antineutrinos via charged current processes involving nucleons and nuclei, neutrino emission by bremsstrahlung and pair processes, and finally elastic scattering off nucleons and nuclei. We choose the SFHo EoS (Steiner, Hempel & Fischer 2013) to calculate spectral emissivities and opacities. sedonu and M1 interpolate them from the NuLib tables (O’Connor 2015). The ASL scheme calculates them following Bruenn (1985), Mezzacappa & Messer (1999), and Hannestad & Raffelt (1998). Three distinct neutrino species are modelled: electron neutrinos ν_e, electron antineutrinos |$\bar{\nu }_\mathrm{ e}$|⁠, and a collective species ν_x for heavy-lepton neutrinos, comprising μ and τ neutrinos and antineutrinos.

2.1 Monte Carlo

Monte Carlo methods can be used to solve the Boltzmann equation. They provide an exact solution to the transport problem in the limit of an infinite particle number. However, Monte Carlo neutrino transport methods tend to be computationally expensive, since a large number of particles is needed to keep the numerical noise in the solution at a low level. These methods have therefore mostly been used on static, spherically symmetric core-collapse supernovae snapshots (Janka & Hillebrandt 1989a, b; Janka 1992; Keil, Raffelt & Janka 2003; Abdikamalov et al. 2012). Similarly, studies of Monte Carlo neutrino transport in neutron star mergers have been carried out by post-processing snapshots, assuming the neutrino field can be approximated as steady-state (Richers et al. 2015). Recently, Foucart et al. (2020) performed the first dynamical hydrodynamic simulation of merging neutron stars up to ∼5 ms post-merger with a computationally particularly efficient Monte Carlo approach. Here, we use the steady-state Monte Carlo code sedonu (Richers et al. 2015) on axisymmetric fluid snapshots to calibrate the parameters of our ASL scheme.

We use a logarithmically spaced energy grid that spans the range E_ν ∈ [0.5, 200] MeV. Neutrinos are simulated as packets passing through and interacting with a static fluid background imported from hydrodynamic simulations. Each packet is specified by the neutrino energy E_ν, the location x, the unit vector |$\hat{{\bf d}}$| describing the propagation direction, and the total number of neutrinos in the packet. A fixed number of neutrino packets are emitted from each grid cell. In particular, sedonu randomly samples the energy-dependent emissivity to determine the neutrino energy, and the packet direction of propagation is drawn from an isotropic distribution. Neutrino packets interact with the fluid by being partially absorbed during the packet propagation, and changing propagation direction when scattered. The distance the packet travels before interacting is randomly drawn from an exponential distribution that depends on the scattering mean free path. If the packet is scattered, it is given an isotropically random new direction. While moving, the packet deposits both energy and lepton number to the traversed grid-cells, from which the corresponding rates of change can be derived. We emit 100 Monte Carlo packets per energy group and grid cell. In this way, integrated quantities such as luminosities and mean neutrino energies are stable in spite of Monte Carlo noise, which is estimated to be a thousand times smaller.

2.2 Two-moment scheme (M1)

A moment scheme is an approximation to the Boltzmann neutrino transport, obtained by evaluating angular integrals of the Boltzmann equation in order to time-evolve the angular moments of the distribution function (Lindquist 1966; Bruenn, Buchler & Yueh 1978; Bruenn 1985; Mezzacappa & Messer 1999). The M1 scheme is a multidimensional, spectral approach where only the zeroth and the first moments of the neutrino distribution function are evolved, which, respectively, describe the spectral energy density and the energy flux density (Castor 2004). The system of equations is closed by an assumed closure relation, which is often analytical (Minerbo 1978; Levermore & Pomraning 1981; Smit, Cernohorsky & Dullemond 1997; Pons, Ibáñez & Miralles 2000; Audit et al. 2002). Generally, the analytical form captures the physics of the transport well in specific regimes, while approximating the exact solution in others. Neutrino–matter interactions are included in appropriate source terms that appear on the right-hand side of the moment equations and which are solved via techniques that are borrowed from finite volume hydrodynamics.

The M1 scheme has been widely applied in simulations that model shock revival in core-collapse supernovae (O’Connor & Ott 2013; Obergaulinger, Janka & Aloy 2014; O’Connor 2015; O’Connor & Couch 2018; Skinner et al. 2019) and recently compared to other transport approaches (Cabezón et al. 2018; O’Connor et al. 2018; Pan et al. 2018), showing good performance and agreement with alternative schemes at the level of |$10\!-\!20{{\ \rm per\ cent}}$| in luminosities, mean energies, and shock radii. However, not so much has been done in the context of binary neutron star mergers. In the work of Foucart et al. (2016), an M1 scheme has been adopted to test the impact of the energy spectrum on the composition of merger outflows. In a subsequent work, Foucart et al. (2018) compared a grey M1 scheme against a Monte Carlo approach, showing the limitations of the M1 analytical closure in recovering the properties of the polar ejecta and therefore in modelling neutrino-driven winds and kilonovae. Another quantitative study of the assumptions employed in analytical closures and their violation has been carried out more recently by post-processing snapshots with Monte Carlo transport (Richers 2020). Despite their limitations, Foucart et al. (2020) have shown that time-dependent neutrino transport in an M1 scheme with an analytical closure that carefully treats the energy spectrum approximates the Monte Carlo solution within |$\sim 10\!-\!20{{\ \rm per\ cent}}$| in the composition and velocity of ejecta, as well as electron neutrino and antineutrino luminosities and mean energies. In Section 5, we explore the performance of the M1 scheme implemented in flash and compare it with our new ASL treatment and the Monte Carlo code sedonu.

2.3 Advanced Spectral Leakage

We start by briefly summarizing the original implementation of the ASL scheme in its main aspects. For more details, the reader is referred to Perego et al. (2014a, 2016). Subsequently, we describe the implementation of Gizzi et al. (2019).

1. The equilibrium-diffusive regime, where |$\tau _{\nu ,\rm {tot}} \gg 1$| and |$\tau _{\nu ,\rm {en}} \gtrsim 1$|⁠. Neutrinos can be treated as a trapped Fermi gas in thermal and weak equilibrium with the matter.

2. The diffusive regime, where |$\tau _{\nu ,\rm {tot}} \gg 1$| and |$\tau _{\nu ,\rm {en}} \lesssim 1$|⁠. Neutrinos still diffuse, but they are not necessarily in thermal equilibrium with the matter.

3. The semitransparent regime, where |$\tau _{\nu ,\rm {tot}} \sim 1$|⁠. Although neutrinos can still be partially absorbed, they begin to decouple from matter, and the surfaces around which this occur are called neutrino-surfaces, which can be identified when τ_ν ∼ 2/3. Low-energy neutrinos decouple earlier than high-energy neutrinos, and therefore the neutrino surfaces for the latter are wider (Endrizzi et al. 2020). In the context of binary neutron star mergers, neutrino absorption occurring in this regime (hereafter referred to as heating) drives mass ejection via so-called neutrino-driven winds (Dessart et al. 2009; Perego et al. 2014a).

4. In the free-streaming regime, where |$\tau _{\nu ,\rm {tot}} \lesssim 1$|⁠, locally produced neutrinos stream out freely, and neutrino absorption is negligible.

1. a blocking parameter, |$\alpha _{\nu ,\rm {blk}}$|⁠, that accounts for Pauli blocking effects due to the large amount of neutrinos locally produced or emitted at the neutrino surface, and for the fact that the interpolation given by equation (3) favours the production rate in the free-streaming regime. Neutrino emission in this regime is assumed to be isotropic, and a fraction of neutrinos is therefore emitted inward and does not contribute to the luminosity.

2. a thermalization parameter, τ_cut, representing the amount of interactions needed to thermalize neutrinos, and which accounts for neutrinos that exchange energy with the fluid while propagating from the innermost to the outermost regions of the remnant, making the neutrino spectrum at the decoupling surface softer.

The computation of the heating rate requires the knowledge of the distribution of the neutrino fluxes coming from the region inside the neutrino surfaces. While the original implementation of the ASL assumes spherically symmetric neutrino fluxes (Perego et al. 2016), the more complex geometry of a binary neutron star merger causes anisotropic neutrino fluxes (Rosswog & Liebendörfer 2003; Dessart et al. 2009; Perego et al. 2014a; Foucart et al. 2016). We have recently implemented an extension of the original ASL scheme that allows the modelling of neutrino-driven winds (Gizzi et al. 2019), which we briefly recap here.

The parameter p has been previously calibrated (Gizzi et al. 2019) by comparing the heating rate maps of electron neutrinos and antineutrinos¹ for a binary neutron star merger snapshot against those obtained with the M1 scheme in flash (Fryxell et al. 2000; O’Connor & Couch 2018). This demonstrated that the new ASL implementation is able to capture the heating rate distribution to better than a factor of 2, with the largest deviations right above the central remnant. In this work, we further improve the analysis of the new ASL scheme and we provide values for |$\alpha _{\nu ,\rm {blk}}$|⁠, τ_cut, and p that are more accurate for neutron star mergers. All the details are described in Section  4.

3 PARTICLE-BASED OPTICAL DEPTH

SPH is a mesh-free, particle-based method to solve the ideal hydrodynamics equations, see Monaghan (2005) and Rosswog (2009, 2015a, 2015b, 2020) for reviews of and recent developments within this method.

1. The SPH particle for which we search the optical depth is labelled a. For the convenience of the subsequent notation we also identify particle a as the ‘zeroth neighbour particle’ and label it b₀.

2. The first task is to find the particle b₁ in a’s neighbour list that has the minimum density among the neighbours and add the contribution, |$\tau _{b_0 \rightarrow b_1}$| to the optical depth of a.

3. Once b₁ has been found, find particle b₂ in b₁’s neighbour list that has the minimum density and so on until a particle |$b_{k_\infty }$| has been found from which neutrinos can escape freely. We use as escape condition ρ ≲ 10⁹g cm⁻³ (Endrizzi et al. 2020).

Unlike in Gizzi et al. (2019) and in most of the literature (see Perego et al. 2014a, George et al. 2020, and Endrizzi et al. 2020 as examples), where the standard approach is to pre-select radial paths on a grid, calculate the integral of equation (1) over these paths, and then take the minimum optical depth, here there is no assumption on the type of path, but it is the matter itself that tells neutrinos how to leak out (see also Perego et al. 2014b for a similar methodology). Fig. 1 shows sketches of both the particle- (top) and grid-based (bottom) methods. The main advantage of our new approach is the affordable computational cost required to get spectral optical depths compared to grid-based integration methods. This is on one side due to the fact that we use directly the particle properties for the calculations, and on the other side because we adopt the fast recursive coordinate bisection tree of Gafton & Rosswog (2011) for neighbour search. Figs 2–3 show on the plane x-z the total and energy optical depths respectively, for case (1) described in Section 4.1. Both of them are calculated on a grid with the algorithm of Gizzi et al. (2019), and with our particle-based approach. We SPH-map on the plane x-z the particle-based optical depths for comparison purposes. We show the case for a low-energy bin (∼10 MeV, upper panels) and for a high-energy bin (∼100 MeV, lower panels) of our energy grid, and for electron neutrinos (top), electron antineutrinos (middle), and heavy-lepton neutrinos (bottom). We also show the map of the relative difference ϵ_τ = (τ_grid − τ_part)/τ_grid. The optical depth appears to be larger along the disc with respect to the grid-based approach, and lower along the poles. Low-density regions in SPH are not as well resolved as the high-density ones, therefore the SPH maps in the former regime are only indicative of what the optical depth could actually be (although we expect it to be rather small anyway). Note that τ_grid is calculated once initial conditions of density, temperature, and electron fraction of the SPH snapshot are mapped on the grid. Therefore, τ_grid suffers from the same uncertainties of τ_part in low-density regions. The higher optical depths along the disc with the particle-based approach leads to more extended neutrino surfaces at high neutrino energies, especially for electron neutrinos, which are the most interacting species given the neutron richness of the material. Nevertheless, our algorithm well captures the expected distribution of the neutrino surfaces according to different neutrino energies and species (see Gizzi et al. 2019 and Perego et al. 2014a for details). Accounting for different resolutions that can be used both in the particle- and grid-based approaches, we estimate the particle-based algorithm to be ≳100 times faster than the grid-based one.

Figure 1.

Simplified sketches of the particle-based (top) and grid-based (bottom) approaches for the optical depth calculation. Shown in grey is a neutron star remnant surrounded by a disc. In the particle-based model, each particle has a smoothing length defining the radius of each circle. The latter encloses the particle neighbours. Unlike the grid-based approach, where outgoing radial paths are a priori selected for calculating the optical depth, in the particle-based case the path along which neutrinos move is defined by the density distribution, and therefore it can be non-straight.

Open in new tabDownload slide

Figure 2.

Total optical depth on the x-z plane for electron neutrinos (top), electron antineutrinos (middle), and heavy-lepton neutrinos (bottom). For each species we show the map for both a low energy (∼10 MeV, top) and a high energy (∼100 MeV, bottom) bin of the spectrum. Moreover, from left to right, we show the maps obtained from the grid-based and particle-based calculations, and a map of the relative difference |$\epsilon _{\tau }= \frac{\tau _{\rm {grid}}-\tau _{\rm {part}}}{\tau _{\rm {grid}}}$|⁠, respectively.

Open in new tabDownload slide

Figure 3.

Same as Fig. 2, but for the energy optical depth.

Open in new tabDownload slide

4 CALIBRATION

4.1 Set-up

We take three snapshots of binary neutron star merger remnants from the SPH simulations of Rosswog et al. (2017), each one representing a specific case for calibrating our scheme:

case (1): an equal mass binary system of 1.4–1.4 M_⊙ at ≈38 ms after merger. This is our late-time configuration, when the matter around the remnant has the highest degree of axial symmetry.

case (2): an equal mass binary system of 1.3–1.3 M_⊙ at ≈18 ms after merger. In this way we include equal mass configurations at earlier post-merger times, when the matter distribution is not fully axially symmetric.

case (3): an unequal mass binary system of 1.2–1.3 M_⊙ at ≈18 ms after merger. We consider this case as representative of unequal mass binaries at early times, with the lowest degree of axial symmetry.

We use cases (2) and (3) to test the assumption of axial symmetry entering equation (5). Table 1 summarizes the properties of each binary configuration.

Since the snapshots are taken from simulations that use a grey leakage scheme (Rosswog & Liebendörfer 2003), we first map the SPH properties on a 3D grid and let the radiation field equilibrate with M1 for 3 ms without changing either the electron fraction or the temperature. Afterwards, we let the electron fraction evolve for 10 ms with M1 to remove grey leakage effects. We do not consider the temperature evolution in M1 since we do not see important differences between the corresponding evolved and non-evolved profiles. The evolution with M1 allows to make a more consistent comparison between neutrino approaches for our parameter calibration. Figs 4–5 show density, temperature, and electron fraction in the equatorial plane and on the plane orthogonal to it after the 10 ms of evolution, and for each case. In order to considerably reduce computational costs when running the Monte Carlo code sedonu we map density, temperature, and the evolved electron fraction from M1 on to a 2D, axially symmetric grid, and run the transport over the obtained profiles. At last, we map the evolved M1 data back to the SPH particles by tri-linear interpolation to run the ASL scheme.

Figure 4.

Density (left column), temperature (middle column), and electron fraction (right column) maps on the equatorial plane at equilibrium, for a 1.4–1.4 M_⊙ binary neutron star merger at ≈38 ms after merger (top row), a 1.3–1.3 M_⊙ binary neutron star merger at ≈18 ms after merger (middle row), and a 1.2–1.3 M_⊙ binary neutron star merger at ≈18 ms after merger (bottom row). From top to bottom the degree of axial symmetry of the remnant decreases. Snapshots are taken from the dynamical simulations of Rosswog et al. (2017).

Open in new tabDownload slide

Figure 5.

Same as in Fig. 4, but on the plane orthogonal to the equatorial one.

Open in new tabDownload slide

4.2 Strategy

As a first step we identify physical quantities that are directly impacted by each ASL parameter:

1. the total luminosities and mean energies, primarily affected by |$\alpha _{\nu ,\rm {blk}}$| and τ_cut.

2. the neutrino flux at different polar angles that an observer located far outside the neutrino surfaces would measure. In this respect, we focus on the trend of the flux with the polar angle for constraining the parameter p in equation (5).

Given that the Monte Carlo approach converges (in the limit of infinite particle numbers) to an exact solution of the transport equation, the best strategy for our parameter calibration would be to extract both quantities (1) and (2) from sedonu and compare with the ASL scheme. However, the assumption of axial symmetry may not be fully appropriate in all cases, especially shortly after merger. We then first test the assumption of azimuthal symmetry by running the M1 transport on both 2D and 3D grid set-ups, and compare the neutrino emission maps for each species. For those snapshots where the impact on the emission is large in at least one of the species, we decide to consider M1 in 3D for the calibration. However, given the limitations of M1 in accurately modelling the distribution of the neutrino flux at different polar angles (Foucart et al. 2018), we make use of M1 purely for calibrating |$\alpha _{\nu ,\rm {blk}}$| and τ_cut, while we keep sedonu to calibrate the heating parameter p. As we shall see, the 2D assumption does not impact the trend of the flux at different polar angles, but only its local values. On the other hand, if the 2D assumption does not affect the emission of any neutrino species sensitively, the calibration of |$\alpha _{\nu ,\rm {blk}}$|⁠, τ_cut, and of the heating parameter p is performed entirely with sedonu.

The parameter space we explore is given by |$\alpha _{\nu ,\rm {blk}} \in [0.35,0.45,0.55,0.65,0.75,0.85]$|⁠, τ_cut ∈ [5, 10, 15, 20], and p ∈ [2, 4, 6, 8, 10, 12]. When calibrating, we assume |$\alpha _{\rm {blk}} \equiv \alpha _{\nu _e,\rm {blk}}= \alpha _{\bar{\nu }_e,\rm {blk}}$|⁠, and |$\alpha _{\nu _x,\rm {blk}}= 0$|⁠, similarly to Perego et al. (2016). We will discuss the accuracy of the assumption |$\alpha _{\nu _e,\rm {blk}}= \alpha _{\bar{\nu }_e,\rm {blk}}$| in Section 5.2. Given that heavy-lepton neutrinos do not contribute to the heating either, we decide to neglect this species entirely when calibrating, and just focus on the other two species. We will anyway report all values of luminosities and mean energies for completeness (see Table 2).

We summarize the values of luminosities and mean energies for the three cases described in Section 4.1 in Table 2. We include the 2D values from sedonu, the 2D and 3D values from M1, and the values from the ASL with the corresponding best parameter set when assuming |$\alpha _{\rm {blk}} \equiv \alpha _{\nu _e,\rm {blk}}= \alpha _{\bar{\nu }_e,\rm {blk}}$|⁠.

5 RESULTS

In Section 5.1, we describe the calibration of each ASL parameter based on a separate analysis for each snapshot. Afterwards, we combine the results in Section 5.2 and discuss the performance of using the same blocking parameter for electron neutrinos and antineutrinos.

5.1 Parameter constraints

5.1.1 Blocking

From Table 2 we can see that by comparing the 2D and 3D luminosities of electron neutrinos and antineutrinos from M1 the impact of the 2D averaging is only of the order of |$\lesssim 7{{\ \rm per\ cent}}$| for both species in case (1). This suggests the usage ofsedonu for the comparison with the ASL scheme. On the other hand, for cases (2) and (3) the 2D averaging implies a reduction in the electron antineutrino luminosity by about a factor of 2 and 10, respectively. We examine the effect of the 2D averaging on the emission maps for each species, and we indeed notice that the 2D assumption affects the neutrino species to a different extent, depending on the respective location of the bulk of the emission. We therefore decide to take the 3D M1 values of luminosities and mean energies as reference for calibrating α_blk and τ_cut in cases (2) and (3). We also notice that the heavy-lepton neutrino luminosity from M1 in 3D and from the ASL is systematically larger by a factor of a few than the one from sedonu in all snapshots. This is different from the results of Foucart et al. (2020). However, the fact that only heavy-lepton neutrinos show large deviations with respect to a Monte Carlo approach is an indication that the cause might be associated with some ingredient in the particular transport approach adopted, and for which heavy-leptons are more affected than the other species. Regarding the ASL, a likely explanation is the treatment of the emission rate, which is calculated as smooth interpolation between production and diffusion rates as shown in equation (2). The diffusion rate depends on the diffusion time-scale, which is estimated via a random-walk argument. As pointed out in Ardevol-Pulpillo et al. (2019), this derivation leads to a steeper decrease of the diffusion time-scale with radius. Combined with the fact that most of the neutrinos escape around the rather small region of the neutrino surface, a lower diffusion time-scale in this region can boost the emission up to more than a factor of 2, depending on the species. For our binary neutron star configurations, this is particularly true for heavy-lepton neutrinos because their sources of production are just pair processes and bremsstrahlung, which are both extremely temperature-dependent. Combined with the fact that heavy-lepton neutrinos decouple at inner and still rather hot regions with respect to the other two species, it is likely that their emission is more affected by the treatment of the diffusion.

Fig. 6 shows line plots of ϵ_L and ϵ_〈E〉 as a function of α_blk, for different values of p and τ_cut. Lines of the same colour are for fixed p, while lines of the same type are for fixed τ_cut. We neglect the cases with τ_cut = 5 and p = 12 to reduce the amount of data to show, but the results do not change. From the top to the bottom row we show the cases from (1) to (3), respectively. Generally, we see that neither ϵ_L nor ϵ_〈E〉 is largely affected by varying p for a given α_blk and τ_cut. The mild dependence can be explained by the fact that p defines the distribution of the heating rather than its intensity. However, because the neutrino absorption is typically more pronounced for electron neutrinos than for antineutrinos, we should expect a somewhat larger dependence of ϵ_L on p for the former species. Nevertheless, since we are evaluating the species-summed ϵ_L, even in case there is such dependence from the electron neutrinos, this must be associated with a rather small |$\epsilon _{L_{\nu _e}}$| (see later Fig. 7), and it is thus not appreciable. On the other hand, the blocking parameter α_blk has a major impact on ϵ_L. This comes directly from the r_ν → (1 − α_blk)r_ν correction to the emission rate r_ν when accounting for blocking, see equation (2). This is different from the impact on ϵ_E, which is basically negligible. The reason is the fact that 〈E_i〉 is computed from the ratio between the luminosity L_i [erg s⁻¹] and the total number of emitted neutrinos of species i per unit time, both affected by blocking for the case of electron neutrinos and antineutrinos. At last, a noticeable dependence to τ_cut can be seen from ϵ_〈E〉, as a consequence of the fact that τ_cut impacts the neutrino spectrum at the decoupling region, and therefore the mean energies of the species. We find that ϵ_〈E〉 has a maximum value of |$\sim 18{{\ \rm per\ cent}}$| and |$\sim 12{{\ \rm per\ cent}}$| at τ_cut = 5 for cases (1) and (2), respectively. Case (3) shows larger values, but limited to |$\lt 30{{\ \rm per\ cent}}$|⁠. On the other hand, ϵ_L varies from |$\sim 10{{\ \rm per\ cent}}$| to more than |$100{{\ \rm per\ cent}}$| for cases (1) and (2), while case (3) shows values ϵ_L ≳ 0.75. The large ϵ_L for any parameter combination for the latter case is a consequence of the |$\alpha _{\rm {blk}}= \alpha _{\nu _e,\rm {blk}}= \alpha _{\bar{\nu }_e,\rm {blk}}$| assumption, which makes it cumbersome to always well catch the luminosities of both electron neutrinos and antineutrinos (see also Fig. 7 and next paragraph). Since ϵ_L is more sensitive than ϵ_〈E〉 to a change in α_blk, we just look at ϵ_L for a first parameter pre-selection. In particular, we find a minimum ϵ_L around α_blk = 0.45 for case (1), α_blk = 0.65 for case (2), and α_blk = 0.75 for case (3).

Figure 6.

Relative difference on the luminosity ϵ_L (left-hand panels) and the mean neutrino energy ϵ_〈E〉 (right-hand panels) as a function of α_blk and for different values of p and τ_cut, calculated by comparing the results from the ASL with the reference solution. From top to bottom, we show cases from (1) to (3). We use sedonu as reference solution for case (1), and M1 for cases (2) and (3). Lines of the same colour are for fixed p, while lines of the same type are for fixed τ_cut. We neglect the cases with τ_cut = 5 to reduce the amount of data to show. The parameter α_blk has a large impact on ϵ_L, but not on ϵ_〈E〉. The thermalization parameter τ_cut has some effect on ϵ_〈E〉, while p mildly impacts both ϵ_L and ϵ_〈E〉.

Open in new tabDownload slide

Figure 7.

Relative difference on the luminosity |$\epsilon _{L_i}$| and the mean energy |$\epsilon _{\langle E_i \rangle }$| for i = ν_e (left-hand panels) and |$i=\bar{\nu }_e$| (right-hand panels), and for case (1) (first and second row), case (2) (third and fourth row), and case (3) (fifth and sixth row), calculated by comparing the results from the ASL with the reference solution. We use sedonu as reference solution for case (1), and M1 for cases (2) and (3). While |$\epsilon _{\langle E_{i} \rangle }$| is always limited to |$\lesssim 23{{\ \rm per\ cent}}$| at the most for electron antineutrinos in case (3), the assumption of a single blocking parameter for both electron neutrinos and antineutrinos leads to inaccurate recoveries of the neutrino luminosities in some cases, with |$\epsilon _{L_{\nu _e}}$| up to |$52{{\ \rm per\ cent}}$|⁠.

Open in new tabDownload slide

5.1.2 Thermalization

We set α_blk to the above values, and look at both |$\epsilon _{L_i}$| and |$\epsilon _{\langle E_i \rangle }$| for each species |$i \in [\nu _e,\bar{\nu }_e]$| in Fig. 7. The first two rows correspond to the results for case (1), the third and fourth row to the results for case (2), and the last two rows for case (3). We show both results for electron neutrinos (left column) and for electron antineutrinos (right column). As anticipated before, |$\epsilon _{L_{\nu _e}}$| shows the largest dependence on p for a given τ_cut. We find that for electron neutrinos |$\epsilon _{L_{\nu _e}} \lesssim 12{{\ \rm per\ cent}}$|⁠, |$\epsilon _{L_{\nu _e}} \lesssim 23{{\ \rm per\ cent}}$|⁠, and |$\epsilon _{L_{\nu _e}} \lesssim 52{{\ \rm per\ cent}}$| for cases (1), (2), and (3), respectively. If we look at electron antineutrinos we have |$\epsilon _{L_{\bar{\nu }_e}}$| up to |$\sim 20{{\ \rm per\ cent}}$|⁠, |$\epsilon _{L_{\bar{\nu }_e}} \lesssim 8{{\ \rm per\ cent}}$|⁠, and |$\epsilon _{L_{\bar{\nu }_e}} \lesssim 31{{\ \rm per\ cent}}$|⁠, respectively. On the other hand, |$\epsilon _{\langle E_i \rangle } \lesssim 12{{\ \rm per\ cent}}$| for each species i in cases (1) and (2), while case (3) shows |$\epsilon _{\langle E_{\bar{\nu }_e} \rangle } \lesssim 23{{\ \rm per\ cent}}$|⁠. While the overall limited values in both |$\epsilon _{L_i}$| and |$\epsilon _{\langle E_i \rangle }$| suggest that any value of τ_cut ∈ [5, 10, 15, 20] seems good enough to describe the thermalization for both cases (1) and (2), in case (3) only either τ_cut = 5 or τ_cut = 10 is able to keep |$\epsilon _{L_{\bar{\nu }_e}} \lt 30{{\ \rm per\ cent}}$|⁠, in spite of the large |$\epsilon _{L_{\nu _e}} \gtrsim 52{{\ \rm per\ cent}}$|⁠. We therefore always take τ_cut = 10 as reference for the ASL luminosities and mean energies shown in Table 2.

We notice that the electron neutrino and antineutrino mean energies are systematically higher for cases (2) and (3) in both the M1 and the ASL by up to ∼1–3 MeV with respect to case (1). We provide an explanation by calculating the average temperatures at which both neutrino species are emitted by means of equation 9 of Rosswog & Liebendörfer (2003). Table 3 shows the results we find. The average temperatures are higher by ≳1 MeV with respect to case (1), implying that most of the emission for both electron neutrinos and antineutrinos comes from hotter regions. Since neutrinos thermalize with matter before free-streaming, the mean neutrino energy is higher if temperatures are higher. The presence of hotter regions is also confirmed by the fact that the maximum temperature seen for case (1) is ∼25 MeV, for case (2) it reaches ∼40 MeV, and for case (3) it is ∼32 MeV. However, the higher electron antineutrino mean energy in case (3) is also a consequence of the rather large |$\epsilon _{\bar{\nu }_e}$| (Fig. 7, right-hand panel on the sixth row).

5.1.3 Heating

To constrain the heating parameter p, we show in the first three rows of Fig. 8 dL_i/dμ as a function of θ for all cases examined, and for i = ν_e (left-hand panels) and |$i=\bar{\nu }_e$| (right-hand panels). The curve from M1 is obtained by taking an average of dL_i/dμ over different azimuthal angles ϕ of the domain for each polar angle θ, while the ASL curve is obtained by means of equation (19), with the best blocking parameter, a reference value of τ_cut = 10,² and the best p from the direct comparison between the polar angle profiles of the luminosities. Unlike case (1), we only show the 2D result from sedonu and the ASL curve of |$\mathrm{ d}L_{\bar{\nu }_e}/\mathrm{ d} \mu$| for cases (2) and (3) because of the impact of the 2D assumption on the electron antineutrino emission. We see that p = 2 nicely fits trends from sedonu for both electron neutrinos and antineutrinos and in all cases. We can also see that M1 overestimates the flux close to the polar axis for both electron neutrinos and antineutrinos in cases (1) and (2), with a maximum disagreement with respect to the result from sedonu by a factor of ≲2 at the pole in case (1). This is in accordance with the previous findings of Foucart et al. (2018, 2020), and it is the result of the approximations introduced by the analytical closure. It is also important to notice that our new estimate of the parameter p = 2 is lower than the one we found in Gizzi et al. (2019), where we estimated p = 8. The reason is the fact that there we chose M1 as source for comparison. As just said, the approximations introduced by the analytical closure inevitably leads to a steeper decrease of the neutrino fluxes with the polar angle, consequently suggesting a larger value of p than the one found here. We justify our decision of using the 2D data from sedonu of electron antineutrinos in cases (2) and (3) for the calibration of p = 2 by noticing that the 2D assumption does not impact the trend of |$\mathrm{ d}L_{\bar{\nu }_e}/\mathrm{ d}\mu$| in M1 (see purple curve in the last row of Fig. 8). The last row shows also the quantity |$\mathrm{ d}L_{\bar{\nu }_e}/\mathrm{ d}\mu$| as a function of the polar angle θ for different azimuthal angles ϕ, obtained with M1 in 3D (green dots). We can clearly see that the variation with ϕ is limited, particularly at angles θ ≲ π/3 ≈ 1.05 where the bulk of neutrino-driven winds is located. This justifies the assumption of axially symmetric fluxes entering equation (5).

Figure 8.

dL_i/dμ as a function of θ. From the first to the third row we show the results for cases from (1) to (3), respectively, for i = ν_e (left) and |$i=\bar{\nu }_\mathrm{ e}$| (right). Dots represent the exact solution from sedonu. The blue curve describes the trend of dL_i/dμ obtained with M1 when performing an azimuthal average over different ϕ at each θ, while the red curve describes the one from the ASL, scaled by some constant, and obtained via equation (19) with the best parameter set of each case. We only show the trends from sedonu and ASL for antineutrinos in cases (2) and (3). The trend from sedonu is overall well described by a |$\rm {cos}^{2}(\theta)$| for both species and in all cases. M1 tends to overestimate the flux close to the polar axis for both electron neutrinos and antineutrinos, as a consequence of the approximations introduced by the analytical closure. The last row shows |$\mathrm{ d}L_{\bar{\nu }_e}/\mathrm{ d} \mu$| as a function of θ for case (2) (left), and case (3) (right). The 2D assumption does not affect the trend of |$\mathrm{ d}L_{\bar{\nu }_e}/\mathrm{ d} \mu$|⁠. Moreover, for a given θ we see a limited variation of |$\mathrm{ d}L_{\bar{\nu }_e}/\mathrm{ d} \mu$| with ϕ at polar angles relevant for neutrino-driven winds (θ ≲ π/3).

Open in new tabDownload slide

5.2 Combination of parameter constraints

Among the three calibrated parameters, both τ_cut and p do not show variations by changing the binary configuration. In particular, τ_cut is always around a value of 10, and p = 2. On the other hand, α_blk is more sensitive than the other parameters when moving from one binary configuration to another. Similar to the results of Perego et al. (2016), the blocking parameter may vary in a range [0.45,0.75], depending on the configuration of the binary and its time after merger. Under the assumption |$\alpha _{\rm {blk}} \equiv \alpha _{\nu _e,\rm {blk}}=\alpha _{\bar{\nu _e},\rm {blk}}$|⁠, we find larger values of α_blk for cases (2) and (3) with respect to case (1). Specifically, α_blk is |$\sim 44{{\ \rm per\ cent}}$| and |$\sim 66{{\ \rm per\ cent}}$| larger, respectively. However, it is important to consider that the values of luminosities and mean energies from M1 that we have taken as reference in cases (2) and (3) might be off by |$\sim 20{{\ \rm per\ cent}}$| with respect to an exact solution to the transport (Foucart et al. 2020), implying that our calibrated α_blk might be slightly affected too.

We find that assuming |$\alpha _{\nu _e,\rm {blk}}= \alpha _{\bar{\nu _e},\rm {blk}}$| is not a good choice at early times after merger. In particular, we notice that higher values of α_blk do not well capture the total luminosity of electron neutrinos (see left-hand panels in the third and fifth row of Fig. 7), pointing to the need for |$\alpha _{\nu _e,\rm {blk}} \lt \alpha _{\bar{\nu }_e,\rm {blk}}$|⁠. This conclusion is consistent with the fact that electron antineutrinos are the most emitted species in merger environments, and therefore likely more affected by blocking effects. Moreover, we find the electron antineutrino gas to have on average fewer energy states available to be populated by new emitted neutrinos of the same species. Indeed, by calculating the difference between the average degeneracy parameter³ of electron neutrinos and antineutrinos, we find positive values for all the cases examined (see eighth column of Table 3). In light of this, we optimize the choice of blocking parameter by exploring in Fig. 9|$\epsilon _{L_{i}}$| for |$i=[\nu _\mathrm{ e},\bar{\nu }_\mathrm{ e}]$| and for the three cases examined in Section 4 by assuming τ_cut = 10, p = 2, and |$\alpha _{\nu _e,\rm {blk}} \ne \alpha _{\bar{\nu _e},\rm {blk}}$|⁠. We confirm that for each binary the value of the blocking corresponding to the minimum |$\epsilon _{L_{i}}$| is lower for i = ν_e than for |$i=\bar{\nu }_\mathrm{ e}$|⁠. For electron neutrinos, the largest blocking parameter at which |$\epsilon _{L_{\nu _\mathrm{ e}}}$| is minimum is found in case (2), with |$\alpha _{\nu _e,\rm {blk}}= 0.55$|⁠. By looking at Table 3, this is on one side due to the larger average emission temperature (⁠|${\langle T \rangle }_{\nu _\mathrm{ e}} = 5.51\:\rm {MeV}$|⁠), for which the emission rate is enhanced considerably due to the ∼T⁵ dependence for charged-current interactions (Rosswog & Liebendörfer 2003), and on the other side to the less neutron-rich environment (⁠|${\langle Y_\mathrm{ e} \rangle }_{\nu _\mathrm{ e}}= 0.104$|⁠) that favours electron captures on protons. Both factors contribute in providing more electron neutrinos, and therefore enhancing Pauli blocking effects. For electron antineutrinos the largest blocking parameter at which |$\epsilon _{L_{\bar{\nu }_\mathrm{ e}}}$| is minimum is found for case (3), with |$\alpha _{\bar{\nu }_\mathrm{ e},\rm {blk}}= 0.75\!-\!0.85$|⁠. Again, this is due to both a rather hot (⁠|${\langle T \rangle }_{\bar{\nu }_e}= 6.41\:\rm {MeV}$|⁠) and neutron-rich material (⁠|${\langle Y_\mathrm{ e} \rangle }_{\nu _\mathrm{ e}}= 0.073$|⁠) that favours positron captures on neutrons. A value of |$\alpha _{\nu _\mathrm{ e},\rm {blk}}= 0.45$| provides |$\epsilon _{L_{\nu _\mathrm{ e}}} \lesssim 20{{\ \rm per\ cent}}$| for all binaries, and we therefore set it as fiducial for this species. Regarding the antineutrinos, the variability of |$\epsilon _{L_{\bar{\nu }_e}}$| with |$\alpha _{\bar{\nu }_e,\rm {blk}}$| is larger and makes the choice of the best blocking parameter more cumbersome. In particular, equal mass binaries prefer |$\alpha _{\bar{\nu }_e,\rm {blk}}= 0.55$| and |$\alpha _{\bar{\nu }_e,\rm {blk}}= 0.65$|⁠, for cases (1) and (2) respectively, while case (3) prefers |$\alpha _{\bar{\nu }_e,\rm {blk}} \gt 0.65$|⁠. We therefore suggest a fiducial value of |$\alpha _{\bar{\nu }_e,\rm {blk}}= 0.55$| for equal mass binaries, such that |$\epsilon _{L_{\bar{\nu }_e}} \lesssim 20{{\ \rm per\ cent}}$|⁠, and |$\alpha _{\bar{\nu }_e,\rm {blk}}= 0.75$| for unequal mass binaries, such that |$\epsilon _{L_{\bar{\nu }_e}} \lt 30{{\ \rm per\ cent}}$|⁠. However, the unequal mass case might need to be explored in other test cases for a more robust gauging of |$\alpha _{\bar{\nu }_e,\rm {blk}}$|⁠. Besides the thermodynamical, compositional, and degeneracy properties of the matter described in Table 3 and determining the extent of Pauli blocking, the increasing value of the blocking parameter for both electron neutrinos and antineutrino species when moving from case (1) to case (3) is also consistent with the fact that generally less massive and/or unequal mass binaries produce larger discs (Rosswog et al. 2000; Bernuzzi et al. 2020; Vincent et al. 2020). As stated in Section 2.3, the blocking parameter takes also into account the reduction of neutrino emission due to inward neutrino fluxes in the semitransparent regime. The presence of larger discs in less massive and/or unequal mass binaries leads to larger neutrino surfaces, and consequently to an overall larger effect of inward neutrino fluxes, therefore contributing to an increase in the size of blocking.

Figure 9.

|$\epsilon _{L_{\nu }}$| as a function of |$\alpha _{\nu ,\rm {blk}} \in [0.35,0.45,0.55,0.65,0.75,0.85]$|⁠, for the three binary set-ups and for both electron neutrinos (green curve) and antineutrinos (orange curve). Plots are obtained with τ_cut = 10 for all binaries, and with p = 2 and p = 4 for equal and unequal mass binaries, respectively. For each binary, the value of the blocking corresponding to the minimum |$\epsilon _{L_{i}}$| is lower for i = ν_e than for |$i=\bar{\nu }_\mathrm{ e}$|⁠. While for |$\alpha _{\nu _e,\rm {blk}}= 0.45$||$\epsilon _{L_{\nu _e}} \lesssim 20{{\ \rm per\ cent}}$| in all binaries, for the antineutrinos the larger variability of |$\epsilon _{L_{\bar{\nu }_e}}$| with |$\alpha _{\bar{\nu }_e,\rm {blk}}$| when exploring different cases leads to a different fiducial value for equal and unequal mass binaries. In particular, |$\alpha _{\bar{\nu }_e,\rm {blk}}=0.55$| for the former case, and |$\alpha _{\bar{\nu }_e,\rm {blk}}=0.75$| for the latter case. In this way, |$\epsilon _{L_{\bar{\nu }_e}} \lesssim 20{{\ \rm per\ cent}}$| and |$\epsilon _{L_{\bar{\nu }_e}} \lt 30{{\ \rm per\ cent}}$|⁠, respectively.

Open in new tabDownload slide

We conclude by showing in Figs 10–11 the distribution of the rate of change of the specific matter internal energy |$\dot{e}$| (units of |$10^{20}\:\rm {erg}\:\rm {g}^{-1}\:\rm {s}^{-1}$|⁠) and of the electron fraction |$\dot{Y}_{\rm {e}}$| in the winds for case (1) and for the different transport approaches, assuming |$\alpha _{\nu _e,\rm {blk}}= 0.45, \alpha _{\bar{\nu }_e,\rm {blk}}= 0.55, \tau _{\rm {cut}}= 10, p=2$|⁠. The upper plots are 3D maps in a box of |$x \times y \times z = 150\:\rm {km} \times 150\:\rm {km} \times 150\:\rm {km}$|⁠, while the lower ones are projections on the x-y plane. We assume a threshold density below which we identify the wind region a value of |$5 \times 10^{9}\:\rm {g\:cm^{-3}}$|⁠. However, this limit in the end also includes the outer regions of the disc, visible as |$\dot{e}\lt 0$| regions around z = 0 km. The plots show the particle distribution, and we recover the values of |$\dot{e}$| and |$\dot{Y}_{\rm {e}}$| of each SPH particle from M1 and sedonu via interpolation from the respective grids. Overall, both the ASL and the M1 |$\dot{e}$| and |$\dot{Y}_{\rm {e}}$| distributions agree well with the solution from sedonu. However, they both show some cooling (⁠|$\dot{e}\lt 0$|⁠) above the remnant which is weaker or absent in sedonu, as well as a few particles with |$\dot{Y}_{\rm {e}}\lt 0$|⁠. More precisely, for the ASL |$\approx 14{{\ \rm per\ cent}}$| of the particles have |$\dot{e}_{\rm {ASL}} \cdot \dot{e}_{\rm {Sed}}\lt 0$|⁠, and |$\approx 3{{\ \rm per\ cent}}$| have |$\dot{Y}_{e,\rm {ASL}} \cdot \dot{Y}_{e,\rm {Sed}}\lt 0$|⁠. For the M1, these numbers are |$\approx 27{{\ \rm per\ cent}}$| and |$\approx 1{{\ \rm per\ cent}}$|⁠, respectively. Considering the overall similarity between the maps of the ASL and M1 we can conclude within the limits of this analysis that the ASL may show similar performances of M1 in dynamical simulations when comparing the wind properties against exact solutions to the transport. Nevertheless, a more robust assessment requires to explore this comparison dynamically and for different binary configurations. The main advantage of our SPH-ASL is in its efficiency. We indeed estimate that, if we had to assume the same time-step as in M1 for a dynamical simulation, and by taking magma2 as our Lagrangian hydrodynamics code (Rosswog 2020), the ratio between the CPU hours spent for the transport and those for the hydrodynamics is about 0.8 per time-step, to be compared with a factor of 10 for the M1 in flash. The number we find is similar to the one in Pan et al. (2018). The ratio for the SPH-ASL could be further cut down if we consider that the optical depth computation may not be required at every time-step, unless the thermodynamical and compositional properties of the matter change considerably. Future dynamical simulation will definitely provide a more robust assessment of the performance.

Figure 10.

3D maps (top row) and projection of the 3D maps on the x-y plane (bottom row) of the rate of change of the internal energy |$\dot{e}$| in the winds (units of |$10^{20}\:\rm {erg}\:\rm {g}^{-1}\:\rm {s}^{-1}$|⁠) for sedonu (left column), the ASL (middle column), and M1 (right column), case (1). The ASL performance is similar to M1, and they both generally reproduce the distribution of |$\dot{e}$| from sedonu. The main difference with respect to sedonu is in somewhat lower |$\dot{e} \gt 0$| and stronger |$\dot{e} \lt 0$| locally. The domain size is |$x \times y \times z = 150\:\rm {km} \times 150\:\rm {km} \times 150\:\rm {km}$|⁠. The map for the ASL is obtained with the parameter set |$\alpha _{\nu _e,\rm {blk}}= 0.45, \alpha _{\bar{\nu }_e,\rm {blk}}= 0.55, \tau _{\rm {cut}}= 10, p=2$|⁠.

Open in new tabDownload slide

Figure 11.

3D maps (top row) and projection of the 3D maps on the x-y plane (bottom row) of the rate of change of the electron fraction |$\dot{Y}_{\rm {e}}$| in the winds for sedonu (left column), the ASL (middle column), and M1 (right column). The ASL performance is very similar to M1, and they both generally reproduce the distribution of |$\dot{Y}_{\rm {e}}$| from sedonu. The latter provides always |$\dot{Y}_{\rm {e}} \gt 0$|⁠, while both the ASL and M1 show a few particles with |$\dot{Y}_{\rm {e}} \lt 0$|⁠. Moreover, |$\dot{Y}_{\rm {e}} \gt 0$| from sedonuappears stronger than in the ASL and M1. The domain size is |$x \times y \times z = 150\:\rm {km} \times 150\:\rm {km} \times 150\:\rm {km}$|⁠. The map for the ASL is obtained with the parameter set |$\alpha _{\nu _e,\rm {blk}}= 0.45, \alpha _{\bar{\nu }_e,\rm {blk}}= 0.55, \tau _{\rm {cut}}= 10, p=2$|⁠.

Open in new tabDownload slide

6 CONCLUSIONS

In this paper, we have provided a detailed calibration analysis of the ASL scheme presented earlier in Gizzi et al. (2019), which is based on the original work of Perego et al. (2016). Our main motivation is the study of neutrino-driven winds emerging from neutron star mergers. The gauging process is performed by post-processing a number of snapshots of a binary neutron star merger remnant. We extract neutrino quantities directly impacted by each parameter, and we compare the ASL results for different parameter combinations with the ones obtained from the Monte Carlo neutrino transport code sedonu (Richers et al. 2015) and from a two-moment scheme (M1) implemented in flash (Fryxell et al. 2000; O’Connor 2015; O’Connor & Couch 2018). In the calibration process we focus on electron-type neutrinos and antineutrinos, since they determine the properties of neutrino-driven winds. We summarize our main findings as follows:

1. Performing neutrino transport in 2D by post-processing initial 3D post-merger configurations is of limited accuracy at early (t ≲ 20 ms) times post-merger. In particular, the 2D averaging can severely impact the total luminosities by more than a factor of 2.

2. The assumption of axially symmetric neutrino fluxes entering the heating rate in the ASL is validated by 3D neutrino transport simulations. In particular, variations of the fluxes with the azimuthal angle ϕ and for a given polar angle θ are limited in the bulk region of neutrino-driven winds.

3. In agreement with Perego et al. (2016), the thermalization parameter τ_cut has an impact mainly on the neutrino mean energies. However, a value of τ_cut = 10 robustly recovers neutrino mean energies within |$\lesssim 25{{\ \rm per\ cent}}$| accuracy.

4. The heating parameter introduced in Gizzi et al. (2019) is recalibrated to a lower value, as a result of the usage of sedonu as reference solution. We find p = 2 to best reproduce the distribution of the neutrino fluxes in all the cases examined. Using an M1 scheme rather than sedonu would lead to the artefact of a larger p than the ones calibrated here because of the approximations introduced by the analytical closure.

5. The blocking parameter α_blk mainly impacts the total neutrino luminosities, in agreement with the results of Perego et al. (2016). Moreover, unlike the other two parameters it is the most sensitive to a change of binary configuration.

6. The assumption |$\alpha _{\rm {blk}} \equiv \alpha _{\nu _e,\rm {blk}}= \alpha _{\bar{\nu }_e,\rm {blk}}$| adopted in Perego et al. (2016) can be rather inaccurate for recovering neutrino luminosities. In particular, electron neutrino luminosities are lower by up to a factor of 2 with respect to the reference solution, suggesting |$\alpha _{\nu _e,\rm {blk}} \lt \alpha _{\bar{\nu }_e,\rm {blk}}$|⁠.

7. Assuming |$\alpha _{\nu _e,\rm {blk}} \ne \alpha _{\bar{\nu }_e,\rm {blk}}$|⁠, we indeed find that |$\alpha _{\nu _e,\rm {blk}}= 0.45$| provides electron neutrino luminosities in agreement with the reference solution at the level of |$\lesssim 20{{\ \rm per\ cent}}$| accuracy, for both equal and unequal mass binaries. Similarly, |$\alpha _{\bar{\nu }_e,\rm {blk}}= 0.55$| results in electron antineutrino luminosities off by |$\lesssim 20{{\ \rm per\ cent}}$| for equal mass binaries.

8. We find |$\alpha _{\bar{\nu }_e,\rm {blk}} = 0.75$| for unequal mass binaries, leading to a relative error in the antineutrino luminosity of |$\lt 30{{\ \rm per\ cent}}$|⁠. However, more test cases might be needed for a more robust evaluation and to reduce systematics.

9. In contrast to Foucart et al. (2020), the heavy-lepton neutrino luminosity is systematically larger by a factor of a few in both the ASL and M1 with respect to an exact solution to the transport. This enhances the overall cooling of the remnant in our snapshot calculations. For the ASL, the most probable explanation is the poor treatment of the diffusion time-scale, which according to Ardevol-Pulpillo et al. (2019) can boost luminosities by more than a factor of 2. We expect this treatment to affect also the electron neutrino and antineutrino luminosities to some extent.

10. The properties of neutrino-driven winds are shaped by the rates of change of internal energy and electron fraction, |$\dot{e}$| and |$\dot{Y}_e$|⁠, respectively. The corresponding maps for a 1.4–1.4 M_⊙ binary demonstrate that for our suggested parameter choice the ASL scheme performs similar to the M1 approach. From the perspective of dynamical simulations, our SPH-ASL comes with the advantage of a better efficiency. In particular, by taking the Lagrangian hydrodynamics code magma2 (Rosswog 2020), we estimate that the ratio between the CPU hours spent for the transport and those spent for the hydrodynamics would be ≲0.8 per time-step, while for the M1 in flash is about 10. In other words, the ASL scheme could be applied in SPH simulations with only a moderate additional computational effort.

Although the geometry of a binary neutron star merger allows neutrinos to escape with more directional freedom than in a spherically symmetric core-collapse supernova, the results we find here show that the blocking parameter can still be quite high in merger remnants. Apart from the different thermodynamics and composition of the matter, the major reason is the disc geometry that increases the effect of inward neutrino fluxes at the neutrino surfaces with respect to a spherically symmetric geometry. In this paper, we have also presented a completely mesh-free, particle-based algorithm to compute spectral, species-dependent optical depths, based on the SPH (Monaghan 1992; Monaghan 2005; Rosswog 2009, 2015a, b, 2020). This algorithm makes our ASL fully grid-independent, and therefore suitable for future SPH dynamical simulations of binary neutron star mergers with neutrino transport.

ACKNOWLEDGEMENTS

We would like to thank Sherwood Richers for developing and supporting our use of sedonu and for a careful reading of the manuscript. This work has been supported by the Swedish Research Council (VR) under grant number 2016- 03657_3, by the Swedish National Space Board under grant number Dnr. 107/16, by the Research Environment grant ‘Gravitational Radiation and Electromagnetic Astrophysical Transients (GREAT)’ funded by the Swedish Research Council (VR) under Dnr 2016-06012, and by the Knut and Alice Wallenberg Foundation under Dnr KAW 2019.0112. EOC is supported by the Swedish Research Council (Project No. 2018-04575 and 2020-00452). We gratefully acknowledge support from COST Action CA16104 ‘Gravitational waves, black holes and fundamental physics’ (GWverse), from COST Action CA16214 ‘The multi-messenger physics and astrophysics of neutron stars’ (PHAROS), and from COST Action MP1304 ‘Exploring fundamental physics with compact stars (NewCompStar)’.

The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC (Centre for High Performance Computing) and on the facilities of the North-German Supercomputing Alliance (HLRN) in Göttingen and Berlin.

DATA AVAILABILITY

The data concerning the initial conditions needed to run the simulations, as well as the extracted neutrino quantities from all neutrino transport codes and neutrino-driven wind maps are available from the authors of this manuscript upon request.

Footnotes

REFERENCES