Neutrino absorption and other physics dependencies in neutrino-cooled black-hole accretion disks
Oliver Just, Stephane Goriely, Hans-Thomas Janka, Shigehiro Nagataki, Andreas Bauswein
MMNRAS , 000–000 (0000) Preprint 18 February 2021 Compiled using MNRAS L A TEX style file v3.0
Neutrino absorption and other physics dependencies inneutrino-cooled black-hole accretion disks
O. Just , ★ , S. Goriely , H.-Th. Janka , S. Nagataki , , & A. Bauswein , GSI Helmholtzzentrum für Schwerionenforschung, Planckstraße 1, 64291 Darmstadt, Germany Astrophysical Big Bang Laboratory, RIKEN Cluster for Pioneering Research, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan Institut d’Astronomie et d’Astrophysique, CP-226, Université Libre de Bruxelles, 1050 Brussels, Belgium Max-Planck-Institut für Astrophysik, Postfach 1317, 85741 Garching, Germany RIKEN Interdisciplinary Theoretical & Mathematical Science Program (iTHEMS), 2-1 Hirosawa, Wako, Saitama, Japan 351-0198 Helmholtz Research Academy Hesse for FAIR (HFHF), GSI Helmholtz Center for Heavy Ion Research,Campus Darmstadt, Planckstraße 1, 64291 Darmstadt, Germany
18 February 2021
ABSTRACT
Black-hole (BH) accretion disks formed in compact-object mergers or collapsars may be majorsites of the rapid-neutron-capture (r-)process, but the conditions determining the electron frac-tion ( 𝑌 𝑒 ) remain uncertain given the complexity of neutrino transfer and angular-momentumtransport. After discussing relevant weak-interaction regimes, we study the role of neutrinoabsorption for shaping 𝑌 𝑒 using an extensive set of simulations performed with two-momentneutrino transport and again without neutrino absorption. We vary the torus mass, BH massand spin, and examine the impact of rest-mass and weak-magnetism corrections in the neu-trino rates. We also test the dependence on the angular-momentum transport treatment bycomparing axisymmetric models using the standard 𝛼 -viscosity with viscous models assum-ing constant viscous length scales ( 𝑙 t ) and three-dimensional magnetohydrodynamic (MHD)simulations. Finally, we discuss the nucleosynthesis yields and basic kilonova properties. Wefind that absorption pushes 𝑌 𝑒 towards ∼ 𝑌 eq 𝑒 by ∼ 𝑌 𝑒 = .
25, leading to a reduced lanthanide fraction and a brighter, earlier, and bluer kilonovathan without absorption. More compact tori with higher neutrino optical depth, 𝜏 , tend to havelower 𝑌 eq 𝑒 up to 𝜏 ∼ 𝑙 t =const. viscosity (MHDtreatment). The solar-like abundance pattern found for our MHD model marginally supportscollapsar disks as major r-process sites, although a strong r-process may be limited to phasesof high mass-infall rates, (cid:164) 𝑀 (cid:38) × − 𝑀 (cid:12) s − . The recent discovery of a binary neutron-star (NS) mergervia gravitational waves and electromagnetic counterparts,GW170817/AT2017gfo/GRB170817 (e.g. Abbott et al. 2017b,a;Chornock et al. 2017; Villar et al. 2017; Kasen et al. 2017; Metzger2019; Tanvir et al. 2017; Waxman et al. 2018; Perego et al. 2017;Gottlieb et al. 2018; Kawaguchi et al. 2018; Mooley et al. 2018),provided long-sought observational support for the idea that NSmergers are prolific sites of the rapid neutron capture (r-) process(e.g. Lattimer et al. 1977; Eichler et al. 1989; Freiburghaus et al.1999; Goriely et al. 2005, 2011; Korobkin et al. 2012; Wanajo et al.2014; Perego et al. 2014; Just et al. 2015a), can be observed as akilo- or macronova in optical frequency bands (e.g. Metzger et al.2010; Roberts et al. 2011; Tanaka & Hotokezaka 2013; Grossmanet al. 2014; Kasen et al. 2015), can produce short gamma-ray burst(GRB) jets (e.g. Eichler et al. 1989; Ruffert & Janka 1998; Rosswoget al. 2003; Nakar 2007; Lee & Ramirez-Ruiz 2007; Rezzolla et al.2011; Paschalidis et al. 2015; Just et al. 2016), and may serve asunique laboratories for exploring the high-density regime of matter (e.g. Bauswein et al. 2017; Margalit & Metzger 2017; Radice et al.2018; Abbott et al. 2018; Rezzolla et al. 2018). One of the mainejecta components encountered in NS mergers (and also in mergersof NSs with black holes, BH) is thought to originate during thefirst few seconds post merger after the central object has formed aBH and the surrounding disk disintegrates as a result of turbulentangular momentum transport and neutrino cooling (Popham et al.1999; Kohri & Mineshige 2002; Beloborodov 2003; Setiawan et al.2004; Shibata et al. 2007; Metzger et al. 2008; Fernández & Met-zger 2013; Just et al. 2015a; Siegel & Metzger 2018; Hossein Nouriet al. 2018; Siegel et al. 2019; Janiuk 2019; Miller et al. 2019a;Fujibayashi et al. 2020a).Apart from NS mergers, neutrino-cooled BH-accretion disksmay also form during the collapse of a strongly rotating massivestar (e.g. MacFadyen & Woosley 1999). Once the innermost corehas collapsed to a BH the surrounding layers of the star with suf-ficient angular momentum settle on circular orbits. These systems,so-called collapsars, are not only potential candidates for poweringlong GRBs (Woosley 1993), but have also been considered as nu-cleosynthesis sites (e.g. Pruet et al. 2004; Surman & McLaughlin © a r X i v : . [ a s t r o - ph . H E ] F e b Just et al. 𝜈 𝑒 ) to electron-antineutrinos( ¯ 𝜈 𝑒 ) characterize the lepton number transport and regulate the elec-tron fraction, 𝑌 𝑒 . In general, a consistent description of neutrinosrequires solving the radiative transfer (i.e. Boltzmann) equation witha total of six (three spatial plus three momentum) degrees of free-dom, which calls for enormous computational efforts (e.g. Mihalas& Mihalas 1984) even without accounting for the possibility of neu-trino flavor oscillations (e.g. Malkus et al. 2012; Wu et al. 2017;Deaton et al. 2018; Richers et al. 2019). While the first simulationsevolving the Boltzmann equations have recently become available(Miller et al. 2019b), their computational demands yet pose strictlimits on the evolution times and prohibit extensive parameter ex-plorations.The typically rather low masses of neutrino-cooled disks andcorrespondingly low optical depths to neutrinos, compared to neu-tron stars formed in core-collapse supernovae (e.g. Janka 2017) oras remnants of NS mergers (e.g. Perego et al. 2014), spur the notionthat the impact of neutrino absorption might be minor or possiblyeven negligible. From the modeling point of view this situationwould be advantageous, because the smaller the impact of neutrinoabsorption is, the more accurate and credible are results obtainedwith schemes incorporating neutrino absorption only approximatelyor not at all, such as purely local trapping schemes (Shibata et al.2007), leakage schemes (e.g. Fernández et al. 2020), leakage pluspost-processing schemes (e.g. Siegel et al. 2019), or M1 schemes(i.e. two-moment transport schemes with a local closure, e.g. Justet al. 2015b). Indirect evidence for a relatively minor relevance ofneutrino absorption might come from the fact that simulations in-cluding neutrino absorption, if only approximately, only found verysmall amounts of neutrino-driven compared to viscously drivenejecta (Just et al. 2015a; Fujibayashi et al. 2020a). Moreover, therecent simulations by Fujibayashi et al. (2020a) of viscous disks,performed in GR using a combination of energy-independent M1and leakage schemes, seem to suggest that neutrino absorption iseven irrelevant for disks more massive than 0 . 𝑀 (cid:12) . These resultsare, however, in stark contrast to the findings of Miller et al. (2019a),whose GRMHD simulations with Boltzmann neutrino transport ad-vocate a substantial sensitivity of 𝑌 𝑒 to absorption-related effectseven for a 0 . 𝑀 (cid:12) disk. Thus, the role of neutrino absorption andits sensitivity to other modeling ingredients still remains unclearand detailed investigations are overdue.The difficulties connected to the neutrino treatment are aggra-vated by the existence of another, comparably challenging modelingingredient, namely angular momentum transport, i.e. the mechanismthat is mainly responsible for accretion, heating, and ejection of diskmaterial. Being a consequence of MHD turbulence driven by the magneto-rotational instability (MRI, e.g. Balbus & Hawley 1991),angular momentum transport in MHD disks requires, in order forit to be modeled properly, that the simulation is performed in threedimensions and with sufficiently high spatial resolution in order toresolve the wavelengths of MRI growth and the relevant scales ofMHD turbulence. Hence, 3D MHD models, of which the first haverecently become available (Siegel & Metzger 2018; Hossein Nouriet al. 2018; Fernández et al. 2019; Christie et al. 2019; Miller et al.2019b) are computationally quite expensive even without neutrinotransport, and many of their properties still remain unexplored orpoorly understood, particularly concerning the 𝑌 𝑒 evolution and itssensitivity to details of the neutrino interactions.Both of the aforementioned requirements for MHD models canbe relaxed by reverting to an approximate mean-field description ofturbulent angular momentum transport, such as embodied by the 𝛼 -viscosity approach (Shakura & Sunyaev 1973) that has been em-ployed in numerous 1D and 2D studies (e.g. Popham et al. 1999; DiMatteo et al. 2002; Chen & Beloborodov 2007; Metzger et al. 2009;Fernández & Metzger 2013; Just et al. 2015a; Fujibayashi et al.2020a). Given their computational efficiency, viscous disk modelshave been studied already for a much broader range of conditionsand input parameters than MHD models. However, although takeninto account by various published results with different degreesof sophistication (e.g. Just et al. 2015a; Fujibayashi et al. 2020a;Miller et al. 2019b), we still lack a comprehensive understandingof the importance of neutrino absorption. Moreover, relatively littleattention has been drawn so far to the sensitivity of the ejecta prop-erties with respect to other components of modeling, such as usinga non-standard prescription for the dynamic viscosity (Fujibayashiet al. 2020a), neglecting rest-mass terms in the 𝛽 -reaction rates orincluding weak magnetism corrections (Horowitz 2002), or usingdifferent initial 𝑌 𝑒 values in the torus.Whether occurring in the course of a NS merger or of a col-lapsar, the possibility that neutrino-cooled disks may be major sitesof r-process elements calls for a profound understanding of all pro-cesses and modeling assumptions that have a leverage on the ejecta 𝑌 𝑒 , first and foremost the interplay between neutrino emission, neu-trino absorption, and angular momentum transport. In this study wetherefore systematically investigate, on the basis of two- and three-dimensional viscous and MHD simulations including M1 neutrinotransport, the impact of neutrino absorption and the sensitivity tothe treatment of angular momentum transport and to the variationof global model parameters. We further test uncertainties connectedto details of the neutrino interaction rates and to the initial electronfraction of the torus. In order to relate the obtained dependencies ofthe hydrodynamical simulations to nucleosynthesis variations andto the kilonova signal, we compute for all models the abundances ofr-process elements and basic properties of the bolometric kilonovalight curve.This paper is organized as follows: In Sect. 2 we first reviewthe equilibrium conditions for weak interactions and corresponding 𝑌 𝑒 values and characteristic timescales and test their sensitivity tocommonly used approximations. Section 3 describes the setup ofour numerical models and of the post-processing steps aiming atevaluating the nucleosynthesis yields and kilonova light curve. InSect. 4 we first summarize basic features of the torus evolution As pointed out by the anti-dynamo theorem (Moffatt 1978), axisymmetricmodels suffer from the inability to efficiently create poloidal magnetic fieldsfrom toroidal fields, a mechanism that is needed to keep the MRI alive andtherefore to sustain angular momentum transport.MNRAS , 000–000 (0000) eutrino absorption in black-hole disks and the neutrino emission, followed by an analysis of the impactof neutrino absorption on the torus evolution and on the outflow.Moreover, we discuss the nucleosynthesis yields and the kilonovaproperties. In Sect. 5 we discuss implications of our results based ona comparison with existing studies. Finally, in Sect. 6 we summarizeand conclude our study. 𝑌 𝐸 Before discussing numerical models we first review the neutrino in-teraction rates, equilibrium conditions, and characteristic timescalesthat are relevant for the evolution of the electron fraction, 𝑌 𝑒 , inneutrino-cooled accretion disks. The interactions mainly responsible for changing 𝑌 𝑒 in neutrino-cooled disks are the nucleonic 𝛽 -processes, namely electron captureon protons, positron capture on neutrons, electron neutrino captureon neutrons, and electron anti-neutrino capture on protons. Theinteraction rates corresponding to these processes are given by (Bruenn 1985; Horowitz 2002): 𝜆 𝑒 − = 𝐾 𝛽 ∫ ∞ 𝜖 𝐹 𝑒 − ( 𝜖 + ) 𝜖 + √︄ − (cid:18) 𝑚 𝑒 𝑐 𝜖 + (cid:19) d 𝜖 , (1a) 𝜆 𝑒 + = 𝐾 𝛽 ∫ ∞ 𝜖 𝜖 𝐹 𝑒 + ( 𝜖 − ) 𝜖 − √︄ − (cid:18) 𝑚 𝑒 𝑐 𝜖 − (cid:19) d 𝜖 , (1b) 𝜆 𝜈 𝑒 = 𝐾 𝛽 ∫ ∞ 𝜖 𝐹 𝜈 𝑒 ( 𝜖 )( − 𝐹 𝑒 − ( 𝜖 + )) 𝜖 + √︄ − (cid:18) 𝑚 𝑒 𝑐 𝜖 + (cid:19) d 𝜖 , (1c) 𝜆 ¯ 𝜈 𝑒 = 𝐾 𝛽 ∫ ∞ 𝜖 𝜖 𝐹 ¯ 𝜈 𝑒 ( 𝜖 )( − 𝐹 𝑒 + ( 𝜖 − )) 𝜖 − √︄ − (cid:18) 𝑚 𝑒 𝑐 𝜖 − (cid:19) d 𝜖 , (1d)where 𝐹 𝑥 ( 𝜖 ) is the distribution function of particle 𝑥 at energy 𝜖 in-tegrated over solid angles in momentum space, 𝑐 the speed of light, 𝑚 𝑒 the electron mass, 𝜖 = 𝑄 𝑛𝑝 + 𝑚 𝑒 𝑐 with 𝑄 𝑛𝑝 being the neutron-proton mass difference, 𝜖 ± = 𝜖 ± 𝑄 𝑛𝑝 , and 𝐾 − 𝛽 = ( )( 𝑚 𝑒 𝑐 ) .The composition of the gas and its thermodynamic properties en-ter the rates through the distribution functions 𝐹 𝑒 ± . The effects ofweak magnetism and nucleon recoil can additionally be taken intoaccount by multiplying the integrands in Eq. (1) by correction fac-tors 𝑅 wm 𝜈 𝑒 / ¯ 𝜈 𝑒 ( 𝜖 ) (see Horowitz 2002 for explicit expressions). Withthe above rates the evolution equation of 𝑌 𝑒 for a Lagrangian fluidelement reads:d 𝑌 𝑒 d 𝑡 = ( 𝜆 𝑒 + + 𝜆 𝜈 𝑒 ) 𝑌 𝑛 − ( 𝜆 𝑒 − + 𝜆 ¯ 𝜈 𝑒 ) 𝑌 𝑝 , (2)where 𝑌 𝑛 / 𝑝 = 𝑛 𝑛 / 𝑝 / 𝑛 𝐵 is the number of free neutrons/protonsrelative to the total number of baryons. In a gas consisting onlyof free neutrons and protons one has 𝑌 𝑝 = 𝑌 𝑒 and 𝑌 𝑛 = ( − 𝑌 𝑒 ) , whereas in a gas composed of nuclei in nuclear statistical We neglect phase space blocking for neutrinos and nucleons as wellas other rate corrections that only play a role at larger densities ( 𝜌 (cid:29) g cm − ) than typically encountered in neutrino-cooled accretion disks. equilibrium (NSE) 𝑌 𝑛 / 𝑝 are functions of density, 𝜌 , temperature, 𝑇 ,and 𝑌 𝑒 . Setting d 𝑌 𝑒 / d 𝑡 = ( 𝜆 𝑒 + + 𝜆 𝜈 𝑒 ) 𝑌 𝑛 − ( 𝜆 𝑒 − + 𝜆 ¯ 𝜈 𝑒 ) 𝑌 𝑝 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 𝜌,𝑇 ,𝑌 eq 𝑒 = , (3)defines an equilibrium value, 𝑌 eq 𝑒 , that would asymptotically bereached by a fluid element with a given density and temperatureand exposed to a given neutrino field. The characteristic timescaleon which 𝑌 𝑒 approaches 𝑌 eq 𝑒 can be estimated as 𝜏 𝛽 = 𝑌 𝑝 ( 𝜆 𝑒 − + 𝜆 ¯ 𝜈 𝑒 ) + 𝑌 𝑛 ( 𝜆 𝑒 + + 𝜆 𝜈 𝑒 ) . (4)Anywhere along a fluid trajectory, weak interactions drive 𝑌 𝑒 to thelocal 𝑌 eq 𝑒 on a local timescale 𝜏 𝛽 . Once 𝜏 𝛽 becomes longer thanthe expansion timescale 𝜏 exp ∼ 𝜌 / (cid:164) 𝜌 ∼ 𝑟 / 𝑣 𝑟 (with 𝑟 and 𝑣 𝑟 beingthe radius and radial velocity of the fluid element) in an expandingoutflow, 𝑌 𝑒 effectively remains constant, i.e. it freezes out. 𝑌 eq 𝑒 In what follows we will briefly discuss three limiting cases of 𝑌 eq 𝑒 and comment on the relevance of each for neutrino-cooled BH-tori. Given the sub-nuclear densities and relatively low neutrino opticaldepths in neutrino-cooled disks, it is reasonable to assume that thebulk 𝑌 𝑒 is to a large extent determined by neutrino emission, i.e.the rates 𝜆 𝑒 ± . In situations when neutrino absorption even becomesnegligible, 𝑌 eq 𝑒 converges to 𝑌 eq , em 𝑒 , which is defined by 𝜆 𝑒 + 𝑌 𝑛 − 𝜆 𝑒 − 𝑌 𝑝 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 𝜌,𝑇 ,𝑌 eq , em 𝑒 = 𝜌 and 𝑇 ). Contours of 𝑌 eq , em 𝑒 ( 𝜌, 𝑇 ) are shown in panel (a) ofFig. 1 overlaid with contours of the electron degeneracy parameter, 𝜂 𝑒 , and the characteristic timescales of neutrino emission, 𝜏 em = 𝑌 𝑝 𝜆 𝑒 − + 𝑌 𝑛 𝜆 𝑒 + ∼ 𝐾 𝛽 ( 𝑘 𝐵 𝑇 ) − (F ( 𝜂 𝑒 ) + F (− 𝜂 𝑒 )) − ∼ (cid:18) 𝑇 (cid:19) − s . (6)As pointed out by Liu (2010) a very good approximation to 𝑌 eq , em 𝑒 ,at least whenever nuclei are absent, can be recovered directly fromthe equation-of-state (EOS) table by exploiting the condition 𝜇 𝜈 ≡ 𝜇 𝑝 − 𝜇 𝑛 + 𝜇 𝑒 = − 𝜇 𝑒 (7)for the chemical potentials 𝜇 𝑖 of species 𝑖 . Contours of 𝑌 𝑒 resultingfrom Eq. (7) are plotted as purple lines in panel (b) of Fig. 1).Several previous studies have discussed the emission equilib-rium defined by Eq. (5) (e.g. Beloborodov 2003; Metzger et al. 2008; The condition 𝜇 𝜈 = − 𝜇 𝑒 follows from the consideration that equal ratesof 𝑝 + 𝑒 − → 𝑛 + 𝜈 𝑒 and 𝑛 + 𝑒 + → 𝑝 + ¯ 𝜈 𝑒 define a kinetic equilibrium,for which 𝜇 𝑝 (cid:164) 𝑛 𝑝 + 𝜇 𝑒 (cid:164) 𝑛 𝑒 − = 𝜇 𝑛 (cid:164) 𝑛 𝑛 − 𝜇 𝑒 (cid:164) 𝑛 𝑒 + and where all (cid:164) 𝑛 𝑖 are equal(see Liu 2010 for more details). We furthermore stress that the quantity 𝜇 𝜈 should only be interpreted as neutrino chemical potential if neutrinos arethermalized, otherwise it is just a placeholder for 𝜇 𝑝 − 𝜇 𝑛 + 𝜇 𝑒 .MNRAS , 000–000 (0000) Just et al. temperature [MeV]10 d e n s i t y [ g / c m ] m01M3A8 em e (a) s s . s . s s s s . . . properties at emission equilibrium Y e q , e m e temperature [MeV]10 d e n s i t y [ g / c m ] = 0 = e em (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . s s . s . s s s s different equilibria Y e q , e m e temperature [MeV]10 d e n s i t y [ g / c m ] with WM corr. Q np = m e =0 em (c) . . . . . . . . . . . . s s . s . s s s s different approximations Y e q , e m e temperature [MeV]10 d e n s i t y [ g / c m ] with WM corr. Q np = m e =0 em (d) different approximations, 20 energy bins . . . . . . . . . . . s s . s . s s s s Y e q , e m e temperature [MeV]10 d e n s i t y [ g / c m ] SFHOonly n,p em (e) . . . . . . . . . . . . . . . . . s s . s . s s s s different NSE ensembles Y e q , e m e temperature [MeV]10 d e n s i t y [ g / c m ] < ( e )>MeV =4 6 8
10 12 14162024 28 < ( e )>MeV =4 6 8
10 12
16 20 24 28 (f) . . . . . . . . .
90 0 . . . . . . . . . mean energy of emitted neutrinos Y e q , e m e Figure 1.
Properties connected to the kinetic emission equilibrium, which is established once the rate of 𝑝 + 𝑒 − → 𝑛 + 𝜈 𝑒 equals that of 𝑛 + 𝑒 + → 𝑝 + ¯ 𝜈 𝑒 .The color map in all panels illustrates 𝑌 eq , em 𝑒 defined by Eqs. (1a), (1b), and (5) and the 4-species NSE composition employed in our numerical simulations. Panel (a): characteristic neutrino emission timescale, 𝜏 em (white lines), electron degeneracy parameter, 𝜂 𝑒 (purple lines), and average density-temperatureevolution of a fiducial numerical model (dashed black line); (b): 𝑌 𝑒 corresponding to 𝜇 𝜈 = − 𝜇 𝑒 (purple lines) and to 𝜇 𝜈 = (c): 𝑌 eq , em 𝑒 computedwith weak-magnetism and recoil corrections (red lines) as well as using the simplification 𝑄 𝑛𝑝 = 𝑚 𝑒 = (d): same as panel (c) but the colormap and lines are obtained using the coarser neutrino energy grid that is employed in numerical simulations of this study; (e): 𝑌 eq , em 𝑒 resulting with the NSEcomposition of the SFHO EOS (red lines) and for a pure neutron-proton gas (purple lines); (f): the mean energies of neutrinos, (cid:104) 𝜖 (cid:105) , emitted from a gas withthe density, temperature, and 𝑌 𝑒 = 𝑌 eq , em 𝑒 given at each point (where some regions less relevant to the freeze out are neglected). All 𝑌 𝑒 contours show valuesof 0.1, 0.2, etc. from top to bottom. Since the dynamical timescales of outflows in neutrino-cooled disks are typically no longer than ∼ 𝜏 em =
100 s contour, where approximately 𝑇 (cid:46) Arcones et al. 2010; Fujibayashi et al. 2020a). The basic notion forthe interpretation of 𝑌 𝑒 in neutrino-cooled disks is that during theirexpansion and cooling fluid elements sitting in and being releasedfrom the torus travel from top right in the 𝜌 − 𝑇 domain (where 𝑌 eq , em 𝑒 (cid:28) . 𝑌 eq , em 𝑒 (cid:38) . 𝑌 𝑒 to 𝑌 eq 𝑒 ( 𝜌, 𝑇 ) until 𝑌 𝑒 freezes out somewhere near the 𝜏 em ∼ . . . . 𝜏 em starts exceeding the expansion timescale. As first pointedout by Chen & Beloborodov (2007), the torus remains mildly degen- erate during its expansion, i.e. with electron-degeneracies 𝜂 𝑒 ∼ 𝑌 𝑒 in sim-ulations of BH-tori can be explained by 𝑌 eq , em 𝑒 , corrections due tothe presence of neutrinos have not been examined so far. MNRAS , 000–000 (0000) eutrino absorption in black-hole disks In the opposite limiting case when neutrino absorption dominatesneutrino emission, such as in neutrino-driven winds, 𝑌 eq 𝑒 will begiven by 𝑌 eq , abs 𝑒 , which fulfills 𝜆 𝜈 𝑒 𝑌 𝑛 − 𝜆 ¯ 𝜈 𝑒 𝑌 𝑝 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) 𝜌,𝑇 ,𝑌 eq , abs 𝑒 = , (8)and the corresponding absorption timescale is 𝜏 abs = 𝑌 𝑛 𝜆 𝜈 𝑒 + 𝑌 𝑝 𝜆 ¯ 𝜈 𝑒 . (9)The electron fraction in absorption equilibrium, 𝑌 eq , abs 𝑒 , dependsmainly (though not solely) on the number densities and mean ener-gies of both neutrino species. Assuming a pure nucleon gas, 𝑌 eq , abs 𝑒 is given by 𝑌 eq , abs 𝑒 = 𝜆 𝜈 𝑒 𝜆 𝜈 𝑒 + 𝜆 ¯ 𝜈 𝑒 , (10)which can further be approximated by 𝑌 eq , abs 𝑒 ∼ (cid:32) + (cid:104) 𝜖 𝜈 𝑒 (cid:105) 𝑛 ¯ 𝜈 𝑒 (cid:104) 𝜖 𝜈 𝑒 (cid:105) 𝑛 𝜈 𝑒 (cid:33) − , ∼ (cid:32) + (cid:104) 𝜖 𝜈 𝑒 (cid:105) 𝐿 𝑁 , ¯ 𝜈 𝑒 (cid:104) 𝜖 𝜈 𝑒 (cid:105) 𝐿 𝑁 ,𝜈 𝑒 (cid:33) − , (11)where 𝑛 𝜈 and 𝐿 𝑁 ,𝜈 are number densities and number fluxes (ornumber luminosities), respectively, for neutrino species 𝜈 and theenergy averages are given by (cid:104) 𝜖 𝜈 (cid:105) = ( ∫ 𝜖 𝐹 𝜈 d 𝜖 )/( ∫ 𝜖 𝐹 𝜈 d 𝜖 ) . Ap-proximate expressions similar to those given in Eq. (11) have beenemployed for the purpose of investigating neutrino-driven windsin numerous studies (e.g. Qian & Woosley 1996; Horowitz & Li1999). The estimate in Eq. (11) neglects mass corrections (i.e. 𝑄 𝑛𝑝 = 𝑚 𝑒 =
0) and ignores Pauli blocking for 𝑒 ± , while the secondline additionally assumes that 𝐿 𝑁 ,𝜈 𝑒 / 𝐿 𝑁 , ¯ 𝜈 𝑒 ≈ 𝑛 𝜈 𝑒 / 𝑛 ¯ 𝜈 𝑒 . In this pa-per we always use Eq. (8) for the computation of 𝑌 eq , abs 𝑒 , because allof the aforementioned assumptions are not entirely justified in thebulk of the torus. During the neutrino-dominated phase the emittedneutrino energies are relatively low (cf. Sect. 2.3 and panel (f) ofFig. 1), electrons are mildly degenerate (cf. panel (a) of Fig. 1), and 𝑛 𝜈 𝑒 > 𝑛 ¯ 𝜈 𝑒 may hold while at the same time 𝐿 𝑁 ,𝜈 𝑒 < 𝐿 𝑁 , ¯ 𝜈 𝑒 (e.g.Wu et al. 2017).In regions surrounding neutrino sources that are approximatelyin emission equilibrium (meaning that d 𝑌 𝑒 / d 𝑡 ≈
0, i.e. the emissiontimescales are short compared to other timescales) roughly the samenumber of 𝜈 𝑒 and ¯ 𝜈 𝑒 neutrinos are emitted per unit of time, suchthat 𝑌 eq , abs 𝑒 is typically close to 0.5. Finally, in the limiting case that the neutrino mean free paths be-come shorter than the hydrodynamic length scales, neutrinos be-come trapped and thermalized by the fluid and attain a Fermi-Diracdistribution that is defined solely by 𝜇 𝜈 and 𝑇 . In the extremecase that no neutrinos diffuse out from local fluid patches (typi-cally for densities above 𝜌 ∼ g cm − ), the total lepton fraction 𝑌 𝑙 = 𝑌 𝑒 + 𝑌 𝜈 𝑒 − 𝑌 ¯ 𝜈 𝑒 (where 𝑌 𝑖 = 𝑛 𝑖 / 𝑛 𝐵 ) remains conserved (i.e.d 𝑌 𝑙 / d 𝑡 = 𝑌 𝑒 becomes an instantaneousfunction of 𝜌, 𝑇 , and 𝑌 𝑙 (see, e.g., Sekiguchi et al. 2012; Perego et al.2016; Ardevol-Pulpillo et al. 2019, for schemes making use of the concept of trapped neutrinos). In neutrino-cooled disks with max-imum densities of only 𝜌 ∼ ... g cm − neutrinos will, if atall, barely reach local thermodynamic equilibrium. However, theneutrino distribution may still be close to thermal, but with vanish-ing chemical potential, 𝜇 𝜈 →
0, because the leakage of neutrinosdrives the number densities down to 𝑛 𝜈 𝑒 − 𝑛 ¯ 𝜈 𝑒 (cid:28) 𝑛 𝑒 − − 𝑛 𝑒 + (see,e.g., Ruffert et al. 1997; Beloborodov 2003). The condition 𝜇 𝜈 (cid:12)(cid:12) 𝜌,𝑇 ,𝑌 eq ,𝜇𝜈 = 𝑒 = equilibrium value, namely 𝑌 eq ,𝜇 𝜈 = 𝑒 . Contours of 𝑌 eq ,𝜇 𝜈 = 𝑒 ( 𝜌, 𝑇 ) are shown in panel (b) of Fig. 1 (red lines), re-vealing that 𝑌 eq ,𝜇 𝜈 = 𝑒 exceeds 𝑌 eq , em 𝑒 by about ∼ . 𝑌 eq ,𝜇 𝜈 = 𝑒 can thus be used as a quantity to roughly estimate(or bracket, together with 𝑌 eq , em 𝑒 representing the opposite limit-ing case) the impact of neutrino absorption. This interpretation issupported by our simulations (cf. Sect. 4.3.2), where we find that 𝑌 eq , em 𝑒 < 𝑌 eq 𝑒 ∼ 𝑌 𝑒 < 𝑌 eq ,𝜇 𝜈 = 𝑒 in the torus at early times duringwhich neutrino absorption is efficient. 𝑌 eq , em 𝑒 to commonly employedapproximations Since 𝑌 eq , em 𝑒 is a proxy for the 𝑌 𝑒 attained in the torus, one can assessthe impact of certain modeling approximation on nucleosynthesisconditions in the ejecta without performing any simulations simplyby checking their influence on 𝑌 eq , em 𝑒 .The first assumption to test is that of neglecting correctionsto 𝜆 𝑒 ± of Eq. (1) associated with the finite electron mass and theneutron-proton mass difference by setting 𝑄 𝑛𝑝 = 𝑚 𝑒 =
0. Testingthis sensitivity is motivated by the fact that many existing disk andmerger models based on grey neutrino leakage schemes employedthis approximation (as originally suggested by Ruffert et al. 1996)in order to reduce the complexity of the integrals and thereforethe computational demands. As the purple lines in panel (c) ofFig. 1 show, this simplification reduces 𝑌 eq , em 𝑒 quite considerablycompared to its original value, namely by about 0 . − . 𝜌 ∼ ... g cm − and 𝑇 ∼ . . . 𝜆 𝑒 ± . This correction has been studied sofar only in the context of neutrino-driven winds, where it was foundto be responsible for increasing 𝑌 𝑒 in the absorption-dominatedregime by about ∼ −
20 % (Horowitz & Li 1999; Pllumbi et al.2015; Goriely et al. 2015). In our case we are instead interested inthe emission equilibrium and find that weak magnetism correctionsshift 𝑌 eq , em 𝑒 towards lower values (because it reduces the absorptioncross section of positrons), but only by a very small amount in theregions relevant to neutrino-cooled disks (cf. red lines in panel (c)of Fig. 1). Opposite to the previously discussed corrections associ-ated with 𝑄 𝑛𝑝 and 𝑚 𝑒 , the impact of the weak-magnetism correc-tion grows with temperature and therefore with the mean energy of The apparent tension with the literature of cold neutron stars (e.g. Yakovlevet al. 2001), where no distinction is being made between 𝑌 eq , em 𝑒 and 𝑌 eq ,𝜇 𝜈 = 𝑒 , can be resolved by realizing that both quantities become iden-tical in the zero-temperature limit.MNRAS , 000–000 (0000) Just et al. emitted neutrinos, which is plotted for 𝜈 𝑒 and ¯ 𝜈 𝑒 assuming emissionequilibrium in panel (f) of Fig. 1.Since we are about to perform energy-dependent neutrinotransport simulations using a limited number of energy zones, weshould ensure that the transport energy grid is able to reproduce thevalues of 𝑌 eq , em 𝑒 as accurately as possible. While in all other panelsquantities related to 𝑌 eq , em 𝑒 have been obtained using a very fineenergy grid for the integrals appearing in Eq. (1), in panel (d) ofFig. 1 we confirm that our transport grid consisting of 20 energybins (see Sect. 3.1.3 for more details of the grid) is able to reproducevery well the quantities plotted in panel (c).The last question we address is whether the equilibrium 𝑌 𝑒 inBH-tori is sensitive to the ensemble of nuclear species taken intoaccount in the EOS. In panel (e) of Fig. 1, we therefore compare 𝑌 eq , em 𝑒 resulting for a pure nucleon gas (purple lines), for the SFHOEOS (red lines, Steiner et al. 2013) that contains a large number ofnuclear species, and a 4-species EOS ( 𝑛, 𝑝, 𝛼 plus one heavy nu-cleus, color map in Fig. 1). The 4-species EOS, used in all remainingpanels of Fig. 1, will be employed also for the numerical modelsin the remainder of this study. Figure 1 reveals that noticeable dif-ferences between the three cases of NSE ensembles appear onlyfor emission timescales 𝜏 em (cid:38) −
10 s, and in particular betweenSFHO and our 4-species EOS only for 𝜏 em >
10 s. Since the freezeout of 𝑌 𝑒 is expected to occur already when 𝜏 em (cid:46) 𝑌 eq , em 𝑒 . After discussing basic aspects of the 𝑌 𝑒 evolution, we now describethe setup of our numerical study of neutrino-cooled BH-tori. We simulate BH-accretion disks with the code AENUS-ALCAR(Obergaulinger 2008; Just et al. 2015b) that handles neutrino trans-port in the M1 approximation and solves the (viscous or magneto-)hydrodynamic and transport equations on a spherical polar grid us-ing Riemann-solver based finite-volume methods. The numericalmethods for all viscous models (cf. Sect. 3.1.2) are exactly the sameas those in Just et al. (2015a) unless stated explicitly otherwise. Wesummarize the evolution equations in Appendix A. For neutrinointeractions, we include the 𝛽 -processes (with rates given by 𝜆 𝑖 ofEqs. (1)) as well as iso-energetic scattering with nucleons and nuclei(Bruenn 1985). Justified by the relatively low temperatures and den-sities encountered in neutrino-cooled disks and by the dominance ofthe included processes under such conditions, we neglect all othertypes of neutrino interactions, such as channels that produce 𝜇 - and 𝜏 -neutrinos. We do not take into account general relativistic effects,except that we use the pseudo-Newtonian gravitational potential byArtemova et al. (1996) to approximately incorporate effects associ-ated with the innermost stable circular orbit (ISCO) and to mimicthe shrinking of the ISCO with increasing BH spin.One of the main goals of this study is to characterize the im-pact of neutrino absorption and quantify its importance relative toneutrino emission in determining 𝑌 𝑒 in the torus and the ejectedmaterial. To this end we conduct for most models two simulations,one with full M1 transport (i.e. including emission and absorption) and one neglecting neutrino absorption as well as any other neutrinoreactions except neutrino emission. In doing so we aim at bracket-ing the maximum error that could be encountered when neglectingabsorption or when treating absorption with more simplistic neu-trino schemes. Although the M1 scheme used in this study is notas accurate as a Boltzmann scheme that solves the six-dimensionalneutrino transfer equations, it is an actual transport scheme – in thesense that it solves time-dependent conservation equations for theneutrino energy density and flux density – and is consistent withthe transfer equations in the diffusion and free streaming limits.The M1 scheme has shown excellent agreement with a ray-by-ray Boltzmann solver (VERTEX, Rampp & Janka 2002; Buras et al.2006) for describing core-collapse supernovae (Just et al. 2018). Itsaccuracy in the context of merger remnants has only been tested fortime independent configurations (Just et al. 2015a; Foucart et al.2018; Sumiyoshi et al. 2021) or short-term simulations (Foucartet al. 2020), where no serious shortcomings were found.Table 1 summarizes the set of investigated models. Our fiducialreference model, m01M3A8, is chosen to have a BH mass of 𝑀 BH = 𝑀 (cid:12) , BH spin parameter of 𝐴 BH = .
8, and a relatively low initialtorus mass, 𝑚 = . 𝑀 (cid:12) . By using such a low torus mass ourresults will bracket the impact of absorption from below, i.e. all diskswith larger masses but otherwise same parameters can be expectedto show an even stronger dependence on absorption. Moreover,another reason for exploring the low-mass regime is that low-masstori may represent, at least approximately, BH-disks formed in thecore of collapsars (Siegel et al. 2019; Miller et al. 2019a).In order to test the sensitivity of the results with respect toglobal parameters and modeling assumptions, we vary the initialdisk mass (indicated by number following “m” in the model name),BH mass and spin (“M” and “A” in model name, respectively),and the treatment of turbulent angular momentum transport (asdescribed in Sect. 3.1.2). Motivated by the discussion in Sect. 2.3,we also test (model m01M3A8-noQM) the impact of neglectingmass corrections, i.e. setting 𝑄 𝑛𝑝 = 𝑚 𝑒 = 𝛽 -processes,Eqs. (1), and of including weak magnetism and recoil corrections (model m01M3A8-wm).We set up the initial torus as equilibrium configuration withconstant specific angular momentum (similarly as in Just et al.2015a; Fernández & Metzger 2013; Fujibayashi et al. 2020a; Siegelet al. 2019). We fix the initial specific entropy at 8 𝑘 𝐵 per baryon,put the location of the initial density maximum at a radius of 40 km ×( 𝑀 BH /( 𝑀 (cid:12) )) , and set the initial electron fraction, 𝑌 𝑒 , to a constantvalue everywhere, namely to 𝑌 𝑒 = . 𝑌 𝑒 = . The ray-by-ray approximation assumes that the radiation field is axisym-metric around each radial coordinate line in a spherical polar coordinatesystem. This approximation was found to be well justified in 3D CCSNsimulations (Glas et al. 2019), whereas in axisymmetry (Skinner et al. 2016;Just et al. 2018) it may artificially promote shock runaway. The reason why our fiducial model does not already include weak mag-netism corrections is not related to any physics argument. After includingweak-magnetism corrections at a later stage of this project, we deemed itunnecessary to repeat a large number of simulations that had already beenperformed without weak magnetism corrections, in particular consideringthat the quantitative impact of this correction is rather small.MNRAS , 000–000 (0000) eutrino absorption in black-hole disks Table 1.
Properties of investigated numerical models. The columns provide from left to right: model name, initial torus mass, BH mass, BH spin, viscous 𝛼 -parameter, the treatment of turbulent angular momentum transport (where std. 𝛼 -vis. or 𝑙 t =const. vis. means that Eq. (13) or Eq. (14) is used, respectively),inclusion of 𝑄 𝑛𝑝 and 𝑚 𝑒 in 𝛽 -interaction rates, inclusion of weak magnetism corrections, initial electron fraction, final simulation time, percentage of initialtorus mass left on the grid at final simulation time, and dimensionality of the simulation. The suffix “-no 𝜈 ” in parenthesis denotes that a separate model wasevolved without neutrino absorption but otherwise similar properties. Entries not in (in) parenthesis refer to models evolved with (without) neutrino absorption.model 𝑚 𝑀 BH 𝐴 BH 𝑌 𝑒 ( 𝑡 =0) 𝛼 vis viscosity mass corr. weak magn. 𝑡 fin 𝑚 fintor / 𝑚 dimensionsname [ 𝑀 (cid:12) ] [ 𝑀 (cid:12) ] treatment included? included? [s] [%]m01M3A8(-no 𝜈 ) 0.01 3 0.8 0.5 0.06 std. 𝛼 -vis. yes no 10 (10) < <
1) 2Dm1M3A8(-no 𝜈 ) 0.1 10 (10) < < 𝜈 ) 0.001 10 (10) < < 𝜈 ) 0.01 5 10 (10) < < 𝜈 ) 10 20 (20) < < 𝜈 ) 3 0.4 10 (10) < < 𝜈 ) 0.9 10 (10) < < 𝛼 𝜈 ) 0.8 0.03 10 (10) < < 𝛼 𝜈 ) 0.1 10 (10) < < 𝜈 ) 0.05 𝑙 t =const. vis. 20 (20) 5.7 (6.81)m01M3A8-mhd(-no 𝜈 ) – MHD 2.1 (2.1) 12.5 (15.0) 3Dm01M3A8-noQm(-no 𝜈 ) 0.06 std. 𝛼 -vis. no 10 (10) < <
1) 2Dm01M3A8-wm yes yes 10 < 𝜈 ) 0.1 no 10 (10) < < and how many r-process elements can be generated self-consistentlyonly by the disk. Thus, the final r-process abundances resulting frommodels with 𝑌 𝑒 = . Most of our models are run in axisymmetry and employ the well-known 𝛼 -viscosity scheme (Shakura & Sunyaev 1973) to mimicturbulent angular momentum transport. We employ for the dynamicviscosity coefficient the prescription 𝜂 vis = 𝛼 vis 𝜌𝑐 𝑖 𝑐 𝑖 Ω 𝐾 (13)(where the isothermal sound speed 𝑐 𝑖 = √︁ 𝑃 / 𝜌 with the gas pressure 𝑃 and the Keplerian angular velocity Ω 𝐾 ), which we call the “stan-dard 𝛼 -viscosity” prescription here, because it was already used ina large number of previous disk studies (e.g. Narayan & Yi 1994;Igumenshchev & Abramowicz 2000; Setiawan et al. 2004; Fernán-dez & Metzger 2013; Just et al. 2015a). We only take into accountthe 𝑟𝜙 - and 𝜃𝜙 -components of the viscous stress tensor, because weonly intend to parametrize the turbulent stresses related to differen-tial rotation. In most cases we use 𝛼 vis = .
06, while we also addedmodels with lower and higher values of 𝛼 vis (cf. models with suffix 𝛼
03 and 𝛼
1, respectively, in Table 1).Since the 𝛼 -viscosity scheme is nothing but a mean-fieldparametrization of unresolved small-scale physics based on heuris-tic dimensional arguments, one is of course left with certain freedomfor defining 𝜂 vis . Fujibayashi et al. (2020a) made use of this libertyand argued that the characteristic length scale 𝑙 t ∼ 𝑐 𝑖 / Ω 𝐾 implic-itly assumed in the prescription of Eq. (13) may be overestimatedat large radii. They instead suggested a different formulation of the 𝛼 -viscosity for which 𝑙 t =const. and 𝜂 vis = 𝛼 vis 𝜌𝑐 𝑠 𝑙 t (14)(with the adiabatic sound speed 𝑐 𝑠 = √︁ 𝛾𝑃 / 𝜌 ). In order to test thesensitivity with respect to the viscosity prescription we implementedthis treatment in model m01M3A8-vis2 and its pendant neglecting neutrino absorption, using 𝑙 t = 𝛼 vis = .
05 (i.e.the same values as used in Fujibayashi et al. 2020a).The third, and most consistent option of treating angular mo-mentum transport is to evolve MHD models, starting with a pre-defined magnetic field configuration; cf. models m01M3A8-mhd(-no 𝜈 ). Moreover, we also switch from a Newtonian axisymmetricdescription to a special relativistic three-dimensional (3D) evolu-tion (see Appendix A for the evolved equations), which is necessaryto overcome time-step limitations (due to unlimited Alfven speedsin Newtonian MHD) and the anti-dynamo theorem (Moffatt 1978),respectively. The initial torus in the MHD models has the same prop-erties as for the viscous models, but is now endowed with a singlepoloidal magnetic-field loop, of which the magnetic field strength ischosen such that the ratio of mass-weighted average of gas pressureto mass-weighted average of magnetic pressure is 𝛽 init ≈ 𝛽 init in our case) as de-scribed in Fernández et al. (2019). We note that the initial magneticfield configuration can have a non-negligible leverage on the torusevolution and ejecta properties (Beckwith et al. 2008; Christie et al.2019), while additionally also the demand on the grid resolution ishigher in MHD models than in viscous disks. However, given thehuge computational expenses to perform the special relativistic 3DMHD simulations, we were unable to address these aspects in thecurrent study. For the same reason, we could only evolve these mod-els until ∼ (cid:62)
10 s for the viscousmodels). Nevertheless, to our knowledge the MHD simulations pre-sented here are the longest performed so far using genuine neutrinotransport (i.e. solving time-dependent conservation equations forneutrinos).
The numerical setup is as follows: For the viscous (MHD) modelsthe radial domain is discretized between 𝑟 min and 2 × km (1 × km) using a logarithmic grid with 𝑁 𝑟 =
400 (320) zones, where 𝑟 min is approximately given by the arithmetic average between theISCO radius and the radius of the event horizon (e.g. 𝑟 min ≈
10 km
MNRAS000
MNRAS000 , 000–000 (0000)
Just et al. for 𝑀 BH = 𝑀 (cid:12) and 𝐴 BH = . 𝜃 = 𝜋 . The 3D MHD models make use of 128 polar zones, of which126 zones are distributed non-uniformly between 𝜋 / < 𝜃 < 𝜋 / Δ 𝜃 ≈ . ◦ at the equator to ≈ ◦ close to thepole, and the two remaining zones attach both ends of the grid tothe polar axis. The azimuthal coordinate, 𝜙 , in the MHD modelsruns from 0 to 𝜋 / 𝜋 / − depending on the model) and beyond that radius decreases roughlyas 𝑟 − . For the MHD models we additionally enforce the conditionthat the magnetic energy should never be larger than 50 times therest-mass energy in a given cell, where both energies are measuredin the comoving frame. In order to analyze the outflow and obtain the nucleosynthesis yieldswe extract fluid trajectories from the hydrodynamical simulationsand post-process the trajectories in a separate step. In all modelsperformed in axisymmetry and employing a viscosity scheme, weuse the available output data from the simulation to numericallyintegrate particle trajectories backward in time starting at a radiusof 𝑟 = cm and at pre-defined times and polar angles for a totalnumber of about 600-1100 trajectories per model. For the 3D MHDmodels we instead distribute 5 × tracer particles of equal massin the initial disk and advect their locations during the simulationwhile recording their thermodynamic properties.The r-process nucleosynthesis is calculated by post-processingthe thus obtained trajectories of the ejected matter. The density andtemperature evolution is taken directly from the hydrodynamicaltrajectories. The initial composition is determined by nuclear sta-tistical equilibrium when the density has dropped below the dripdensity and matter has cooled below 𝑇 = K during the expan-sion. As soon as the temperature has fallen below 10 GK, furtherchanges of the composition are followed by a full network calcula-tion including all 5000 nuclear species from protons up to 𝑍 = 𝛽 -stability and the neutron-drip line.All charged-particle fusion reactions as well as their reverse reac-tions on light elements up to Th isotopes are included in addition toradiative neutron captures and photodisintegrations on all speciesup to 𝑍 =
110 isotopes. The reaction rates on light species aretaken from the NETGEN library, which includes all the latest com-pilations of experimentally determined reaction rates (Xu et al.2013). Experimentally unknown reactions are estimated with theTALYS code (Koning et al. 2005; Goriely et al. 2008) on the ba-sis of the Skyrme Hartree-Fock-Bogolyubov (HFB) nuclear massmodel, HFB-21 (Goriely et al. 2010). On top of these reactions, 𝛽 -decays as well as 𝛽 -delayed neutron emission probabilities arealso included, the corresponding rates being taken from the rela-tivistic mean-field plus ramdom-phase-approximation calculation(Marketin et al. 2016).All fission rates, i.e. the neutron-induced, photo-induced, 𝛽 -delayed and spontaneous fission rates, are estimated on the basis ofthe HFB-14 fission paths (Goriely et al. 2007; Goriely 2015) and thenuclear level densities within the combinatorial approach (Gorielyet al. 2008) are obtained with the same single-particle scheme and pairing strength. The neutron-induced fission rates are calculatedwith the TALYS code for all nuclei with 90 (cid:54) 𝑍 (cid:54)
110 (Gorielyet al. 2009). Similarly, the 𝛽 -delayed and spontaneous fission ratesare estimated with the same TALYS fission barrier penetration cal-culation. The 𝛽 -delayed fission rate takes into account the full com-petition between the fission, neutron and photon channels, weightedby the population probability given by the 𝛽 -decay strength function(Kodama & Takahashi 1975). The fission fragment yield distribu-tion is estimated with the renewed statistical scission-point modelbased on microscopic ingredients, the so-called SPY model, as de-scribed in Lemaître et al. (2019). The possibility to observe the torus ejecta and r-process yieldsdirectly by means of the emitted kilonova opens up prospects ofquantitatively determining the global parameters of the merging bi-nary, e.g. the post-merger torus mass, its composition etc. However,many dependencies still remain unexplored that establish the linkbetween the input physics of hydrodynamic models and the resultingkilonova produced by the ejected material. In order to explore somebasic sensitivities we employ a simple and approximate kilonovacomputation scheme to translate the variations of 𝑌 𝑒 and of nuclearmass fractions found for our hydrodynamic models into variationsof observable features, such as the time, photospheric temperature,and bolometric luminosity at peak emission.Depending on the intended level of consistency, various moreand less sophisticated approaches have previously been used forcomputing the kilonova light curve based on the ejecta proper-ties (e.g. Li & Paczyński 1998; Kulkarni 2005; Metzger et al.2010; Grossman et al. 2014; Kasen et al. 2015; Perego et al.2017; Kawaguchi et al. 2018; Metzger 2019; Hotokezaka & Nakar2020; Korobkin et al. 2020). In our study we employ an approx-imate and to-our-knowledge new method that makes direct use ofthe mass and velocity of each outflow trajectory as well as thecorresponding nucleosynthesis data (lanthanide plus actinide frac-tion, 𝑋 LA ≡ 𝑋 lan + 𝑋 act , and radioactive heating rate per mass, 𝑞 rad ). While not being as accurate as sophisticated radiative transferschemes, our method has the advantage over several existing ap-proximate schemes that it captures the detailed dependence of massdensity, heating rate, and 𝑋 LA on the velocity coordinate, 𝑣 , insteadof assuming idealized functions. In the following we outline themethod only briefly, while explicit equations are provided in Ap-pendix B. We assume spherical symmetry, homology (i.e. 𝑟 ∝ 𝑣 ),and employ the grey approximation for the light curve, such thatthe only coordinate degrees of freedom are the expansion time andejecta velocity. The velocity space is discretized by about 𝑁 𝑣 ∼ km. A mass-weighted inte-gration of the corresponding quantities then yields 𝑚 ( 𝑣 ) , 𝑋 LA ( 𝑣, 𝑡 ) , and 𝑞 ( 𝑣, 𝑡 ) , where 𝑚 ( 𝑣 ) equals the ejecta mass with velocity greaterthan 𝑣 and 𝑞 ≡ 𝜖 therm 𝑞 rad is the effective heating rate including thethermalization factor 𝜖 therm that is obtained from interpolation ofvalues provided for case “Random” in Table 1 of Barnes et al.(2016). We finally need opacities, 𝜅 , in order to compute a bolo-metric light curve from these data. Lacking opacity data based ondetailed atomic models (e.g. Kasen et al. 2013; Tanaka et al. 2020),we approximate 𝜅 as a function of 𝑋 LA and 𝑇 , which is calibratedby fitting kilonova light curves from Kasen et al. (2017); see Ap-pendix B for more details. Although this simplified prescription is acrude approximation to atomic opacities, it captures the main effectthat we are interested in, namely the correlation between 𝜅 and the MNRAS , 000–000 (0000) eutrino absorption in black-hole disks lanthanide content of the ejecta. After preparation of the data in thedescribed manner, we solve a two-moment system of conservationequations for the energy density and flux density of photons usingthe M1 closure. We first review, mostly on the basis of the fiducial modelsm01M3A8(-no 𝜈 ), some generic features of the neutrino emissionand the weak interaction freeze out. Although many of the followingaspects are not entirely new (e.g. Metzger et al. 2008; Fernández &Metzger 2013; Just et al. 2015a), we briefly summarize them hereand introduce quantities, which will subsequently serve as diagnos-tic tools to assess the impact of neutrino absorption.Figures 2 and 3 show several global quantities as functionsof time, and Figs. 4, 5, and 6 provide contours of density, 𝑌 𝑒 ,and temperature, as well as of neutrino emission timescales andabsorption timescales for different times of evolution. Here and inthe following we consider as the torus all material that is locatedbetween the inner radial boundary of the computational domainand the radius 𝑟 ≡ km. The torus evolution can roughly bedivided into two phases (Metzger et al. 2008): In the first, neutrino-dominated phase neutrino cooling is efficient in removing heat fromviscous processes, i.e. the neutrino emission timescale, estimatedas (cid:104) 𝜏 em (cid:105) = ∫ 𝑟<𝑟 𝑛 𝐵 d 𝑉 ∫ 𝑟<𝑟 (cid:164) 𝑛 em d 𝑉 (15)(where (cid:164) 𝑛 em = 𝜆 𝑒 − 𝑛 𝑝 + 𝜆 𝑒 + 𝑛 𝑛 with 𝜆 𝑖 from Eq. (1) and proton-/neutron number densities 𝑛 𝑝 / 𝑛 ) is shorter than the viscoustimescale, computed as (cid:104) 𝜏 vis (cid:105) = 𝑚 tor (cid:164) 𝑚 tor (16)(where 𝑚 tor is the total baryonic mass within the sphere of radius 𝑟 ); see panel (d) in Fig. 2. Moreover, during the neutrino-dominatedphase fluid elements accreting onto the BH typically radiate away asizable fraction of their rest-mass energy in the form of neutrinos,i.e. 𝜂 𝜈 = 𝐿 𝜈 /( (cid:164) 𝑀 BH 𝑐 ) (cid:38) − 𝜂 𝜈 plunges and viscosity becomesa powerful agent of mass ejection (e.g. Fernández & Metzger 2013;Just et al. 2015a).The evolution of 𝑌 𝑒 as well as of the equilibrium values in-troduced in Sect. 2, mass averaged up to a radius of 10 km, isdepicted in the top rows of Figs. 7 (for models with different torusmasses) and 8 (for models with different viscosity treatments). Dur-ing the neutrino-dominated phase a large fraction of the torus isclose to weak equilibrium, meaning that neutrino emission is ef-ficient enough for 𝑌 𝑒 to roughly track its respective equilibriumvalue 𝑌 eq 𝑒 . Therefore, in all models that start with an initial valueof 𝑌 𝑒 = .
5, the torus average, (cid:104) 𝑌 𝑒 (cid:105) (black lines), first drops fromits initial value until approximately reaching (cid:104) 𝑌 eq 𝑒 (cid:105) (green or redlines for models with or without neutrino absorption, respectively),and then rises again in an attempt to match (cid:104) 𝑌 eq 𝑒 (cid:105) , which in turn in-creases owing to viscous expansion and a subsiding level of electron Recall that 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 for models without neutrino absorption. degeneracy in the torus. However, the decreasing rates of neutrinoproduction (i.e. growing emission timescales 𝜏 em ) in the more andmore diluting torus thwart and ultimately terminate the evolutionof 𝑌 𝑒 in the torus.Closely connected to the aforementioned features is the behav-ior of the ratio of 𝜈 𝑒 to ¯ 𝜈 𝑒 number luminosities, 𝐿 𝑁 ,𝜈 𝑒 / 𝐿 𝑁 , ¯ 𝜈 𝑒 ,which is plotted in panel (f) of Fig. 2. After a short initialphase of 𝐿 𝑁 ,𝜈 𝑒 / 𝐿 𝑁 , ¯ 𝜈 𝑒 (cid:29)
1, during which the torus with ini-tially 𝑌 𝑒 = . 𝐿 𝑁 ,𝜈 𝑒 / 𝐿 𝑁 , ¯ 𝜈 𝑒 remains fairly close to, but marginally be-low, unity during the neutrino-dominated phase. The condition 𝐿 𝑁 ,𝜈 𝑒 ∼ 𝐿 𝑁 , ¯ 𝜈 𝑒 is basically equivalent to d (cid:104) 𝑌 𝑒 (cid:105)/ d 𝑡 ∼ 𝑌 𝑒 ∼ 𝑌 eq 𝑒 and that neutrino-emissiontimescales are short or comparable to the dynamical (i.e. viscous)timescales. The fact that 𝐿 𝑁 ,𝜈 𝑒 < 𝐿 𝑁 , ¯ 𝜈 𝑒 derives from the tendencythat 𝑌 𝑒 keeps running behind the increasing values of 𝑌 eq 𝑒 . Once thetorus becomes non-radiative, the number luminosities drop, 𝑌 𝑒 de-couples from 𝑌 eq 𝑒 , and the result is a growing disparity between 𝑌 𝑒 and 𝑌 eq 𝑒 , which is counteracted by the system with boosting theproduction rates of ¯ 𝜈 𝑒 relative to those of 𝜈 𝑒 . This explains the dropof 𝐿 𝑁 ,𝜈 𝑒 / 𝐿 𝑁 , ¯ 𝜈 𝑒 at late times.The efficiency by which neutrino emission can change theaverage electron fraction of the torus can be measured by the fraction 𝑚 em / 𝑚 tor of torus material that produces neutrinos on timescales 𝜏 em shorter than a certain threshold value; here we use 𝜏 em < 𝑡 ∼
100 ms). We thereforeuse this quantity to measure the freeze-out time, 𝑡 em , as the timewhen 𝑚 em / 𝑚 tor drops below 10 %. We caution, however, that 𝑡 em is only a crude estimate and that fluid elements freeze out overan extended range of time before and after 𝑡 em (see, e.g., plotsof 𝑚 FO in the bottom panels of Figs. 7 and 8). For our set ofmodels, 𝑡 em lies between 50 and 300 ms (cf. Table 2), while it isprolonged for disks with greater masses and for more massive andfaster rotating BHs. In general, torus configurations with overalllower temperatures freeze out earlier. In the extreme case of lowdisk compactness (i.e. low values of disk mass over disk size) andtherefore low temperature, the bulk of the torus may not even beable to achieve weak equilibrium (in the sense that (cid:104) 𝑌 𝑒 (cid:105) ≈ (cid:104) 𝑌 eq 𝑒 (cid:105) )to begin with. This is the case for model m01M10A8 with a 10 𝑀 (cid:12) central BH and for model m001M3A8 having a 0.001 𝑀 (cid:12) torus(cf. left panel of Fig. 7, where (cid:104) 𝑌 𝑒 (cid:105) is barely changing from itsinitial value of 0.5). Table 2 also reveals that the freeze-out time isroughly inversely proportional to 𝛼 vis , which is not surprising as 𝛼 vis essentially regulates the accretion timescale.Table 2 also lists for all models the mass accretion rates into theBH measured at 𝑡 em . Most values, albeit with large scatter, lie within (cid:164) 𝑀 BH ( 𝑡 em ) ∼ − . . . − 𝑀 (cid:12) s − , which is broadly consistentwith analytical estimates (see, e.g., Siegel et al. 2019; De & Siegel2020, who call this value the ignition accretion rate). Finally, inorder to facilitate comparison of the neutrino properties observed inour models with studies using more or less sophisticated neutrinoschemes, we provide in Table 2 the total percentage of accretedrest-mass energy radiated away by neutrinos and the mean energiesof all emitted neutrinos. The late-time variations of (cid:104) 𝑌 𝑒 (cid:105) visible in Figs. 7 and 8 are not caused byweak interactions but by material leaving the control volume.MNRAS000
100 ms). We thereforeuse this quantity to measure the freeze-out time, 𝑡 em , as the timewhen 𝑚 em / 𝑚 tor drops below 10 %. We caution, however, that 𝑡 em is only a crude estimate and that fluid elements freeze out overan extended range of time before and after 𝑡 em (see, e.g., plotsof 𝑚 FO in the bottom panels of Figs. 7 and 8). For our set ofmodels, 𝑡 em lies between 50 and 300 ms (cf. Table 2), while it isprolonged for disks with greater masses and for more massive andfaster rotating BHs. In general, torus configurations with overalllower temperatures freeze out earlier. In the extreme case of lowdisk compactness (i.e. low values of disk mass over disk size) andtherefore low temperature, the bulk of the torus may not even beable to achieve weak equilibrium (in the sense that (cid:104) 𝑌 𝑒 (cid:105) ≈ (cid:104) 𝑌 eq 𝑒 (cid:105) )to begin with. This is the case for model m01M10A8 with a 10 𝑀 (cid:12) central BH and for model m001M3A8 having a 0.001 𝑀 (cid:12) torus(cf. left panel of Fig. 7, where (cid:104) 𝑌 𝑒 (cid:105) is barely changing from itsinitial value of 0.5). Table 2 also reveals that the freeze-out time isroughly inversely proportional to 𝛼 vis , which is not surprising as 𝛼 vis essentially regulates the accretion timescale.Table 2 also lists for all models the mass accretion rates into theBH measured at 𝑡 em . Most values, albeit with large scatter, lie within (cid:164) 𝑀 BH ( 𝑡 em ) ∼ − . . . − 𝑀 (cid:12) s − , which is broadly consistentwith analytical estimates (see, e.g., Siegel et al. 2019; De & Siegel2020, who call this value the ignition accretion rate). Finally, inorder to facilitate comparison of the neutrino properties observed inour models with studies using more or less sophisticated neutrinoschemes, we provide in Table 2 the total percentage of accretedrest-mass energy radiated away by neutrinos and the mean energiesof all emitted neutrinos. The late-time variations of (cid:104) 𝑌 𝑒 (cid:105) visible in Figs. 7 and 8 are not caused byweak interactions but by material leaving the control volume.MNRAS000 , 000–000 (0000) Just et al. m [ M / s ] (a) accretion and outflow rates accretedejected 10 m [ m t o r ] (b) torus and ejected mass torus ( r <10 km)ejected ( r >10 km) 10 ( L e + L e ) / ( M B H c ) (c) neutr. emission efficiency [ s ] (d) characteristic timescales < em >< abs >< vis > 10 m / m t o r ( t ) (e) mass affected by em./abs. m ( em < 1s) m ( abs < 1s) 0.51.01.52.02.53.0 L N , e / L N , e (f) ratio of e / e number loss rates < r > [ k m ] (g) avg. torus radius <> [ g / c m ] (h) avg. torus density < T > [ M e V ] (i) avg. torus temperature <> (j) mean energy of e releasedabsorbed 10 o p t (k) optical depth along equator @ = released @ = absorbed a b s (l) absorption factor ee time [s]0.10.20.30.40.50.6 Y e (m) ejecta Y e at r = 10 km time [s]10 v / c (n) ejecta vel. at r = 10 km time [s]10 s [ k B / nu c l e o n ] (o) ejecta entropy at r = 10 km m = 0.001 Mm = 0.01 Mm = 0.1 M Figure 2.
Time evolution of global properties for models m001M3A8 (blue lines), m01M3A8 (black lines), and m1M3A8 (orange lines). Thin lines ofsame color denote the corresponding “no 𝜈 ” models ignoring neutrino absorption. Panel (a): mass-accretion rates onto the BH and ejecta mass-flux rates at 𝑟 = km; (b): mass of torus (i.e. at 𝑟 < km) and ejected material (i.e. at 𝑟 > km) ; (c): total neutrino luminosity (measured at 𝑟 =
500 km inthe laboratory frame) divided by rest-mass energy accretion rate onto the BH; (d): characteristic timescales of emission, Eq. (15), absorption, Eq. (21), andviscous expansion, Eq. (16); (e): fraction of the torus mass for which 𝜏 em < 𝜏 abs < (f): ratio of 𝜈 𝑒 over ¯ 𝜈 𝑒 numberfluxes measured at 𝑟 =
500 km in the laboratory frame; (g), (h), (i): mass-weighted averages of radial coordinate, density, and temperature, respectively; (j): mean energy of released 𝜈 𝑒 as measured at 𝑟 =
500 km (solid lines) and of absorbed 𝜈 𝑒 (dashed lines); (k): the optical depth along the equator for the energiesshown in panel (j); (l): absorption factor, Eq. (20); (m), (n), (o): mass-weighted angular averages of outflow material crossing the sphere at 𝑟 = km of 𝑌 𝑒 ,velocity, and entropy per baryon, respectively. MNRAS , 000–000 (0000) eutrino absorption in black-hole disks Table 2.
Properties characterizing weak interactions in the disk: Minimum of average 𝑌 𝑒 , maximum optical depth at a neutrino energy of 20 MeV along equator,estimated times when emission ( 𝑡 em ) and absorption ( 𝑡 abs ) become irrelevant, mass accretion rates onto the BH at times 𝑡 em and 𝑡 abs , percentage of accretedrest-mass energy radiated in the form of neutrinos until 𝑡 fin , mean energy of radiated electron neutrinos ( 𝜈 𝑒 ) and electron antineutrinos ( ¯ 𝜈 𝑒 ) computed as (cid:104) 𝜖 (cid:105) 𝜈 rel = 𝐿 𝜈 / 𝐿 𝜈,𝑁 where 𝐿 𝜈 and 𝐿 𝜈,𝑁 are measured at 𝑟 =
500 km in the laboratory frame, mean energy of all absorbed neutrinos (i.e. neutrinos capturedby protons or neutrons) defined by Eq. (19). Entries not in (in) parenthesis refer to models evolved with (without) neutrino absorption.model (cid:104) 𝑌 𝑒 (cid:105) mintor 𝜏 maxopt (
20 MeV ) 𝑡 em / abs (cid:164) 𝑀 BH ( 𝑡 em / abs ) ¯ 𝜂 𝜈 (cid:104) 𝜖 (cid:105) 𝜈 𝑒 / ¯ 𝜈 𝑒 rel (cid:104) 𝜖 (cid:105) 𝜈 𝑒 / ¯ 𝜈 𝑒 abs name w/ (w/o) abs. w/ (w/o) abs. [ms] [10 − 𝑀 (cid:12) 𝑠 − ]) [%] [MeV] [MeV]m01M3A8 0.227 (0.153) 14.0 (11.1) 105/ 60 1.0/3.2 4.2 14.2/17.5 23.6/28.1m1M3A8 0.247 (0.073) 47.1 (106 ) 175/135 5.2/8.6 5.2 13.7/17.1 35.0/38.1m001M3A8 0.346 (0.336) 1.44 (1.0 ) 55/ 17 0.3/1.6 1.2 12.5/15.3 18.9/23.5m01M5A8 0.261 (0.228) 6.08 (7.10) 120/ 60 1.2/4.2 2.8 12.6/15.3 19.3/23.7m01M10A8 0.370 (0.366) 1.47 (1.50) 145/ 32 1.4/7.3 0.8 9.50/11.7 14.0/17.9m01M3A4 0.240 (0.177) 12.0 (17.1) 85/ 35 1.0/4.9 1.7 13.0/16.3 21.6/26.4m01M3A9 0.224 (0.148) 15.4 (23.4) 110/ 72 1.1/2.7 6.3 14.8/18.1 24.9/29.4m01M3A8- 𝛼
03 0.176 (0.118) 17.2 (11.9) 235/130 0.5/1.6 4.6 12.9/15.6 20.7/25.0m01M3A8- 𝛼 While the observations described in the previous section remain ap-plicable also for models m01M3A8-vis2(-no 𝜈 ), several quantitativedifferences appear when using Eq. (14) instead of Eq. (13) to ex-press the dynamic viscosity coefficient 𝜂 vis . As shown in Fig. 9, thequantity 𝑐 𝑖 / Ω 𝐾 , which is employed in the conventional 𝛼 -viscosityprescription (cf. Eq. (13)) as a proxy for the length scale of tur-bulent eddies, 𝑙 t , grows roughly linearly with radius. In contrast, 𝑙 t = const. = 𝜈 ). This mis-match between length scales implies that the impact of viscosity inthe 𝑙 t =const. scheme is comparable to the conventional 𝛼 -viscosityscheme only at small radii, whereas it becomes relatively weaker atlarger radii. Since the torus is expanding with time (cf. the radialcoordinate of the center-of-mass, (cid:104) 𝑟 (cid:105) 𝜌 , shown in panel (g) of Fig. 3),the aforementioned circumstance explains why most evolutionaryfeatures agree well at early times between models m01M3A8 (blacklines in Fig. 3) and m01M3A8-vis2 (green lines), but start to divergeat later times. Approximately once the torus in model m01M3A8-vis2 reaches (cid:104) 𝑟 (cid:105) 𝜌 (cid:38)
100 km, the torus expansion decelerates, anddensities and temperatures drop more slowly (see (cid:104) 𝜌 (cid:105) 𝜌 and (cid:104) 𝑇 (cid:105) 𝜌 in panels (h) and (i) of Fig. 3, respectively). This causes neutrinoemission to remain efficient for a longer time (cf. panels (c) and (e)in Fig. 3) and to impact a greater fraction of the expanding ejectarelative to tori with a standard 𝛼 -viscosity. This is further discussedin Sect. 4.3.3. The basic evolution of the MHD models, m01M3A8-mhd(-no 𝜈 ), iswell in agreement with that reported in previous studies of MHDdisks (e.g. De Villiers et al. 2003; McKinney et al. 2014; Siegel& Metzger 2018; Fernández et al. 2019): Right after the start ofthe simulation the poloidal fields are wound up into toroidal fields,which soon become the dominant field component in the mainbody of the disk. Meanwhile, the MRI starts to grow and quicklygenerates turbulence that fuels angular momentum transport. In thesurroundings of the disk a magnetized corona forms, while the polarregions develop a funnel of strong radial 𝐵 -field, which dominates all baryonic energies and persists until the end of the simulation.Importantly, even though the flow pattern is turbulent right fromthe beginning in the MHD models – whereas it remains ratherlaminar in the viscous models during the neutrino-dominated phase– the features of neutrino emission and weak freeze out describedin Sect. 4.1 are similarly observed for the MHD models (see, e.g.,Fig. 3). Specific features of the MHD models in comparison to theviscous models will be discussed in Sect. 4.3.3.Comparing the MHD models with and without neutrino ab-sorption (e.g. panel (h) in Fig. 3 showing the average torus density)we observe that the evolution of the optically thin torus notice-ably lags behind that of the full transport model. The most likelyexplanation for this is that the almost perfect neutrino cooling inthe no-absorption model until about 𝑡 ∼
600 ms reduces the diskthickness so dramatically that the dominant MRI modes becomenumerically under-resolved. This explanation is supported by thesudden transition to a faster disk expansion once neutrino coolinghas become inefficient (see, e.g., strong decline of torus density inpanel (h) at about 𝑡 ∼
600 ms). A higher grid resolution could there-fore solve this issue, however, the necessary grid refinement factormay be large, and we were unable to run additional high-resolutionsimulations given our limited budget of computing time. The sup-pression of the MRI in model m01M3A8-mhd-no 𝜈 to some extentdelays the ejection of material and, by doing so, facilitates ejectawith relatively high 𝑌 𝑒 (see Sect. 4.3.3 for a deeper discussion ofthis aspect). Before investigating the quantitative impact of neutrino absorp-tion, we first have a look at some qualitative features. In Figs. 4–6,contour maps are shown of the neutrino emission timescales, 𝜏 em ,and absorption timescales, 𝜏 abs , for models m01M3A8, m1M3A8,m01M3A8-mhd, and their corresponding “no 𝜈 ” pendants in whichneutrino absorption is neglected. The conditions for weak interac-tions observed in our models motivate the definition of four char-acteristic regions in the torus, which are schematically depicted in MNRAS000
600 ms). A higher grid resolution could there-fore solve this issue, however, the necessary grid refinement factormay be large, and we were unable to run additional high-resolutionsimulations given our limited budget of computing time. The sup-pression of the MRI in model m01M3A8-mhd-no 𝜈 to some extentdelays the ejection of material and, by doing so, facilitates ejectawith relatively high 𝑌 𝑒 (see Sect. 4.3.3 for a deeper discussion ofthis aspect). Before investigating the quantitative impact of neutrino absorp-tion, we first have a look at some qualitative features. In Figs. 4–6,contour maps are shown of the neutrino emission timescales, 𝜏 em ,and absorption timescales, 𝜏 abs , for models m01M3A8, m1M3A8,m01M3A8-mhd, and their corresponding “no 𝜈 ” pendants in whichneutrino absorption is neglected. The conditions for weak interac-tions observed in our models motivate the definition of four char-acteristic regions in the torus, which are schematically depicted in MNRAS000 , 000–000 (0000) Just et al. m [ M / s ] (a) accretion and outflow rates accretedejected 10 m [ m t o r ] (b) torus and ejected mass torus ( r <10 km)ejected ( r >10 km) 10 ( L e + L e ) / ( M B H c ) (c) neutr. emission efficiency [ s ] (d) characteristic timescales < em >< abs >< vis > 10 m / m t o r ( t ) (e) mass affected by em./abs. m ( em < 1s) m ( abs < 1s) 0.51.01.52.02.53.0 L N , e / L N , e (f) ratio of e / e number loss rates < r > [ k m ] (g) avg. torus radius <> [ g / c m ] (h) avg. torus density < T > [ M e V ] (i) avg. torus temperature <> (j) mean energy of e releasedabsorbed 10 o p t (k) optical depth along equator @ = released @ = absorbed a b s (l) absorption factor ee time [s]0.10.20.30.40.50.6 Y e (m) ejecta Y e at r = 10 km time [s]10 v / c (n) ejecta vel. at r = 10 km time [s]10 s [ k B / nu c l e o n ] (o) ejecta entropy at r = 10 km std. -viscosity l t =const. viscosityMHD Figure 3.
Same as Fig. 2 but for the series of models with different treatments of angular momentum transport, m01M3A8 (black lines), m01M3A8-vis2 (bluelines), and m01M3A8-mhd (orange lines).
Fig. 10. The B region (green) in Fig. 10 can be identified in Figs. 4–6by the overlap region, where both neutrino emission (white solidlines) and absorption (white dashed lines) are efficient. This region,if present at a given time, is located around the densest parts of thetorus and typically extends to radii of 𝑟 ∼ −
250 km. In thisregion the electron fraction is close to 𝑌 eq 𝑒 (cp. Eq. (3)). During early evolution times and for high mass accretion rates, the northernand southern surface layers of the torus are surrounded by a region(denoted A in Fig. 10), in which the temperatures are relatively low,such that absorption dominates neutrino emission, i.e. 𝜏 abs < 𝜏 em .Here, 𝑌 𝑒 is pushed towards 𝑌 eq , abs 𝑒 . Moreover, neutrino heating inthis region may even be dynamically important and launch a thermal MNRAS , 000–000 (0000) eutrino absorption in black-hole disks Figure 4.
Snapshots of model m01M3A8 (top row) and the corresponding model without neutrino absorption, m01M3A8-no 𝜈 (bottom row). Each panel showsthe density, 𝜌 , electron fraction, 𝑌 𝑒 , emission timescale, 𝜏 em , and absorption timescale, 𝜏 abs , at evolution times 𝑡 = , 𝑇 = . , , and 3 MeV, and the white lines on the right-hand side mark the surfaces ofregions where the emission timescale (solid) and absorption timescale (dashed) are shorter than 1 s as well as the location where both timescales are equal(dotted). Figure 5.
Same as Fig. 4 but for models m1M3A8 (upper panels) and m1M3A8-no 𝜈 (lower panels) with initial torus mass of 0.1 𝑀 (cid:12) instead of 0.01 𝑀 (cid:12) . wind . Next, since the torus optical depth is largest in the equatorialdirection, the edge of the torus around the midplane experiences rel- Given the low torus masses for most of our models, a neutrino-drivenwind is observed only in model m1M3A8 during a short initial phase, as atively low neutrino absorption rates, while neutrino emission ratesmay still be significant. In this region (C in Fig. 10) 𝑌 𝑒 → 𝑌 eq , em 𝑒 . can be seen when comparing the mass ejection rates in panel (a) of Fig. 2between models m1M3A8 and m1M3A8-no 𝜈 .MNRAS000
Same as Fig. 4 but for models m1M3A8 (upper panels) and m1M3A8-no 𝜈 (lower panels) with initial torus mass of 0.1 𝑀 (cid:12) instead of 0.01 𝑀 (cid:12) . wind . Next, since the torus optical depth is largest in the equatorialdirection, the edge of the torus around the midplane experiences rel- Given the low torus masses for most of our models, a neutrino-drivenwind is observed only in model m1M3A8 during a short initial phase, as atively low neutrino absorption rates, while neutrino emission ratesmay still be significant. In this region (C in Fig. 10) 𝑌 𝑒 → 𝑌 eq , em 𝑒 . can be seen when comparing the mass ejection rates in panel (a) of Fig. 2between models m1M3A8 and m1M3A8-no 𝜈 .MNRAS000 , 000–000 (0000) Just et al.
Figure 6.
Same as Fig. 4 but for models m01M3A8-mhd (upper panels) and m01M3A8-mhd-no 𝜈 (lower panels), where the viscosity is replaced by magneticfields, and showing data that has been obtained by averaging along the azimuthal direction. Y e m = 0.001 M , std. -vis. Y e m = 0.01 M , std. -vis. < Y eq,em e >< Y eq,abs e >< Y eq e >< Y eq, =0 e >< Y e > Y e m = 0.1 M , std. -vis. time [s]0.00.20.40.60.81.0 Y F O e , m F O / m e j time [s]0.00.20.40.60.81.0 Y F O e , m F O / m e j m FO / m ej Y FO e time [s]0.00.20.40.60.81.0 Y F O e , m F O / m e j Figure 7.
Top row:
Mass-weighted averages of the electron fraction, 𝑌 𝑒 (black lines), as well as the equilibrium quantities 𝑌 eq , em 𝑒 (Eq. (5), red lines), 𝑌 eq , abs 𝑒 (Eq. (8), magenta lines), 𝑌 eq 𝑒 (Eq. (3, green lines)), and 𝑌 eq ,𝜇 𝜈 = 𝑒 (Eq. (12), blue lines) for models m001M3A8(-no 𝜈 ) (left panel), m01M3A8(-no 𝜈 ) (middlepanel), and m1M3A8(-no 𝜈 ) (right panel). Thick lines (thin lines) belong to models including (neglecting) neutrino absorption. The averages include all materialwithin radii of 10 km. Note that the reason for (cid:104) 𝑌 𝑒 (cid:105) not remaining exactly constant after weak freeze out is that torus material leaves the considered volumedue to mass ejection at 10 km and mass accretion onto the BH. Bottom row:
Mass fraction of outflow material for which 𝑌 𝑒 is already frozen out, 𝑚 FO / 𝑚 ej ,i.e. with weak interaction timescales longer than 10 s (golden lines), and average 𝑌 𝑒 of all ejecta trajectories that are freezing out during a short time window Δ 𝑡 centered around the current time 𝑡 , 𝑌 FO 𝑒 . MNRAS , 000–000 (0000) eutrino absorption in black-hole disks Y e m = 0.01 M , std. -vis. Y e m = 0.01 M , l t =const. vis. < Y eq,em e >< Y eq,abs e >< Y eq e >< Y eq, =0 e >< Y e > Y e m = 0.01 M , MHD time [s]0.00.20.40.60.81.0 Y F O e , m F O / m e j time [s]0.00.20.40.60.81.0 Y F O e , m F O / m e j m FO / m ej Y FO e time [s]0.00.20.40.60.81.0 Y F O e , m F O / m e j Figure 8.
Same as Fig. 7 but for models m01M3A8-vis2(-no 𝜈 ) and m01M3A8-mhd(-no 𝜈 ). radius [cm]10 l t u r b characteristic length scale of -viscosity c i / K at t = 10ms c i / K at t = 100ms c i / K at t = 1s c i / K at t = 10s l turb = 9km Figure 9.
Comparison of 𝑙 t = 𝑐 𝑖 / Ω 𝐾 (where 𝑐 𝑖 = √︁ 𝑃 / 𝜌 and Ω 𝐾 is theKeplerian angular velocity) used as characteristic length scale in the standard 𝛼 -viscosity approach with 𝑙 t = 𝜈 ).The profiles of 𝑙 t = 𝑐 𝑖 / Ω 𝐾 are computed along the equator for the indicatedtimes based on the data of model m01M3A8. Finally, region D subsumes all torus material where weak interac-tions take place on timescales much longer than the viscous (orexpansion) timescales and 𝑌 𝑒 is basically frozen out.The four different cases sketched in Fig. 10 represent dif-ferent regimes of mass accretion rates, (cid:164) 𝑀 BH . During the evolu-tion of a sufficiently massive torus, region A will disappear first,when (cid:164) 𝑀 BH (cid:46) . 𝑀 (cid:12) s − , followed by region B once the op-tical depth for neutrinos falls below unity (approximately once (cid:164) 𝑀 BH (cid:46) . 𝑀 (cid:12) s − ). Finally, region C vanishes around the freeze-out time, 𝑡 em , when (cid:164) 𝑀 BH drops below ∼ . 𝑀 (cid:12) s − . As a result,a fixed-mass torus (such as formed in a NS merger) may never ex- hibit regions A, B, or C if its initial mass is too low to achieve thecorresponding mass-accretion rates.As can be seen in Fig. 6, the turbulent behavior of the MHDmodels strongly distorts instantaneous iso-contours and sometimeseven disrupts them into multiple patches. However, these perturba-tions are stochastic in nature, hence, a classification based on Fig. 10and described above still holds for the MHD models but should thenbe interpreted in the time-averaged sense, i.e. considering propertiesthat are averaged over a few dynamical timescales. Next, since to our knowledge this quantity has rarely been system-atically reported so far for neutrino-cooled disks, we have a brieflook at the optical depth, 𝜏 opt ( 𝜖 ) = ∫ 𝜅 abs ( 𝜖 ) d 𝑟 , (17)where 𝜅 abs ( 𝜖 ) is the total opacity as function of neutrino energy 𝜖 . Although the optical depth, and therefore the neutrino emissionrates, depend on the direction in which neutrinos leave the torus, weare only interested in the maximum optical depth, which is obtainedwhen performing the integral in Eq. (17) along the equatorial di-rection through the entire computational domain. A straightforwardestimate of 𝜏 opt can be obtained by adopting for the neutrino energythe mean energy of released neutrinos measured far away from thetorus, given by 𝐿 𝜈 / 𝐿 𝜈,𝑁 . We use the approximate formula (e.g.Bruenn 1985) 𝜅 abs ( 𝜖 ) ≈ 𝜎 𝑚 𝑒 𝑐 ( 𝑔 𝐴 + ) 𝑛 𝐵 𝜖 , (18)for calculating all optical depths in this paper (where 𝜎 = . × − cm and 𝑔 𝐴 = . 𝜖 and 𝜏 opt ( 𝜖 ) computed in MNRAS000
Comparison of 𝑙 t = 𝑐 𝑖 / Ω 𝐾 (where 𝑐 𝑖 = √︁ 𝑃 / 𝜌 and Ω 𝐾 is theKeplerian angular velocity) used as characteristic length scale in the standard 𝛼 -viscosity approach with 𝑙 t = 𝜈 ).The profiles of 𝑙 t = 𝑐 𝑖 / Ω 𝐾 are computed along the equator for the indicatedtimes based on the data of model m01M3A8. Finally, region D subsumes all torus material where weak interac-tions take place on timescales much longer than the viscous (orexpansion) timescales and 𝑌 𝑒 is basically frozen out.The four different cases sketched in Fig. 10 represent dif-ferent regimes of mass accretion rates, (cid:164) 𝑀 BH . During the evolu-tion of a sufficiently massive torus, region A will disappear first,when (cid:164) 𝑀 BH (cid:46) . 𝑀 (cid:12) s − , followed by region B once the op-tical depth for neutrinos falls below unity (approximately once (cid:164) 𝑀 BH (cid:46) . 𝑀 (cid:12) s − ). Finally, region C vanishes around the freeze-out time, 𝑡 em , when (cid:164) 𝑀 BH drops below ∼ . 𝑀 (cid:12) s − . As a result,a fixed-mass torus (such as formed in a NS merger) may never ex- hibit regions A, B, or C if its initial mass is too low to achieve thecorresponding mass-accretion rates.As can be seen in Fig. 6, the turbulent behavior of the MHDmodels strongly distorts instantaneous iso-contours and sometimeseven disrupts them into multiple patches. However, these perturba-tions are stochastic in nature, hence, a classification based on Fig. 10and described above still holds for the MHD models but should thenbe interpreted in the time-averaged sense, i.e. considering propertiesthat are averaged over a few dynamical timescales. Next, since to our knowledge this quantity has rarely been system-atically reported so far for neutrino-cooled disks, we have a brieflook at the optical depth, 𝜏 opt ( 𝜖 ) = ∫ 𝜅 abs ( 𝜖 ) d 𝑟 , (17)where 𝜅 abs ( 𝜖 ) is the total opacity as function of neutrino energy 𝜖 . Although the optical depth, and therefore the neutrino emissionrates, depend on the direction in which neutrinos leave the torus, weare only interested in the maximum optical depth, which is obtainedwhen performing the integral in Eq. (17) along the equatorial di-rection through the entire computational domain. A straightforwardestimate of 𝜏 opt can be obtained by adopting for the neutrino energythe mean energy of released neutrinos measured far away from thetorus, given by 𝐿 𝜈 / 𝐿 𝜈,𝑁 . We use the approximate formula (e.g.Bruenn 1985) 𝜅 abs ( 𝜖 ) ≈ 𝜎 𝑚 𝑒 𝑐 ( 𝑔 𝐴 + ) 𝑛 𝐵 𝜖 , (18)for calculating all optical depths in this paper (where 𝜎 = . × − cm and 𝑔 𝐴 = . 𝜖 and 𝜏 opt ( 𝜖 ) computed in MNRAS000 , 000–000 (0000) Just et al. Y e → Y eq,abs e τ abs < τ em A Y e → Y eq,em e τ em < τ abs C d Y e /d t ≈ 0 τ dyn ≪ τ em , τ abs D BH torus · M BH ∼ 0.1…1 M ⊙ s −1 > Y e → Y eq e τ abs ∼ τ em B Y e ∼ Y eq,abs e τ abs < τ em B C DB · M BH ∼ 0.01…0.1 M ⊙ s −1 Y e ∼ Y eq,abs e τ abs < τ em C D · M BH ∼ 0.001…0.01 M ⊙ s −1 Y e ∼ Y eq,abs e τ abs < τ em D · M BH ∼ 0.001 M ⊙ s −1 < Figure 10.
Schematic illustration of the characteristic weak-interaction regimes encountered in neutrino-cooled disks and their corresponding equilibriumelectron fractions in dependence of the mass accretion rate onto the central BH, (cid:164) 𝑀 BH . The torus can roughly be divided into a region where only neutrinoabsorption is relevant (A), where both neutrino emission and absorption are relevant (B), where only neutrino emission is relevant (C), and where all weakinteractions are inefficient (D). The cases 1 to 4 indicate different regimes of mass accretion rates onto the BH. See Sect. 2 for the definition of the corresponding 𝑌 𝑒 -equilibria and the emission/absorption timescales, as well as Sect. 4.2.1 for a discussion of the regions A, B, C, and D. this way for various models. However, the numbers resulting in thiscase for 𝜖 and 𝜏 opt ( 𝜖 ) are systematically underrated, because theydisregard the fact that preferrably neutrinos of higher energy areabsorbed, owing to the 𝜖 dependence of 𝜅 abs . Therefore, a moreappropriate energy for measuring the impact of absorption is givenby the average energy of all neutrinos captured by nucleons per unitof time, i.e. 𝜖 abs = ∫ 𝑟<𝑟 (cid:164) 𝑒 abs d 𝑉 ∫ 𝑟<𝑟 (cid:164) 𝑛 abs d 𝑉 , (19)where (cid:164) 𝑛 abs = 𝜆 ¯ 𝜈 𝑒 𝑛 𝑝 + 𝜆 𝜈 𝑒 𝑛 𝑛 and (cid:164) 𝑒 abs is the corresponding energy-absorption rate that results after replacing 𝜖 by 𝜖 in the rates 𝜆 𝜈 𝑒 / ¯ 𝜈 𝑒 , cf. Eq. (1). The resulting changes in 𝜖 and 𝜏 opt ( 𝜖 ) ∝ 𝜖 are quite significant, approximately a factor of 2 and 4, respec-tively, during the neutrino-dominated phase. For the fiducial model,m01M3A8, with a relatively low torus mass of 0 . 𝑀 (cid:12) the opticaldepth computed in this way exceeds 10 during the first ∼
20 msand drops below 𝜏 opt = 𝑡 ≈
60 ms. The sharp depen-dence of the optical depth on the detailed neutrino energy spectrumhighlights the importance of using an energy-dependent neutrinotransport scheme for investigating optical-depth related effects inneutrino-cooled disks.We provide for each model in Table 2 the time-integrated meanenergies of released and absorbed neutrinos as well as the maxi-mum value of the optical depth attained by each model duringits evolution, 𝜏 maxopt . We find higher values of 𝜏 maxopt for models that The reason why in Table 2 we employ a fixed neutrino energy of 𝜖 =
20 MeV to compute maximum optical depths instead of 𝜖 abs from Eq. (19) lead to more compact torus configurations, namely for larger diskmasses, smaller BH masses, higher BH spins, and lower values ofthe viscous 𝛼 -parameter. An enhanced role of absorption for fasterspinning BHs has also been reported in Fernández et al. (2015). Itturns out that 𝜏 maxopt is not only a useful measure for the importance ofneutrino absorption in a given torus configuration, but it also corre-lates, though only approximately, with the average electron fractionof the torus and, therefore, of the ejecta, as can be seen in Fig. 11.The reason for this correlation is simple: High optical depths tendto be found in tori with high densities and therefore more degener-ate electron distributions and correspondingly low values of 𝑌 eq , em 𝑒 (see, e.g., Fig. 1).Lastly, to enable the comparison with popular neutrino leak-age schemes we compute another measure for the importance ofneutrino absorption, namely the quantity 𝜒 abs = ∫ 𝑟<𝑟 ( (cid:164) 𝑛 em − (cid:164) 𝑛 abs ) d 𝑉 ∫ 𝑟<𝑟 (cid:164) 𝑛 em d 𝑉 (20)that we call absorption factor here and which measures the fractionof neutrinos that are released by the system with respect to the totalnumber of produced neutrinos (see panel (l) in Figs. 2 and 3). Thisquantity is an integral version of the quenching factor commonlyused in neutrino leakage schemes (see, e.g., Ruffert et al. 1996 andthe discussion in Sect. 5.1) to approximate the effective reductionof neutrino emission rates due to optical-depths effects. The valueof 𝜒 abs vanishes if neutrino absorption exactly balances neutrino is simply to enable a straightforward comparison between all models, eventhose neglecting neutrino absorption. MNRAS , 000–000 (0000) eutrino absorption in black-hole disks < Y e > m i n t o r u s m01M3A8m1M3A8m001M3A8m01M5A8m01M10A8m01M3A4m01M3A9m01M3A8- 03m01M3A8- 1m01M3A8-vis2m01M3A8-mhdm01M3A8-noQmm01M3A8-wmm01_ye01 Y e , e j ( m e j + m f i n t o r ) / m t o r (20 MeV)0510152025303540 m i n e r t e j / m e j [ % ] (20 MeV)10 X l a n + X a c t (20 MeV)10 [ X r d / X n d ] Figure 11.
Global comparison of models visualizing data of Tables 2, 3, and 4. For each torus model we show as function of its maximum optical depththe average electron fraction of the ejected material, 𝑌 ej 𝑒 , minimum of the average electron fraction of the torus, (cid:104) 𝑌 𝑒 (cid:105) mintor , mass fraction of ejected fluid thatremained unaffected to weak interactions, 𝑚 inertej / 𝑚 ej , the sum of lanthanide and actinide mass fractions, 𝑋 LA = 𝑋 lan + 𝑋 act , and the ratio of 3rd-peak to2nd-peak abundances normalized to the solar value, [ 𝑋 ( )/ 𝑋 ( ] (cid:12) . The mass fractions are defined in the caption of Table 4. Filled (open) symbols referto models including (ignoring) neutrino absorption. emission, while it approaches unity if neutrino absorption becomesnegligible. For our disk models we find values of 𝜒 abs ∼ . − . The relative importance of neutrino absorption with respect to neu-trino emission can be assessed in a straightforward manner by adapt-ing the quantities of Sect. 4.1.1 to the analysis of absorption, suchas the average absorption timescale, (cid:104) 𝜏 abs (cid:105) = ∫ 𝑟<𝑟 𝑛 𝐵 d 𝑉 ∫ 𝑟<𝑟 (cid:164) 𝑛 abs d 𝑉 , (21)and the fraction of the torus in which significant absorption ratesare measured (i.e. where 𝜏 abs < (cid:104) 𝜏 abs (cid:105) , is an approximatemeasure of the time needed for absorption reactions to drive (cid:104) 𝑌 𝑒 (cid:105) to its equilibrium value, (cid:104) 𝑌 eq 𝑒 (cid:105) , recalling that 𝑌 eq 𝑒 is the equilibrium 𝑌 𝑒 for given thermodynamic state under the presence of neutrinos(cf. Eq. (3)). Remarkably, even for the fiducial model, m01M3A8,with a rather low initial mass of 0.01 𝑀 (cid:12) the absorption timescale is (cid:46)
10 ms initially and shorter than the viscous timescales during thefirst 20 −
30 ms of evolution. This implies that even for this low-masssystem absorption is at least at early times non-negligible in regulat-ing 𝑌 eq 𝑒 and therefore 𝑌 𝑒 . Absorption reactions become eventuallyunimportant around 𝑡 ∼ 𝑡 abs =
60 ms, where 𝑡 abs is defined as thetime when only 10 % of nucleons in the torus experience absorptionrates larger than 1 s − . When varying the global model parameters (cf. Table 2) 𝑡 abs scales about linearly with the emission freeze-outtimes, 𝑡 em , such that roughly 𝑡 abs ∼ . × 𝑡 em and 𝑡 abs (cid:46)
100 ms inmost models. The corresponding mass accretion rates measured attimes 𝑡 abs are ∼ − × − 𝑀 (cid:12) s − (cf. Table 2). 𝑌 𝑒 We now discuss the impact of neutrino absorption on the evolutionof 𝑌 𝑒 in the torus. We again take a look at Fig. 7, where torus-averaged versions are plotted of 𝑌 𝑒 (black lines), 𝑌 eq 𝑒 (green line), 𝑌 eq , em 𝑒 (red line), 𝑌 eq , abs 𝑒 (magenta line), and 𝑌 eq ,𝜇 𝜈 = 𝑒 (blue line)and models with neutrino absorption (thick lines) are comparedto models without (thin lines). In all models except the low-massmodel, m001M3A8, neutrino absorption is responsible for a signif-icant increment of (cid:104) 𝑌 eq 𝑒 (cid:105) (recalling that (cid:104) 𝑌 eq 𝑒 (cid:105) = (cid:104) 𝑌 eq , em 𝑒 (cid:105) for theno 𝜈 models), as a result of two effects.The first effect is that captures of 𝜈 𝑒 / ¯ 𝜈 𝑒 on 𝑛 / 𝑝 shift thekinetic 𝛽 -equilibrium towards less neutron-rich conditions, i.e. (cid:104) 𝑌 eq 𝑒 (cid:105) > (cid:104) 𝑌 eq , em 𝑒 (cid:105) for models with neutrino absorption. The ex-planation is not far to seek and generic for neutrino-cooled disks:Since the emission timescales are sufficiently short in the neutrino-dominated phase, the torus emits roughly the same number of 𝜈 𝑒 and ¯ 𝜈 𝑒 per unit of time. Neutrino irradiation by itself would thussaturate at 𝑌 𝑒 → 𝑌 eq , abs 𝑒 ≈ . 𝑌 eq , abs 𝑒 > 𝑌 eq , em 𝑒 and there-fore 𝑌 eq 𝑒 > 𝑌 eq , em 𝑒 . At later times, well after neutrino absorptionhas become irrelevant, this situation changes, first because of thedecline of 𝐿 𝑁 ,𝜈 𝑒 / 𝐿 𝑁 , ¯ 𝜈 𝑒 (which without recombination would leadto very low values of 𝑌 eq , abs 𝑒 ), and second because of nuclear re- MNRAS , 000–000 (0000) Just et al. combination, which ultimately drives both equilibrium quantitiesto the same value of 𝑌 eq , em 𝑒 ≈ 𝑌 eq , abs 𝑒 ≈ 𝑍 ℎ / 𝐴 ℎ ∼ .
5, where 𝑍 ℎ =
25 and 𝐴 ℎ =
54 are the charge number and mass number ofthe representative nucleus in our 4-species EOS.The second effect on 𝑌 𝑒 related to the (non-)inclusion of neu-trino absorption is of dynamical nature and is the reason why (cid:104) 𝑌 eq , em 𝑒 (cid:105) is reduced in models neglecting neutrino absorptions com-pared to models including them. This reduction stems from the factthat cooling is enhanced in those models, because neutrinos canstream out freely from the torus without experiencing scatteringsand re-absorption by the fluid. The boosted energy release ratesimply higher densities (see panel (h) in Fig. 2) and therefore moredegenerate electron distributions and ultimately lower values of 𝑌 eq , em 𝑒 . This finding illustrates that a proper treatment of the energytransport can have just as important consequences for the evolutionof 𝑌 𝑒 as the lepton transport in that an underestimation (overestima-tion) of net cooling rates leads to higher (lower) values of 𝑌 𝑒 in thetorus.Both aforementioned effects together lead to a reduction of theminimum value of the average torus 𝑌 𝑒 , (cid:104) 𝑌 𝑒 (cid:105) mintor , which is given inTable 2 and plotted in Fig. 11 for models not accounting for neutrinobackreactions. The size of this reduction grows for models with over-all higher optical depth, and it is most dramatic ( Δ (cid:104) 𝑌 𝑒 (cid:105) mintor ≈ . 𝑀 (cid:12) . Further-more, (cid:104) 𝑌 𝑒 (cid:105) mintor appears to decrease monotonically with the torusmass – or probably with any other global parameter that leads tohigher densities in the torus – only if absorption reactions are absent.If taken into account, neutrino absorptions may limit or even reversethis trend for sufficiently high optical depths, as observed for thesequence of models with increasing torus mass (compare (cid:104) 𝑌 𝑒 (cid:105) mintor of models m001M3A8, m01M3A8, and m1M3A8 with that of thecorresponding no 𝜈 models).On a final note, we point out that the blue lines in Fig. 7 thatrepresent (cid:104) 𝑌 eq ,𝜇 𝜈 = 𝑒 (cid:105) , i.e. 𝑌 𝑒 resulting if neutrinos were in thermody-namic equilibrium with vanishing chemical potential, bracket (cid:104) 𝑌 eq 𝑒 (cid:105) from above during almost the entire evolution. After having examined the evolution of 𝑌 𝑒 in the disks, we nowinvestigate how neutrino absorptions influence the properties of theejected material. Different types of outflows may be encountered inneutrino-cooled disks, e.g. neutrino-driven or viscous, or magne-tohydrodynamically launched outflows (e.g. Fernández & Metzger2013; Just et al. 2015a; Siegel & Metzger 2018; Miller et al. 2019a).In our set of rather low-mass models, viscous ejecta by far dom-inate neutrino-driven outflows (e.g. Just et al. 2015a). For MHDmodels one might be able to distinguish material ejected throughMHD-driven turbulent viscosity from ejecta expelled by means ofother MHD related effects (see, e.g., Fernández et al. 2019). How-ever, finding suitable diagnostic criteria is not trivial and beyond thescope of our study, hence we do not discriminate between differentejecta components here. We note that since our special relativistic MHD models cannot de-scribe the general relativistic Blandford-Znajek process (Blandford & Znajek1977), they do not exhibit electromagnetic jets along the polar axis.
Before considering the thermodynamic properties of the ejecta,we summarize our results concerning the ejecta masses: The totalejection efficiency, 𝑚 ej / 𝑚 , in the considered set of models (seeTable 3) is always in the ballpark of ∼ −
30 %. These valuesare reduced by (at most) a few per cent in models where neutrinofeedback is ignored, because more efficient neutrino cooling tendsto counteract the ejection processes. Furthermore, our results arein broad agreement with previous studies (e.g. Fernández 2015;Just et al. 2015a; Fernández et al. 2020; Fujibayashi et al. 2020a)concerning variations of 𝑚 , 𝑀 BH , and 𝐴 BH : 𝑚 ej / 𝑚 decreasesfor more massive tori (mostly because those evolve longer in theneutrino-cooled phase), increases with BH spin (thanks to a reduc-tion of the ISCO that favors mass ejection; see Fernández 2015), andincreases with BH mass (mainly because tori are less compact andtherefore transition earlier into the adiabatic phase). We note thatFernández et al. (2020) report an opposite trend, namely decreas-ing ejecta masses with more massive BHs, but they also increase 𝑀 BH / 𝑟 𝑑 , i.e. the disk compactness (where 𝑟 𝑑 is the radius of max-imum density in the initial torus) together with 𝑀 BH , whereas wekeep this quantity constant.In the MHD models, the mass of material ejected until thefinal simulation times, 𝑡 𝑓 , is not finalized owing to computationallimitations that have prevented us from continuing the simulation.The final ejecta mass can be estimated by assuming that all materialat radii smaller than 𝑟 = km at 𝑡 = 𝑡 𝑓 will ultimately becomeejected (cf. Tables 1 and 3), because BH accretion has basicallyceased by 𝑡 = 𝑡 𝑓 . This gives 𝑚 ej / 𝑚 ≈ . ≈
40 % reported by Fernández et al. (2019)may to some extent be explained by our weaker initial magnetic field(such a tendency is supported by the results of Christie et al. 2019),but could also be connected to our simplified treatment of GR inthe form of the pseudo-Newtonian gravitational potential.
Next we take a look at the electron fraction of the ejecta. Althoughmost features discussed in this section are generic for all models,we focus first on the model sequence with increasing torus mass,i.e. m001M3A8, m01M3A8, and m1M3A8, representing cases withnegligible, modest, and significant impact of 𝜈 -absorption, respec-tively. The dependence on other model parameters and modelingingredients is addressed in the next section. Panels (a) and (m) inFig. 2 show the mass-loss rate and the spherically averaged 𝑌 𝑒 ,respectively, measured at a fixed radius of 𝑟 = km. For the se-quence of models with increasing torus mass and employing thestandard 𝛼 -viscosity a reduction of the ejecta 𝑌 𝑒 in models withoutcompared to models with absorption is found analogous to whatwas found for the torus 𝑌 𝑒 in Sect. 4.2.4.Since the extraction radius of 𝑟 = km is far away fromthe torus, properties measured at this radius cannot inform aboutthe times when outflow material actually attained its final 𝑌 𝑒 . Thisinformation is provided in the bottom panels of Fig. 7, where thefraction of ejected material with already finalized 𝑌 𝑒 values is plot-ted as function of time ( 𝑚 FO / 𝑚 ej ) as well as the instantaneous 𝑌 𝑒 value of material freezing out at time 𝑡 ( 𝑌 FO 𝑒 ). In these plots we iden-tify material for which 𝑌 𝑒 has frozen out by the criterion 𝜏 𝛽 >
10 sfor the weak interaction timescale 𝜏 𝛽 . The figures reveal that mate-rial freezes out over an extended period of time, which is roughlycentered around 𝑡 em , and that freeze-out material is more neutron- MNRAS , 000–000 (0000) eutrino absorption in black-hole disks Table 3.
Basic properties of the material ejected until the end of each simulation: ejecta mass, ejecta mass having 𝑌 𝑒 < .
25, fraction of ejecta barely affectedby weak interactions, average electron fraction and velocity of ejecta measured at 𝑟 = km, average entropy of ejecta measured for each trajectory attemperature 𝑇 = 𝑚 ej / 𝑚 𝑚 𝑌 𝑒 < . / 𝑚 𝑚 inertej / 𝑚 ej 𝑌 𝑒, ej ( 𝑟 = km ) 𝑣 ej ( 𝑟 = km ) 𝑠 ej ( 𝑇 = ) name [%] [%] [%] [ 𝑐 ] [ 𝑘 𝑏 /baryon]m01M3A8(-no 𝜈 ) 22 (21) 0.5 (10) 8.4 (9.2) 0.322 (0.276) 0.048 (0.043) 26.8 (27.8)m1M3A8(-no 𝜈 ) 22 (17) 0.1 (11) 5.3 (7.7) 0.333 (0.241) 0.045 (0.040) 19.6 (22.1)m001M3A8(-no 𝜈 ) 28 (28) 0 (0) 24 (24) 0.405 (0.403) 0.050 (0.048) 37.3 (36.5)m01M5A8(-no 𝜈 ) 24 (23) 0 (0) 14 (19) 0.358 (0.342) 0.045 (0.045) 32.0 (30.2)m01M10A8(-no 𝜈 ) 27 (28) 0 (0) 32 (37) 0.420 (0.421) 0.046 (0.044) 39.1 (35.7)m01M3A4(-no 𝜈 ) 12 (11) 1.5 (3.3) 21 (24) 0.325 (0.303) 0.037 (0.038) 23.6 (24.9)m01M3A9(-no 𝜈 ) 28 (25) 0 (9.2) 12 (9.1) 0.331 (0.280) 0.051 (0.051) 27.4 (28.2)m01M3A8- 𝛼 𝜈 ) 19 (18) 0 (7.2) 5.1 (5.8) 0.306 (0.276) 0.033 (0.032) 27.2 (27.2)m01M3A8- 𝛼 𝜈 ) 24 (22) 3 (9.3) 13 (14) 0.338 (0.290) 0.057 (0.054) 27.4 (28.5)m01M3A8-vis2(-no 𝜈 ) 31 (26) 0 (0) 2.5 (2.9) 0.347 (0.337) 0.040 (0.048) 30.5 (33.3)m01M3A8-mhd(-no 𝜈 ) 15 (14) 1.7 (3.0) 2.5 (10.5) 0.304 (0.279) 0.125 (0.139) 23.7 (23.6)m01M3A8-noQm(-no 𝜈 ) 21 (20) 8.6 (17) 5.8 (11) 0.271 (0.241) 0.046 (0.042) 28.4 (28.9)m01M3A8-wm 23 0.6 8.2 0.318 0.044 26.4m01M3A8-ye01(-no 𝜈 ) 22 (20) 15 (16) 6.9 (14) 0.234 (0.203) 0.046 (0.047) 26.7 (30.3) Table 4.
Quantities characterizing the nucleosynthesis yields in the ejecta and kilonova light curve: Mass fractions of lanthanides (with atomic charge numberfulfilling 58 (cid:54) 𝑍 (cid:54) (cid:54) 𝐴 (cid:54) (cid:54) 𝐴 (cid:54) 𝐴 (cid:62) km at the end of the simulation will become ejected. Entries not in (in) parenthesisrefer to models evolved with (without) neutrino absorption.model 𝑋 lan 𝑋 act 𝑋 𝑋 [ 𝑋 / 𝑋 ] (cid:12) 𝑡 peak 𝐿 peak 𝑇 peak name [day] [10 erg/s] [10 K]m01M3A8(-no 𝜈 ) 0.032 (0.124) 2.9e-5 (2.1e-3) 0.280 (0.326) 8e-3 (0.118) 0.051 (0.628) 2.00 (2.69) 5.30 (3.12) 2.70 (2.39)m1M3A8(-no 𝜈 ) 3e-3 (0.093) 9.0e-8 (3.5e-3) 0.192 (0.287) 3e-4 (0.144) 0.003 (0.872) 4.23 (7.32) 24.3 (9.68) 3.02 (2.15)m001M3A8(-no 𝜈 ) 3e-3 (7e-4 ) 2.4e-5 (1.8e-5) 0.025 (0.032) 8e-4 (4e-4 ) 0.054 (0.021) 0.54 (0.49) 9.02 (10.6) 5.81 (6.03)m01M5A8(-no 𝜈 ) 8e-3 (0.021) 2.5e-5 (1.1e-5) 0.182 (0.237) 4e-3 (5e-3 ) 0.039 (0.037) 1.48 (1.72) 8.64 (6.57) 3.72 (3.20)m01M10A8(-no 𝜈 ) 2e-4 (1e-4 ) 4.8e-6 (4.8e-6) 0.011 (0.011) 2e-4 (0.000) 0.026 (8e-4 ) 0.70 (0.57) 35.5 (50.8) 6.67 (8.22)m01M3A4(-no 𝜈 ) 0.055 (0.138) 4.5e-5 (2.3e-3) 0.310 (0.294) 0.017 (0.119) 0.096 (0.699) 2.00 (2.44) 2.37 (1.72) 2.61 (2.33)m01M3A9(-no 𝜈 ) 0.028 (0.098) 1.7e-5 (1.4e-3) 0.210 (0.270) 7e-3 (0.093) 0.056 (0.592) 1.72 (2.56) 8.24 (4.16) 3.39 (2.46)m01M3A8- 𝛼 𝜈 ) 0.030 (0.074) 4.5e-5 (7.4e-4) 0.248 (0.346) 0.011 (0.062) 0.075 (0.310) 2.21 (2.98) 3.99 (2.69) 2.58 (2.55)m01M3A8- 𝛼 𝜈 ) 0.043 (0.111) 4.3e-5 (1.4e-3) 0.223 (0.267) 0.015 (0.114) 0.118 (0.736) 2.00 (2.32) 5.33 (3.78) 2.66 (2.45)m01M3A8-vis2(-no 𝜈 ) 8e-4 (3e-3 ) 1.4e-5 (3.3e-5) 0.026 (0.078) 6e-4 (2e-3 ) 0.038 (0.041) 0.90 (1.15) 28.4 (17.6) 4.78 (3.93)m01M3A8-mhd(-no 𝜈 ) 0.053 (0.080) 1.3e-3 (3.1e-3) 0.176 (0.146) 0.066 (0.125) 0.648 (1.485) 1.15 (1.10) 15.6 (17.4) 3.46 (3.23)m01M3A8-noQm(-no 𝜈 ) 0.157 (0.186) 1.7e-3 (7.9e-3) 0.337 (0.257) 0.144 (0.277) 0.739 (1.866) 2.69 (2.98) 3.61 (2.89) 2.54 (2.29)m01M3A8-wm 0.043 3.9e-5 0.307 0.016 0.089 2.10 4.75 2.83m01M3A8-ye01(-no 𝜈 ) 0.115 (0.199) 3.0e-3 (6.6e-3) 0.359 0.140 (0.244) 0.674 (1.343) 2.44 (2.56) 3.88 (3.58) 2.51 (2.44) rich for models without neutrino absorption even well beyond thetime 𝑡 abs at which neutrino absorption becomes irrelevant. This lastobservation implies that ejecta trajectories can carry to a certain de-gree memory of the conditions prevalent when neutrino absorptionswere still active, even if this material is expelled later and via viscousmechanisms that are completely unrelated to neutrino-absorption.We remark that Miller et al. (2019a) come to a similar conclusionin their investigation of neutrino-related effects in MHD disks.More insight about the 𝑌 𝑒 evolution of outflow particles alongtheir journey from deep inside the torus until free expansion can begained by considering Fig. 12, where, as functions of temperature(with inverted 𝑥 -axis to resemble the time evolution), the average 𝑌 𝑒 is plotted together with its equilibria (top row) as well as the char-acteristic timescales (bottom row) of neutrino absorption (magentalines), emission (red lines), and expansion (blue lines). This analysisis inspired by a related one performed recently by Fujibayashi et al.(2020b) to investigate NS remnants of NS-NS mergers. In order to obtain an average ¯ 𝑋 of a quantity 𝑋 at fixed temperature, we firstmap 𝑋 for each trajectory onto a temperature grid and then computea mass-weighted average of 𝑋 if, and only if, for a given temper-ature more than 50 % of ejecta (in terms of mass) can be found.For the temperature mapping we only take into account the pointconnected to the latest evolution time if a certain temperature isreached multiple times (e.g. due to turnover motions). The charac-teristic timescales plotted in the lower panels of each row in Fig. 12are obtained by plugging previously computed quantities definedat fixed temperature into the defining equations, Eqs. (6), (9), and 𝜏 expansion = ¯ 𝑟 / ¯ 𝑣 for emission, absorption, and expansion, respec-tively.The results in Fig. 12 corroborate our previous interpretationconcerning the (equilibrium) conditions of the torus 𝑌 𝑒 : At hightemperatures, i.e. deep inside the torus (region B in Fig. 10), thesoon-to-be ejected outflow material is close to the respective equi-librium, 𝑌 𝑒 → 𝑌 eq 𝑒 > 𝑌 eq , em 𝑒 ( 𝑌 𝑒 → 𝑌 eq , em 𝑒 ) for models with MNRAS000
Quantities characterizing the nucleosynthesis yields in the ejecta and kilonova light curve: Mass fractions of lanthanides (with atomic charge numberfulfilling 58 (cid:54) 𝑍 (cid:54) (cid:54) 𝐴 (cid:54) (cid:54) 𝐴 (cid:54) 𝐴 (cid:62) km at the end of the simulation will become ejected. Entries not in (in) parenthesisrefer to models evolved with (without) neutrino absorption.model 𝑋 lan 𝑋 act 𝑋 𝑋 [ 𝑋 / 𝑋 ] (cid:12) 𝑡 peak 𝐿 peak 𝑇 peak name [day] [10 erg/s] [10 K]m01M3A8(-no 𝜈 ) 0.032 (0.124) 2.9e-5 (2.1e-3) 0.280 (0.326) 8e-3 (0.118) 0.051 (0.628) 2.00 (2.69) 5.30 (3.12) 2.70 (2.39)m1M3A8(-no 𝜈 ) 3e-3 (0.093) 9.0e-8 (3.5e-3) 0.192 (0.287) 3e-4 (0.144) 0.003 (0.872) 4.23 (7.32) 24.3 (9.68) 3.02 (2.15)m001M3A8(-no 𝜈 ) 3e-3 (7e-4 ) 2.4e-5 (1.8e-5) 0.025 (0.032) 8e-4 (4e-4 ) 0.054 (0.021) 0.54 (0.49) 9.02 (10.6) 5.81 (6.03)m01M5A8(-no 𝜈 ) 8e-3 (0.021) 2.5e-5 (1.1e-5) 0.182 (0.237) 4e-3 (5e-3 ) 0.039 (0.037) 1.48 (1.72) 8.64 (6.57) 3.72 (3.20)m01M10A8(-no 𝜈 ) 2e-4 (1e-4 ) 4.8e-6 (4.8e-6) 0.011 (0.011) 2e-4 (0.000) 0.026 (8e-4 ) 0.70 (0.57) 35.5 (50.8) 6.67 (8.22)m01M3A4(-no 𝜈 ) 0.055 (0.138) 4.5e-5 (2.3e-3) 0.310 (0.294) 0.017 (0.119) 0.096 (0.699) 2.00 (2.44) 2.37 (1.72) 2.61 (2.33)m01M3A9(-no 𝜈 ) 0.028 (0.098) 1.7e-5 (1.4e-3) 0.210 (0.270) 7e-3 (0.093) 0.056 (0.592) 1.72 (2.56) 8.24 (4.16) 3.39 (2.46)m01M3A8- 𝛼 𝜈 ) 0.030 (0.074) 4.5e-5 (7.4e-4) 0.248 (0.346) 0.011 (0.062) 0.075 (0.310) 2.21 (2.98) 3.99 (2.69) 2.58 (2.55)m01M3A8- 𝛼 𝜈 ) 0.043 (0.111) 4.3e-5 (1.4e-3) 0.223 (0.267) 0.015 (0.114) 0.118 (0.736) 2.00 (2.32) 5.33 (3.78) 2.66 (2.45)m01M3A8-vis2(-no 𝜈 ) 8e-4 (3e-3 ) 1.4e-5 (3.3e-5) 0.026 (0.078) 6e-4 (2e-3 ) 0.038 (0.041) 0.90 (1.15) 28.4 (17.6) 4.78 (3.93)m01M3A8-mhd(-no 𝜈 ) 0.053 (0.080) 1.3e-3 (3.1e-3) 0.176 (0.146) 0.066 (0.125) 0.648 (1.485) 1.15 (1.10) 15.6 (17.4) 3.46 (3.23)m01M3A8-noQm(-no 𝜈 ) 0.157 (0.186) 1.7e-3 (7.9e-3) 0.337 (0.257) 0.144 (0.277) 0.739 (1.866) 2.69 (2.98) 3.61 (2.89) 2.54 (2.29)m01M3A8-wm 0.043 3.9e-5 0.307 0.016 0.089 2.10 4.75 2.83m01M3A8-ye01(-no 𝜈 ) 0.115 (0.199) 3.0e-3 (6.6e-3) 0.359 0.140 (0.244) 0.674 (1.343) 2.44 (2.56) 3.88 (3.58) 2.51 (2.44) rich for models without neutrino absorption even well beyond thetime 𝑡 abs at which neutrino absorption becomes irrelevant. This lastobservation implies that ejecta trajectories can carry to a certain de-gree memory of the conditions prevalent when neutrino absorptionswere still active, even if this material is expelled later and via viscousmechanisms that are completely unrelated to neutrino-absorption.We remark that Miller et al. (2019a) come to a similar conclusionin their investigation of neutrino-related effects in MHD disks.More insight about the 𝑌 𝑒 evolution of outflow particles alongtheir journey from deep inside the torus until free expansion can begained by considering Fig. 12, where, as functions of temperature(with inverted 𝑥 -axis to resemble the time evolution), the average 𝑌 𝑒 is plotted together with its equilibria (top row) as well as the char-acteristic timescales (bottom row) of neutrino absorption (magentalines), emission (red lines), and expansion (blue lines). This analysisis inspired by a related one performed recently by Fujibayashi et al.(2020b) to investigate NS remnants of NS-NS mergers. In order to obtain an average ¯ 𝑋 of a quantity 𝑋 at fixed temperature, we firstmap 𝑋 for each trajectory onto a temperature grid and then computea mass-weighted average of 𝑋 if, and only if, for a given temper-ature more than 50 % of ejecta (in terms of mass) can be found.For the temperature mapping we only take into account the pointconnected to the latest evolution time if a certain temperature isreached multiple times (e.g. due to turnover motions). The charac-teristic timescales plotted in the lower panels of each row in Fig. 12are obtained by plugging previously computed quantities definedat fixed temperature into the defining equations, Eqs. (6), (9), and 𝜏 expansion = ¯ 𝑟 / ¯ 𝑣 for emission, absorption, and expansion, respec-tively.The results in Fig. 12 corroborate our previous interpretationconcerning the (equilibrium) conditions of the torus 𝑌 𝑒 : At hightemperatures, i.e. deep inside the torus (region B in Fig. 10), thesoon-to-be ejected outflow material is close to the respective equi-librium, 𝑌 𝑒 → 𝑌 eq 𝑒 > 𝑌 eq , em 𝑒 ( 𝑌 𝑒 → 𝑌 eq , em 𝑒 ) for models with MNRAS000 , 000–000 (0000) Just et al. Y e m = 0.001 M , std. -viscosity Y eq,em e Y eq,abs e Y eq e Y e m = 0.01 M , std. -viscosity m = 0.1 M , std. -viscosity T [MeV]10 t i m e s c a l e s [ s ] emissionabsorptionexpansion 10 T [MeV] 10 T [MeV]0.00.10.20.30.40.50.6 Y e m = 0.01 M , l t =const. viscosity Y eq,em e Y eq,abs e Y eq e Y e m = 0.01 M , MHD T [MeV]10 t i m e s c a l e s [ s ] emissionabsorptionexpansion 10 T [MeV] Figure 12.
Averages of various quantities across the outflow trajectories at given temperatures for models m001M3A8, m01M3A8, and m1M3A8 (top row) aswell as models m01M3A8-vis2 m01M3A8-mhd (bottom row). The top panels of each column depict averages of 𝑌 𝑒 (black lines) as well as its different limitingcases (see Sect. 2 for the definitions). The bottom panels provide averages of the characteristic timescales of emission, ¯ 𝜏 em , absorption, ¯ 𝜏 abs , and expansion,¯ 𝜏 exp ≡ ¯ 𝑟 / ¯ 𝑣 . On average, the effect of neutrino emission (absorption) starts to freeze out once ¯ 𝜏 em ( ¯ 𝜏 abs ) becomes longer than ¯ 𝜏 exp . Averages are only computedfor a given temperature if in terms of mass more than 50 % of the ejecta reach that temperature, hence ¯ 𝑌 𝑒 does not start at the initial value of 𝑌 𝑒 = . (without) neutrino feedback. Soon after the ejection sets in and ma-terial leaves the densest layers of the torus, the expansion timescalebecomes the shortest of the three timescales, which on average hap-pens at temperatures 𝑇 ≈ 𝑌 𝑒 , starts to divergefrom its equilibrium value, which continues to rise as a result ofdecreasing densities and temperatures. After a few more neutrinointeractions and at slightly lower temperatures, ¯ 𝑌 𝑒 then finally levelsoff to remain constant.While the qualitative behavior of ¯ 𝑌 𝑒 is similar for both typesof models, with and without absorption, quantitative differencesappear that are all the more pronounced for more massive tori. Thereasons are in agreement with the findings of Sect. 4.2.4: First,the starting values of 𝑌 𝑒 right before ejection are higher to beginwith for more massive tori, namely about ¯ 𝑌 eq 𝑒 − ¯ 𝑌 eq , em 𝑒 ∼ . − .
1. Second, as can be deduced from the plotted timescales ofemission and absorption, neutrino absorptions provide an additional boost in raising ¯ 𝑌 𝑒 also after ¯ 𝑌 𝑒 starts to decouple from ¯ 𝑌 eq 𝑒 . Inmodel m1M3A8, the absorption rates even dominate emission ratesduring this phase, indicating that a significant fraction of the ejectahave experienced conditions as in region A of Fig. 10 during theirexpansion. The third and last reason for elevated 𝑌 𝑒 values in modelswith absorption is that 𝑌 eq , em 𝑒 itself is already higher as a result offinite diffusion timescales out of the torus and therefore less efficientneutrino cooling compared to the no-absorption case.The top row of Fig. 13 shows the ejecta mass distribution in 𝑌 𝑒 ,measured at 𝑟 = km, as well as the abundance patterns of nucle-osynthesis yields produced in the ejecta for the three models withincreasing torus mass. In agreement with many existing results forviscous models of neutrino-cooled disks (e.g. Just et al. 2015a; Wuet al. 2017; Fujibayashi et al. 2020a), most models show a peak closeto 𝑌 𝑒 ∼ . − .
3, whereas the left boundary of the 𝑌 𝑒 distribution issystematically shifted to lower values for models without neutrino MNRAS , 000–000 (0000) eutrino absorption in black-hole disks absorption, more so for higher torus masses. Since the productionof lanthanides (and heavier elements) is typically activated for 𝑌 𝑒 near 0 .
25 (e.g. Kasen et al. 2015), the nucleosynthesis pattern ofthe heavier elements (with mass numbers 𝐴 (cid:38) 𝑌 𝑒 = . − .
3. It is therefore not surprising to observe amodest (substantial) boost in heavy-element and, in particular, lan-thanide production when ignoring neutrino absorption for a torusmass of 0.01 𝑀 (cid:12) (0.1 𝑀 (cid:12) ). Remarkably, in both cases a nearly per-fect solar-like abundance pattern (shown as circles in Fig. 13) isobtained only when neglecting neutrino absorption. Another inter-esting observation is that when including neutrino absorption, andonly when doing so, the efficiency of heavy-element productiongrows non-monotonically with torus mass: The 0.01 𝑀 (cid:12) torus inthe fiducial model m01M3A8 produces more heavy elements rel-ative to its initial mass than both the lighter and heavier tori inmodels m001M3A8 and m1M3A8, respectively. We note that Justet al. (2015a); Fernández et al. (2020); Fujibayashi et al. (2020a)report a consistent trend for growing torus masses. We now analyze the sensitivity of the ejecta properties and the im-pact of absorption to the variation of model parameters or modelingassumptions. Global properties of the ejecta can be found for allmodels in Table 3, namely the total mass, 𝑚 ej , average electronfraction, 𝑌 𝑒, ej , velocity, 𝑣 ej , and entropy per baryon, 𝑠 ej . Table 4provides mass fractions of characteristic classes of elements, e.g.lanthanides and actinides. Some of the results of Tables 2, 3, and 4are visualized in Fig. 11 as function of torus optical depth and inFig. 14 as differences relative to the fiducial model. Black hole properties.
When increasing the BH mass while keep-ing 𝑟 𝑑 ∝ 𝑀 BH (where 𝑟 𝑑 is the radius at maximum density of theinitial disk) the initial disk becomes less compact. As a consequence,the densities and temperatures are lower for more massive BHs thanfor less massive BHs, and therefore neutrino emission favors anequilibrium with higher electron fractions 𝑌 eq , em 𝑒 (cf. Fig. 1). More-over, due to the reduced optical depth the overall impact of neutrinoabsorption becomes less significant. Roughly speaking, increasing 𝑀 BH has a similar effect on the 𝑌 𝑒 evolution as decreasing the torusmass.On the other hand, varying the spin parameter of the BH be-tween 0 . < 𝐴 BH < . 𝑟 𝑑 =const.) appears to resultin relatively small differences concerning the 𝑌 𝑒 evolution, whichis likely related to the fact that the spin parameter does not lead todramatically different thermodynamic conditions and optical depthsin the torus. We note, however, that the sensitivity could be moresignificant for more massive tori, in which a larger fraction of theejecta is launched as neutrino-driven winds (Fernández 2015).We caution the reader that since our approximate initial diskmodels are constructed by hand in a way to reproduce the mostbasic global parameters (i.e. 𝑚 , 𝑀 BH , and 𝐴 BH , as well as ini-tial electron fraction and entropy), the local configuration does notperfectly agree with actual post-merger disks or disks in collapsarengines that exhibit the same values of 𝑚 , 𝑀 BH , and 𝐴 BH . Forinstance, as pointed out by Fernández et al. (2020), the initial sizeof the disk, approximately measured by the radius of the densitymaximum 𝑟 𝑑 , represents an additional relevant degree of freedomcharacterizing the disk distribution. Hence, some of the tendenciesas function of 𝑀 BH and 𝐴 BH may depart from those found herewhen considering different initial disk distributions. Finite-mass corrections and weak magnetism.
In modelsm01M3A8-noQm(-no 𝜈 ) we set the electron mass 𝑚 𝑒 = 𝑄 𝑛𝑝 = 𝛽 -reaction rates,cf. Eqs. (1). In Sect. 2.3 we already found that this simplificationcan underestimate the electron fractions corresponding to emissionequilibrium, 𝑌 eq , em 𝑒 , by up to ∼ . 𝜌 − 𝑇 regimecovered by the torus during the neutrino-dominated phase of evolu-tion. It is therefore not surprising to find more neutron-rich ejectafor these models, in which 𝑌 𝑒, ej is reduced by about 0.05. The nucle-osynthesis pattern matches very well the solar abundance pattern.Incidentally, for the fiducial model the noQm simplification has asimilar net effect as ignoring neutrino absorption.Considering now model m01M3A8-wm, in which weak mag-netism is taken into account in contrast to the fiducial modelm01M3A8, the differences appear to be minor and the reductionof the ejecta 𝑌 𝑒 is only 0.004. The small size of the impact of weakmagnetism was anticipated already in Sect. 2.3 and can be ascribedto the relatively low neutrino energies involved in the disk evolution. 𝑙 𝑡 =const. viscosity treatment. Probably the most challengingmodeling ingredient, apart from neutrino transport, is connected tothe description of turbulent angular momentum transport. For thepurpose of assessing how the 𝑌 𝑒 evolution and the impact of absorp-tion depends on the treatment of angular momentum transport, wedid not just vary the viscous 𝛼 -parameter, but also included modelswith two different viscosity prescriptions, Eqs. (13) and (14), as wellas a 3D MHD model. The first and quite remarkable observationwhen comparing the 𝑌 𝑒 distributions and abundance patterns is thatthe ejecta properties appear to be less sensitive to varying the viscous 𝛼 -parameter than they are to the choice of the turbulent viscosityscheme (i.e. standard 𝛼 -viscosity vs. 𝑙 t =const. viscosity vs. MHD).This result supports the notion that each viscosity scheme generatesits own characteristic flow pattern with correspondingly differentejecta properties. In general, this notion implies that it might beimpossible to reproduce all ejecta features of MHD models by justvarying the 𝛼 -parameter in a certain viscosity scheme. Fernándezet al. (2019) come to a similar conclusion when comparing viscousmodels with an MHD model.We first discuss the alternative viscosity prescription expressedby Eq. (14). First of all, not only the average 𝑌 𝑒 is enhanced com-pared to the standard 𝛼 -viscosity model, but, probably even moreimportantly, also the bottom value in the 𝑌 𝑒 –mass distribution,which lies at about 𝑌 𝑒 ≈ .
3. As a consequence, for this modelthe abundances of elements with
𝐴 >
130 are dramatically re-duced. Moreover, the overall difference between models with andwithout absorption is relatively small compared to the standard vis-cosity. Why is that? As already mentioned in Sect. 4.1.2 the effectivestrength of the 𝑙 t =const. viscosity employed in model m01M3A8-vis2 is reduced at later times compared to the case m01M3A8 usingthe standard 𝛼 -viscosity. This entails a slower and more gradualfreeze out and emission of outflow material (cp. middle panelsin Fig. 7 with left panels in Fig. 8). However, since neutrino ab-sorptions become inefficient roughly around the same time in bothmodels ( 𝑡 abs =
75 ms and 60 ms in the vis2 model and the referencemodel, respectively) the overall relevance of absorptions is natu-rally expected to be smaller in the vis2 model. This aspect, namelythat more time is spent in a state where only neutrino emission isefficient, is one reason for the low sensitivity of neutrino absorptionin this model. Another, related reason is connected by the circum-stance that the equilibrium value for neutrino emission, 𝑌 eq , em 𝑒 , istime-dependent and grows from 𝑌 eq , em 𝑒 ∼ . 𝑌 eq , em 𝑒 ∼ 𝑌 eq , abs 𝑒 during the disk evolution, whereas neutrino absorption pushes 𝑌 𝑒 MNRAS , 000–000 (0000) Just et al. Y e M [ M ] m =0.1 Mm =0.01 Mm =0.001 M -6 -5 -4 -3 -2 -1 A M a ss f r a c t i on Y e M [ M ] M BH = 3 MM BH = 5 MM BH = 10 M -6 -5 -4 -3 -2 -1 A M a ss f r a c t i on Y e M [ M ] A BH = 0.4 A BH = 0.8 A BH = 0.9 -6 -5 -4 -3 -2 -1 A M a ss f r a c t i on Y e M [ M ] reference model Q np = m e = 0weak. magn. included -6 -5 -4 -3 -2 -1 A M a ss f r a c t i on Y e M [ M ] vis =0.03 vis =0.06 vis =0.1 -6 -5 -4 -3 -2 -1 A M a ss f r a c t i on Y e M [ M ] std. -vis. l t =const. -vis.MHD -6 -5 -4 -3 -2 -1 A M a ss f r a c t i on Figure 13.
Left column:
Histograms for the mass distribution versus electron fraction 𝑌 𝑒 as measured at radii of 10 km in the hydrodynamic simulations. Right column:
Corresponding abundance distributions of nuclei synthesized in the ejecta as function of mass number, 𝐴 . The colors refer to the same modelsthat are plotted on the left. Mass fractions corresponding to the models are normalized to sum up to unity, while the solar abundance pattern (depicted by opencircles) is normalized to the 𝐴 =
130 mass fraction of model m01m3A8. In all panels the thick (thin) lines are used for models including (neglecting) neutrinoabsorption. The black lines always refer to the same, fiducial model, m01M3A8(-no 𝜈 ). From top to bottom (only) the following ingredients are varied withrespect to those of the fiducial model: Initial torus mass, black hole mass, black hole spin, viscous 𝛼 parameter, neutrino interaction physics ( 𝑄 𝑛𝑝 and 𝑚 𝑒 corrections (green lines) and weak magnetism correction (red lines)), treatment of turbulent viscosity ( 𝑙 t =const. viscosity (green lines), and MHD (red lines)).MNRAS , 000–000 (0000) eutrino absorption in black-hole disks Y e , ej no -absorption m / M : 0.01 0.1 m / M : 0.01 0.001 M BH / M : 3 10 A BH : 0.4 0.9 vis : 0.03 0.1std. -vis. l t =const. vis.std. -vis. MHD Q np = m e = 0weak magn. Y e : 0.5 0.1 change of ejecta Y e
15 10 5 0 5 10 15( m ej + m fintor )/ m [%] change of ejecta mass Figure 14.
Bar diagram visualizing the absolute change of the average ejecta 𝑌 𝑒 (left panel) and the final ejecta mass (right panel) with variation of theindicated model parameters and modeling assumptions. Differences are measured with respect to the fiducial model, m01M3A8 (for which { 𝑌 𝑒, ej , ( 𝑚 ej + 𝑚 fintor )/ 𝑚 } = { . ,
22 % } ), except for the variation of the BH spin and 𝛼 vis , where the reference models are m01M3A4 ( { . ,
12 % } ) and m01M3A8- 𝛼 { . ,
19 % } ), respectively. Y eq,em e > ( em < 1 s)0.00.20.40.60.81.0 m F O / m e j Y e at freeze out vs. torus Y e std. -viscosity l t =const. viscosityMHD Figure 15.
Total mass of all trajectories for which 𝑌 𝑒 is frozen out (i.e. withweak interaction timescales 𝜏 𝛽 >
10 s) as function of the equilibrium 𝑌 𝑒 averaged over regions with efficient neutrino emission (i.e. where 𝜏 em < 𝑌 𝑒 , whereas in the model with 𝑙 t =const.viscosity efficient freeze out only takes place at a higher torus 𝑌 𝑒 . almost at all times towards 𝑌 eq , abs 𝑒 ∼ .
5. In other words, for fluidelements that freeze out at a very late stage of evolution it does notmatter whether neutrino absorption is still efficient or not, because 𝑌 𝑒 in any case saturates at some high value (cid:38) . − . 𝑌 𝑒 evolution for modelswith different viscosity may be misleading because of the differ- ent viscous accretion timescales of each model. For that reason weplot in Fig. 15 the cumulative freeze-out masses, 𝑚 FO , of materialthat becomes ejected as a function of the average equilibrium 𝑌 𝑒 of the torus, (cid:104) 𝑌 eq , em 𝑒 (cid:105) , i.e. 𝑌 eq , em 𝑒 averaged over the region where 𝜏 em < (cid:104) 𝑌 eq , em 𝑒 (cid:105) grows monotonically, except at very earlyand very late times, this quantity is a more suitable measure ofthe weak-interaction regime of the torus than the time coordinate.Looking at Fig. 15 we find additional support for the arguments ofthe previous paragraph: In the 𝑙 t =const. model the freeze out takesplace at systematically higher values of (cid:104) 𝑌 eq , em 𝑒 (cid:105) compared to bothalternative descriptions of angular momentum transport. MHD treatment.
Before looking at the results of the MHD mod-els, we first comment on some technical difficulties and caveats.First, we point out that in some previous studies of neutrino-cooledMHD disks, outflows crossing a certain radius before a certaintime were removed from the analysis (e.g. Siegel & Metzger 2018;Miller et al. 2019a) in order to minimize spurious effects relatedto the early transient, during which the disk establishes a new dy-namical equilibrium. While such a procedure is well motivated, itentails the risk of cutting out too much or too little material be-cause of the difficulty to unambiguously demarcate spurious fromproper ejecta. Unfortunately, this risk is particularly high for a low-mass torus that we consider herein, because the time during which 𝑌 eq 𝑒 (cid:46) .
25 is relatively short and largely overlaps with the initialtransient phase (cf. top right panel in Fig. 8), which probably lastsuntil about 𝑡 (cid:46) −
20 ms. Moreover, it is a priori not clear whethersuch a cut-off truly captures all spurious ejecta, because it mayhappen that during the initial transient fluid elements with low 𝑌 𝑒 are artificially elevated to a higher, but still gravitationally boundorbit, on which weak interaction rates are low and 𝑌 𝑒 thus remainsalmost constant. Since those fluid elements are most likely ejectedduring the subsequent evolution, the early transient could thus havea non-negligible impact also on the late-time ejecta. Thus, we donot disregard ejecta based on the time of ejection in this study. Asecond point to mention is, however, that we remove from the out- MNRAS000
20 ms. Moreover, it is a priori not clear whethersuch a cut-off truly captures all spurious ejecta, because it mayhappen that during the initial transient fluid elements with low 𝑌 𝑒 are artificially elevated to a higher, but still gravitationally boundorbit, on which weak interaction rates are low and 𝑌 𝑒 thus remainsalmost constant. Since those fluid elements are most likely ejectedduring the subsequent evolution, the early transient could thus havea non-negligible impact also on the late-time ejecta. Thus, we donot disregard ejecta based on the time of ejection in this study. Asecond point to mention is, however, that we remove from the out- MNRAS000 , 000–000 (0000) Just et al. flow analysis all ejecta crossing the sphere at 𝑟 = km at polarangles 𝜃 / 𝜋 < . 𝜃 / 𝜋 > .
9. Owing to our coarse 𝜃 -resolutionnear the poles we noticed an artificial drag of outflow material to-wards the poles, and since the mass- and energy-density have tobe reset near the poles whenever reaching low values, the hydrody-namic properties of axis-near material cannot be trusted. Therefore,all diagnostic quantities reported in this paper, except the ejectamasses ( 𝑚 ej and 𝑚 𝑌 𝑒 < . , cf. Table 3), neglect material ejected inthe aforementioned polar-angle intervals. The third issue arising inthe MHD models is that, as already explained in Sect. 4.1.3, theno-absorption model shows a throttled MRI activity, most likelybecause of its systematically reduced disk thickness. Hence, thecomparison to the case including absorption is somewhat distorted.Based on similar arguments as raised in the previous comparisonbetween different viscosity prescriptions, it can be expected thata fully developed MRI and therefore more efficient angular mo-mentum transport in the “no 𝜈 ” model would have facilitated earliermass ejection and more neutron-rich conditions in the ejecta. Ifthis assessment is correct, then the consequence would have been amore significant difference between the two cases with and withoutneutrino absorption than currently observed.We now discuss the obtained results, keeping the aforemen-tioned issues in mind. Compared to both types of viscous models,the ejecta 𝑌 𝑒 –mass distribution for the MHD models (red lines inbottom panel of Fig. 13) is broader and reaches down to significantlylower values of 𝑌 𝑒 even though the average 𝑌 𝑒 of the ejecta is not toodifferent (cf. Table 3). Even when including neutrino absorptions,the 𝑌 𝑒 distribution seems to provide sufficiently neutron-rich con-ditions to produce heavy elements with a final abundance patternthat is almost resembling the solar one. The bottom right panels ofFig. 12 as well as Fig. 15 suggest that the low end of the 𝑌 𝑒 distribu-tion is connected to fluid elements that freeze out rather early, withina few tens of milliseconds when the equilibrium conditions in thetorus are still more neutron rich. Matter ejection at early times thusappears to be more powerful in the MHD model than in the viscousmodels. This tendency, which was also reported by Fernández et al.(2019), is a likely consequence of the fact that viscous models inthe neutrino-dominated phase are, by design, barely turbulent at all.Most of the material ejected in viscous models is launched only onceor after neutrino cooling has become inefficient at about 𝑡 ∼ 𝑡 em (e.g. lower panels of Figs. 7 and 8). Since the neutrino-dominatedphase is typically characterized by more neutron-rich equilibriumconditions, turbulent MHD tori thus seem to provide systematicallymore favorable conditions for the ejection of low- 𝑌 𝑒 material thanviscous tori. Overall beneficial for the ejection of low- 𝑌 𝑒 materialare also the shorter expansion timescales (cf. bottom right panel inFig. 12), i.e. higher ejecta velocities compared to the viscous mod-els. Apart from these differences, the global characteristics of ourMHD models are remarkably similar to the viscous models, e.g. thebehavior of (cid:104) 𝑌 eq 𝑒 (cid:105) and the freeze-out 𝑌 𝑒 (cf. Fig. 8). Initial electron fraction.
Most of our models start with a ratherlarge electron fraction of 𝑌 𝑒 = .
5, which was chosen intentionallyin order to study precisely the r-process material that is producedself-consistently during and due to the disk evolution. While 𝑌 𝑒 = . 𝑌 𝑒 ∼ .
1. For a direct comparison,we have also set up a model with 𝑌 𝑒 = . 𝑚 = . 𝑀 (cid:12) , 𝑀 BH = 𝑀 (cid:12) , and 𝐴 BH = . 𝑌 𝑒 of 0.23 (compared to 0.32 for 𝑌 𝑒 = . 𝑌 eq 𝑒 . In order to estimate thefraction of the ejecta, 𝑚 inertej / 𝑚 ej , that effectively remains unaffectedby (or inert to) weak interactions we compute along each outflowtrajectory the weak interaction timescale, 𝜏 𝛽 = ( 𝜏 − + 𝜏 − ) − , aswell as the expansion timescale, 𝜏 exp = 𝑟 / 𝑣 𝑟 , and we count theoutflow particle as unaffected if during its evolution from 𝑡 = 𝑟 = km 𝜏 𝛽 never becomes shorter than thetime Δ 𝑡 that the particle spent in the torus before its ejection. Weapproximate Δ 𝑡 by the total time during which 𝜏 𝛽 < 𝜏 exp . We findvalues of 𝑚 inertej / 𝑚 ej (see Table 3) of about O(
10 % ) for most ofour models. This implies that the nucleosynthesis yields for ten(or tens of) percent of the final ejecta are more determined by theinitial conditions than by the conditions during the disk evolution.As expected, the inert fraction of the ejecta turns out to be larger fordisks with lower optical depth (e.g. for lower disk mass, higher BHmass). In the last part of this study we briefly examine the bolometric kilo-nova light curves associated with the disk models. We employ themethods outlined in Sect. 3.3 and Appendix B. The primary point ofthis exercise is not to provide detailed predictions for observers, butrather for modelers to get some basic idea about the yet poorly ex-plored sensitivity of the three main kilonova properties (bolometricluminosity, photospheric temperature, photospheric velocity) withrespect to modeling variations considered in this study. In orderto remedy the fact that in some models the ejecta masses are notconverged yet at the times when the simulation had to be stopped, 𝑡 𝑓 , and the heating rate would thus be underestimated, we rescaleall trajectory masses by such that their sum equals the ejecta plustorus mass at 𝑡 = 𝑡 𝑓 .In Fig. 16 we plot for five models the mass distribution, 𝑚 ( 𝑣 ) ,electron fraction, 𝑌 𝑒 , lanthanide plus actinide fraction, 𝑋 LA , andphoton specific opacity, 𝜅 , as function of the expansion velocity, 𝑣 / 𝑐 , as well as the photospheric temperatures, 𝑇 phot and the photo-spheric velocities, 𝑣 phot / 𝑐 , as functions of time. The quantities 𝑇 phot and 𝑣 phot are given by the temperatures and velocities at the locationwhere the instantaneous optical depth of the ejecta equals unity. Thebolometric luminosities, computed as 𝐿 bol ( 𝑡 ) = 𝜋 ( 𝑣𝑡 ) 𝐹 ( 𝑡 ) withradiation flux density 𝐹 ( 𝑡 ) and at constant velocity 𝑣 = . 𝑐 , as wellas the integral heating rates, 𝑄 ≡ ∫ 𝜌𝑞 d 𝑉 , are depicted for the samemodels in Fig. 17. Thin lines in Figs. 16 and 17 again denote the cor-responding models that were evolved without neutrino absorption.The ejecta clearly exhibit a radial structure, with large variations in 𝑋 LA and 𝜅 . As expected, the opacities for models ignoring absorp-tion tend to be enhanced compared to models including absorption,all the more for models with increasing optical depth and thereforehigher sensitivity to neutrino absorption.Some of our light curves exhibit maximum values for the bolo-metric luminosities well before the ejecta become optically thin,namely at 𝑡 (cid:46) . 𝐿 bol , solid lines in Fig. 17) becomes equal to theinstantaneous heating rate ( 𝑄 , dashed lines in Fig. 17). At that timea large fraction of the ejecta becomes optically thin and radiation, MNRAS , 000–000 (0000) eutrino absorption in black-hole disks m =0.001 M std. -viscosity m ( v ) [ M ] Y e log X LA (1day) [cm /g] T phot [K] v phot / c m =0.01 M std. -viscosity m =0.1 M std. -viscosity m =0.01 Ml t =const.viscosity v / c m =0.01 M MHD v / c v / c v / c time [day]10 time [day]10 Figure 16.
Spherically averaged ejecta properties (first four columns from left), namely mass distribution, electron fraction, lanthanide plus actinide fraction,and opacity measured at 𝑡 = 𝑣 are obtained by assuming that 𝑣 remains constant for each trajectory after crossing 𝑟 = km. The high-velocity tailthat is more pronounced for more massive, viscous tori (i.e. up to 𝑣 ≈ . 𝑐 in model m1M3A8) represents a neutrino-driven wind, which owing to its briefappearance and relatively small mass is sampled rather poorly, explaining the step-like features in 𝑚 ( 𝑣 ) . The vertical arrows mark the times when the ejectastart to become optically thin, i.e. when the bolometric luminosity equals the instantaneous heating rate. which was produced earlier but remained trapped, becomes releasedat once, creating an excess of 𝐿 bol with respect to 𝑄 . This suddenrelease of radiation, i.e. the so-called diffusion wave Waxman et al.(e.g. 2018); Kasen & Barnes (e.g. 2019); Hotokezaka & Nakar (e.g.2020), has been hypothesized by Hotokezaka & Nakar (2020) toexplain the light curve steepening observed after 𝑡 ∼ 𝐿 bol = 𝑄 , and we provide the lu-minosity, temperature, and photospheric velocity at this time for allmodels in Table 4. Using this definition of the peak, we measure typi-cal peak times of a few days, peak luminosities of 10 − erg s − ,and peak temperatures of 1 − × K. When comparing models with and without neutrino absorption, we observe relative differ-ences for the peak luminosities and times that are scattered around ∼ −
80 %, where models including absorption throughout ex-hibit earlier and brighter peaks due to their systematically reducedlanthanide content. Moreover, models without absorption producea peak at lower temperatures, i.e. in a red-shifted color regime.The relative differences are similarly large when comparing modelswith a different treatment of angular momentum transport (see rightpanel in Fig. 17): Compared to the standard 𝛼 -viscosity treatment,the 𝑙 t =const. viscosity produces an earlier and more luminous peakwith photospheric temperatures enhanced by a factor of about two,which is explained by the systematically reduced fraction of lan- MNRAS000
80 %, where models including absorption throughout ex-hibit earlier and brighter peaks due to their systematically reducedlanthanide content. Moreover, models without absorption producea peak at lower temperatures, i.e. in a red-shifted color regime.The relative differences are similarly large when comparing modelswith a different treatment of angular momentum transport (see rightpanel in Fig. 17): Compared to the standard 𝛼 -viscosity treatment,the 𝑙 t =const. viscosity produces an earlier and more luminous peakwith photospheric temperatures enhanced by a factor of about two,which is explained by the systematically reduced fraction of lan- MNRAS000 , 000–000 (0000) Just et al. time [day]10 L b o l , Q [ e r g / s ] bol. luminosity, different torus masses L bol Q m = 0.001 Mm = 0.01 Mm = 0.1 M time [day]10 L b o l , Q [ e r g / s ] bol. luminosity, different ang. mom. transport L bol Q m = 0.01 M , std. -vis. m = 0.01 M , l t =const. vis. m = 0.01 M , MHD Figure 17.
Bolometric kilonova luminosity, 𝐿 bol (solid lines), and total radioactive heating rate (including the effect of thermalization), 𝑄 = ∫ 𝜌𝑞 d 𝑉 (dashedlines), of ejecta computed using post-processed trajectory data of the hydrodynamic models. Left panel:
For models with increasing initial torus mass, namelym001M3A8(-no 𝜈 ) (blue lines), m01M3A8(-no 𝜈 ) (black lines), and m1M3A8(-no 𝜈 ) (orange lines). Right panel:
For models with different treatment of angularmomentum transport, namely m01M3A8(-no 𝜈 ) (black lines), m01M3A8-vis2(-no 𝜈 ) (blue lines), and m01M3A8-mhd(-no 𝜈 ) (orange lines). Thick (thin) linesbelong to models including (neglecting) neutrino absorption. time [day]10 L b o l , Q [ e r g / s ] bol. luminosity, different simplifications L bol Q full v -dependence X LA ( v )= < X LA > , q ( v )= < q > X LA ( v )= < X LA > , q ( v )= < q > , m ( v ) v Figure 18.
Same as Fig. 17 but only for model m1M3A8 and additionallyshowing the results for two simplifications made in the computation ofthe light curve: The blue line depicts the case where the composition andradioactive heating rate are taken to be homogeneous throughout the ejecta.The case represented by the orange line additionally replaces the originalmass distribution by a power law such that 𝑚 ( 𝑣 ) ∝ 𝑣 − while keeping thetotal mass unchanged. thanide and actinide elements. The results for the MHD model endup somewhat in between those for the two different viscosities. Adistinct feature of the MHD models is that the ejecta span a larger velocity range, up to 𝑣 ∼ . 𝑐 , whereas ejecta in the viscous modelsdo not exceed 𝑣 ∼ . 𝑐 . This fast ejecta tail is probably the reasonfor the earlier and brighter peak in the MHD model compared tothe fiducial model using the standard 𝛼 viscosity despite the (for themost part) higher opacities.A particular advantage of our kilonova model compared toothers using manually constructed ejecta properties is the fact thatthe distribution of mass, heating rates and lanthanide fractions along 𝑣 comes directly from hydrodynamical simulations. For this reasonit is worth exploring in Fig. 18 for model m1M3A8 the differencesthat would result when averaging 𝑋 LA ( 𝑣 ) and 𝑞 ( 𝑣 ) over the ejecta(blue line) and when additionally replacing the original 𝑚 ( 𝑣 ) bya power-law as 𝑚 ( 𝑣 ) ∝ 𝑣 − (orange line). Both approximationsinduce considerable changes on the ∼
50 % level in the luminosity,at least at times 𝑡 (cid:54)
10 days, i.e. before the optically thin phase ofemission.Even though the present examination of the kilonova is keptrather brief and simple and does not address spectroscopic proper-ties, the results clearly advocate the importance of accurate neutrino-transport modeling as well as a careful treatment of turbulent angularmomentum transport for reliable predictions of the kilonova lightcurve.
MNRAS , 000–000 (0000) eutrino absorption in black-hole disks Table 5.
Compilation of average electron fraction of the ejected material, 𝑌 𝑒, ej , for models of neutrino-cooled disks available in the literature that consider thesame astrophysical setup and that differ in modeling aspects investigated in this study, namely the initial electron fraction, 𝑌 𝑒 , turbulent viscosity treatment,type of neutrino treatment, and inclusion or omission of 𝑄 𝑛𝑝 and 𝑚 𝑒 terms. The names used for the neutrino treatment have the following meaning: “leakage” → classical leakage scheme as in Ruffert et al. (1996). “leak.+abs.” → same as before but augmented with a scheme to describe net 𝜈 -absorption in opticallythin regions. “leak. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 )” → simplified leakage with 𝜒 𝜈 𝑒 = 𝜒 ¯ 𝜈 𝑒 and therefore 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 (see Sect. 5.1.1). “grey M1+leak.” → combination ofschemes in which leakage (M1) is employed in regions of high (low) optical depth. “no abs. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 )” → neutrino absorption is completely ignored. SeeSect. 5.1 for more detailed explanations. The ’*’ indicates that the given value of 𝑌 𝑒, ej refers only to early ejecta produced within 𝑡 (cid:46) −
200 ms. The ’**’means the same as ’*’, but additionally the value of 𝑌 𝑒, ej , as it was not provided, had to be estimated based on other information found in that reference.Reference 𝑚 𝑀 BH 𝐴 BH 𝑌 𝑒 viscosity neutrino 𝑄 𝑛𝑝 and 𝑚 𝑒 𝑌 𝑒, ej [ 𝑀 (cid:12) ] [ 𝑀 (cid:12) ] treatment treatment included? [1]Fernández et al. (2020) 0.03 3 0.8 0.1 std. 𝛼 -vis. leak.+abs. yes 0.28Fernández et al. (2019) 0.03 3 0.8 0.1 std. 𝛼 -vis. leak. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 ) no 0.20Just et al. (2015a) 0.03 3 0.8 0.1 std. 𝛼 -vis. spectral M1 yes 0.27m01M3A8 (this work) 𝛼 -vis. spectral M1 yes 0.32m01M3A8-noQm-no 𝜈 (this work) 𝛼 -vis. no abs. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 ) no 0.24Fujibayashi et al. (2020a) 0.1 3 0.8 0.07-0.5 𝑙 t =const. vis. grey M1+leak. yes 0.310.1 3 0.8 0.07-0.5 𝑙 t =const. vis. no abs. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 ) yes 0.30m01M3A8-vis2 (this work) 𝑙 t =const. vis. spectral M1 yes 0.35m01M3A8-vis2-no 𝜈 (this work) 𝑙 t =const. vis. no abs. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 ) yes 0.34Siegel & Metzger (2018) 0.03 3 0.8 0.1 MHD leakage no 0.18Siegel et al. (2019) 0.016 3 0.8 0.5 MHD leakage no (cid:46) . ∗∗ Fernández et al. (2019) 0.03 0.8 0.1 MHD leak. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 ) no 0.16Miller et al. (2019b) 0.12 2.58 0.69 0.1 MHD Boltzmann yes ∼ . − . ∗∗ Miller et al. (2019a) 0.02 3 0.8 0.5 MHD Boltzmann yes 0 . ∗ m01M3A8-mhd (this work) Given the relevance for heavy-element production, the availablenumber and degree of sophistication of models for neutrino-cooledBH accretion disks is growing quickly. In the following we brieflyput the results of our study into context by comparison with a varietyof models available in the literature. The basic model setup usedhere (equilibrium torus, 𝛼 -viscosity or single-loop initial magneticfield configuration) is similar to what was used in most previousstudies of viscous or MHD disks, and several features are in broadagreement with the literature. One difference to several existingmodels is that we initiated most of our simulations with 𝑌 𝑒 = . 𝑌 𝑒 , 𝑌 𝑒, ej , is overall shifted compared to models starting with a 𝑌 𝑒 = . 𝑌 𝑒 should be similar in our simulations compared tosimulations from the literature. Note that in this case, since 𝑚 ( 𝑣 ) diverges at small velocities, a minimumvelocity, 𝑣 , must be imposed, below which d 𝑚 / d 𝑣 = 𝑣 is fixed by thecondition that the average velocity of the power-law distribution is the sameas for the original distribution. Before considering individual studies we first clarify some impor-tant points regarding neutrino leakage schemes, which otherwisemay become sources of confusion. Classical leakage schemes, suchas introduced by Ruffert et al. (1996) and Rosswog & Liebendör-fer (2003), are often mentioned to be unable to describe neutrinoabsorption. This is however only true for regions with positive net absorption rates (i.e. absorption minus emission rates, such as inregion A of the idealized model sketched in Fig. 10). The impact ofabsorption in all other regions (such as region B in Fig. 10) can inprinciple be captured by leakage schemes, which can be understoodby the following considerations: Defining emission and absorptionrates of 𝜈 𝑒 and ¯ 𝜈 𝑒 , respectively, as 𝑅 em 𝜈 𝑒 ≡ 𝜆 𝑒 − 𝑌 𝑝 , (22a) 𝑅 abs 𝜈 𝑒 ≡ 𝜆 𝜈 𝑒 𝑌 𝑛 , (22b) 𝑅 em¯ 𝜈 𝑒 ≡ 𝜆 𝑒 + 𝑌 𝑛 , (22c) 𝑅 abs¯ 𝜈 𝑒 ≡ 𝜆 ¯ 𝜈 𝑒 𝑌 𝑝 , (22d)the rate of change of 𝑌 𝑒 can be written asd 𝑌 𝑒 d 𝑡 = −( 𝑅 em 𝜈 𝑒 − 𝑅 abs 𝜈 𝑒 ) + ( 𝑅 em¯ 𝜈 𝑒 − 𝑅 abs¯ 𝜈 𝑒 ) = − 𝑅 eff 𝜈 𝑒 + 𝑅 eff¯ 𝜈 𝑒 = − 𝑅 em 𝜈 𝑒 𝜒 𝜈 𝑒 + 𝑅 em¯ 𝜈 𝑒 𝜒 ¯ 𝜈 𝑒 , (23)where in the second line 𝑅 eff 𝜈 ≡ 𝑅 em 𝜈 − 𝑅 abs 𝜈 and in the third line 𝑅 eff 𝜈 = 𝑅 em 𝜈 𝜒 𝜈 are used for 𝜈 = 𝜈 𝑒 , ¯ 𝜈 𝑒 . Leakage schemes nowcompute d 𝑌 𝑒 / d 𝑡 in the form provided by the third line of Eq. (23) MNRAS000
200 ms. The ’**’means the same as ’*’, but additionally the value of 𝑌 𝑒, ej , as it was not provided, had to be estimated based on other information found in that reference.Reference 𝑚 𝑀 BH 𝐴 BH 𝑌 𝑒 viscosity neutrino 𝑄 𝑛𝑝 and 𝑚 𝑒 𝑌 𝑒, ej [ 𝑀 (cid:12) ] [ 𝑀 (cid:12) ] treatment treatment included? [1]Fernández et al. (2020) 0.03 3 0.8 0.1 std. 𝛼 -vis. leak.+abs. yes 0.28Fernández et al. (2019) 0.03 3 0.8 0.1 std. 𝛼 -vis. leak. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 ) no 0.20Just et al. (2015a) 0.03 3 0.8 0.1 std. 𝛼 -vis. spectral M1 yes 0.27m01M3A8 (this work) 𝛼 -vis. spectral M1 yes 0.32m01M3A8-noQm-no 𝜈 (this work) 𝛼 -vis. no abs. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 ) no 0.24Fujibayashi et al. (2020a) 0.1 3 0.8 0.07-0.5 𝑙 t =const. vis. grey M1+leak. yes 0.310.1 3 0.8 0.07-0.5 𝑙 t =const. vis. no abs. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 ) yes 0.30m01M3A8-vis2 (this work) 𝑙 t =const. vis. spectral M1 yes 0.35m01M3A8-vis2-no 𝜈 (this work) 𝑙 t =const. vis. no abs. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 ) yes 0.34Siegel & Metzger (2018) 0.03 3 0.8 0.1 MHD leakage no 0.18Siegel et al. (2019) 0.016 3 0.8 0.5 MHD leakage no (cid:46) . ∗∗ Fernández et al. (2019) 0.03 0.8 0.1 MHD leak. ( 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 ) no 0.16Miller et al. (2019b) 0.12 2.58 0.69 0.1 MHD Boltzmann yes ∼ . − . ∗∗ Miller et al. (2019a) 0.02 3 0.8 0.5 MHD Boltzmann yes 0 . ∗ m01M3A8-mhd (this work) Given the relevance for heavy-element production, the availablenumber and degree of sophistication of models for neutrino-cooledBH accretion disks is growing quickly. In the following we brieflyput the results of our study into context by comparison with a varietyof models available in the literature. The basic model setup usedhere (equilibrium torus, 𝛼 -viscosity or single-loop initial magneticfield configuration) is similar to what was used in most previousstudies of viscous or MHD disks, and several features are in broadagreement with the literature. One difference to several existingmodels is that we initiated most of our simulations with 𝑌 𝑒 = . 𝑌 𝑒 , 𝑌 𝑒, ej , is overall shifted compared to models starting with a 𝑌 𝑒 = . 𝑌 𝑒 should be similar in our simulations compared tosimulations from the literature. Note that in this case, since 𝑚 ( 𝑣 ) diverges at small velocities, a minimumvelocity, 𝑣 , must be imposed, below which d 𝑚 / d 𝑣 = 𝑣 is fixed by thecondition that the average velocity of the power-law distribution is the sameas for the original distribution. Before considering individual studies we first clarify some impor-tant points regarding neutrino leakage schemes, which otherwisemay become sources of confusion. Classical leakage schemes, suchas introduced by Ruffert et al. (1996) and Rosswog & Liebendör-fer (2003), are often mentioned to be unable to describe neutrinoabsorption. This is however only true for regions with positive net absorption rates (i.e. absorption minus emission rates, such as inregion A of the idealized model sketched in Fig. 10). The impact ofabsorption in all other regions (such as region B in Fig. 10) can inprinciple be captured by leakage schemes, which can be understoodby the following considerations: Defining emission and absorptionrates of 𝜈 𝑒 and ¯ 𝜈 𝑒 , respectively, as 𝑅 em 𝜈 𝑒 ≡ 𝜆 𝑒 − 𝑌 𝑝 , (22a) 𝑅 abs 𝜈 𝑒 ≡ 𝜆 𝜈 𝑒 𝑌 𝑛 , (22b) 𝑅 em¯ 𝜈 𝑒 ≡ 𝜆 𝑒 + 𝑌 𝑛 , (22c) 𝑅 abs¯ 𝜈 𝑒 ≡ 𝜆 ¯ 𝜈 𝑒 𝑌 𝑝 , (22d)the rate of change of 𝑌 𝑒 can be written asd 𝑌 𝑒 d 𝑡 = −( 𝑅 em 𝜈 𝑒 − 𝑅 abs 𝜈 𝑒 ) + ( 𝑅 em¯ 𝜈 𝑒 − 𝑅 abs¯ 𝜈 𝑒 ) = − 𝑅 eff 𝜈 𝑒 + 𝑅 eff¯ 𝜈 𝑒 = − 𝑅 em 𝜈 𝑒 𝜒 𝜈 𝑒 + 𝑅 em¯ 𝜈 𝑒 𝜒 ¯ 𝜈 𝑒 , (23)where in the second line 𝑅 eff 𝜈 ≡ 𝑅 em 𝜈 − 𝑅 abs 𝜈 and in the third line 𝑅 eff 𝜈 = 𝑅 em 𝜈 𝜒 𝜈 are used for 𝜈 = 𝜈 𝑒 , ¯ 𝜈 𝑒 . Leakage schemes nowcompute d 𝑌 𝑒 / d 𝑡 in the form provided by the third line of Eq. (23) MNRAS000 , 000–000 (0000) Just et al. and approximate the local quenching factors 𝜒 𝜈 based on (oftenproblem-dependent) conditions related to the neutrino optical depth.Since all right-hand sides of Eq. (23) are identical, the quenchingfactors are nothing but the difference between emission rate andabsorption rate normalized to the emission rate, i.e. 𝜒 𝜈 = ( 𝑅 em 𝜈 − 𝑅 abs 𝜈 )/ 𝑅 em 𝜈 , (24)which is a local version of the factor 𝜒 abs of Eq. (20) that is plottedin Figs. 2 and 3. Therefore, classical leakage schemes are alreadycapable of describing both the asymptotic 𝑌 𝑒 in neutrino-less opti-cally thin conditions, 𝑌 eq , em 𝑒 (since 𝜒 𝜈 = 𝑅 em 𝜈 𝑒 = 𝑅 em¯ 𝜈 𝑒 ,which is equal to Eq. (5)), and its absorption-modified version, 𝑌 eq 𝑒 (since the third line of Eq. (23) is equal to Eq. (3) up to errorsentering the computation of 𝜒 𝜈 ). From Eq. (23) it becomes clearthat 𝑌 eq 𝑒 depends sensitively on the assumptions entering the com-putation of the 𝜒 𝜈 factors and, in particular, on the resulting ratio 𝜒 𝜈 𝑒 / 𝜒 ¯ 𝜈 𝑒 . If, for instance, the same quenching factors are chosen forboth neutrino species, i.e. 𝜒 𝜈 𝑒 = 𝜒 ¯ 𝜈 𝑒 , then again 𝑅 em 𝜈 𝑒 = 𝑅 em¯ 𝜈 𝑒 andtherefore 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 , which is equivalent to disregarding neu-trino absorptions entirely concerning the equilibrium value of 𝑌 𝑒 .To our knowledge 𝑌 eq 𝑒 has not been discussed so far in the contextof leakage schemes, which is why the accuracy of 𝑌 eq 𝑒 in mod-els of neutrino-cooled disks using leakage schemes is unknown atthis point. The few comparisons of leakage schemes with transportmethods that exist so far (Richers et al. 2015; Foucart et al. 2016;Ardevol-Pulpillo et al. 2019; Gizzi et al. 2019) – and performedin most cases only in the context of proto-neutron stars – reportrather modest agreement for the total luminosities and the ratio of 𝜈 𝑒 to ¯ 𝜈 𝑒 luminosities when using the methodology of Ruffert et al.1996 or Rosswog & Liebendörfer 2003 (in contrast to the goodperformance of more recent leakage variants by Perego et al. 2016;Ardevol-Pulpillo et al. 2019). However, in disks the agreement isunknown and might be both better or worse.Various extensions to classical leakage schemes have been sug-gested in order to approximate net neutrino absorption in opticallythin regions (Fernández & Metzger 2013; Perego et al. 2016; Siegel& Metzger 2018; Ardevol-Pulpillo et al. 2019). The basic strategyoften consists of irradiating (subdomains of) the fluid configurationwith a neutrino field of constant luminosity and mean energy, sim-ilar to what is done in light-bulb schemes in the context of CCSNe(e.g. Bethe & Wilson 1985; Janka & Mueller 1996). While someof these extensions might be more accurate than others, rarely anysuch scheme has been benchmarked against a Boltzmann solution,which at this point makes it basically impossible to assess whetheravailable simulations using such schemes typically under- or overes-timate 𝑌 𝑒 in region A of neutrino-cooled disks (cf. Fig. 10). Frankly,in such regions the accuracy of the M1 scheme (employed in ourstudy) is not well known either. The results by Miller et al. (2019a)(cf. next section for a discussion) suggest that 𝑌 𝑒 , and therefore thequantitative impact of absorption, may be underrated in leakage andM1 schemes.Finally, although this aspect is probably less relevant for semi-transparent disks, for the sake of completeness we mention thatclassical neutrino leakage schemes cannot properly describe dy-namic conditions of very high optical depth, in which neutrinosare trapped and become advected by the fluid. This shortcomingis tackled in implementations of Sekiguchi et al. (2012); Peregoet al. (2016); Ardevol-Pulpillo et al. (2019) by additionally solv-ing an advection equation for trapped neutrinos and assuming thesetrapped neutrinos to equilibrate with the fluid under conditions ofhigh optical depth. Using a leakage scheme in combination with neutrino irradiationfrom a manually constructed neutrino field, Fernández et al. (2020)investigated the impact of absorption based on post-processingmethods, namely by comparing the net contributions of each ofthe four 𝛽 -reactions (cf. Eq. (1)) along outflow trajectories. Lack-ing at this point detailed information about their equilibrium values 𝑌 eq , em 𝑒 or 𝑌 eq 𝑒 , we cannot directly compare the impact of absorptionin region B (cf. Fig. 10) of their simulations to that found for ourmodels. Nevertheless, they obtain remarkably similar ejecta 𝑌 𝑒 aswe do and draw the same conclusion regarding the impact of neu-trino absorption, namely that neutrino absorption becomes relevantfor disk masses above ∼ . 𝑀 (cid:12) . In Fernández et al. (2019) similarviscous models are presented, which in contrast to those of Fernán-dez et al. (2020) assume 𝜒 𝜈 𝑒 = 𝜒 ¯ 𝜈 𝑒 (and therefore 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 )and ignore terms containing 𝑚 𝑒 and 𝑄 𝑛𝑝 in the 𝛽 -rates. In agree-ment with our findings that both of these measures reduce 𝑌 𝑒 in theejecta, they obtain an average electron fraction in the ejecta, 𝑌 𝑒, ej ,that is reduced by 0.08 (cf. Table 5).Fujibayashi et al. (2020a), employing the 𝑙 t =const. scheme todescribe viscous angular momentum transport, also switch off neu-trino irradiation for one of their models using a combination of greyleakage and M1 transport (cf. Sekiguchi et al. 2012 for details ofthis neutrino scheme). They find a very small impact of absorption,which is fully compatible to the minor difference that we find be-tween our models m01M3A8-vis2 and m01M3A8-vis2-no 𝜈 . Thisagreement suggests that the small impact of absorption observedin Fujibayashi et al. (2020a) is to a lesser extent a consequence ofGR effects (which in that study were newly included in contrast toprevious simulations of viscous disks) but mainly a ramification ofthe 𝑙 t =const. scheme, which compared to the standard 𝛼 -viscosityscheme underrates the impact of neutrino absorption by decelerat-ing the evolution during the neutrino-dominated phase as discussedin Sects. 4.1.2 and 4.3.3. As a result of the extended neutrino-dominated phase, the torus 𝑌 𝑒 follows 𝑌 eq 𝑒 until higher values arereached, which shifts the final 𝑌 𝑒 pattern of the ejecta and corre-spondingly reduces the abundances of elements with 𝐴 > 𝑌 𝑒 distributions peaking between 𝑌 𝑒 = . − . ∼ . 𝑀 (cid:12) ). Our MHD modelsconfirm the basic tendency that more neutron-rich material can beejected compared to viscous models, owing mainly to the morevigorous turbulence during the neutrino-dominated phase and theshorter expansion timescales. However, the average 𝑌 𝑒 in our ejecta(cf. 𝑌 𝑒, ej in Table 3) is still considerably higher than theirs. A part ofthis discrepancy can be accounted to the initial 𝑌 𝑒 , which is 𝑌 𝑒 = . 𝑌 𝑒 in thosestudies can be attributed to the fact that they employ the simpli-fication of 𝑄 𝑛𝑝 = 𝑚 𝑒 = 𝑌 eq , em 𝑒 in the torus by about 0.05-0.1. Moreover, the enhancement of 𝑌 eq 𝑒 with respect to 𝑌 eq , em 𝑒 due to neutrino absorptions is not capturedin Fernández et al. (2019) and Christie et al. (2019), because, like MNRAS , 000–000 (0000) eutrino absorption in black-hole disks for the viscous models mentioned above, they use the same opticaldepth for 𝜈 𝑒 and ¯ 𝜈 𝑒 .At the time of this writing the most complete simulations ofneutrino-cooled disks – including GRMHD and a Monte-CarloBoltzmann-solver for neutrino transfer – have been presented byMiller et al. (2019b,a), while Miller et al. (2019a) particularly fo-cused on the role of neutrino absorption in raising 𝑌 𝑒 in the outflow.Using a similar initial torus as Siegel et al. (2019) with a mass of0 . 𝑀 (cid:12) , Miller et al. (2019a) find ejecta with significantly reducedneutron densities, which is in line with the notion suggested byour results, namely that neutrino absorption as well as a correcttreatment of the terms including 𝑄 𝑛𝑝 and 𝑚 𝑒 in the 𝛽 -rates drive 𝑌 𝑒 to considerably larger values. The absence of any lanthanidesor heavier elements in the ejecta reported in Miller et al. (2019a)points to an even more significant impact of neutrino absorptionthan we find in the present study. It cannot be ruled out that our M1approximation for the neutrino transport, or neglecting GR metriceffects, or something else could attenuate the impact of absorptionin our models. Nevertheless, the differences between Miller et al.(2019a) and our models might not be so dramatic after all, becauseMiller et al. (2019a) only follow ejected material for evolution timesbetween 𝑡 ∼
70 and 150 ms. The fact that the average torus 𝑌 𝑒 isstill as low as 𝑌 𝑒 ∼ .
25 at the end of their simulation (cf. theirFig. 1) leaves open the possibility that neutron-rich material with 𝑌 𝑒 (cid:46) .
25 may still be expelled at later times. Moreover, in contrastto our study the impact of absorption is not tested by conducting twoseparate simulations with and without neutrino absorption. Instead,Miller et al. (2019a) estimate the impact of absorption by compar-ing the 𝑌 𝑒 pattern of the full transport simulation with a secondpattern that is obtained by integrating d 𝑌 𝑒 / d 𝑡 along the same fluidtrajectories using the same neutrino emission rates but neglectingthe absorption rates (assuming that we correctly understood theirprocedure). By doing so, however, the emission rates – which them-selves depend on 𝑌 𝑒 – become inconsistent with 𝑌 𝑒 , most likelyleading to an overestimated difference between the two cases withand without neutrino absorption. An important question is whether neutrino-cooled disks, either asremnants of NS-NS or NS-BH mergers or as collapsar engines,could represent major sources of the heaviest r-process elements,namely lanthanides and 3rd-peak nuclei. A first, interesting ob-servation in this regard is that for otherwise similar global modelparameters (i.e. 𝑀 BH , 𝐴 BH , 𝑚 ) a disk formed after a compact-object merger ( 𝑌 𝑒 ∼ .
1) would produce heavy r-process elementsmore efficiently than a collapsar disk ( 𝑌 𝑒 ∼ . 𝑌 𝑒 . While this difference, attributed entirely tothe initial 𝑌 𝑒 , needs to be kept in mind when comparing the entireelement production in disk ejecta, in the following we are only inter-ested in the ejecta fraction with conditions shaped self-consistentlyduring the disk evolution.The abundance distributions plotted in Fig. 13 suggest thatneutrino-cooled disks can be prolific sources of heavy (i.e. 𝐴 > [ 𝑋 / 𝑋 ] (cid:12) (cid:62)
1; see Table 4and Fig. 11). Basically all of the viscous models with neutrino ab-sorption, 𝑌 𝑒 = .
5, and without the 𝑄 𝑛𝑝 = 𝑚 𝑒 = 𝑌 𝑒 .The situation seems to look more favorable in the case of MHDdisks, because lower values of 𝑌 𝑒 can be reached by a larger fractionof outflow material than in viscous disks. Indeed, we find highervalues of [ 𝑋 / 𝑋 ] (cid:12) than for the corresponding viscous modelsand an abundance pattern that is only mildly sub-solar in the 𝐴 >
130 regime, mild enough for the gap to be explained by other, forinstance nuclear uncertainties. Hence, the results indicate that MHDdisks could indeed produce 3rd-peak material prolifically enoughto represent a major source, may it be in collapsars or NS mergers.However, we refrain from drawing strong conclusions at this point,on the one hand because of the poorly constrained uncertaintiesconnected to the initial conditions and the initial transient. On theother hand, at this point we have not investigated the dependenceof the thermodynamic conditions in the disk – and therefore of 𝑌 eq 𝑒 and of the ejecta 𝑌 𝑒 – on the numerical resolution, simplybecause of computational limitations. Another reason to be carefulnot to overinterpret our results is the apparent tension between ourrather optimistic results with the pessimistic results by Miller et al.(2019a), which may have its origin in our approximate M1 neutrinotreatment or the crude approximation of GR effects.For the case of a collapsar disk an additional major uncertainty,apart from those related to neutrino transport and MHD turbulence,is connected to the formation and resulting structure of the collapsardisk. In contrast to the isolated, finite-mass disks studied in thiswork (as well as in Siegel et al. 2019 and Miller et al. 2019a), themass reservoir in collapsar disks is continuously fed from infallingand circularizing stellar material, such that the evolution close tothe BH is quasi-stationary and mostly a function of the infallingmass flow rate, (cid:164) 𝑀 infall , which is presumably comparable to theBH accretion rate, i.e. (cid:164) 𝑀 infall ≈ (cid:164) 𝑀 BH . In the most optimistic caseneutrino absorption could, for whatever reason, be barely relevantand the mass ejection processes so swift that also neutrino emissionis unable to raise the electron fraction while ejecta travel from deepinside the torus to large radii. Assuming such an idealized situation,and supposing that our tori can be interpreted as collapsar disks,one can estimate the lower limit of the ejecta 𝑌 𝑒 by considering thevalues 𝑌 eq , em 𝑒 corresponding to emission equilibrium (cf. Eq. (5))attained in the torus for given mass accretion rates onto the centralBH. To this end we plot in Fig. 19 for all models of this study 𝑌 eq , em 𝑒 (as well as 𝑌 eq 𝑒 ) averaged over regions where weak interactions areefficient (i.e. where 𝜏 em < 𝑌 eq , em 𝑒 ( (cid:164) 𝑀 BH ) between different models provides an idea aboutthe level of uncertainty that is connected to the unknown structureof collapsar disks. If we assume that our set of models is somewhatrepresentative (which is by no means guaranteed) to bracket thespectrum of possible 𝑌 eq , em 𝑒 ( (cid:164) 𝑀 infall ) curves realized in collapsars,then we might infer that the infall rate needs to be at least as highas (cid:164) 𝑀 infall (cid:38) × − 𝑀 (cid:12) s − in order for 𝑌 𝑒 to attain a sufficientlylow value that can enable strong r-processing, 𝑌 𝑒 ≈ .
25. For lowermass-accretion rates, the r-process could operate with sufficientefficiency only up to the 2nd peak or, for (cid:164) 𝑀 infall (cid:46) − 𝑀 (cid:12) s − , MNRAS000
25. For lowermass-accretion rates, the r-process could operate with sufficientefficiency only up to the 2nd peak or, for (cid:164) 𝑀 infall (cid:46) − 𝑀 (cid:12) s − , MNRAS000 , 000–000 (0000) Just et al. M BH [ M /s]0.00.10.20.30.40.50.60.70.8 Y e strong r-processweak r-process equilibrium Y e vs mass accretion rate m01M3A8m1M3A8m001M3A8m01M5A8m01M10A8m01M3A4m01M3A9< Y eq,em e >( em < 1 s)< Y eq e >( em < 1 s) 10 M BH [ M /s]0.00.10.20.30.40.50.60.70.8 Y e strong r-processweak r-process equilibrium Y e vs mass accretion rate m01M3A8- 03m01M3A8- 1m01M3A8-vis2m01M3A8-mhdm01M3A8-noQmm01M3A8-wmm01_ye01 Figure 19.
Equilibrium electron fractions, averaged over emission-efficient regions where 𝜏 em < ∼ × − 𝑀 (cid:12) s − in order for the torus to have (cid:104) 𝑌 eq , em 𝑒 (cid:105) < .
25 and thus the possibility toenable a strong r-process in the ejecta. not at all. We stress again that this condition on 𝑌 eq 𝑒 in the torus is arather conservative limit concerning 𝑌 𝑒 in the ejecta. The latter willmost likely be driven to higher values along the outflow trajectorydue to the effects of neutrino emission and absorption.By means of Fig. 19, our models thus provide an additional,indirect constraint on the progenitor structure, which needs to besupplemented to the two (obvious) conditions that a BH must beformed and that material carries sufficient angular momentum tostay on circular orbits and form a disk. While the model E15 ofHeger et al. (2000) was estimated in Fig. 1 of Siegel et al. (2019)to fulfill these constraints, Aloy & Obergaulinger (2020) expectrather high BH masses of 𝑀 BH ≈ . .
3, and thereforeprobably too low mass accretion rates, for their models 35OC and350B (which were taken from Woosley & Heger 2006), respectively.Moreover, the sensitivity of 𝑌 eq 𝑒 with respect to the mass-accretionrate could imply a natural dependency of the produced r-processpattern on the progenitor structure. For instance, one might speculatethat progenitors with sufficient angular momentum to create disksearly during the collapse might tend to generate a stronger r-processthan low angular momentum progenitors, for which disks only formlate and with low accretion rates. Neutrino-cooled disks are likely to be important production sites ofheavy elements, but our understanding of the detailed compositionof the ejecta and its sensitivity with respect to the physics of theneutrino treatment is still incomplete. The aim of this study was toobtain a better quantitative understanding of the neutrino emissionand absorption effects during the long-term evolution of a disk. Tothis end, we investigated, for the first time systematically over a widerange of global parameters, the impact of including or not includingneutrino absorption in regulating the composition in the disk andthe ejecta. We varied the disk mass, BH mass and spin, details ofthe neutrino interaction rates, and the initial electron fraction of the disk. We furthermore tested the sensitivity to varying the viscous 𝛼 -parameter in the standard 𝛼 -viscosity formalism, to using an entirelydifferent viscosity prescription, and to replacing the axisymmetricviscous evolution by a 3D MHD evolution. In a post-processing step,we computed the nucleosynthesis yields of the ejecta and estimatedbasic properties of the expected kilonova.The main results of our study are:(i) We identify four characteristic regions for weak interactions in typ-ical neutrino-cooled disks (sketched in Fig. 10): At high mass accre-tion rates, (cid:164) 𝑀 BH (cid:38) − 𝑀 (cid:12) s − , neutrino absorption competes withneutrino emission in the bulk of the torus and it effectively raisesthe 𝑌 𝑒 equilibrium value from the pure-emission case, 𝑌 eq , em 𝑒 (cf.Eq. (5) and Fig. 1) to 𝑌 eq 𝑒 = 𝑌 eq , em 𝑒 + Δ 𝑌 𝑒 (cf. Eq. (3)), where typically Δ 𝑌 𝑒 ≈ . − .
1. In the surface layers close to the symmetry axisneutrino absorption dominates and drives 𝑌 𝑒 towards 𝑌 eq , abs 𝑒 ∼ . (cid:164) 𝑀 BH (cid:38) − 𝑀 (cid:12) s − neutrino absorption becomes irrelevant compared to emission and 𝑌 eq 𝑒 ≈ 𝑌 eq , em 𝑒 . Finally, below (cid:164) 𝑀 BH ∼ − 𝑀 (cid:12) s − also neutrinoemission ceases and 𝑌 𝑒 freezes out.(ii) Neutrino absorption also has an indirect leverage on the values of 𝑌 eq , em 𝑒 , because absorption attenuates the rates of neutrino cool-ing and therefore keeps the electron degeneracy lower than withoutabsorption (compare thick and thin red lines in Fig. 7). Moreover, 𝑌 eq , em 𝑒 is artificially reduced by ∼ . − . 𝑄 𝑛𝑝 , and electron mass in the neu-trino emission rates, 𝜆 𝑒 ± (cf. Eq. (1)), e.g. in models that employneutrino emission rates as formulated in Ruffert et al. (1996) orRosswog & Liebendörfer (2003) (see Table 5). In contrast, 𝑌 eq , em 𝑒 is almost insensitive to the inclusion of weak-magnetism effects (cf.Fig. 1 and Table 2).(iii) The electron fraction of the torus (approximately given by (cid:104) 𝑌 𝑒 (cid:105) mintor in Table 2) roughly anti-correlates with the torus optical depth toneutrinos, 𝜏 maxopt (cf. Fig. 11), for 𝜏 maxopt (cid:46) −
10, because a higheroptical depth tends to be reached by more compact tori exhibit-ing a higher level of electron degeneracy. Increasing the BH mass,
MNRAS , 000–000 (0000) eutrino absorption in black-hole disks 𝑀 BH , while keeping the initial disk size proportional to 𝑀 BH , low-ers the optical depth and neutron density of the torus and of theejecta, whereas varying the BH spin only has a weak impact. Forlarger optical depths, 𝜏 maxopt (cid:38) −
10, the anti-correlation betweenoptical depth and 𝑌 𝑒 saturates and tends to be reversed due to thecounteracting effects of neutrino absorption. As a result of this non-monotonic behavior, we observe the highest efficiency of heavyelement production for the model with intermediate torus mass of0 . 𝑀 (cid:12) within the sequence of models varying the torus mass (see,e.g., Fig. 14).(iv) Assuming a constant length 𝑙 t (parametrizing the scale of turnovermotions) instead of the conventional prescription for the 𝛼 -viscosityeffectively slows down angular momentum transport at late timeswhen the torus has expanded to radii 𝑟 (cid:29) 𝑙 t . This late-time decel-eration extends the phase of efficient neutrino emission and allows 𝑌 𝑒 to trace its secularly increasing equilibrium value until it hasreached higher values. As a consequence, the ejecta are less neutronrich and their 𝑌 𝑒 is less sensitive to the effects of neutrino absorp-tion (as in Fujibayashi et al. 2020a). The fact that the two viscosityprescriptions lead to qualitatively different torus dynamics and nu-cleosynthesis results highlights the importance of a self-consistentMHD description of angular momentum transport.(v) Replacing the viscous 2D treatment by a 3D MHD description, weobserve a broader 𝑌 𝑒 –mass distribution reaching down to lower 𝑌 𝑒 values. This trend, which is in agreement with Siegel & Metzger(2018); Fernández et al. (2019) and Miller et al. (2019b), is con-nected to the circumstance that the flow pattern of the MHD modelsis violently turbulent already during the neutrino-dominated phase,whereas the flow pattern of viscous models remains rather laminarduring that phase. Our results are overall less optimistic for pro-ducing low 𝑌 𝑒 material than those of Siegel & Metzger (2018) andFernández et al. (2019) (who employ a more approximate treat-ment of neutrinos as well as 𝑄 𝑛𝑝 = 𝑚 𝑒 =
0) but more optimisticthan those of Miller et al. (2019a) (who apply a more sophisticatedneutrino solver than we do but only evolve until about 150 ms).Compared to the models using the standard 𝛼 -viscosity we observea relatively weak impact of neutrino absorption, but as a reason wesuspect insufficient grid resolution to capture the MRI in the modelignoring neutrino absorption owing to a geometrically thinner diskcompared to the case including neutrino absorption.(vi) For our set of models with disk masses up to 0 . 𝑀 (cid:12) , neutrino ab-sorption results in a rise of the average ejecta 𝑌 𝑒 by about 0 . − . −
30. However, despite this reduction the abun-dance pattern in the
𝐴 >
130 mass range is still rather close tothe solar pattern for sufficiently compact and degenerate disks (e.g.as in our fiducial model with 𝑚 = . 𝑀 (cid:12) , 𝑀 BH = 𝑀 (cid:12) ,and 𝐴 BH = . 𝛼 -viscosity. Thus, our results support, if only marginally, the possi-bility of neutrino-cooled disks in NS mergers or collapsars beingmajor sources of 𝐴 >
130 elements.(vii) Since 𝑌 eq , em 𝑒 < 𝑌 eq 𝑒 ∼ 𝑌 𝑒 one can obtain an estimate of the lowestpossible 𝑌 𝑒 that can be reached in ejecta of collapsar disks forgiven mass accretion rates by mapping the equilibrium electronfractions of all our models to the corresponding BH accretion rates(cf. Fig. 19). If our models were representative for conditions incollapsars, the resulting 𝑌 eq , em 𝑒 ( (cid:164) 𝑀 infall ) relations would imply that 𝑌 𝑒 < .
25 can only be reached for infall mass fluxes of (cid:164) 𝑀 infall > × − 𝑀 (cid:12) s − and, correspondingly, in progenitors providingsufficient angular momentum for a disk to be formed at such highmass infall rates.(viii) Using an approximate kilonova model, we find variations of 40 −
80 % in the peak times and luminosities when comparing modelswith neutrino transport to models without, and variations of similarorder of magnitude when changing the prescription for angularmomentum transport. We observe a considerable sensitivity alsoto the detailed ejecta structure, 𝑚 ( 𝑣 ) , as well as when averagingthe composition and heating rates over the entire ejecta in velocityspace. These modeling uncertainties must be taken into accountwhen decyphering future kilonova observations.The results of our study demonstrate that the nucleosynthesisin outflows from neutrino-cooled disks is by far not universal andmay delicately depend both on the astrophysical parameters as wellas on the modeling assumptions. We stress again that the impact ofabsorption is even more significant for heavier disks that typicallyoccur in NS mergers than found here for the 0 . 𝑀 (cid:12) disk models.Considering that in most cases the ejected lanthanide fractions dif-fer by a factor of a few and more between models with and withoutneutrino absorption, we conclude that reliable nucleosynthesis andkilonova predictions require continued efforts to improve the quan-titative understanding of weak interactions in neutrino-cooled disks.Having said that, a profound knowledge of the MHD processes thatdetermine the thermodynamic state of the fluid, and by that regulatethe equilibrium values 𝑌 eq , em 𝑒 and 𝑌 eq 𝑒 in the disk, is every bit asimportant in order to faithfully predict 𝑌 𝑒 in the ejecta.Our study, although comprehensive in a great variety of mod-eling aspects, still contains a number of limitations. Due to itsapproximate nature, the M1 method used here might over- or under-estimate the impact of absorption, and the omission of GR effects(apart from the ones captured by the Artemova-potential) representsan additional source of uncertainty. Moreover, future investigationswill have to overcome the shortcomings connected to the manuallyconstructed initial configuration of the gas and, even more impor-tantly, of the magnetic field (see Christie et al. 2019, for a first suchstudy not accounting for neutrino absorption). Last but not least,resolution studies will have to elaborate the conditions on resolu-tion and/or grid configuration required in order to obtain convergedequilibrium values 𝑌 eq , em 𝑒 and 𝑌 eq 𝑒 in neutrino-MHD disks. ACKNOWLEDGMENTS
We thank Miguel Aloy, Hirotaka Ito, Martin Obergaulinger, andGabriel Martínez-Pinedo for stimulating discussions. OJ was sup-ported by the Special Postdoctoral Researchers (SPDR) programRIKEN. OJ and AB acknowledge support by the European Re-search Council (ERC) under the European Union’s Horizon 2020research and innovation programme under grant agreement No.759253. SG is an FRS-FNRS research associate and his workwas supported by the Fonds de la Recherche Scientifique (FNRS)and the Fonds Wetenschappelijk Onderzoek-Vlaanderen (FWO) un-der the EOS Project No O022818F. HTJ acknowledges fundingby the Deutsche Forschungsgemeinschaft (DFG, German ResearchFoundation) through grants SFB-1258 “Neutrinos and Dark Mat-ter in Astro- and Particle Physics (NDM)” and under Germany’sExcellence Strategy through Excellence Cluster ORIGINS (EXC2094)—390783311. SN is partially supported by JSPS Grants-in-Aid for Scientific Research KAKENHI (A) 19H00693”, Pio-neering Program of RIKEN for Evolution of Matter in the Uni-verse (r-EMU), and Interdisciplinary Theoretical and MathematicalSciences Program (iTHEMS) of RIKEN. AB acknowleges sup-port by Deutsche Forschungsgemeinschaft (DFG, German ResearchFoundation) - Project-ID 279384907 - SFB 1245 and - Project-ID138713538 - SFB 881 (“The Milky Way System”, subproject A10).
MNRAS000
MNRAS000 , 000–000 (0000) Just et al.
We acknowledge computational resources by the HOKUSAI super-computer at RIKEN and by the Max Planck Computing and Data Fa-cility (MPCDF). Nucleosynthesis calculations benefited from com-putational resources made available on the Tier-1 supercomputerof the Fédération Wallonie-Bruxelles infrastructure funded by theWalloon Region under the grant agreement no. 1117545.
APPENDIX A: EQUATIONS FOR THE HYDRODYNAMICMODELS
In this section we briefly summarize the equations that are solvedin our hydrodynamic simulations. For the viscous models we solvethe Newtonian Navier-Stokes equations in axisymmetry: 𝜕 𝑡 𝜌 + ∇ 𝑗 ( 𝜌𝑣 𝑗 ) = , (A1a) 𝜕 𝑡 ( 𝜌𝑌 𝑒 ) + ∇ 𝑗 ( 𝜌𝑌 𝑒 𝑣 𝑗 ) = 𝑄 N , (A1b) 𝜕 𝑡 ( 𝜌𝑣 𝑖 ) + ∇ 𝑗 ( 𝜌𝑣 𝑖 𝑣 𝑗 + 𝑃 g − 𝑇 𝑖 𝑗 vis ) = − 𝜌 ∇ 𝑖 Φ + 𝑄 𝑖 M , (A1c) 𝜕 𝑡 𝑒 t + ∇ 𝑗 ( 𝑣 𝑗 𝑒 t + 𝑣 𝑗 𝑃 g − 𝑣 𝑖 𝑇 𝑖 𝑗 vis ) = − 𝜌𝑣 𝑗 ∇ 𝑗 Φ + 𝑄 E + 𝑣 𝑗 𝑄 𝑗 M , (A1d)where 𝜌, 𝑣 𝑖 , 𝑌 𝑒 , 𝑃 g , 𝑒 t , Φ , 𝑇 𝑖 𝑗 vis are the baryonic mass density, veloc-ity, electron fraction, gas pressure, total (i.e. internal plus kinetic)energy density, gravitational potential, and viscous stress tensor, re-spectively, while 𝑄 N , 𝑄 𝑖 M , 𝑄 E stand for the source terms related tothe exchange of lepton number, momentum, and energy between thegas and neutrinos, respectively (cf. Eqs. (A6)). The viscous stresstensor is generally given by 𝑇 𝑖 𝑗 vis = 𝜂 vis (∇ 𝑖 𝑣 𝑗 + ∇ 𝑗 𝑣 𝑖 − 𝛿 𝑖 𝑗 ∇ 𝑘 𝑣 𝑘 ) , (A2)where we take into account only the 𝑇 𝑟 𝜙 vis and 𝑇 𝜃 𝜙 vis components.For the MHD models we solve in all three spherical polarcoordinates the special relativistic MHD counterpart of Eqs. (A1),namely 𝜕 𝑡 𝐷 + ∇ 𝑗 ( 𝐷𝑣 𝑗 ) = , (A3a) 𝜕 𝑡 ( 𝐷𝑌 𝑒 ) + ∇ 𝑗 ( 𝐷𝑌 𝑒 𝑣 𝑗 ) = 𝑄 N , (A3b) 𝜕 𝑡 ( 𝑆 𝑖 ) + ∇ 𝑗 ( 𝑆 𝑖 𝑣 𝑗 + 𝑃 ∗ g − 𝑏 𝑖 𝐵 𝑗 / 𝑊 ) = − 𝐷 ∇ 𝑖 Φ + 𝑄 𝑖 M , (A3c) 𝜕 𝑡 𝜏 + ∇ 𝑗 ( 𝑣 𝑗 𝜏 + 𝑣 𝑗 𝑃 ∗ g − 𝑏 𝐵 𝑗 / 𝑊 ) = − 𝐷𝑣 𝑗 ∇ 𝑗 Φ + 𝑄 E + 𝑣 𝑗 𝑄 𝑗 M , (A3d) 𝜕 𝑡 𝐵 𝑖 + ∇ 𝑗 ( 𝑣 𝑗 𝐵 𝑖 − 𝑣 𝑖 𝐵 𝑗 ) = , (A3e) ∇ 𝑖 𝐵 𝑖 = , (A3f)in which 𝑊 = ( − 𝑣 𝑗 𝑣 𝑗 / 𝑐 ) − / is the Lorentz factor, 𝑃 ∗ g = 𝑃 g + 𝑏 / 𝐵 𝑖 ( 𝑏 𝑖 ) the three-vector (four-vector) of the magnetic field in the laboratory (comoving) frame.The relationship between conserved and primitive variables is givenby: 𝐷 = 𝜌𝑊 , (A4a) 𝑆 𝑖 = 𝜌ℎ ∗ 𝑊 𝑣 𝑖 − 𝑏 𝑏 𝑖 / 𝑐 , (A4b) 𝜏 = 𝜌ℎ ∗ 𝑊 𝑐 − 𝑃 ∗ g − 𝑏 𝑏 − 𝜌𝑊𝑐 , (A4c) 𝑏 = 𝑊 ( 𝑣 𝑗 𝐵 𝑗 )/ 𝑐 , (A4d) 𝑏 𝑖 = 𝑊𝑣 𝑖 ( 𝑣 𝑗 𝐵 𝑗 )/ 𝑐 + 𝐵 𝑖 / 𝑊 , (A4e) with ℎ ∗ = ( 𝜌𝑐 + 𝑒 i + 𝑃 g + 𝑏 )/( 𝜌𝑐 ) being the total specific enthalpyand 𝑒 i the gas internal energy density.The above equations are solved along with the two-momentequations of neutrino transport, which read: 𝜕 𝑡 𝐸 + ∇ 𝑗 ( 𝐹 𝑗 + 𝑣 𝑗 𝐸 )+ 𝑃 𝑖 𝑗 ∇ 𝑖 𝑣 𝑗 − 𝜕 𝜖 ( 𝜖 𝑃 𝑖 𝑗 ∇ 𝑖 𝑣 𝑗 ) = 𝑆 ( ) , (A5a) 𝜕 𝑡 𝐹 𝑖 + ∇ 𝑗 ( 𝑐 𝑃 𝑖 𝑗 + 𝑣 𝑗 𝐹 𝑖 )+ 𝐹 𝑗 ∇ 𝑗 𝑣 𝑖 − 𝜕 𝜖 ( 𝜖𝑄 𝑖 𝑗𝑘 ∇ 𝑗 𝑣 𝑘 ) = 𝑆 ( ) ,𝑖 , (A5b)where 𝐸, 𝐹 𝑖 , 𝑃 𝑖 𝑗 , 𝑄 𝑖 𝑗𝑘 are the angular moments of rank 0 , , , 𝜖 . The higher moments, 𝑃 𝑖 𝑗 , 𝑄 𝑖 𝑗𝑘 , are estimated by expressing them as local functions of 𝐸, 𝐹 𝑖 , in the form suggested by Minerbo (1978) with the correspond-ing 3rd-order moment given in Just et al. (2015b). Since Eqs. (A5)are Newtonian and we want to avoid unphysical effects related tothe frame-dependent terms (i.e. terms containing the fluid velocities 𝑣 𝑗 ), we limit the fluid velocities entering Eqs. (A5) to not exceed0.1c by absolute value and neglect the azimuthal components. Thecoupling between hydrodynamics and neutrino transport is accom-plished by means of the source terms, which are related by: 𝑄 N = − 𝑚 B ∫ ∞ ( 𝑆 ( ) 𝜈 𝑒 − 𝑆 ( ) ¯ 𝜈 𝑒 ) d 𝜖𝜖 , (A6a) 𝑄 𝑖 M = − 𝑐 ∑︁ 𝜈 ∫ ∞ 𝑆 ( ) ,𝑖𝜈 d 𝜖 , (A6b) 𝑄 E = − ∑︁ 𝜈 ∫ ∞ 𝑆 ( ) 𝜈 d 𝜖 , (A6c)where 𝑚 B = . × − g.Both the (magneto-/viscous-)hydrodynamics equations as wellas the transport equations are evolved using PPM reconstruction(Colella & Woodward 1984) in the formulation of Mignone (2014),the HLLE Riemann solver, and 2nd-order Runge-Kutta time inte-gration. In order to satisfy the divergence constraint, Eq. (A3f), atall times, the magnetic fields are defined on a staggered grid and theinduction equation, Eq. (A3e), is evolved using the UCT scheme asin Del Zanna et al. (2007). For the recovery of the primitive vari-ables from the conserved variables in the MHD models we employa hybrid scheme, which for the first few iterations tries to find asolution with the 3D Newton-Raphson method as in Cerdá-Duránet al. (2008) and, if not successful, reverts to the more robust 1Dbracketing scheme recently presented by Kastaun et al. (2021) incombination with the 1D root-finding procedure of Chandrupatla(1997). APPENDIX B: METHOD OF LIGHT CURVECOMPUTATION
In Sect. 3.3 we already outlined the basic assumptions and theprocedure for obtaining the distribution of mass, 𝑚 ( 𝑣 ) , specific nu-clear heating rate, 𝑞 ( 𝑣 ) , and lanthanide plus actinide mass fraction, 𝑋 LA ( 𝑣 ) , along the homologous velocity coordinate 𝑣 . In this sec-tion we explain the procedure to kilonova light curve. We start withthe same assumptions as were invoked in Pinto & Eastman (2000);Wollaeger et al. (2018); Rosswog et al. (2018) for analytic schemesto construct the evolution equation for the radiation energy density, 𝐸 . However, we relax the Eddington approximation, 𝑃 → 𝐸 / 𝑃 / 𝐸 as function of 𝐸 andenergy-flux density, 𝐹 , using the closure relation of Minerbo (1978). MNRAS , 000–000 (0000) eutrino absorption in black-hole disks time [d]10 L [ e r g / s ] toy model lightcurves Wollaeger 2018, SuperNuWollaeger 2018, analyticour scheme= 10cm /g= 100cm /g Figure B1.
Comparison of light curves obtained using our two-momentscheme, Eqs. (B3), with reference solutions from Wollaeger et al. (2018)(data taken from their Fig. 5) for a toy model of expanding ejecta with aconstant opacity of 𝜅 =
10 cm g − and 100 cm g − . Black lines refer toour scheme, while red lines and blue lines refer to the Monte-Carlo and theanalytic scheme of Wollaeger et al. (2018), respectively. Hence, we end up with a set of two-moment equations for
𝐸, 𝐹 thatis very similar to the neutrino-transport system in Eqs. (A5). In con-trast to the latter, however, we evolve the energy-integrated equationswith a velocity explicitly given by 𝑣 / 𝑐 = 𝑟 /( 𝑐𝑡 ) ≡ 𝑥 as a functionof a single coordinate, 𝑥 . This leads to:d 𝐸 d 𝑡 + 𝑐𝑡𝑥 𝜕𝜕𝑥 ( 𝑥 𝐹 ) + 𝐸𝑡 = 𝑞𝜌 , (B1a)1 𝑐 d 𝐹 d 𝑡 + 𝑡𝑥 𝜕𝜕𝑥 ( 𝑥 𝑃 ) + 𝐹𝑐𝑡 = − 𝜅𝜌𝐹 , (B1b)where d / d 𝑡 = 𝜕 / 𝜕𝑡 + 𝑣𝜕 / 𝜕𝑟 is the Lagrangian time derivative andall quantities only depend on 𝑥 and 𝑡 . After multiplying by 𝑡 andsubstituting 𝑡 by the dimensionless coordinate˜ 𝑡 = ln 𝑡𝑡 , (B2)where 𝑡 is a fiducial timescale, we end up with:d 𝐸 d˜ 𝑡 + 𝑐𝑥 𝜕𝜕𝑥 ( 𝑥 𝐹 ) = 𝑞𝜌𝑡 𝑒 ˜ 𝑡 − 𝐸 , (B3a)1 𝑐 d 𝐹 d˜ 𝑡 + 𝑥 𝜕𝜕𝑥 ( 𝑥 𝑃 ) = − 𝜅𝜌𝐹𝑡 𝑒 ˜ 𝑡 − 𝐹𝑐 , (B3b)The motivation for the above manipulation was to bring Eqs. (B3)into a form where the left-hand side is formally equivalent to thetwo-moment equations in the laboratory frame but with the radialcoordinate, 𝑟 , replaced by the velocity coordinate, 𝑥 . This allowsto employ the same numerical methods as used for the neutrinotransport and described in Just et al. (2015b). We verify the methodand corresponding numerical scheme by comparing the result fora simple toy model of expanding ejecta with a constant opacity( 𝜅 =
10 cm g − or 𝜅 =
100 cm g − ) with reference solutionsby Wollaeger et al. (2018); cf. Fig. B1 and see Wollaeger et al.(2018) for the detailed model specifications. The agreement withthe transfer scheme SuperNu of Wollaeger et al. (2018) is excellentand slightly better than for the analytic method used in Wollaegeret al. (2018).The final ingredient needed to compute light curves is the translation between ejecta properties and photon opacities, 𝜅 . Inreality, 𝜅 is a function of the detailed composition as well as thethermodynamic and ionization state (e.g. Kasen et al. 2013; Tanakaet al. 2020). In this study we employ a simplified parametrization of 𝜅 in terms of only 𝑋 LA considering that lanthanides, if present in theejecta, tend to dominate the opacity. Additionally, we approximatelytake into account the reduction of possible line-transitions, andtherefore of 𝜅 , due to electron recombination below temperatures 𝑇 (cid:46) 𝐾 . We employ the following functional dependence: 𝜅 ( 𝑋 LA , 𝑇 ) = max { 𝜅 , min { 𝜅 , . × ( 𝑋 LA / − ) 𝛼 ( 𝑋 LA ) }× min {( 𝑇 / ) , } (B4)where 𝜅 = . g − , 𝜅 = .
65 cm g − , and 𝛼 = . = . 𝑋 LA < − ( > − ). The temperature is computed assumingthat radiation is in thermal equilibrium with the gas in regions ofsignificant optical depth. The form of 𝜅 as in Eq. (B4) was chosenwith the motivation to reproduce a set of light curves by Kasenet al. (2017), in which the lanthanide content of the radiating ejectawas varied from 𝑋 LA = . 𝑋 LA = − . The ejecta massand bulk velocity for this reference case are 0 . 𝑀 (cid:12) and 0 . 𝑐 ,respectively; see Kasen et al. (2017) for the precise form of theanalytic mass distribution. With the aim to reproduce the light curvesby Kasen et al. (2017), we use for this test (and only for this test) aradioactive heating rate that is presumably similar, though possiblynot identical, to the one used in Kasen et al. (2017). It is takenfrom Lippuner & Roberts (2015) and given by 𝑞 rad [erg g − s − ] = ( . × 𝑡 − . + . × 𝑒 − 𝑡 day / . ) with 𝑡 day beingthe time in units of days. Thermalization is included as describedin Sect. 3.3. As can be seen in Figure B2, our radiation solver,Eqs. (B3), combined with the opacity prescription of Eq. (B4), canreasonably well reproduce the original light curves by Kasen et al.(2017) and their dependence on 𝑋 LA . For comparison, Fig. B2 alsoshows the bolometric light curve of GW170817 (Waxman et al.2018). REFERENCES
Abbott B. P., et al., 2017a, ApJ, 848, L12Abbott B. P., et al., 2017b, ApJ, 848, L13Abbott B. P., et al., 2018, Physical Review Letters, 121, 161101Aloy M.-Á., Obergaulinger M., 2020, arXiv e-prints, p. arXiv:2008.03779Arcones A., Martínez-Pinedo G., Roberts L. F., Woosley S. E., 2010, A&A,522, A25+Ardevol-Pulpillo R., Janka H.-T., Just O., Bauswein A., 2019, MNRAS, 485,4754Artemova I. V., Bjoernsson G., Novikov I. D., 1996, ApJ, 461, 565Balbus S. A., Hawley J. F., 1991, ApJ, 376, 214Banerjee S., Tanaka M., Kawaguchi K., Kato D., Gaigalas G., 2020, ApJ,901, 29Barnes J., Kasen D., Wu M.-R., Martínez-Pinedo G., 2016, ApJ, 829, 110Bauswein A., Ardevol Pulpillo R., Janka H.-T., Goriely S., 2014, ApJ, 795,L9Bauswein A., Just O., Janka H.-T., Stergioulas N., 2017, ApJ, 850, L34Beckwith K., Hawley J. F., Krolik J. H., 2008, ApJ, 678, 1180Beloborodov A. M., 2003, ApJ, 588, 931Bethe H. A., Wilson J. R., 1985, ApJ, 295, 14Blandford R. D., Znajek R. L., 1977, MNRAS, 179, 433Bruenn S. W., 1985, ApJS, 58, 771Buras R., Rampp M., Janka H.-T., Kifonidis K., 2006, A&A, 447, 1049Cerdá-Durán P., Font J. A., Antón L., Müller E., 2008, A&A, 492, 937Chandrupatla T. R., 1997, Advances in Engineering Software, 28, 145Chen W., Beloborodov A. M., 2007, ApJ, 657, 383Chornock R., et al., 2017, ApJ, 848, L19MNRAS , 000–000 (0000) Just et al.
Figure B2.
Comparison of bolometric luminosity obtained using our two-moment scheme, Eqs. (B3), with the simplified opacity prescription of Eq. (B4)with radiative transfer solutions including detailed atomic opacities (Kasen et al. 2017) for a power-law ejecta distribution bearing different lanthanide massfractions, 𝑋 LA . Left and right panels show the same data but use different scaling of the time axis in order to facilitate the comparison between the data sets.The blue dotted line shows the heating rate, 𝑞 , powering all light curves. Consistent with the reference solutions, our light curves exhibit earlier and brighterpeaks for decreasing values of 𝑋 LA . Dots with error bars show for comparison the bolometric luminosity observed for GW170817 (Waxman et al. 2018).Christie I. M., Lalakos A., Tchekhovskoy A., Fernández R., Foucart F.,Quataert E., Kasen D., 2019, MNRAS, 490, 4811Colella P., Woodward P. R., 1984, Journal of Computational Physics, 54,174De S., Siegel D., 2020, arXiv e-prints, p. arXiv:2011.07176De Villiers J.-P., Hawley J. F., Krolik J. H., 2003, ApJ, 599, 1238Deaton M. B., O’Connor E., Zhu Y. L., Bohn A., Jesse J., Foucart F., DuezM. D., McLaughlin G. C., 2018, Phys. Rev. D, 98, 103014Del Zanna L., Zanotti O., Bucciantini N., Londrillo P., 2007, A&A, 473, 11Di Matteo T., Perna R., Narayan R., 2002, ApJ, 579, 706Eichler D., Livio M., Piran T., Schramm D. N., 1989, Nature, 340, 126Fernández R., 2015, MNRAS, 452, 2071Fernández R., Metzger B. D., 2013, MNRAS, 435, 502Fernández R., Kasen D., Metzger B. D., Quataert E., 2015, MNRAS, 446,750Fernández R., Tchekhovskoy A., Quataert E., Foucart F., Kasen D., 2019,MNRAS, 482, 3373Fernández R., Foucart F., Lippuner J., 2020, MNRAS, 497, 3221Foucart F., et al., 2015, Phys. Rev. D, 91, 124021Foucart 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., ScheelM. 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, L121Fujibayashi S., Shibata M., Wanajo S., Kiuchi K., Kyutoku K., SekiguchiY., 2020a, Phys. Rev. D, 101, 083029Fujibayashi S., Wanajo S., Kiuchi K., Kyutoku K., Sekiguchi Y., Shibata M.,2020b, ApJ, 901, 122Gizzi D., O’Connor E., Rosswog S., Perego A., Cabezón R. M., Nativi L.,2019, MNRAS, 490, 4211Glas R., Just O., Janka H. T., Obergaulinger M., 2019, ApJ, 873, 45Goriely S., 2015, The European Physical Journal A, 51, 22Goriely S., Demetriou P., Janka H., Pearson J. M., Samyn M., 2005, NuclearPhysics A, 758, 587Goriely S., Samyn M., Pearson J. M., 2007, Phys. Rev. C, 75, 064312Goriely S., Hilaire S., Koning A. J., 2008, A&A, 487, 767Goriely S., Hilaire S., Koning A. J., Sin M., Capote R., 2009, Phys. Rev. C,79, 024612Goriely S., Chamel N., Pearson J. M., 2010, Phys. Rev. C, 82, 035804 Goriely S., Bauswein A., Janka H.-T., 2011, ApJ, 738, L32+Goriely S., Bauswein A., Just O., Pllumbi E., Janka H.-T., 2015, MNRAS,452, 3894Gottlieb O., Nakar E., Piran T., Hotokezaka K., 2018, MNRAS, 479, 588Grossman D., Korobkin O., Rosswog S., Piran T., 2014, MNRAS, 439, 757Heger A., Langer N., Woosley S. E., 2000, ApJ, 528, 368Horowitz C. J., 2002, Phys. Rev. D, 65, 043001Horowitz C. J., Li G., 1999, Phys. Rev. Lett., 82, 5198Hossein Nouri F., et al., 2018, Phys. Rev. D, 97, 083014Hotokezaka K., Nakar E., 2020, ApJ, 891, 152Hotokezaka K., Kiuchi K., Kyutoku K., Okawa H., Sekiguchi Y.-i., ShibataM., Taniguchi K., 2013, Phys. Rev. D, 87, 024001Igumenshchev I. V., Abramowicz M. A., 2000, ApJS, 130, 463Janiuk A., 2019, ApJ, 882, 163Janka H.-T., 2017, Neutrino Emission from Supernovae. Springer, Cham,p. 1575, doi:10.1007/978-3-319-21846-5_4Janka H. T., Mueller E., 1996, A&A, 306, 167Just O., Bauswein A., Pulpillo R. A., Goriely S., Janka H.-T., 2015a, MN-RAS, 448, 541Just O., Obergaulinger M., Janka H.-T., 2015b, MNRAS, 453, 3386Just O., Obergaulinger M., Janka H.-T., Bauswein A., Schwarz N., 2016,ApJ, 816, L30Just O., Bollig R., Janka H.-T., Obergaulinger M., Glas R., Nagataki S.,2018, MNRAS,Kasen D., Barnes J., 2019, ApJ, 876, 128Kasen D., Badnell N. R., Barnes J., 2013, ApJ, 774, 25Kasen D., Fernández R., Metzger B. D., 2015, MNRAS, 450, 1777Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Na-ture, 551, 80Kastaun W., Kalinani J. V., Ciolfi R., 2021, Phys. Rev. D, 103, 023018Kawaguchi K., Shibata M., Tanaka M., 2018, ApJ, 865, L21Kodama T., Takahashi K., 1975, Nuclear Physics A, 239, 489Kohri K., Mineshige S., 2002, ApJ, 577, 311Koning A. J., Hilaire S., Duijvestijn M. C., 2005, in Haight R. C., ChadwickM. B., Kawano T., Talou P., eds, American Institute of Physics Con-ference Series Vol. 769, International Conference on Nuclear Data forScience and Technology. pp 1154–1159, doi:10.1063/1.1945212Korobkin O., Rosswog S., Arcones A., Winteler C., 2012, MNRAS, 426,1940Korobkin O., et al., 2020, arXiv e-prints, p. arXiv:2004.00102Kulkarni S. R., 2005, ArXiv Astrophysics e-prints,MNRAS , 000–000 (0000) eutrino absorption in black-hole disks Lattimer J. M., Mackie F., Ravenhall D. G., Schramm D. N., 1977, ApJ,213, 225Lee W. H., Ramirez-Ruiz E., 2007, New Journal of Physics, 9, 17Lemaître J.-F., Goriely S., Hilaire S., Sida J.-L., 2019, Phys. Rev. C, 99,034612Li L.-X., Paczyński B., 1998, ApJ, 507, L59Lippuner J., Roberts L. F., 2015, ApJ, 815, 82Liu M.-Q., 2010, Research in Astronomy and Astrophysics, 11, 91MacFadyen A. I., Woosley S. E., 1999, ApJ, 524, 262Malkus A., Kneller J. P., McLaughlin G. C., Surman R., 2012, Phys. Rev. D,86, 085015Margalit B., Metzger B. D., 2017, ApJ, 850, L19Marketin T., Huther L., Martínez-Pinedo G., 2016, Phys. Rev. C, 93, 025805McKinney J. C., Tchekhovskoy A., Sadowski A., Narayan R., 2014, MN-RAS, 441, 3177Metzger B. D., 2019, Living Reviews in Relativity, 23, 1Metzger B. D., Piro A. L., Quataert E., 2008, MNRAS, 390, 781Metzger B. D., Piro A. L., Quataert E., 2009, MNRAS, 396, 304Metzger B. D., et al., 2010, MNRAS, 406, 2650Mignone A., 2014, Journal of Computational Physics, 270, 784Mihalas D., Mihalas B. W., 1984, Foundations of radiation hydrodynamics.New York, Oxford University PressMiller J. M., Sprouse T. M., Fryer C. L., Ryan B. R., DolenceJ. C., Mumpower M. R., Surman R., 2019a, arXiv e-prints, p.arXiv:1912.03378Miller J. M., et al., 2019b, Phys. Rev. D, 100, 023008Minerbo G. N., 1978, J. Quant. Spectrosc. Radiative Transfer, 20, 541Moffatt H. K., 1978, Magnetic field generation in electrically conductingfluids. Cambridge, England, Cambridge University Press, 1978. 353 p.Mooley K. P., et al., 2018, ApJ, 868, L11Nagataki S., Mizuta A., Sato K., 2006, ApJ, 647, 1255Nakamura K., Kajino T., Mathews G. J., Sato S., Harikae S., 2015, A&A,582, A34Nakar E., 2007, Phys. Rep., 442, 166Narayan R., Yi I., 1994, ApJ, 428, L13Obergaulinger M., 2008, Dissertation, Technische Universität München,MünchenPaschalidis V., Ruiz M., Shapiro S. L., 2015, ApJ, 806, L14Perego A., Rosswog S., Cabezón R. M., Korobkin O., Käppeli R., ArconesA., Liebendörfer M., 2014, MNRAS, 443, 3134Perego A., Cabezón R. M., Käppeli R., 2016, ApJS, 223, 22Perego A., Radice D., Bernuzzi S., 2017, ApJ, 850, L37Pinto P. A., Eastman R. G., 2000, The Astrophysical Journal, 530, 744Pllumbi E., Tamborra I., Wanajo S., Janka H.-T., Hüdepohl L., 2015, ApJ,808, 188Popham R., Woosley S. E., Fryer C., 1999, ApJ, 518, 356Pruet J., Thompson T. A., Hoffman R. D., 2004, ApJ, 606, 1006Qian Y., Woosley S. E., 1996, ApJ, 471, 331Radice D., Perego A., Hotokezaka K., Fromm S. A., Bernuzzi S., RobertsL. F., 2018, ApJ, 869, 130Rampp M., Janka H., 2002, A&A, 396, 361Rezzolla L., Giacomazzo B., Baiotti L., Granot J., Kouveliotou C., AloyM. A., 2011, ApJ, 732, L6+Rezzolla L., Most E. R., Weih L. R., 2018, ApJ, 852, L25Richers S., Kasen D., O’Connor E., Fernández R., Ott C. D., 2015, ApJ,813, 38Richers S. A., McLaughlin G. C., Kneller J. P., Vlasenko A., 2019, Phys.Rev. D, 99, 123014Roberts L. F., Kasen D., Lee W. H., Ramirez-Ruiz E., 2011, ApJ, 736, L21+Rosswog S., Liebendörfer M., 2003, MNRAS, 342, 673Rosswog S., Ramirez-Ruiz E., Davies M. B., 2003, MNRAS, 345, 1077Rosswog S., Sollerman J., Feindt U., Goobar A., Korobkin O., WollaegerR., Fremling C., Kasliwal M. M., 2018, A&A, 615, A132Ruffert M., Janka H. T., 1998, A&A, 338, 535Ruffert M., Janka H.-T., Schaefer G., 1996, A&A, 311, 532Ruffert M., Janka H.-T., Takahashi K., Schaefer G., 1997, A&A, 319, 122Sa¸dowski A., Narayan R., Tchekhovskoy A., Abarca D., Zhu Y., McKinneyJ. C., 2015, MNRAS, 447, 49 Sekiguchi Y., Kiuchi K., Kyutoku K., Shibata M., 2012, Progress of Theo-retical and Experimental Physics, 2012, 01A304Setiawan S., Ruffert M., Janka H.-T., 2004, MNRAS, 352, 753Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337Shibata M., Sekiguchi Y., Takahashi R., 2007, Progress of TheoreticalPhysics, 118, 257Siegel D. M., Metzger B. D., 2018, ApJ, 858, 52Siegel D. M., Barnes J., Metzger B. D., 2018, preprint,( arXiv:1810.00098 )Siegel D. M., Barnes J., Metzger B. D., 2019, Nature, 569, 241Skinner M. A., Burrows A., Dolence J. C., 2016, ApJ, 831, 81Steiner A. W., Hempel M., Fischer T., 2013, ApJ, 774, 17Sumiyoshi K., Fujibayashi S., Sekiguchi Y., Shibata M., 2021, ApJ, 907, 92Surman R., McLaughlin G. C., 2005, ApJ, 618, 397Tanaka M., Hotokezaka K., 2013, ApJ, 775, 113Tanaka M., Kato D., Gaigalas G., Kawaguchi K., 2020, MNRAS, 496, 1369Tanvir N. R., et al., 2017, ApJ, 848, L27Villar V. A., et al., 2017, The Astrophysical Journal, 851, L21Wanajo S., Sekiguchi Y., Nishimura N., Kiuchi K., Kyutoku K., Shibata M.,2014, ApJ, 789, L39Waxman E., Ofek E. O., Kushnir D., Gal-Yam A., 2018, MNRAS, 481, 3423Wollaeger R. T., et al., 2018, MNRAS, 478, 3298Woosley S. E., 1993, ApJ, 405, 273Woosley S. E., Heger A., 2006, ApJ, 637, 914Wu M.-R., Tamborra I., Just O., Janka H.-T., 2017, Phys. Rev. D, 96, 123015Xu Y., Goriely S., Jorissen A., Chen G. L., Arnould M., 2013, A&A, 549,A106Yakovlev D., Kaminker A., Gnedin O., Haensel P., 2001, Physics Reports,354, 1MNRAS000