Calibration of the Advanced Spectral Leakage scheme for neutron star merger simulations, and extension to smoothed-particle hydrodynamics
Davide Gizzi, Christoffer Lundman, Evan O'Connor, Stephan Rosswog, Albino Perego
MMNRAS , 1–18 (2020) Preprint 18 February 2021 Compiled using MNRAS L A TEX style file v3.0
Extension of the Advanced Spectral Leakage scheme for neutron starmerger simulations
D. Gizzi, ★ C. Lundman, E. O’Connor, S. Rosswog, A. Perego , The Oskar Klein Centre, Department of Astronomy, Albanova, Stockholm University, SE-106 91 Stockholm, Sweden Dipartimento di Fisica, Università degli Studi di Trento, via Sommarive 14, 38123 Trento, Italy INFN-TIFPA, Trento Institute for Fundamental Physics and Applications, via Sommarive 14, I-38123 Trento, Italy
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We extend a neutrino transport approximation, called Advanced Spectral Leakage (ASL), with the purpose of modeling neutrino-driven winds in neutron star mergers. Based on a number of snapshots we gauge the ASL parameters by comparing against boththe two-moment (M1) scheme implemented in the
FLASH code and the Monte Carlo neutrino code
Sedonu . The ASL schemecontains three parameters, the least robust of which results to be a blocking parameter for electron neutrinos and anti-neutrinos.The parameter steering the angular distribution of neutrino heating is re-calibrated compared to the earlier work of Gizzi et al.(2019). We also present a new, fast and mesh-free algorithm for calculating spectral optical depths, which, when using SmoothedParticle Hydrodynamics (SPH), makes the neutrino transport completely particle-based. We estimate a speed-up of a factorof (cid:38)
100 in the optical depth calculation when comparing to a grid-based approach. In the suggested calibration we recoverluminosities and mean energies within 25%. A comparison of the rates of change of internal energy and electron fraction in theneutrino-driven wind suggests comparable accuracies of ASL and M1, but a higher computational efficiency of the ASL scheme.We estimate that the ratio between the CPU hours spent on the ASL neutrino scheme and those spent on the hydrodynamicsis (cid:46) . FLASH . Key words: neutrinos, radiative transfer, hydrodynamics, star: neutron, stars: supernovae: general
The combined detection of gravitational and electromagnetic wavesfrom the neutron star merger event GW170817 (Abbott et al. 2017a,c)has solved several decades-old puzzles. For example, the detectionof the short Gamma-Ray Burst (sGRB) GRB170817A ∼ .
74 safter the merger (Savchenko et al. 2017; Goldstein et al. 2017) con-firmed the long suspected association between sGRBs and compactbinary mergers (Eichler et al. 1989; Paczynski 1986; Narayan et al.1992). Moreover, the detection of the macronova/kilonova transientAT2017gfo (Chornock et al. 2017; Kasliwal et al. 2017; Kilpatricket al. 2017; Pian et al. 2017; Coulter et al. 2017; Soares-Santos et al.2017; Lipunov et al. 2017; Arcavi et al. 2017; Evans et al. 2017;Smartt et al. 2017) confirmed binary neutron star mergers as a major,and possibly dominant, r-process nucleosynthesis site (Rosswog et al.2018; Kasen et al. 2017; Drout et al. 2017), thus confirming earliertheoretical predictions on the subject (Lattimer & Schramm 1974;Eichler et al. 1989; Rosswog et al. 1999; Freiburghaus et al. 1999).The event GW170817 also allowed to place important constraints onthe Equation of State (EoS) of nuclear matter (Most et al. 2018; Deet al. 2018; Abbott et al. 2018; Radice et al. 2018a; Radice & Dai2019; Coughlin et al. 2019; Kiuchi et al. 2019; Bauswein et al. 2017;Jiang et al. 2019, 2020) and provided an independent measurement ★ E-mail: [email protected] of the Hubble parameter (Abbott et al. 2017b; Dhawan et al. 2020).The luminous blue component of the macronova transient AT2017gfohas also emphasized the role played by neutrinos in changing the elec-tron fraction 𝑌 𝑒 of the ejected matter. Weak interactions modify theelectron fraction in a fair portion of the ejecta from initial values of 𝑌 𝑒 ∼ .
05 to values above 𝑌 𝑒 = .
25, where the resulting nucleosyn-thesis changes abruptly (Korobkin et al. 2012; Lippuner & Roberts2015) and no more lanthanides and heavier nuclei are produced. Thisincrease in 𝑌 𝑒 also drastically reduces the optical opacities (Kasenet al. 2013; Tanaka & Hotokezaka 2013) and leads to blue (ratherthan red) transients after about a day rather than a week. Fits to theblue component of the light curve from macronova models suggestejecta of several 10 − 𝑀 (cid:12) with velocity ∼ . 𝑐 . There is a generalagreement that secular mass ejection from the remnant disk is neededto achieve such large ejecta amounts (e.g. Ciolfi & Kalinani (2020);Nedora et al. (2020); Radice et al. (2018b)), but it is still a matter ofdebate to which extent the different proposed secular ejection chan-nels contribute to the blue macronova. The ejection channels includewinds driven by magnetic fields (Ciolfi et al. 2017; Ciolfi & Kalinani2020), viscosity (Fernández & Metzger 2013; Fernández et al. 2019;Fujibayashi et al. 2018, 2020; Siegel & Metzger 2017), neutrinos(Dessart et al. 2009; Perego et al. 2014a; Martin et al. 2015) andspiral-waves (Nedora et al. 2019, 2020). Our motivation here areneutrino-driven winds, for which no systematic study in dynamicalsimulations of binary neutron star mergers exists to date. For this © a r X i v : . [ a s t r o - ph . H E ] F e b Gizzi D. et al. reason, we have recently implemented an extension to the originalAdvanced Spectral Leakage (ASL) scheme (Perego et al. 2016) withthe purpose of modelling the neutrino absorption responsible forlaunching the winds in merger simulations (Gizzi et al. 2019).In this paper we extend the analysis of the ASL scheme by per-forming a careful calibration for the three parameters that enter thescheme. These parameters have an impact on both neutrino emis-sion and absorption, and must therefore be gauged carefully. Two ofthese parameters were already introduced in Perego et al. (2016), butcalibrated for core-collapse supernovae simulations. They define theneutrino emission and represent two physical effects. The first is ablocking parameter, accounting for both Pauli blocking effects dueto the fermionic nature of neutrinos and for inward neutrino fluxeswhen calculating the neutrino emission above the decoupling region,while the second is a thermalization parameter, describing the num-ber of weak interactions needed to thermalize neutrinos. We added athird parameter in Gizzi et al. (2019) specifically for modelling theneutrino absorption responsible for launching the neutrino-drivenwinds, and we gauged it by comparing the outcome of the simu-lations with the two-moment scheme (M1) approach implementedin
FLASH (Fryxell et al. 2000; O’Connor 2015; O’Connor & Couch2018). Here, we use both the M1 and the Monte Carlo neutrino trans-port code
Sedonu (Richers et al. 2015) as sources of comparison.We use neutrino quantities directly impacted by each parameter toperform our calibration and we extract them by post-processing anumber of snapshots of binary neutron star remnants.When using the ASL scheme, we present a new, mesh-free im-plementation for calculating spectral, multi-group optical depths forSmoothed-Particle Hydrodynamics (SPH) (Monaghan 1992; Mon-aghan 2005; Rosswog 2009, 2015a,b, 2020) simulations. This makesour ASL fully mesh-free and ideal for dynamical SPH simulationsof merging neutron stars with neutrino effects.The paper is structured as follows. We begin with a description ofthe neutrino transport codes in Sec. 2. We then describe the set ofsnapshots that we use, as well as the calibration strategy in Sec. 3.The results of the calibration are discussed in Sec. 4. We explainand test the new optical depth algorithm in Appendix A, includinga comparison with grid-based calculations. In Sec. 5 we summariseand draw our conclusions.
The neutrino transport approaches we use for our analysis are: theMonte Carlo code
Sedonu , the two-moment scheme (M1) developedin the Eulerian hydrodynamics code
FLASH , and the ASL scheme. Allof them are spectral, (i.e., neutrino energy is treated as an independentvariable and neutrino-matter interactions are calculated for differentneutrino energies), an important ingredient given the energy-squareddependence entering neutrino cross sections (see e.g. Burrows &Thompson (2004) and references therein). Moreover, we choose thesame set of weak interactions: production and absorption of electronneutrinos and anti-neutrinos via charged current processes involvingnucleons and nuclei, neutrino emission by bremsstrahlung and pairprocesses, and finally elastic scattering off nucleons and nuclei. Wechoose the SFHo EoS (Steiner et al. 2013) to calculate spectral emis-sivities and opacities.
Sedonu and M1 interpolate them from the
NuLib tables (O’Connor 2015). The ASL scheme calculates themfollowing Bruenn (1985); Mezzacappa & Messer (1999); Hannestad& Raffelt (1998). Three distinct neutrino species are modelled: elec-tron neutrinos 𝜈 𝑒 , electron anti-neutrinos ¯ 𝜈 𝑒 and a collective species 𝜈 𝑥 for heavy-lepton neutrinos, comprising 𝜇 and 𝜏 neutrinos andanti-neutrinos. Monte Carlo methods can be used to solve the Boltzmann equation.They provide an exact solution to the transport problem in the limit ofan infinite particle number. However, Monte Carlo neutrino transportmethods tend to be computationally expensive, since a large numberof particles is needed to keep the numerical noise in the solution at alow level. These methods have therefore mostly been used on static,spherically symmetric core-collapse supernovae snapshots (Janka& Hillebrandt 1989a,b; Janka 1992; Keil et al. 2003; Abdikamalovet al. 2012). Similarly, studies of Monte Carlo neutrino transport inneutron star mergers have been carried out by post-processing snap-shots, assuming the neutrino field can be approximated as steady-state(Richers et al. 2015). Recently, Foucart et al. (2020) performed thefirst dynamical hydrodynamic simulation of merging neutron stars upto ∼ Sedonu (Richers et al. 2015) on axisymmetric fluid snapshotsto calibrate the parameters of our ASL scheme.We use a logarithmically spaced energy grid that spans the range 𝐸 𝜈 ∈ [ . , ] MeV. Neutrinos are simulated as packets passingthrough and interacting with a static fluid background imported fromhydrodynamic simulations. Each packet is specified by the neutrinoenergy 𝐸 𝜈 , the location x , the unit vector ˆ d describing the propa-gation direction, and the total number of neutrinos in the packet. Afixed number of neutrino packets are emitted from each grid cell. Inparticular, Sedonu randomly samples the energy-dependent emis-sivity to determine the neutrino energy, and the packet direction ofpropagation is drawn from an isotropic distribution. Neutrino packetsinteract with the fluid by being partially absorbed during the packetpropagation, and changing propagation direction when scattered. Thedistance the packet travels before interacting is randomly drawn froman exponential distribution that depends on the scattering mean freepath. If the packet is scattered, it is given an isotropically random newdirection. While moving, the packet deposits both energy and lep-ton number to the traversed grid-cells, from which the correspondingrates of change can be derived. We emit 100 Monte Carlo packets perenergy group and grid cell. In this way, integrated quantities such asluminosities and mean neutrino energies are stable in spite of MonteCarlo noise, which is estimated to be a thousand times smaller.
A moment scheme is an approximation to the Boltzmann neutrinotransport, obtained by evaluating angular integrals of the Boltzmannequation in order to time-evolve the angular moments of the distribu-tion function (Lindquist 1966; Bruenn et al. 1978; Bruenn 1985; Mez-zacappa & Messer 1999). The M1 scheme is a multi-dimensional,spectral approach where only the 0th and the 1st moments of the neu-trino distribution function are evolved, which respectively describethe 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 (Audit et al. 2002; Levermore & Pomraning1981; Minerbo 1978; Smit et al. 1997; Pons et al. 2000). Gener-ally, the analytical form captures the physics of the transport wellin specific regimes, while approximating the exact solution in oth-ers. Neutrino-matter interactions are included in appropriate sourceterms that appear on the right-hand side of the moment equations
MNRAS , 1–18 (2020)
SL calibration for BNS merger simulations and which are solved via techniques that are borrowed from finitevolume hydrodynamics.The M1 scheme has been widely applied in simulations thatmodel shock revival in core-collapse supernovae (O’Connor 2015;O’Connor & Couch 2018; O’Connor & Ott 2013; Obergaulinger et al.2014; Skinner et al. 2019) and recently compared to other transportapproaches (Pan et al. 2018; Cabezón et al. 2018; O’Connor et al.2018), showing good performance and agreement with alternativeschemes at the level of 10 −
20% in luminosities, mean energies andshock radii. However, not so much has been done in the context ofbinary neutron star mergers. In the work of Foucart et al. (2016)an M1 scheme has been adopted to test the impact of the energyspectrum on the composition of merger outflows. In a subsequentwork Foucart et al. (2018) compared a grey M1 scheme against aMonte Carlo approach, showing the limitations of the M1 analyticalclosure in recovering the properties of the polar ejecta and thereforein modelling neutrino-driven winds and kilonovae. Another quanti-tative study of the assumptions employed in analytical closures andtheir violation has been carried out more recently by post-processingsnapshots with Monte Carlo transport (Richers 2020). Despite theirlimitations, Foucart et al. (2020) have shown that time-dependentneutrino transport in an M1 scheme with an analytical closure thatcarefully treats the energy spectrum approximates the Monte Carlosolution within ∼ −
20% in the composition and velocity of ejecta,as well as electron neutrino and anti-neutrino luminosities and meanenergies. In Sec. 4 we explore the performance of the M1 schemeimplemented in
FLASH and compare it with our new ASL treatmentand the Monte Carlo code
Sedonu . We start by briefly summarizing the original implementation of theASL scheme in its main aspects. For more details, the reader isreferred to Perego et al. (2014a, 2016). Subsequently, we describethe implementation of Gizzi et al. (2019).The ASL is a spectral neutrino scheme, ideal for computationallyinexpensive, multi-dimensional hydrodynamic simulations of core-collapse supernovae (Perego et al. 2016; Pan et al. 2018; Ebinger et al.2020a; Curtis et al. 2019; Ebinger et al. 2020b) and binary neutronstar mergers (Perego et al. 2014a) with neutrino transport effects.We assume a logarithmically-spaced energy grid spanning the range[3,300] MeV. The change in the internal energy and electron fractionof the matter due to weak interactions is estimated by means ofcalculating local neutrino emission and absorption rates. However,their computation depends on an integrated quantity, the spectral,flavour-dependent optical depth , defined at position x as: 𝜏 𝜈 ( 𝐸, x ) = ∫ 𝛾 : x →+∞ 𝜆 𝜈 ( 𝐸, x’ ( 𝑠 )) d 𝑠, (1)where 𝜆 𝜈 ( 𝐸, x’ ) is the neutrino mean free path for neutrinos ofspecies 𝜈 at energy 𝐸 and position x’ , calculated over a path 𝛾 . Herewe introduce a new, mesh-free algorithm to compute 𝜏 𝜈 ( 𝐸, x ) ina SPH context, so that neither the hydrodynamics nor the neutrinotransport depend on any grid implementation, see Appendix A fora detailed description. Depending on the interactions that enter thecomputation of 𝜆 𝜈 ( 𝐸, x’ ) , we can distinguish between the total optical depth 𝜏 𝜈, tot ( 𝐸, x ) , accounting equally for both absorptionand elastic scattering interactions, and the energy optical depth 𝜏 𝜈, en ( 𝐸, x ) , giving more emphasis to energy-exchange interactionsbetween neutrinos and matter. From the above definitions of 𝜏 𝜈 different regimes can be defined: • The equilibrium-diffusive regime, where 𝜏 𝜈, tot (cid:29) 𝜏 𝜈, en (cid:38)
1. Neutrinos can be treated as a trapped Fermi gas inthermal and weak equilibrium with the matter. • The diffusive regime, where 𝜏 𝜈, tot (cid:29) 𝜏 𝜈, en (cid:46) • The semi-transparent regime, where 𝜏 𝜈, tot ∼
1. Althoughneutrinos can still be partially absorbed, they begin to decouplefrom matter, and the surfaces around which this occur are calledneutrino-surfaces, which can be identified when 𝜏 𝜈 ∼ /
3. Low-energy neutrinos decouple earlier than high-energy neutrinos, andtherefore the neutrino surfaces for the latters are wider (Endrizziet al. 2020). In the context of binary neutron star mergers, neutrinoabsorption occurring in this regime (hereafter referred to as heating )drives mass ejection via so-called neutrino-driven winds (Peregoet al. 2014a; Dessart et al. 2009). • In the free-streaming regime , where 𝜏 𝜈, tot (cid:46)
1, locally pro-duced neutrinos stream out freely, and neutrino absorption isnegligible.The local emission rate is given by 𝑟 𝜈 ( 𝐸, x ) = ( − 𝛼 𝜈, blk ) ˜ 𝑟 𝜈 ( 𝐸, x ) Ψ 𝜈 ( x ) exp (− 𝜏 𝜈, en ( 𝐸, x )/ 𝜏 cut ) , (2)where ˜ 𝑟 𝜈 ( 𝐸, x ) is given as smooth interpolation between the produc-tion rate 𝑟 𝜈, prod ( 𝐸, x ) and the diffusion rate 𝑟 𝜈, diff ( 𝐸, x ) :˜ 𝑟 𝜈 ( 𝐸, x ) = 𝑟 𝜈, prod ( 𝐸, x ) 𝑟 𝜈, diff ( 𝐸, x ) 𝑟 𝜈, prod ( 𝐸, x ) + 𝑟 𝜈, diff ( 𝐸, x ) , (3)and Ψ 𝜈 ( x ) = ∫ +∞ ˜ 𝑟 𝜈 ( 𝐸, x ) 𝑒 − 𝜏 𝜈, en ( 𝐸, x )/ 𝜏 cut 𝐸 d 𝐸 ∫ +∞ ˜ 𝑟 𝜈 ( 𝐸, x ) 𝐸 d 𝐸 . (4)Eq. 2 depends on two parameters, so far calibrated in the context ofcore-collapse supernovae simulations:1) a blocking parameter 𝛼 𝜈, blk , that accounts for Pauli block-ing effects due to the large amount of neutrinos locally produced oremitted at the neutrino surface, and for the fact that the interpolationgiven by Eq. 3 favours the production rate in the free-streamingregime. Neutrino emission in this regime is assumed to be isotropic,and a fraction of neutrinos is therefore emitted inward and does notcontribute to the luminosity.2) a thermalization parameter 𝜏 cut , representing the amountof interactions needed to thermalize neutrinos, and which accountsfor neutrinos that exchange energy with the fluid while propagatingfrom the innermost to the outermost regions of the remnant, makingthe neutrino spectrum at the decoupling surface softer.The computation of the heating rate requires the knowledgeof the distribution of the neutrino fluxes coming from the regioninside the neutrino surfaces. While the original implementation ofthe ASL assumes spherically symmetric neutrino fluxes (Peregoet al. 2016), the more complex geometry of a binary neutron starmerger causes anisotropic neutrino fluxes (Perego et al. 2014a;Dessart et al. 2009; Foucart et al. 2016; Rosswog & Liebendörfer2003). We have recently implemented an extension of the originalASL scheme that allows the modelling of neutrino-driven winds MNRAS , 1–18 (2020)
Gizzi D. et al. (Gizzi et al. 2019), which we briefly recap here.The presence of an opaque disk around the central remnant ofa binary neutron star merger makes it difficult for neutrinos toescape along the equatorial region. Therefore an observer locatedfar from the decoupling region is expected to receive the maximumneutrino flux at angles near the poles, and the minimum fluxalong the equatorial region itself (Rosswog & Liebendörfer 2003).Accordingly, we expect most of the neutrino heating to occur alongthe polar regions. We abbreviate as 𝑛 𝜈,𝜏 (cid:46) ( 𝐸, x ) the number ofneutrinos with energy between 𝐸 and 𝐸 + 𝑑𝐸 available for absorptionat position x in the semi-transparent and free streaming regimes. Itcan be approximated as: 𝑛 𝜈,𝜏 (cid:46) ( 𝐸, x ) = ( + 𝑝 )( + 𝛽 𝜈 cos 𝑝 ( 𝜃 )) + 𝑝 + 𝛽 𝜈 𝑙 𝜈 ( 𝐸, | x |) 𝜋 | x | 𝑐𝜇 𝜈 ( 𝐸, x ) , (5)where the first term depends of the polar angle 𝜃 defined with respectto the polar axis of the remnant, and the second term is similar tothe one in Perego et al. (2016), with a spherically symmetric partdivided by the flux factor 𝜇 𝜈 ( 𝐸, x ) . The angular term is assumed tobe axially-symmetric around the polar axis, since some increasingdegree of axial symmetry is generally expected from a few tens ofms irrespective of the mass ratio (Perego et al. 2019). Moreover, itcontains the quantity 𝛽 𝜈 , which is defined as: 𝛽 𝜈 = Λ 𝜈 ( 𝜃 ≈ ° ) Λ 𝜈 ( 𝜃 ≈ ° ) − , (6) Λ 𝜈 ( 𝜃 ) being the luminosity per solid angle at 𝜃 . The computation of Λ 𝜈 ( 𝜃 ) is done with a spectral version of the prescription of Rosswog& Liebendörfer (2003). In the context of merger simulations 𝛽 𝜈 > 𝛽 𝜈 = 𝑝 , providing the flux modula-tion as a function of 𝜃 . The spherically symmetric term contains thespectral number of neutrinos 𝑙 𝜈 per unit time of energy 𝐸 at radius 𝑅 = | x | , calculated by following (Gizzi et al. 2019). The flux factoris computed analytically, similarly to O’Connor & Ott (2010):1 𝜇 𝜈 ( 𝐸, x ) = (cid:40) . 𝜏 𝜈, tot ( 𝐸, x ) + 𝜏 𝜈, tot ( 𝐸, x ) ≤ /
32 otherwise . (7)The flux factor describes the average cosine of the propagation angleof the streaming neutrinos. For an observer far from the neutrinosurfaces the neutrino distribution function peaks radially, therefore 𝜇 𝜈 →
1. On the other hand, for an observer close to the neutrinosurface 𝜇 𝜈 → / 𝑝 has been previously calibrated (Gizzi et al. 2019)by comparing the heating rate maps of electron neutrinos and anti-neutrinos for a binary neutron star merger snapshot against those ob-tained with the M1 scheme in FLASH (Fryxell et al. 2000; O’Connor& Couch 2018). This demonstrated that the new ASL implementationis able to capture the heating rate distribution to better than a factorof 2, with the largest deviations right above the central remnant. Inthis work we further improve the analysis of the new ASL schemeand we provide values for 𝛼 𝜈, blk , 𝜏 cut and 𝑝 that are more accuratefor neutron star mergers. All the details are described in Sec. 3. Re-absorption by heavy-lepton neutrinos is negligible, see Gizzi et al.(2019) for details.
Table 1.
Summary of the binary properties for each case examined. We reportthe masses 𝑚 and 𝑚 of each star in the binary, the mass ratio 𝑞 = 𝑚 / 𝑚 and the time after merger.Case 𝑚 ( 𝑀 (cid:12) ) 𝑚 ( 𝑀 (cid:12) ) 𝑞 time after merger (ms)1 1.4 1.4 1.00 ≈
382 1.3 1.3 1.00 ≈
183 1.3 1.2 0.92 ≈ We take three snapshots of binary neutron star merger remnants fromthe SPH simulations of Rosswog et al. (2017), each one representinga specific case for calibrating our scheme:case 1) an equal mass binary system of 1.4-1.4 𝑀 (cid:12) at ≈ 𝑀 (cid:12) at ≈ 𝑀 (cid:12) at ≈
18 ms after merger. We consider this case as representative ofunequal mass binaries at early times, with the lowest degree of axialsymmetry.We use cases 2) and 3) to test the assumption of axial sym-metry entering Eq. 5. Table 1 summarises the properties of eachbinary configuration. Since the snapshots are taken from simulationsthat use a grey leakage scheme (Rosswog & Liebendörfer 2003),we first map the SPH properties on a 3D grid and let the radiationfield equilibrate with M1 for 3 ms without changing either theelectron fraction or the temperature. Afterwards, we let the electronfraction evolve for 10 ms with M1 to remove grey leakage effects.We do not consider the temperature evolution in M1 since we donot see important differences between the corresponding evolvedand non-evolved profiles. The evolution with M1 allows to makea more consistent comparison between neutrino approaches forour parameter calibration. Figs. 1-2 show density, temperature andelectron fraction in the equatorial plane and on the plane orthogonalto it after the 10 ms of evolution, and for each case. In order toconsiderably reduce computational costs when running the MonteCarlo code
Sedonu we map density, temperature and the evolvedelectron fraction from M1 onto a 2D, axially symmetric grid, andrun the transport over the obtained profiles. At last, we map theevolved M1 data back to the SPH particles by tri-linear interpolationto run the ASL scheme.
As a first step we identify physical quantities that are directlyimpacted by each ASL parameter:1) the total luminosities and mean energies, primarily affected by 𝛼 𝜈, blk and 𝜏 cut .2) the neutrino flux at different polar angles that an observer MNRAS000
As a first step we identify physical quantities that are directlyimpacted by each ASL parameter:1) the total luminosities and mean energies, primarily affected by 𝛼 𝜈, blk and 𝜏 cut .2) the neutrino flux at different polar angles that an observer MNRAS000 , 1–18 (2020)
SL calibration for BNS merger simulations
100 0 100x [km]1000100 y [ k m ] density(g cm ) 100 0 100x [km]temperature(MeV) 100 0 100x [km]electron fraction
100 0 100x [km]1000100 y [ k m ] density(g cm ) 100 0 100x [km]temperature(MeV) 100 0 100x [km]electron fraction
100 0 100x [km]1000100 y [ k m ] density(g cm ) 100 0 100x [km]temperature(MeV) 100 0 100x [km]electron fraction Figure 1. Left: density, middle: temperature, and right: electron fraction maps on the equatorial plane at equilibrium, for top row : a 1.4-1.4 𝑀 (cid:12) binary neutronstar merger at ≈
38 ms after merger, middle row : a 1.3-1.3 𝑀 (cid:12) binary neutron star merger at ≈
18 ms after merger, bottom row : a 1.2-1.3 𝑀 (cid:12) binary neutronstar merger at ≈
18 ms after merger. From top to bottom the degree of axial symmetry of the remnant decreases. Snapshots are taken from the dynamicalsimulations of Rosswog et al. (2017). located far outside the neutrino surfaces would measure. In thisrespect, we focus on the trend of the flux with the polar angle forconstraining the parameter 𝑝 in Eq. 5.Given that the Monte Carlo approach converges (in the limitof infinite particle numbers) to an exact solution of the transport equation, the best strategy for our parameter calibration would be toextract both quantities 1) and 2) from Sedonu and compare withthe ASL scheme. However, the assumption of axial symmetry maynot be fully appropriate in all cases, especially shortly after merger.We then first test the assumption of azimuthal symmetry by runningthe M1 transport on both 2D and 3D grid setups, and compare the
MNRAS , 1–18 (2020)
Gizzi D. et al.
100 0 100x [km]1000100 z [ k m ] density(g cm ) 100 0 100x [km]temperature(MeV) 100 0 100x [km]electron fraction
100 0 100x [km]1000100 z [ k m ] density(g cm ) 100 0 100x [km]temperature(MeV) 100 0 100x [km]electron fraction
100 0 100x [km]1000100 z [ k m ] density(g cm ) 100 0 100x [km]temperature(MeV) 100 0 100x [km]electron fraction Figure 2.
Same as in Fig. 1, but on the plane orthogonal to the equatorial one. neutrino emission maps for each species. For those snapshots wherethe impact on the emission is large in at least one of the species,we decide to consider M1 in 3D for the calibration. However, giventhe limitations of M1 in accurately modelling the distribution ofthe neutrino flux at different polar angles (Foucart et al. 2018), wemake use of M1 purely for calibrating 𝛼 𝜈, blk and 𝜏 cut , while wekeep Sedonu to calibrate the heating parameter 𝑝 . As we shallsee, the 2D assumption does not impact the trend of the flux atdifferent polar angles, but only its local values. On the other hand,if the 2D assumption does not affect the emission of any neutrinospecies sensitively, the calibration of 𝛼 𝜈, blk , 𝜏 cut and of the heatingparameter 𝑝 is performed entirely with Sedonu .The parameter space we explore is given by 𝛼 𝜈, blk ∈[ . , . , . , . , . , . ] , 𝜏 cut ∈ [ , , , ] ,and 𝑝 ∈ [ , , , , , ] . When calibrating, we assume 𝛼 blk ≡ 𝛼 𝜈 𝑒 , blk = 𝛼 ¯ 𝜈 𝑒 , blk , and 𝛼 𝜈 𝑥 , blk =
0, similarly to Peregoet al. (2016). We will discuss the accuracy of the assumption 𝛼 𝜈 𝑒 , blk = 𝛼 ¯ 𝜈 𝑒 , blk in Sec. 4.2. Given that heavy-lepton neutrinos do not contribute to the heating either, we decide to neglect this speciesentirely when calibrating, and just focus on the other two species.We will anyway report all values of luminosities and mean energiesfor completeness (see Table 2).During the calibration process, when exploring luminosities andmean energies, we first compute: 𝜖 𝐿 = 𝜖 𝐿 𝜈𝑒 + 𝜖 𝐿 ¯ 𝜈𝑒 , (8) 𝜖 (cid:104) 𝐸 (cid:105) = 𝜖 (cid:104) 𝐸 𝜈𝑒 (cid:105) + 𝜖 (cid:104) 𝐸 ¯ 𝜈𝑒 (cid:105) , (9)where 𝜖 𝐿 𝑖 = | 𝐿 𝑖, ref − 𝐿 𝑖, ASL | 𝐿 𝑖, ref (10)and 𝜖 (cid:104) 𝐸 𝑖 (cid:105) = |(cid:104) 𝐸 𝑖, ref (cid:105) − (cid:104) 𝐸 𝑖, ASL (cid:105)|(cid:104) 𝐸 𝑖, ref (cid:105) , (11) 𝑖 labeling neutrinos of species 𝑖 ∈ [ 𝜈 𝑒 , ¯ 𝜈 𝑒 ] , while 𝐿 𝑖, ref and (cid:104) 𝐸 𝑖, ref (cid:105) MNRAS000
0, similarly to Peregoet al. (2016). We will discuss the accuracy of the assumption 𝛼 𝜈 𝑒 , blk = 𝛼 ¯ 𝜈 𝑒 , blk in Sec. 4.2. Given that heavy-lepton neutrinos do not contribute to the heating either, we decide to neglect this speciesentirely when calibrating, and just focus on the other two species.We will anyway report all values of luminosities and mean energiesfor completeness (see Table 2).During the calibration process, when exploring luminosities andmean energies, we first compute: 𝜖 𝐿 = 𝜖 𝐿 𝜈𝑒 + 𝜖 𝐿 ¯ 𝜈𝑒 , (8) 𝜖 (cid:104) 𝐸 (cid:105) = 𝜖 (cid:104) 𝐸 𝜈𝑒 (cid:105) + 𝜖 (cid:104) 𝐸 ¯ 𝜈𝑒 (cid:105) , (9)where 𝜖 𝐿 𝑖 = | 𝐿 𝑖, ref − 𝐿 𝑖, ASL | 𝐿 𝑖, ref (10)and 𝜖 (cid:104) 𝐸 𝑖 (cid:105) = |(cid:104) 𝐸 𝑖, ref (cid:105) − (cid:104) 𝐸 𝑖, ASL (cid:105)|(cid:104) 𝐸 𝑖, ref (cid:105) , (11) 𝑖 labeling neutrinos of species 𝑖 ∈ [ 𝜈 𝑒 , ¯ 𝜈 𝑒 ] , while 𝐿 𝑖, ref and (cid:104) 𝐸 𝑖, ref (cid:105) MNRAS000 , 1–18 (2020)
SL calibration for BNS merger simulations are the luminosity and mean energy of species 𝑖 from the referencesolution (either Sedonu or M1). We compare 𝜖 𝐿 and 𝜖 (cid:104) 𝐸 (cid:105) and even-tually consider the error among the two more sensitive to a variationof the parameters. We look for regions in the parameter space wherethis error is minimal for a first pre-selection, and subsequently ex-plore 𝜖 𝐿 𝑖 and 𝜖 (cid:104) 𝐸 𝑖 (cid:105) individually for a more detailed analysis.When calibrating 𝑝 with Sedonu we assume an analytical form ofthe neutrino fluxes as function of the polar angle in the ASL scheme.In particular, the flux of neutrinos of species 𝑖 leaving the source atsome distance 𝑅 is: F 𝑖 = 𝑅 𝑑𝐿 𝑖 𝑑 Ω n , (12) n being the unit vector orthogonal to the surface 𝑑𝑆 = 𝑅 𝑑 Ω , and 𝑑 Ω = 𝑑 cos ( 𝜃 ) 𝑑𝜙 = 𝑑𝜇𝑑𝜙 being the solid angle of polar angle 𝜃 andazimuthal angle 𝜙 . For axial symmetry with 𝑑 Ω = 𝜋𝑑𝜇 we have: F 𝑖 = 𝜋𝑅 (cid:16) 𝑑𝐿 𝑖 𝑑𝜇 (cid:17) n . (13)Given the analytical angular term in Eq. 5, we assume that: 𝑑𝐿 𝑖 𝑑 Ω = 𝜋 ( + 𝑝 )( + 𝛽 𝑖 cos 𝑝 ( 𝜃 )) + 𝑝 + 𝛽 𝑖 𝐿 𝑖 (14)and 𝑑𝐿 𝑖 𝑑𝜇 = ( + 𝑝 )( + 𝛽 𝑖 cos 𝑝 ( 𝜃 )) + 𝑝 + 𝛽 𝑖 𝐿 𝑖 . (15)We caution the reader that Eq. 15 is purely to understand which valueof 𝑝 best reproduces the trend of the flux as function of the polarangle from Sedonu . We do not attempt to estimate the magnitude ofthe flux at each polar angle via Eq. 15 because it is not possible toextract this information with a leakage scheme.We summarise the values of luminosities and mean energies for thethree cases described in Sec. 3.1 in Table 2. We include the 2D valuesfrom
Sedonu , the 2D and 3D values from M1, and the values fromthe ASL with the corresponding best parameter set when assuming 𝛼 blk ≡ 𝛼 𝜈 𝑒 , blk = 𝛼 ¯ 𝜈 𝑒 , blk . We start by providing in Sec. 4.1 constraints on the parameters foreach case described in Sec. 3.1. Afterwards, we combine them inSec. 4.2 and discuss the performance of using the same blockingparameter for electron neutrinos and anti-neutrinos.
From Table 2 we can see that by comparing the 2D and 3D lumiositiesof electron neutrinos and anti-neutrinos from M1 the impact of the 2Daveraging is only of the order of (cid:46)
7% for both species. This suggeststhe usage of
Sedonu for the comparison with the ASL scheme. Wealso notice that the heavy-lepton neutrino luminosity from M1 in 3Dis ∼ . Sedonu . When comparingusing the ASL with the best parameter set the leakage value is about ∼ . 𝜖 𝐿 and 𝜖 (cid:104) 𝐸 (cid:105) as a function of 𝛼 blk , fordifferent values of 𝑝 and 𝜏 cut . Lines of the same color are for fixed 𝑝 , while lines of the same type are for fixed 𝜏 cut . We neglect thecases with 𝜏 cut = 𝑝 =
12 to reduce the amount of data toshow, but the results do not change. We see that neither 𝜖 𝐿 nor 𝜖 (cid:104) 𝐸 (cid:105) is largely affected by varying 𝑝 for a given 𝛼 blk and 𝜏 cut . The milddependence can be explained by the fact that 𝑝 defines the distributionof the heating rather than its intensity. However, because the neutrinoabsorption is typically more pronounced for electron neutrinos thanfor anti-neutrinos, we should expect a somewhat larger dependence of 𝜖 𝐿 on 𝑝 for the former species. Nevertheless, since we are evaluatingthe species-summed 𝜖 𝐿 , even in case there is such dependence fromthe electron neutrinos, this must be associated to a rather small 𝜖 𝐿 𝜈𝑒 (see later Fig. 4), and it is thus not appreciable. On the other hand,the blocking parameter 𝛼 blk has a major impact on 𝜖 𝐿 . This comesdirectly from the 𝑟 𝜈 → ( − 𝛼 blk ) 𝑟 𝜈 correction to the emission rate 𝑟 𝜈 when accounting for blocking, see Eq. 2. This is different fromthe impact on 𝜖 𝐸 , which is basically negligible. The reason is thefact that (cid:104) 𝐸 𝑖 (cid:105) is computed from the ratio between the luminosity 𝐿 𝑖 [erg/s] and the total number of emitted neutrinos of species 𝑖 per unittime, both affected by blocking for the case of electron neutrinos andanti-neutrinos. At last, a noticeable dependence to 𝜏 cut can be seenfrom 𝜖 (cid:104) 𝐸 (cid:105) , as a consequence of the fact that 𝜏 cut impacts the neutrinospectrum at the decoupling region, and therefore the mean energiesof the species. We find that 𝜖 (cid:104) 𝐸 (cid:105) has a maximum value of ∼ 𝜏 cut =
5, while 𝜖 𝐿 varies from ∼
10% to more than 100%. Forthis reason, we look at 𝜖 𝐿 for a first pre-selection of the parameters.In particular, we find a minimum 𝜖 𝐿 around 𝛼 blk = . − . 𝛼 blk = .
45 and look at both 𝜖 𝐿 𝑖 and 𝜖 (cid:104) 𝐸 𝑖 (cid:105) for eachspecies 𝑖 ∈ [ 𝜈 𝑒 , ¯ 𝜈 𝑒 ] . Fig. 4 shows 𝜖 𝐿 𝑖 (top panels) and 𝜖 (cid:104) 𝐸 𝑖 (cid:105) (bottompanels) in the 𝑝 × 𝜏 cut parameter space. As anticipated before, 𝜖 𝐿 𝜈𝑒 shows the largest dependence on 𝑝 for a given 𝜏 cut , but 𝜖 𝐿 𝜈𝑒 (cid:46) 𝜖 (cid:104) 𝐸 𝑖 (cid:105) (cid:46)
12% for each species 𝑖 . Electron anti-neutrinosshow a somewhat larger 𝜖 𝐿 ¯ 𝜈𝑒 up to ∼ 𝜖 𝐿 𝑖 and 𝜖 (cid:104) 𝐸 𝑖 (cid:105) , any value of 𝜏 cut ∈[ , , , ] seems good enough to describe the thermalization inthis case. We take 𝜏 cut =
10 as reference for the ASL luminosities andmean energies in Table 2. To constrain 𝑝 we show in Fig. 5 𝑑𝐿 𝑖 / 𝑑𝜇 asa function of 𝜃 for 𝑖 = 𝜈 𝑒 (left) and 𝑖 = ¯ 𝜈 𝑒 (right). The curve from M1is obtained by taking an average of 𝑑𝐿 𝑖 / 𝑑𝜇 over different azimuthalangles 𝜙 of the domain for each polar angle 𝜃 , while the ASL curveis obtained by means of Eq. 15, and by taking 𝛼 blk = .
45, a given
MNRAS , 1–18 (2020)
Gizzi D. et al.
Table 2.
Summary of the values of luminosities and mean energies for each species from
Sedonu (2nd column), M1 in 2D (3rd colum), M1 in 3D (4th column)and the ASL (5th column), for each of the three cases (1st column) examined. For the ASL, we also specify the best parameter set resulting from the calibration(6th column) when assuming 𝛼 blk ≡ 𝛼 𝜈 𝑒 , blk = 𝛼 ¯ 𝜈 𝑒 , blk . The reported numbers from Sedonu are accurate in spite of the Monte Carlo noise, which is estimatedto be a thousand times smaller.Case Sedonu (2D) M1(2D) M1(3D) ASL [ 𝛼 blk , 𝜏 cut , 𝑝 ]1 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − [0.45,10,2] (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
24 MeV (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
37 MeV (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
38 MeV (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
70 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
05 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
00 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
10 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
88 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
87 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
27 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
28 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
63 MeV2 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − [0.65,10,2] (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
65 MeV (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
89 MeV (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
11 MeV (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
43 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
74 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
69 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
23 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
61 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
01 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
98 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
05 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
54 MeV3 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑒 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − 𝐿 𝜈 𝑥 = . · erg s − [0.75,10,2] (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
88 MeV (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
60 MeV (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
90 MeV (cid:104) 𝐸 𝜈 𝑒 (cid:105) = .
61 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
17 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
10 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
61 MeV (cid:104) 𝐸 ¯ 𝜈 𝑒 (cid:105) = .
85 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
42 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
39 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
56 MeV (cid:104) 𝐸 𝜈 𝑥 (cid:105) = .
62 MeV blk L blk E p = 2, cut = 10p = 4, cut = 10p = 6, cut = 10p = 8, cut = 10p = 10, cut = 10p = 2, cut = 15p = 4, cut = 15p = 6, cut = 15p = 8, cut = 15p = 10, cut = 15p = 2, cut = 20p = 4, cut = 20p = 6, cut = 20p = 8, cut = 20p = 10, cut = 20 Figure 3. left: 𝜖 𝐿 and right: 𝜖 (cid:104) 𝐸 (cid:105) as a function of 𝛼 blk , for different values of 𝑝 and 𝜏 cut , case 1). Lines of the same color are for fixed 𝑝 , while lines of thesame type are for fixed 𝜏 cut . We neglect the cases with 𝜏 cut = 𝛼 blk has a large impact on 𝜖 𝐿 , but not on 𝜖 (cid:104) 𝐸 (cid:105) . Thethermalization parameter 𝜏 cut has some effect on 𝜖 (cid:104) 𝐸 (cid:105) , while 𝑝 mildly impacts both 𝜖 𝐿 and 𝜖 (cid:104) 𝐸 (cid:105) . 𝜏 cut = and the best parameter 𝑝 from the direct comparisonbetween the polar angle profiles of the luminosities. In particular,we see that 𝑝 = Sedonu for both electron We just use this value of 𝜏 cut to illustrate the calibration of 𝑝 . Other choicesof 𝜏 cut would have been equivalent. neutrinos and anti-neutrinos. We can also see that M1 overestimatesthe flux in the regions 𝜃 < 𝜋 / 𝜃 > 𝜋 / Sedonu by a factor of (cid:46)
MNRAS000
MNRAS000 , 1–18 (2020)
SL calibration for BNS merger simulations p c u t blk = 0.45, electron neutrinos L e p c u t blk = 0.45, electron anti neutrinos L e p c u t blk = 0.45, electron neutrinos E e p c u t blk = 0.45, electron anti neutrinos E e Figure 4. top: 𝜖 𝐿 𝑖 and bottom: 𝜖 (cid:104) 𝐸 𝑖 (cid:105) for 𝑖 = 𝜈 𝑒 (left panels) and 𝑖 = ¯ 𝜈 𝑒 (right panels), case 1). Apart from 𝜖 𝐿 ¯ 𝜈𝑒 which shows values up to ∼ 𝜖 𝐿 𝜈𝑒 aswell as both 𝜖 (cid:104) 𝐸 𝑖 (cid:105) show a limited relative error of (cid:46)
12% at the most. estimate of the parameter 𝑝 = 𝑝 =
8. The reason is the factthat there we chose M1 as source for comparison. As just said, theapproximations introduced by the analytical closure inevitably leadsto a steeper decrease of the neutrino fluxes with the polar angle,consequently suggesting a larger value of 𝑝 than the one found here. Unlike case 1), by comparing the 2D and 3D luminosities of electronneutrinos and anti-neutrinos from M1 in Table 2 we see that the 2Daveraging impacts the electron anti-neutrino luminosity by a factorof 2, while the difference for electron neutrinos is limited to ∼ 𝛼 blk and 𝜏 cut .Fig. 6 shows line plots of 𝜖 𝐿 and 𝜖 (cid:104) 𝐸 (cid:105) . Similarly to case 1), we cali-brate 𝛼 blk by looking at 𝜖 𝐿 only, given that we find a maximum valueof 𝜖 (cid:104) 𝐸 (cid:105) ∼
12% at 𝜏 cut =
5. We observe a remarkable difference inin the best 𝛼 blk with respect to case 1), with 𝛼 blk = .
65 providing 𝜖 𝐿 (cid:46)
20% in most of the cases. 𝛼 blk = .
55 provides a similar per-formance for the cases 𝑝 = 𝛼 blk = .
65 we look at 𝜖 𝐿 𝑖 and 𝜖 (cid:104) 𝐸 𝑖 (cid:105) individually in Fig. 7.Contrary to case 1), the difference in the maximum relative errorbetween 𝜖 𝐿 𝜈𝑒 and 𝜖 𝐿 ¯ 𝜈𝑒 is larger, with 𝜖 𝐿 𝜈𝑒 (cid:46)
23% for some param-eter combination, while 𝜖 𝐿 ¯ 𝜈𝑒 (cid:46) 𝛼 blk to always accommodate accurate luminosities for bothneutrino species. On the other hand, 𝜖 (cid:104) 𝐸 𝑖 (cid:105) (cid:46)
10% for any combi-nation. Nevertheless, the rather limited 𝜖 𝐿 𝑖 and 𝜖 (cid:104) 𝐸 𝑖 (cid:105) point again to any 𝜏 cut ∈ [ , , , ] as rather suitable thermalization parameter.The left and middle panels of Fig. 8 shows 𝑑𝐿 𝑖 / 𝑑𝜇 for 𝛼 blk = . 𝜏 cut =
10. While the plot on the left shows alsoM1 in 3D because of the limited impact of the 2D assumption onthe electron neutrino emission, the plot in the middle shows only the2D result from
Sedonu , and compare it with the ASL for electronanti-neutrinos. We again find 𝑝 = Sedonu , while the performance of M1 on the left is similar to case 1)because of the usage of an analytical closure. We justify our decisionof using the 2D data from
Sedonu of electron anti-neutrinos for thecalibration of 𝑝 = 𝑑𝐿 ¯ 𝜈 𝑒 / 𝑑𝜇 in M1, see plot on the right. By choosing 𝑝 =
2, from the top-left panel in Fig. 7 we see that either 𝜏 cut = 𝜏 cut =
15 or 𝜏 cut =
20 provide 𝜖 𝐿 𝜈𝑒 < 𝜖 𝐿 ¯ 𝜈𝑒 (cid:46) 𝜖 (cid:104) 𝐸 (cid:105) 𝑖 (cid:46) 𝜏 cut =
10 as reference for the values ofluminosities and mean energies in Table 2. The plot on the right ofFig. 8 shows also the quantity 𝑑𝐿 ¯ 𝜈 𝑒 / 𝑑𝜇 as function of the polar angle 𝜃 for different azimuthal angles 𝜙 , obtained with M1 in 3D. We canclearly see that the variation with 𝜙 is limited, particularly at angles 𝜃 (cid:46) 𝜋 / ≈ .
05 where the bulk of neutrino-driven winds is located.This justifies the assumption of axially symmetric fluxes enteringEq. 5.Comparing the parameter set of this case with case 1), only theblocking parameter appears to be significantly different, and largerby ∼ ∼ − MNRAS , 1–18 (2020) Gizzi D. et al. d L e / d ( e r g s ) blk = 0.45, cut = 10, p = 2 M1(3D, avg)ASL,scaledSedonu d L e / d ( e r g s ) blk = 0.45, cut = 10, p = 2 M1(3D, avg)ASL,scaledSedonu
Figure 5. 𝑑𝐿 𝑖 / 𝑑𝜇 as a function of 𝜃 for Left: 𝑖 = 𝜈 𝑒 and Right: 𝑖 = ¯ 𝜈 𝑒 , case 1). Dots represent the exact solution from Sedonu . The blue curve describesthe trend of 𝑑𝐿 𝑖 / 𝑑𝜇 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 Eq. 15 with 𝛼 blk = . 𝜏 cut =
10. The trend from
Sedonu is overall well described by a cos ( 𝜃 ) for both species.M1 tends to overestimate the flux at 𝜃 < 𝜋 / 𝜃 > 𝜋 / from hotter regions with respect to case 1). Since neutrinos thermal-ize with matter before free-streaming, the mean neutrino energy ishigher if temperatures are higher. The presence of hotter regions isalso confirmed by the fact that the maximum temperature seen forcase 1) is ∼
25 MeV, while for case 2) reaches ∼
40 MeV. It is im-portant to consider that the values of luminosities and mean energiesfrom M1 that we have taken as reference might be off by ∼
20% withrespect to an exact solution to the transport (Foucart et al. 2020),implying that our calibrated 𝛼 blk might be slightly affected too.Similarly to case 1), the heavy-lepton neutrino luminosities are sys-tematically higher for both M1 and the ASL, with the latter providinga luminosity off by a factor of 2. Similarly to case 2), by comparing the 2D and 3D luminosities ofelectron neutrinos and anti-neutrinos from M1 in Table 2 we see thatthe 2D averaging implies a reduction in the anti-neutrino luminosityby about a factor of 10. We therefore adopt the values of luminositiesand mean energies from M1 in 3D as fiducial for the calibration.Fig.9 shows line plots of 𝜖 𝐿 and 𝜖 (cid:104) 𝐸 (cid:105) . Clearly, this case is the hardestone for our calibration process. Indeed, values of 𝜖 𝐿 are significantlylarger than for the other two cases, with a minimum of 𝜖 𝐿 ≈ . 𝛼 blk = . 𝜖 (cid:104) 𝐸 (cid:105) is also larger overall, but its values are < 𝛼 blk = .
75 and look at 𝜖 𝐿 𝑖 and 𝜖 (cid:104) 𝐸 𝑖 (cid:105) individually in Fig. 10.The choice of a single 𝛼 blk makes it cumbersome to well catch theluminosities of both electron neutrinos and anti-neutrinos. While 𝜖 𝐿 ¯ 𝜈𝑒 (cid:46) 𝜖 𝐿 𝜈𝑒 is above 40% and reaches ∼
52% in some cases.Mean energies are better recovered, with 𝜖 (cid:104) 𝐸 𝜈𝑒 (cid:105) <
10% and 𝜖 (cid:104) 𝐸 ¯ 𝜈𝑒 (cid:105) (cid:46) 𝜏 cut is best by first constraining 𝑝 .Fig. 11 shows our results. The left and middle panels show 𝑑𝐿 / 𝑑𝜇 for 𝛼 blk = .
75 and a given 𝜏 cut =
10 taken as reference. Similarly tocase 2), we do not show the curve from M1 in 3D for electron anti-neutrinos in the middle plot. In addition, the M1 curve of the plot onthe left is below the
Sedonu data as a result of the difference by abouta factor of 2 in the electron neutrino luminosity (see Table 2), contraryto what is seen for cases 1) and 2). Similarly to cases 1) and 2) we find 𝑝 = 𝑑𝐿 ¯ 𝜈 𝑒 / 𝑑𝜇 as function of 𝜃 , and that the deviation with 𝜙 of 𝑑𝐿 ¯ 𝜈 𝑒 / 𝑑𝜇 at each 𝜃 is rather small, especially at 𝜃 (cid:46) 𝜋 /
3. Setting 𝑝 = 𝜏 cut = 𝜖 𝐿 ¯ 𝜈𝑒 < 𝜖 𝐿 𝜈𝑒 (cid:38) 𝜖 (cid:104) 𝐸 ¯ 𝜈𝑒 (cid:105) < 𝜏 cut =
10 is also reasonable, providing 𝜖 𝐿 ¯ 𝜈𝑒 < 𝜖 𝐿 𝜈𝑒 ≈ 𝜖 (cid:104) 𝐸 ¯ 𝜈𝑒 (cid:105) ∼ 𝜖 (cid:104) 𝐸 𝜈𝑒 (cid:105) (cid:46)
10% regardless of 𝜏 cut . The blocking parameter is nowlarger by ∼
66% with respect to the value of case 1). Moreover,similarly to case 2) the electron neutrino and anti-neutrino meanenergies are higher in both the M1 and the ASL (see Table 2) withrespect to case 1), as a result of higher average emission temperatures(Table 3), confirmed also by the higher maximum temperature seenfor case 3) of ∼
32 MeV. However, the higher electron anti-neutrinomean energy is also a consequence of the rather large 𝜖 ¯ 𝜈 𝑒 (Fig. 10,bottom right panel). As already seen for cases 1) and 2), the heavy-lepton neutrino luminosities are systematically higher than Sedonu for both M1 and the ASL, with the latter being off by a factor of 3.
Among the three calibrated parameters, both 𝜏 cut and 𝑝 do not showvariations by changing the binary configuration. In particular, 𝜏 cut is always around a value of 10, and 𝑝 = 𝛼 blk is more sensitivethan the other parameters when moving from one binary configu-ration to another. Similar to the results of Perego et al. (2016), theblocking parameter may vary in a range [0.45,0.75], depending onthe configuration of the binary and its time after merger. Under theassumption 𝛼 blk ≡ 𝛼 𝜈 𝑒 , blk = 𝛼 ¯ 𝜈 𝑒 , blk , we find larger values of 𝛼 blk for cases 2) and 3). However, assuming 𝛼 𝜈 𝑒 , blk = 𝛼 ¯ 𝜈 𝑒 , blk is not agood choice at early times after merger. In particular, we notice thathigher values of 𝛼 blk do not well capture the total luminosity ofelectron neutrinos (see top-left panels in Figs. 7- 10), pointing tothe need for 𝛼 𝜈 𝑒 , blk < 𝛼 ¯ 𝜈 𝑒 , blk . This conclusion is consistent withthe fact that electron anti-neutrinos are the most emitted species inmerger environments, and therefore likely more affected by blockingeffects. Moreover, we find the electron anti-neutrino gas to have onaverage fewer energy states available to be populated by new emittedneutrinos of the same species. Indeed, by calculating the differencebetween the average degeneracy parameter of electron neutrinosand anti-neutrinos, we find positive values for all the cases examined(see 8th column of Table 3). In light of this, we optimize the choiceof blocking parameter by exploring in Fig 12 𝜖 𝐿 𝑖 for 𝑖 = [ 𝜈 𝑒 , ¯ 𝜈 𝑒 ] and for the three cases examined in Sec. 3 by assuming 𝜏 cut = 𝑝 = 𝛼 𝜈 𝑒 , blk ≠ 𝛼 ¯ 𝜈 𝑒 , blk . We confirm that for each binary the The degeneracy parameter is locally defined as 𝜂 𝑖 = 𝜇 𝑖 / 𝑇 , where 𝜇 𝑖 isthe chemical potential of neutrino of species 𝑖 and 𝑇 is the temperature.MNRAS , 1–18 (2020) SL calibration for BNS merger simulations blk L blk E p = 2, cut = 10p = 4, cut = 10p = 6, cut = 10p = 8, cut = 10p = 10, cut = 10p = 2, cut = 15p = 4, cut = 15p = 6, cut = 15p = 8, cut = 15p = 10, cut = 15p = 2, cut = 20p = 4, cut = 20p = 6, cut = 20p = 8, cut = 20p = 10, cut = 20 Figure 6.
Same as in Fig. 3 but for case 2). Similarly to case 1), 𝛼 blk has a large impact on 𝜖 𝐿 , but not on 𝜖 (cid:104) 𝐸 (cid:105) . The thermalization parameter 𝜏 cut has someeffect on 𝜖 (cid:104) 𝐸 (cid:105) , while 𝑝 mildly impacts both 𝜖 𝐿 and 𝜖 (cid:104) 𝐸 (cid:105) . p c u t blk = 0.65, electron neutrinos L e p c u t blk = 0.65, electron anti neutrinos L e p c u t blk = 0.65, electron neutrinos E e p c u t blk = 0.65, electron anti neutrinos E e Figure 7. top: 𝜖 𝐿 𝑖 and bottom: 𝜖 (cid:104) 𝐸 𝑖 (cid:105) for 𝑖 = 𝜈 𝑒 (left panels) and 𝑖 = ¯ 𝜈 𝑒 (right panels), case 2). Unlike case 1), 𝜖 𝐿 𝜈𝑒 shows now the largest values up to ∼ 𝜖 𝐿 ¯ 𝜈𝑒 (cid:46)
8% and 𝜖 (cid:104) 𝐸 𝑖 (cid:105) (cid:46) Table 3.
Average densities (2nd and 3rd columns), temperatures (4th and 5th columns) and electron fractions (6th and 7th columns) at which electron neutrinosand anti-neutrinos are emitted, and difference between average electron anti-neutrino and electron neutrino degeneracy parameters (8th column), for each of thethree cases examined (1st column). The computation of the averages follows equation 9 of Rosswog & Liebendörfer (2003). For the computation of the neutrinochemical potentials entering the degeneracy parameters, we assume weak equilibrium.Case (cid:104) 𝜌 (cid:105) 𝜈 𝑒 [10 g · cm − ] (cid:104) 𝜌 (cid:105) ¯ 𝜈 𝑒 [10 g · cm − ] (cid:104) 𝑇 (cid:105) 𝜈 𝑒 [MeV] (cid:104) 𝑇 (cid:105) ¯ 𝜈 𝑒 [MeV] (cid:104) 𝑌 𝑒 (cid:105) 𝜈 𝑒 (cid:104) 𝑌 𝑒 (cid:105) ¯ 𝜈 𝑒 (cid:104) 𝜂 (cid:105) ¯ 𝜈 𝑒 − (cid:104) 𝜂 (cid:105) 𝜈 𝑒 , 1–18 (2020) Gizzi D. et al. d L e / d ( e r g s ) blk = 0.65, cut = 10, p = 2 M1(3D, avg)ASL,scaledSedonu 0.5 1.0 1.5 2.0 2.50.20.40.60.81.0 d L e / d ( e r g s ) blk = 0.65, cut = 10, p = 2 ASL,scaledSedonu d L e / d ( e r g s ) M1(3D, avg)M1(2D)M1(3D)
Figure 8. Left: 𝑑𝐿 𝜈 𝑒 / 𝑑𝜇 and Middle: 𝑑𝐿 ¯ 𝜈 𝑒 / 𝑑𝜇 as a function of 𝜃 , case 2). Dots represent the exact solution from Sedonu . The blue curve describes thetrend obtained with M1 when performing an azimuthal average over different 𝜙 at each 𝜃 , while the red curve describes the one from the ASL, scaled by someconstant, and obtained via Eq. 15 with 𝛼 blk = . 𝜏 cut =
10. We only show the trends from
Sedonu and the ASL for anti-neutrinos. The trend from
Sedonu isagain described by a cos ( 𝜃 ) dependence. The performance of M1 is the same as for case 1), as a result of the usage of an analytical closure. Right: 𝑑𝐿 ¯ 𝜈 𝑒 / 𝑑𝜇 as a function of 𝜃 , when running M1 in 3D and azimuthally-averaging at each 𝜃 (blue curve), and in 2D (purple curve), case 2). The 2D assumption does notaffect the trend with 𝜃 . In the same plot we show 𝑑𝐿 ¯ 𝜈 𝑒 / 𝑑𝜇 from M1 in 3D as function of 𝜃 for different angles 𝜙 (green dots). For a given 𝜃 , the limitedvariation of 𝑑𝐿 ¯ 𝜈 𝑒 / 𝑑𝜇 with 𝜙 especially at polar angles relevant for neutrino-driven winds ( 𝜃 (cid:46) 𝜋 /
3) justifies the assumption of axially symmetric fluxes whencalculating the heating rate in the ASL. blk L blk E p = 2, cut = 10p = 4, cut = 10p = 6, cut = 10p = 8, cut = 10p = 10, cut = 10p = 2, cut = 15p = 4, cut = 15p = 6, cut = 15p = 8, cut = 15p = 10, cut = 15p = 2, cut = 20p = 4, cut = 20p = 6, cut = 20p = 8, cut = 20p = 10, cut = 20 Figure 9.
Same as in Figs. 3-6 but for case 3). Unlike cases 1) and 2), values of 𝜖 𝐿 are larger, with a minimum of ≈ .
75 at 𝛼 blk = .
75. Similarly, values of 𝜖 (cid:104) 𝐸 (cid:105) are between 21% and 29%. value of the blocking corresponding to the minimum 𝜖 𝐿 𝑖 is lowerfor 𝑖 = 𝜈 𝑒 than for 𝑖 = ¯ 𝜈 𝑒 . For electron neutrinos, the largest block-ing parameter at which 𝜖 𝐿 𝜈𝑒 is minimum is found in case 2), with 𝛼 𝜈 𝑒 , blk = .
55. By looking at Table 3, this is on one side due tothe larger average emission temperature ( (cid:104) 𝑇 (cid:105) 𝜈 𝑒 = .
51 MeV), forwhich the emission rate is enhanced considerably due to the ∼ 𝑇 dependence for charged-current interactions (Rosswog & Liebendör-fer 2003), and on the other side to the less neutron-rich environment( (cid:104) 𝑌 𝑒 (cid:105) 𝜈 𝑒 = . 𝜖 𝐿 ¯ 𝜈𝑒 is minimum is found forcase 3), with 𝛼 ¯ 𝜈 𝑒 , blk = . − .
85. Again, this is due to both a ratherhot ( (cid:104) 𝑇 (cid:105) ¯ 𝜈 𝑒 = .
41 MeV) and neutron-rich material ( (cid:104) 𝑌 𝑒 (cid:105) 𝜈 𝑒 = . 𝛼 𝜈 𝑒 , blk = . 𝜖 𝐿 𝜈𝑒 (cid:46)
20% for all binaries, and we therefore set it as fidu-cial for this species. Regarding the anti-neutrinos, the variability of 𝜖 𝐿 ¯ 𝜈𝑒 with 𝛼 ¯ 𝜈 𝑒 , blk is larger and makes the choice of the best block-ing parameter more cumbersome. In particular, equal mass binariesprefer 𝛼 ¯ 𝜈 𝑒 , blk = .
55 and 𝛼 ¯ 𝜈 𝑒 , blk = .
65, for cases 1) and 2) respec-tively, while case 3) prefers 𝛼 ¯ 𝜈 𝑒 , blk > .
65. We therefore suggest afiducial value of 𝛼 ¯ 𝜈 𝑒 , blk = .
55 for equal mass binaries, such that 𝜖 𝐿 ¯ 𝜈𝑒 (cid:46) 𝛼 ¯ 𝜈 𝑒 , blk = .
75 for unequal mass binaries, suchthat 𝜖 𝐿 ¯ 𝜈𝑒 < 𝛼 ¯ 𝜈 𝑒 , blk .Beside the thermodynamical, compositional and degeneracy prop-erties of the matter described in Table 3 and determining the extentof Pauli blocking, the increasing value of the blocking parameterfor both electron neutrinos and anti-neutrino species when movingfrom case 1) to case 3) is also consistent with the fact that gener-ally less massive and/or unequal mass binaries produce larger disks(Rosswog et al. 2000; Vincent et al. 2020; Bernuzzi et al. 2020). Asstated in Sec. 2.3, the blocking parameter takes also into account thereduction of neutrino emission due to inward neutrino fluxes in the MNRAS000
75 for unequal mass binaries, suchthat 𝜖 𝐿 ¯ 𝜈𝑒 < 𝛼 ¯ 𝜈 𝑒 , blk .Beside the thermodynamical, compositional and degeneracy prop-erties of the matter described in Table 3 and determining the extentof Pauli blocking, the increasing value of the blocking parameterfor both electron neutrinos and anti-neutrino species when movingfrom case 1) to case 3) is also consistent with the fact that gener-ally less massive and/or unequal mass binaries produce larger disks(Rosswog et al. 2000; Vincent et al. 2020; Bernuzzi et al. 2020). Asstated in Sec. 2.3, the blocking parameter takes also into account thereduction of neutrino emission due to inward neutrino fluxes in the MNRAS000 , 1–18 (2020)
SL calibration for BNS merger simulations p c u t blk = 0.75, electron neutrinos L e p c u t blk = 0.75, electron anti neutrinos L e p c u t blk = 0.75, electron neutrinos E e p c u t blk = 0.75, electron anti neutrinos E e Figure 10. top: 𝜖 𝐿 𝑖 and bottom: 𝜖 (cid:104) 𝐸 𝑖 (cid:105) for 𝑖 = 𝜈 𝑒 (left panels) and 𝑖 = ¯ 𝜈 𝑒 (right panels), case 3). Unlike the previous cases, 𝜖 𝐿 𝑖 for both species show largervalues, with 𝜖 𝐿 𝜈𝑒 reaching ∼
52% for some parameter combination. 𝜖 (cid:104) 𝐸 ¯ 𝜈𝑒 (cid:105) is also larger, with values up to ∼ 𝜖 (cid:104) 𝐸 𝜈𝑒 (cid:105) <
10% for anycombination. d L e / d ( e r g s ) blk = 0.75, cut = 10, p = 2 M1(3D, avg)ASL,scaledSedonu 0.5 1.0 1.5 2.0 2.50.51.01.52.02.53.0 d L e / d ( e r g s ) blk = 0.75, cut = 10, p = 2 ASL,scaledSedonu d L e / d ( e r g s ) M1(3D, avg)M1(2D)M1(3D)
Figure 11. Left: 𝑑𝐿 𝜈 𝑒 / 𝑑𝜇 and Middle: 𝑑𝐿 ¯ 𝜈 𝑒 / 𝑑𝜇 as a function of 𝜃 , case 3). Dots represent the exact solution from Sedonu . The blue curve describes thetrend obtained with M1 when performing an azimuthal average over different 𝜙 at each 𝜃 , while the red curve describes the one from the ASL, scaled by someconstant, and obtained via Eq. 15 with 𝛼 blk = . 𝜏 cut =
10. We only show the trend from
Sedonu and the ASL for anti-neutrinos. The trend from
Sedonu is again described by a cos ( 𝜃 ) dependence. The M1 curve for electron neutrinos shows lower values than Sedonu because of a factor of about 2 of differencein the total emission (see Table 2).
Right: 𝑑𝐿 ¯ 𝜈 𝑒 / 𝑑𝜇 as a function of 𝜃 , when running M1 in 3D and azimuthally-averaging at each 𝜃 (blue curve), and in 2D(purple curve), case 3). The 2D assumption does not affect the trend of the curve, as already seen for case 2). Similarly, for a given 𝜃 we see a limited variationof 𝑑𝐿 ¯ 𝜈 𝑒 / 𝑑𝜇 with 𝜙 at polar angles relevant for neutrino-driven winds ( 𝜃 (cid:46) 𝜋 / semi-transparent regime. The presence of larger disks in less massiveand/or unequal mass binaries leads to larger neutrino surfaces, andconsequently to an overall larger effect of inward neutrino fluxes,therefore contributing to an increase in the size of blocking.We conclude by showing in Figs. 13-14 the distribution of therate of change of the specific matter internal energy (cid:164) 𝑒 (units of10 erg g − s − ) and of the electron fraction (cid:164) 𝑌 e in the windsfor case 1) and for the different transport approaches, assuming 𝛼 𝜈 𝑒 , blk = . , 𝛼 ¯ 𝜈 𝑒 , blk = . , 𝜏 cut = , 𝑝 =
2. The upper plotsare 3D maps in a box of 𝑥 × 𝑦 × 𝑧 =
150 km ×
150 km ×
150 km,while the lower ones are projections on the x-y plane. We assume as threshold density below which we identify the wind region a value of5 · g cm − . However, this limit in the end also includes the outerregions of the disk, visible as (cid:164) 𝑒 < (cid:164) 𝑒 and (cid:164) 𝑌 e of each SPH particle from M1 and Sedonu via interpolation fromthe respective grids. Overall, both the ASL and the M1 (cid:164) 𝑒 and (cid:164) 𝑌 e dis-tributions agree well with the solution from Sedonu . However, theyboth show some cooling ( (cid:164) 𝑒 <
0) above the remnant which is weakeror absent in
Sedonu , as well as a few particles with (cid:164) 𝑌 e <
0. Moreprecisely, for the ASL ≈
14% of the particles have (cid:164) 𝑒 ASL · (cid:164) 𝑒 Sed < ≈
3% have (cid:164) 𝑌 𝑒, ASL · (cid:164) 𝑌 𝑒, Sed <
0. For the M1, these numbers are
MNRAS , 1–18 (2020) Gizzi D. et al. ≈
27% and ≈
1% respectively. Considering the overall similaritybetween the maps of the ASL and M1 we can conclude within thelimits of this analysis that the ASL may show similar performancesof M1 in dynamical simulations when comparing the wind propertiesagainst exact solutions to the transport. Nevertheless, a more robustassessment requires to explore this comparison dynamically and fordifferent binary configurations. The main advantage of our SPH-ASLis in its efficiency. We indeed estimate that, if we had to assume thesame timestep as in M1 for a dynamical simulation, and by takingMAGMA2 as our Lagrangian hydrodynamics code (Rosswog 2020),the ratio between the CPU hours spent for the transport and those forthe hydrodynamics is about 0.8 per timestep, to be compared witha factor of 10 for the M1 in
FLASH . The number we find is similarto the one in Pan et al. (2018). The ratio for the SPH-ASL could befurther cut down if we consider that the optical depth computationmay not be required at every timestep, unless the thermodynami-cal and compositional properties of the matter change considerably.Future dynamical simulation will definitely provide a more robustassessment of the performance.
In this paper we have provided a detailed calibration analysis ofthe Advanced Spectral Leakage (ASL) scheme presented earlier inGizzi et al. (2019), which is based on the original work of Peregoet al. (2016). Our main motivation is the study of neutrino-drivenwinds emerging from neutron star mergers. The gauging processis performed by post-processing a number of snapshots of abinary neutron star merger remnant. We extract neutrino quantitiesdirectly impacted by each parameter, and we compare the ASLresults for different parameter combinations with the ones obtainedfrom the Monte Carlo neutrino transport code
Sedonu (Richerset al. 2015) and from a two-moment scheme (M1) implementedin
FLASH (Fryxell et al. 2000; O’Connor 2015; O’Connor &Couch 2018). In the calibration process we focus on electron-typeneutrinos and anti-neutrinos, since they determine the properties ofneutrino-driven winds. We summarise our main findings as follows: • Performing neutrino transport in 2D by post-processing ini-tial 3D post-merger configurations is of limited accuracy at early( 𝑡 (cid:46)
20 ms) times post-merger. In particular, the 2D averaging canseverely impact the total luminosities by more than a factor of 2. • The assumption of axially symmetric neutrino fluxes enter-ing the heating rate in the ASL is validated by 3D neutrino transportsimulations. In particular, variations of the fluxes with the azimuthalangle 𝜙 and for a given polar angle 𝜃 are limited in the bulk regionof neutrino-driven winds. • In agreement with Perego et al. (2016), the thermalizationparameter 𝜏 cut has an impact mainly on the neutrino mean energies.However, a value of 𝜏 cut =
10 robustly recovers neutrino meanenergies within (cid:46)
25% accuracy. • The heating parameter introduced in Gizzi et al. (2019) isre-calibrated to a lower value, as a result of the usage of
Sedonu asreference solution. We find 𝑝 = Sedonu would lead to the artifact of a larger 𝑝 than theones calibrated here because of the approximations introduced bythe analytical closure. • The blocking parameter 𝛼 blk mainly impacts the total neu-trino luminosities, in agreement with the results of Perego et al.(2016). Moreover, unlike the other two parameters it is the mostsensitive to a change of binary configuration. • The assumption 𝛼 blk ≡ 𝛼 𝜈 𝑒 , blk = 𝛼 ¯ 𝜈 𝑒 , blk adopted in Perego et al.(2016) can be rather inaccurate for recovering neutrino luminosities.In particular, electron neutrino luminosities are lower by up toa factor of 2 with respect to the reference solution, suggesting 𝛼 𝜈 𝑒 , blk < 𝛼 ¯ 𝜈 𝑒 , blk . • Assuming 𝛼 𝜈 𝑒 , blk ≠ 𝛼 ¯ 𝜈 𝑒 , blk , we indeed find that 𝛼 𝜈 𝑒 , blk = . (cid:46)
20% accuracy, for both equaland unequal mass binaries. Similarly, 𝛼 ¯ 𝜈 𝑒 , blk = .
55 results in elec-tron anti-neutrino luminosities off by (cid:46)
20% for equal mass binaries. • We find 𝛼 ¯ 𝜈 𝑒 , blk = .
75 for unequal mass binaries, leadingto a relative error in the anti-neutrino luminosity of < • In contrast to Foucart et al. (2020), the heavy-lepton neu-trino luminosity is systematically larger by a factor of a few inboth the ASL and M1 with respect to an exact solution to thetransport. This enhances the overall cooling of the remnant in oursnapshot calculations. For the ASL, the most probable explanationis the poor treatment of the diffusion timescale, which according toArdevol-Pulpillo et al. (2019) can boost luminosities by more thana factor of 2. We expect this treatment to affect also the electronneutrino and anti-neutrino luminosities to some extent. • The properties of neutrino-driven winds are shaped by therates of change of internal energy and electron fraction, (cid:164) 𝑒 and (cid:164) 𝑌 𝑒 ,respectively. The corresponding maps for a 1.4-1.4 𝑀 (cid:12) binarydemonstrate that for our suggested parameter choice the ASLscheme performs similar to the M1 approach. From the perspectiveof dynamical simulations, our SPH-ASL comes with the advantageof a better efficiency. In particular, by taking the Lagrangianhydrodynamics code MAGMA2 (Rosswog 2020), we estimate thatthe ratio between the CPU hours spent for the transport and thosespent for the hydrodynamics would be (cid:46) . FLASH is about 10. In other words, the ASL schemecould be applied in SPH simulations with only a moderate additionalcomputational effort.Although the geometry of a binary neutron star merger al-lows neutrinos to escape with more directional freedom than in aspherically symmetric core-collapse supernova, the results we findhere show that the blocking parameter can still be quite high inmerger remnants. Apart from the different thermodynamics andcomposition of the matter, the major reason is the disk geometrythat increases the effect of inward neutrino fluxes at the neutrinosurfaces with respect to a spherically symmetric geometry. In thispaper we have also presented a completely mesh-free, particle-basedalgorithm to compute spectral, species-dependent optical depths,based on the Smoothed-Particle Hydrodynamics (SPH) (Monaghan1992; Monaghan 2005; Rosswog 2009, 2015a,b, 2020). Thisalgorithm makes our ASL fully grid-independent, and thereforesuitable for future SPH dynamical simulations of binary neutron starmergers with neutrino transport.
MNRAS000
MNRAS000 , 1–18 (2020)
SL calibration for BNS merger simulations , blk L ,blk ,blk Figure 12. 𝜖 𝐿 𝜈 as a function of 𝛼 𝜈, blk ∈ [ . , . , . , . , . , . ] , for the three binary setups and for both electron neutrinos (green curve) andanti-neutrinos (orange curve). Plots are obtained with 𝜏 cut =
10 for all binaries, and with 𝑝 = 𝑝 = 𝜖 𝐿 𝑖 is lower for 𝑖 = 𝜈 𝑒 than for 𝑖 = ¯ 𝜈 𝑒 . While for 𝛼 𝜈 𝑒 , blk = . 𝜖 𝐿 𝜈𝑒 (cid:46)
20% in allbinaries, for the anti-neutrinos the larger variability of 𝜖 𝐿 ¯ 𝜈𝑒 with 𝛼 ¯ 𝜈 𝑒 , blk when exploring different cases leads to a different fiducial value for equal and unequalmass binaries. In particular, 𝛼 ¯ 𝜈 𝑒 , blk = .
55 for the former case, and 𝛼 ¯ 𝜈 𝑒 , blk = .
75 for the latter case. In this way, 𝜖 𝐿 ¯ 𝜈𝑒 (cid:46)
20% and 𝜖 𝐿 ¯ 𝜈𝑒 <
30% respectively. x ( k m ) y ( k m ) z ( k m ) e Sed x ( k m ) y ( k m ) z ( k m ) e ASL x ( k m ) y ( k m ) z ( k m ) e M1
150 100 50 0 50 100 150 x[km] y [ k m ] e Sed
150 100 50 0 50 100 150 x[km] y [ k m ] e ASL
150 100 50 0 50 100 150 x[km] y [ k m ] e M1 Figure 13. Top:
3D maps bottom: projection of the 3D maps on the x-y plane of the rate of change of the internal energy (cid:164) 𝑒 in the winds (units of 10 ergg − s − ) for left: Sedonu , middle: the ASL, right: M1, case 1). The ASL performance is similar to M1, and they both generally reproduce the distribution of (cid:164) 𝑒 from Sedonu .The main difference with respect to
Sedonu is in somewhat lower (cid:164) 𝑒 > (cid:164) 𝑒 < 𝑥 × 𝑦 × 𝑧 =
150 km ×
150 km ×
150 km.The map for the ASL is obtained with the parameter set 𝛼 𝜈 𝑒 , blk = . , 𝛼 ¯ 𝜈 𝑒 , blk = . , 𝜏 cut = , 𝑝 = ACKNOWLEDGEMENTS
We would like to thank Sherwood Richers for developing andsupporting our use of
Sedonu and for a careful reading of themanuscript. This work has been supported by the Swedish ResearchCouncil (VR) under grant number 2016- 03657_3, by the SwedishNational Space Board under grant number Dnr. 107/16, by the Re-search Environment grant "Gravitational Radiation and Electromag-netic Astrophysical Transients (GREAT)" funded by the SwedishResearch council (VR) under Dnr 2016-06012 and by the Knut andAlice Wallenberg Foundation under Dnr KAW 2019.0112. EOC is supported by the Swedish Research Council (Project No. 2018-04575and 2020-00452). We gratefully acknowledge support from COSTAction CA16104 "Gravitational waves, black holes and fundamen-tal physics" (GWverse), from COST Action CA16214 "The multi-messenger physics and astrophysics of neutron stars" (PHAROS),and from COST Action MP1304 "Exploring fundamental physicswith compact stars (NewCompStar)".The simulations were performed on resources provided by theSwedish National Infrastructure for Computing (SNIC) at PDC (Cen-tre for High Performance Computing) and on the facilities of the
MNRAS , 1–18 (2020) Gizzi D. et al. x ( k m ) y ( k m ) z ( k m ) Y e,Sed x ( k m ) y ( k m ) z ( k m ) Y e,ASL x ( k m ) y ( k m ) z ( k m ) Y e,M1
150 100 50 0 50 100 150 x[km] y [ k m ] Y e,Sed
150 100 50 0 50 100 150 x[km] y [ k m ] Y e,ASL
150 100 50 0 50 100 150 x[km] y [ k m ] Y e,M1 Figure 14. Top:
3D maps and bottom: projection of the 3D maps on the x-y plane of the rate of change of the electron fraction (cid:164) 𝑌 e in the winds for left: Sedonu , middle: the ASL, right: M1. The ASL performance is very similar to M1, and they both generally reproduce the distribution of (cid:164) 𝑌 e from Sedonu . The latterprovides always (cid:164) 𝑌 e >
0, while both the ASL and M1 show a few particles with (cid:164) 𝑌 e <
0. Moreover, (cid:164) 𝑌 e > Sedonu appears stronger than in the ASL and M1.The domain size is 𝑥 × 𝑦 × 𝑧 = × × 𝛼 𝜈 𝑒 , blk = . , 𝛼 ¯ 𝜈 𝑒 , blk = . , 𝜏 cut = , 𝑝 = North-German Supercomputing Alliance (HLRN) in Göttingen andBerlin.
DATA AVAILABILITY
The data concerning the initial conditions needed to run the simula-tions, as well as the extracted neutrino quantities from all neutrinotransport codes and neutrino-driven wind maps are available fromthe authors of this manuscript upon request.
REFERENCES
Abbott B. P., et al., 2017a, Physical Review Letters, 119, 161101Abbott B. P., et al., 2017b, Nature, 551, 85Abbott B. P., et al., 2017c, ApJ, 848, L12Abbott B. P., et al., 2018, Physical Review Letters, 121, 161101Abdikamalov E., et al., 2012, The Astrophysical Journal, 755, 111Arcavi I., et al., 2017, Nature, 551, 64Ardevol-Pulpillo R., Janka H. T., Just O., Bauswein A., 2019, MNRAS, 485,4754Audit E., et al., 2002, arXiv e-prints, pp astro–ph/0206281Bauswein A., Just O., Janka H.-T., Stergioulas N., 2017, The AstrophysicalJournal, 850, L34Bernuzzi S., et al., 2020, MNRAS, 497, 1488Bruenn S. W., 1985, ApJS, 58, 771Bruenn S. W., Buchler J. R., Yueh W. R., 1978, Ap&SS, 59, 261Burrows A., Thompson T. A., 2004, in Fryer C. L., ed., Astrophysics and SpaceScience Library Vol. 302, Astrophysics and Space Science Library. pp133–174, doi:10.1007/978-0-306-48599-2_5 Cabezón R. M., Pan K.-C., Liebendörfer M., Kuroda T., Ebinger K., Heini-mann O., Perego A., Thielemann F.-K., 2018, A&A, 619, A118Castor J. I., 2004, Radiation Hydrodynamics. Cambridge University Press,doi:10.1017/CBO9780511536182Chornock R., et al., 2017, ApJL, 848, L19Ciolfi R., Kalinani J. V., 2020, ApJ, 900, L35Ciolfi R., Kastaun W., Giacomazzo B., Endrizzi A., Siegel D. M., Perna R.,2017, Phys. Rev. D, 95, 063016Coughlin M. W., Dietrich T., Margalit B., Metzger B. D., 2019, MNRAS,489, L91Coulter D. A., et al., 2017, Science, 358, 1556Curtis S., et al., 2019, ApJ, 870, 2De S., Finstad D., Lattimer J., et al., 2018, Phys. Rev. Lett., 121, 091102Dessart L., Ott C. D., Burrows A., Rosswog S., Livne E., 2009, ApJ, 690,1681Dhawan S., Bulla M., Goobar A., et al., 2020, The Astrophysical Journal,888, 67Drout M. R., et al., 2017, Science, 358, 1570Ebinger K., et al., 2020a, ApJ, 888, 91Ebinger K., et al., 2020b, ApJ, 888, 91Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126Endrizzi A., et al., 2020, European Physical Journal A, 56, 15Evans P. A., et al., 2017, Science, 358, 1565Fernández R., Metzger B. D., 2013, MNRAS, 435, 502Fernández R., Tchekhovskoy A., Quataert E., Foucart F., Kasen D., 2019,MNRAS, 482, 3373Foucart F., O’Connor E., Roberts L., Kidder L. E., Pfeiffer H. P., ScheelM. A., 2016, Phys. Rev. D, 94, 123016Foucart F., Duez M. D., Kidder L. E., Nguyen R., Pfeiffer H. P., Scheel M. A.,2018, Phys. Rev. D, 98, 063007Foucart F., Duez M. D., Hebert F., Kidder L. E., Pfeiffer H. P., Scheel M. A.,2020, ApJ, 902, L27Freiburghaus C., Rosswog S., Thielemann F.-K., 1999, ApJ, 525, L121MNRAS000
Abbott B. P., et al., 2017a, Physical Review Letters, 119, 161101Abbott B. P., et al., 2017b, Nature, 551, 85Abbott B. P., et al., 2017c, ApJ, 848, L12Abbott B. P., et al., 2018, Physical Review Letters, 121, 161101Abdikamalov E., et al., 2012, The Astrophysical Journal, 755, 111Arcavi I., et al., 2017, Nature, 551, 64Ardevol-Pulpillo R., Janka H. T., Just O., Bauswein A., 2019, MNRAS, 485,4754Audit E., et al., 2002, arXiv e-prints, pp astro–ph/0206281Bauswein A., Just O., Janka H.-T., Stergioulas N., 2017, The AstrophysicalJournal, 850, L34Bernuzzi S., et al., 2020, MNRAS, 497, 1488Bruenn S. W., 1985, ApJS, 58, 771Bruenn S. W., Buchler J. R., Yueh W. R., 1978, Ap&SS, 59, 261Burrows A., Thompson T. A., 2004, in Fryer C. L., ed., Astrophysics and SpaceScience Library Vol. 302, Astrophysics and Space Science Library. pp133–174, doi:10.1007/978-0-306-48599-2_5 Cabezón R. M., Pan K.-C., Liebendörfer M., Kuroda T., Ebinger K., Heini-mann O., Perego A., Thielemann F.-K., 2018, A&A, 619, A118Castor J. I., 2004, Radiation Hydrodynamics. Cambridge University Press,doi:10.1017/CBO9780511536182Chornock R., et al., 2017, ApJL, 848, L19Ciolfi R., Kalinani J. V., 2020, ApJ, 900, L35Ciolfi R., Kastaun W., Giacomazzo B., Endrizzi A., Siegel D. M., Perna R.,2017, Phys. Rev. D, 95, 063016Coughlin M. W., Dietrich T., Margalit B., Metzger B. D., 2019, MNRAS,489, L91Coulter D. A., et al., 2017, Science, 358, 1556Curtis S., et al., 2019, ApJ, 870, 2De S., Finstad D., Lattimer J., et al., 2018, Phys. Rev. Lett., 121, 091102Dessart L., Ott C. D., Burrows A., Rosswog S., Livne E., 2009, ApJ, 690,1681Dhawan S., Bulla M., Goobar A., et al., 2020, The Astrophysical Journal,888, 67Drout M. R., et al., 2017, Science, 358, 1570Ebinger K., et al., 2020a, ApJ, 888, 91Ebinger K., et al., 2020b, ApJ, 888, 91Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126Endrizzi A., et al., 2020, European Physical Journal A, 56, 15Evans P. A., et al., 2017, Science, 358, 1565Fernández R., Metzger B. D., 2013, MNRAS, 435, 502Fernández R., Tchekhovskoy A., Quataert E., Foucart F., Kasen D., 2019,MNRAS, 482, 3373Foucart F., O’Connor E., Roberts L., Kidder L. E., Pfeiffer H. P., ScheelM. A., 2016, Phys. Rev. D, 94, 123016Foucart F., Duez M. D., Kidder L. E., Nguyen R., Pfeiffer H. P., Scheel M. A.,2018, Phys. Rev. D, 98, 063007Foucart F., Duez M. D., Hebert F., Kidder L. E., Pfeiffer H. P., Scheel M. A.,2020, ApJ, 902, L27Freiburghaus C., Rosswog S., Thielemann F.-K., 1999, ApJ, 525, L121MNRAS000 , 1–18 (2020)
SL calibration for BNS merger simulations Fryxell B., et al., 2000, ApJS, 131, 273Fujibayashi S., et al., 2018, The Astrophysical Journal, 860, 64Fujibayashi S., et al., 2020, Phys. Rev. D, 101, 083029Gafton E., Rosswog S., 2011, MNRAS, 418, 770George M., et al., 2020, arXiv e-prints, p. arXiv:2009.04046Gizzi D., O’Connor E., Rosswog S., Perego A., Cabezón R. M., Nativi L.,2019, MNRAS, 490, 4211Goldstein A., et al., 2017, ApJ, 848, L14Hannestad S., Raffelt G., 1998, The Astrophysical Journal, 507, 339Janka H. T., 1992, A&A, 256, 452Janka H. T., Hillebrandt W., 1989a, A&AS, 78, 375Janka H. T., Hillebrandt W., 1989b, A&A, 224, 49Jiang J.-L., et al., 2019, The Astrophysical Journal, 885, 39Jiang J.-L., et al., 2020, The Astrophysical Journal, 892, 55Kasen D., Badnell N. R., Barnes J., 2013, ApJ, 774, 25Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Nature,551, 80Kasliwal M. M., Nakar E., Singer L. P., Kaplan D. L., et al. 2017, Science,358, 1559Keil M. T., Raffelt G. G., Janka H.-T., 2003, ApJ, 590, 971Kilpatrick C. D., et al., 2017, Science, 358, 1583Kiuchi K., Kyutoku K., Shibata M., Taniguchi K., 2019, The AstrophysicalJournal, 876, L31Korobkin O., Rosswog S., Arcones A., Winteler C., 2012, MNRAS, 426,1940Lattimer J. M., Schramm D. N., 1974, ApJ, (Letters), 192, L145Levermore C. D., Pomraning G. C., 1981, ApJ, 248, 321Liebendörfer M., Whitehouse S. C., Fischer T., 2009, The AstrophysicalJournal, 698, 1174Lindquist R. W., 1966, Annals of Physics, 37, 487Lippuner J., Roberts L. F., 2015, ApJ, 815, 82Lipunov V. M., et al., 2017, ApJ, 850, L1Martin D., et al., 2015, The Astrophysical Journal, 813, 2Mezzacappa A., Messer O. E. B., 1999, Journal of Computational and AppliedMathematics, 109, 281Minerbo G. N., 1978, J. Quant. Spectrosc. Radiative Transfer, 20, 541Monaghan J. J., 1992, Ann. Rev. Astron. Astrophys., 30, 543Monaghan J. J., 2005, Reports on Progress in Physics, 68, 1703Most E. R., Weih L. R., Rezzolla L., Schaffner-Bielich J., 2018, PhysicalReview Letters, 120, 261103Narayan R., Paczynski B., Piran T., 1992, The Astrophysical Journal, 395,L83Nedora V., Bernuzzi S., Radice D., Perego A., Endrizzi A., Ortiz N., 2019,ApJ, 886, L30Nedora V., et al., 2020, arXiv e-prints, p. arXiv:2008.04333O’Connor E., 2015, ApJS, 219, 24O’Connor E. P., Couch S. M., 2018, The Astrophysical Journal, 854, 63O’Connor E., Ott C. D., 2010, Classical and Quantum Gravity, 27, 114103O’Connor E., Ott C. D., 2013, ApJ, 762, 126O’Connor E., et al., 2018, Journal of Physics G Nuclear Physics, 45, 104001Obergaulinger M., Janka H. T., Aloy M. A., 2014, MNRAS, 445, 3169Paczynski B., 1986, The Astrophysical Journal, 308, L43Pan K.-C., Mattes C., O’Connor E. P., Couch S. M., Perego A., Arcones A.,2018, Journal of Physics G: Nuclear and Particle Physics, 46, 014001Perego A., Rosswog S., Cabezón R. M., Korobkin O., Käppeli R., ArconesA., Liebendörfer M., 2014a, MNRAS, 443, 3134Perego A., Gafton E., Cabezón R., Rosswog S., Liebendörfer M., 2014b,A&A, 568, A11Perego A., Cabezón R. M., Käppeli R., 2016, ApJS, 223, 22Perego A., Bernuzzi S., Radice D., 2019, European Physical Journal A, 55,124Pian E., et al., 2017, Nature, 551, 67Pons J. A., Ibáñez J. M., Miralles J. A., 2000, MNRAS, 317, 550Radice D., Dai L., 2019, European Physical Journal A, 55, 50Radice D., Perego A., Zappa F., Bernuzzi S., 2018a, ApJ, 852, L29Radice D., Perego A., Hotokezaka K., Fromm S. A., Bernuzzi S., RobertsL. F., 2018b, The Astrophysical Journal, 869, 130Richers S., 2020, Phys. Rev. D, 102, 083017 Richers S., Kasen D., O’Connor E., Fernández R., Ott C. D., 2015, ApJ, 813,38Rosswog S., 2009, New Astron. Rev., 53, 78Rosswog S., 2015a, Living Reviews of Computational Astrophysics (2015),1Rosswog S., 2015b, MNRAS, 448, 3628Rosswog S., 2020, MNRAS, 498, 4230Rosswog S., Liebendörfer M., 2003, MNRAS, 342, 673Rosswog S., Liebendörfer M., Thielemann F.-K., Davies M., Benz W., PiranT., 1999, A & A, 341, 499Rosswog S., et al., 2000, A&A, 360, 171Rosswog S., Feindt U., Korobkin O., Wu M. R., Sollerman J., Goobar A.,Martinez-Pinedo G., 2017, Classical and Quantum Gravity, 34, 104001Rosswog S., Sollerman J., Feindt U., Goobar A., Korobkin O., Wollaeger R.,Fremling C., Kasliwal M. M., 2018, A&A, 615, A132Savchenko V., et al., 2017, ApJ, 848, L15Siegel D. M., Metzger B. D., 2017, Phys. Rev. Lett., 119, 231102Skinner M. A., Dolence J. C., Burrows A., Radice D., Vartanyan D., 2019,ApJS, 241, 7Smartt S. J., et al., 2017, Nature, 551, 75Smit J. M., Cernohorsky J., Dullemond C. P., 1997, A&A, 325, 203Soares-Santos M., et al., 2017, The Astrophysical Journal, 848, L16Steiner A. W., Hempel M., Fischer T., 2013, ApJ, 774, 17Tanaka M., Hotokezaka K., 2013, ApJ, 775, 113Vincent T., et al., 2020, Phys. Rev. D, 101, 044053
APPENDIX A: PARTICLE-BASED OPTICAL DEPTH
Smoothed Particle Hydrodynamics (SPH) is a mesh-free, particle-based method to solve the ideal hydrodynamics equations, see Mon-aghan (2005); Rosswog (2009, 2015a,b, 2020) for reviews of andrecent developments within this method. In the following, we labeleach SPH particle and related physical quantities by means of label a .Each particle 𝑎 has a smoothing length ℎ 𝑎 , which defines the parti-cle’s local interaction radius. In particular, given two particles a and b ,the latter is a neighbour of the former if the distance Δ 𝑟 𝑎𝑏 = | x 𝑎 − x 𝑏 | between a and b satisfies Δ 𝑟 𝑎𝑏 ≤ ℎ 𝑎 . From Eq. 1 we can define theoptical depth at x 𝑎 as: 𝜏 𝜈 ( 𝐸, x 𝑎 ) = ∫ 𝛾 : x 𝑎 →+∞ 𝜆 𝜈 ( 𝐸, x’ ( 𝑠 )) d 𝑠, (A1)or equivalently: 𝜏 𝜈 ( 𝐸, x 𝑎 ) = ∫ 𝛾 : x 𝑎 →+∞ 𝜅 ( 𝐸, x’ ( 𝑠 )) 𝜌 ( x’ ( 𝑠 )) d 𝑠, (A2)where 1 / 𝜆 𝜈 ( 𝐸, x’ ( 𝑠 )) = 𝜅 ( 𝐸, x’ ( 𝑠 )) 𝜌 ( x’ ( 𝑠 )) , being 𝜅 ( 𝐸, x’ ) theopacity of energy 𝐸 at position x’ and 𝜌 ( x’ ) the density at posi-tion x’ . Obviously, the optical depth at a given location depends onthe exact path 𝛾 along which the neutrino escapes. Neutrinos that areproduced inside a neutron star merger remnant will typically escapefrom the remnant after a potentially large number of interactions.While each radiation particle follows an individual path, particleshave larger probabilities to escape if they are locally moving in adirection of increasing mean free path. In our meshfree algorithm toefficiently calculate optical depths we try to mimick such a behaviour.We approximate the integral in Eq. (A2) as a sum of contributionsbetween pairs of SPH-particles so that we move from SPH-particleto SPH-particle until a transparent matter region has been reached.The choice of the next SPH particle on the route to escape, is guidedby the aim to pick the direction that maximizes the local mean freepath. In practice we use the mass density as a proxy, i.e. each SPHparticle chooses as the next particle the one in its neighbour list that MNRAS , 1–18 (2020) Gizzi D. et al. has the smallest density. In other words we are following the localnegative density gradient until neutrinos can escape freely.More specifically, we proceed via the following steps:(i) The SPH particle for which we search the optical depth islabelled 𝑎 . For the convenience of the subsequent notation we alsoidentify particle 𝑎 as the "zeroth neighbour particle" and label it 𝑏 .(ii) The first task is to find the particle 𝑏 in 𝑎 ’s neighbour listthat has the minimum density among the neighbours and add thecontribution, 𝜏 𝑏 → 𝑏 to the optical depth of 𝑎 .(iii) Once 𝑏 has been found, find particle 𝑏 in 𝑏 ’s neighbourlist that has the minimum density and so on until a particle 𝑏 𝑘 ∞ hasbeen found from which neutrinos can escape freely. We use as escapecondition 𝜌 (cid:46) g / cm (Endrizzi et al. 2020).In other words, we discretize the integral of Eq. A1 as: 𝜏 𝜈,𝑎 ( 𝐸 ) = 𝑘 esc ∑︁ 𝑘 = 𝜏 𝑏 𝑘 → 𝑏 𝑘 + ( 𝐸 ) , (A3)where 𝜏 𝑏 𝑘 → 𝑏 𝑘 + ( 𝐸 ) = ∫ (cid:174) 𝑥 𝑏𝑘 + (cid:174) 𝑥 𝑏𝑘 𝑑𝑠𝜆 𝜈 ( 𝐸, x’ ( 𝑠 ))≈ |(cid:174) 𝑥 𝑏 𝑘 + − (cid:174) 𝑥 𝑏 𝑘 | (cid:18) 𝜆 𝑏 𝑘 + ( 𝐸 ) + 𝜆 𝑏 𝑘 ( 𝐸 ) (cid:19) (A4)and we have applied the trapezoidal rule in the last step. The localmean free paths at the particle positions are calculated followingequations 1 and 2 of Perego et al. (2016), for the total and the energymean free paths respectively.Unlike in Gizzi et al. (2019) and in most of the literature (see Peregoet al. (2014a); George et al. (2020); Endrizzi et al. (2020) as exam-ples), where the standard approach is to pre-select radial paths on agrid, calculate the integral of Eq. 1 over these paths, and then takethe minimum optical depth, here there is no assumption on the typeof path, but is the matter itself that tells neutrinos how to leak out(see also Perego et al. (2014b) for a similar methodology). Fig. A1shows sketches of both the particle (top) and grid-based (bottom)methods. The main advantage of our new approach is the affordablecomputational cost required to get spectral optical depths comparedto grid-based integration methods. This is due on one side to the factthat we use directly the particle properties for the calculations, andon the other side because we adopt the fast recursive bisection treeof Gafton & Rosswog (2011) for neighbour search.Figs. A2-A3 show on the plane x-z the total and energy optical depthsrespectively, for case 1) described in Sec. 3.1. Both of them are calcu-lated on a grid with the algorithm of Gizzi et al. (2019), and with ourparticle-based approach. We SPH-map on the plane x-z the particle-based optical depths for comparison purposes. We show the case fora low-energy bin ( ∼
10 MeV, upper panels) and for a high-energybin ( ∼
100 MeV, lower panels) of our energy grid, and for electronneutrinos (top), electron anti-neutrinos (middle), and heavy-leptonneutrinos (bottom). We also show the map of the relative difference 𝜖 𝜏 = ( 𝜏 grid − 𝜏 part )/ 𝜏 grid . The optical depth appears to be larger alongthe disk with respect to the grid-based approach, and lower along thepoles. Low density regions in SPH are not as well resolved as the highdensity ones, therefore the SPH maps in the former regime are onlyindicative of what the optical depth could actually be (although weexpect it to be rather small anyway). The higher optical depths alongthe disk with the particle-based approach leads to more extendedneutrino surfaces at high neutrino energies, especially for electronneutrinos, which are the most interacting species given the neutronrichness of the material. Nevertheless, our algorithm well captures Figure A1.
Simplified sketches of the particle-based (top) and grid-based(bottom) approaches for the optical depth calculation. Shown in grey is aneutron star remnant surrounded by a disk. In the particle-based model,each particle has a smoothing length defining the radius of each circle. Thelatter encloses the particle neighbours. Unlike the grid-based approach, whereoutgoing radial paths are a priori selected for calculating the optical depth,in the particle-based case the path along which neutrinos move is defined bythe density distribution, and therefore it can be non-straight. the expected distribution of the neutrino surfaces according to dif-ferent neutrino energies and species (see Gizzi et al. (2019); Peregoet al. (2014a) for details). Accounting for different resolutions thatcan be used both in the particle- and grid-based approaches, we esti-mate the particle-based algorithm to be (cid:38)
100 times faster than thegrid-based one.
This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000
This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000 , 1–18 (2020)
SL calibration for BNS merger simulations
200 0 200 x [km] z [ k m ] grid
200 0 200 x [km] part
200 0 200 x [km] x [km] z [ k m ]
200 0 200 x [km]
200 0 200 x [km] x [km] z [ k m ]
200 0 200 x [km]
200 0 200 x [km] x [km] z [ k m ]
200 0 200 x [km]
200 0 200 x [km] x [km] z [ k m ]
200 0 200 x [km]
200 0 200 x [km] x [km] z [ k m ]
200 0 200 x [km]
200 0 200 x [km] Figure A2.
Total optical depth on the x-z plane for top: electron neutrino, middle: electron anti-neutrinos and bottom: heavy-lepton neutrinos. For each specieswe 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 themaps obtained from the grid-based and particle-based calculations, and a map of the relative difference 𝜖 𝜏 = 𝜏 grid − 𝜏 part 𝜏 grid , respectively.MNRAS , 1–18 (2020) Gizzi D. et al.
200 0 200 x [km] z [ k m ] grid
200 0 200 x [km] part
200 0 200 x [km] x [km] z [ k m ]
200 0 200 x [km]
200 0 200 x [km] x [km] z [ k m ]
200 0 200 x [km]
200 0 200 x [km] x [km] z [ k m ]
200 0 200 x [km]
200 0 200 x [km] x [km] z [ k m ]
200 0 200 x [km]
200 0 200 x [km] x [km] z [ k m ]
200 0 200 x [km]
200 0 200 x [km] Figure A3.
Same as Fig. A2, but for the energy optical depth.MNRAS000