Impact of neutrino pair-production rates in Core-Collapse Supernovae

In this paper, we present a careful study on the impact of neutrino pair-production on core-collapse supernovae via spherically-symmetric, general-relativistic simulations of two different massive star progenitors with energy-dependent neutrino transport. We explore the impact and consequences of both the underlying microphysics and the implementation in the radiation transport algorithms on the supernova evolution, neutrino signal properties, and the explosion dynamics. We consider the two dominant neutrino pair-production processes found in supernovae, electron-positron annihilation as well as nucleon-nucleon bremsstrahlung in combination with both a simplified and a complete treatment of the processes in the radiation transport algorithms. We find that the use of the simplified prescription quantitatively impacts the neutrino signal at the 10\% level and potentially the supernova dynamics, as we show for the case of a 9.6M$_\odot$ progenitor. We also show that the choice of nucleon-nucleon bremsstrahlung interaction can also have a quantitative impact on the neutrino signal. A self-consistent treatment with state-of-the-art microphysics is suggested for precision simulations of core collapse, however the simplified treatment explored here is both computationally less demanding and results in a qualitatively similar evolution.


I. INTRODUCTION
Core-collapse supernovae (CCSNe) represent the last stage of massive star evolution for stars more massive than 8M and, along with neutron star-neutron star mergers and Type-Ia supernovae, are one of the main channels of galactic nucleosynthesis [1,2]. Not only do CCSNe contribute to the production of heavy elements but they are also the main birth site of neutron stars and stellar mass black holes.
Supernovae are also true multimessenger events, producing neutrinos, gravitational waves, as well as photons. The most readily available observable is the electromagnetic signal, for example, the Zwicky Transient Facility observed over 800 CCSNe in 2018 [3]. In the fortunate case of a galactic supernova, the two other channels, gravitational waves and neutrinos, become possible [4,5]. At the onset of the explosion, the outside layers of the progenitor star shroud the core and prevent photons from carrying direct information from the core. Neutrinos and gravitational waves are the only direct channels helping us deciphering the physics of the early explosion. The supernova mechanism is thus still largely observationally unconstrained. Regardless, numerical simulations performed by different groups allow us to test and refine the theories we have [6][7][8][9][10][11][12].
The main theory is the neutrino-driven supernova mechanism. Once the fusion reactions in the core stop and gravity overcomes the electron degeneracy pressure, the collapse begins. At nuclear density the core stiffens and the collapse stops and rebounds outwards. The information about this bounce propagates through the in-falling mater, reaching supersonic velocities and creating a shock. This shock propagates out flowing against the ram pressure of the infalling layers and dissociating the nuclei accreting through. In doing so, the shock loses energy and ends up stalling. The neutrino-driven mechanism is the idea [13] that the neutrinos can re-energize the shock by transferring energy from the cooling protoneutron star (PNS) to the material behind the shock through absorption in the so-called gain layer. Studies have shown that this heating is very sensitive to the neutrino spectrum, which in turn is sensitive to the emission and absorption processes [12,[14][15][16][17]. While much progress has been made in 3D [9,[18][19][20][21][22][23][24][25][26][27][28][29], the theory is not yet completed. Progress on all fronts is needed to further constrain CCSN theories and the underlying physics.
Neutrino transport is one of the most difficult aspects of modeling supernova simulations. A completely selfconsistent treatment of neutrinos would involve solving the 6D Boltzmann equation over the course of the simulation, along with capturing all of the important interactions with the medium. This is too computationally expensive, especially with the required resolution, though see [30]. Therefore, we often resort to approximate schemes such as flux limited diffusion or truncated moment schemes. In this study, we use a moment scheme [31,32]. We essentially treat the neutrinos like a fluid, evolving the energy density and momentum density. Neutrinos are produced inside of the PNS, an optically thick medium, and diffuse out into semitransparent, optically-thin matter in the gain region behind the shock. This change in the qualitative nature of the environment makes the neutrino transport complex as neutrinos transition from being strongly coupled to the matter to a free-flowing behavior, and therefore neither assumption can be used globally to simplify the problem. The main type of neutrinos interacting in the gain region are the electron type neutrinos and antineutrinos through charged-current interactions. However, heavy lepton neutrinos (ν µ ,ν µ , ν τ , andν τ ), which mainly cool the PNS, also play a major role. Hence, their interactions with matter need to be treated as accurately as possible. The main production channel for heavy-lepton neutrinos is via pair-production, where a pair consisting of a neutrino and an antineutrino is formed. The dominate production processes for these pairs in CC-SNe include electron-positron annihilation and nucleonnucleon bremsstrahlung. Charged-current interactions (either emission or absorption) of single heavy-lepton neutrinos with muons or taus are suppressed due to those charged lepton's large mass, although see [33][34][35].
In this paper, different treatments for the thermal pairproduction processes are tested. These interactions are challenging to treat as they involve not just one neutrino, but two, which necessitates the coupling of the species and of the energy bins. This had often lead to approximations for their inclusion in neutrino transport algorithms [36]. As part of this paper, we assess one such approximation with the goal of reducing computational expense while maintaining the fidelity of the solution. Not only are the neutrino pair-production processes computationally complex, another problem inherent in the nucleon-nucleon bremsstrahlung interaction is its uncertain nuclear physics. For this reason, we consider two different ways of treating the interaction. First we consider the commonly used one pion exchange (OPE) formalism by Hannestad and Raffelt [37]. We also consider a recent T-matrix formalism formulated by Guo et al. [38] based on chiral effective field theory fitted to experimental phase shifts. We test these two different formalisms as well as a simplified version for the nucleon-nucleon bremsstrahlung based on [39]. For electron-positron annihilation we follow the formalism described in Bruenn et al. [40] as well as a simplified version [36]. We perform 1D simulations for each of the six combinations of different treatments for two different progenitors. We use a 20-M progenitor, a model studied across many CCSN codes in [6] and a 9.6-M progenitor, which has the property of exploding in 1D simulations. We explore the impact of all the different treatments on the early supernova evolution, the explosion parameters, and the neutrino luminosities and mean energies. Furthermore, we scrutinize the validity of our simplified approximation used in order to inform future multidimensional simulations on the impact of heavy-lepton neutrino pair-production treatments.
The outline of this paper is as follows. In § II, we overview the simulation code we use, GR1D, and we describe the different interactions involved in our study and their implementation in GR1D and NuLib. We also present the two progenitors and their history of use in CCSN simulations. In § III A, we describe the results on the 20M progenitor and § III B the ones of the 9.6M progenitor. We finally conclude in § IV.

A. GR1D
For all the simulations presented in this paper we use the general-relativistic radiation-hydrodynamic code GR1D [36,41]. For the neutrino transport, GR1D uses a moment scheme [32,42]. It evolves the 0th and 1st moment of the neutrino distribution function for multiple neutrino species and multiple neutrino energies. The neutrino-matter interaction terms (completely local) are solved implicitly while the non-local spatial fluxes are solved explicitly. The evolution is done in the coordinate (or laboratory) frame but full velocity dependence is included in the neutrino-matter interactions and to order v/c in the spatial transport terms. We present the model moment evolution equations here, highlighting the neutrino-matter interaction source terms and refer the reader to [36] for full details, and where E and F r are the zeroth and first moments of the species and energy-dependent neutrino distribution functions, P rr is the 2nd moment, and in the M1 approximation is taken as an analytic expression involving the first two moments. Here, and in the following, we suppress the energy and species dependence of these moments and source terms, unless needed. α and X are metric functions, ∂ [...] refers to the energy-space fluxes, and G t/r and C t/r are the geometric and neutrino-matter source terms, respectively. For the full expression for ∂ [...] and G t/r we refer the reader to [36], since in this paper we focus on the neutrino-matter interactions, we explicitly write C t/r here and describe each term below, In GR1D, neutrino-matter interactions fall into four categories. (i) [S α e/a ] Charged-current neutrino-matter interactions, where electron type neutrinos and antineutrinos are absorbed or emitted from the matter. (ii) [S α iso ] Elastic scattering interactions, where neutrinos of all types scatter on nucleons and nuclei. These scatters change the neutrino direction but maintain their energy. For the emission, absorption and the elastic scattering interactions, we treat the source terms in the following way: where η is the emissivity, κ a and κ s are the absorption and scattering opacities respectively, u α is the fluid fourvelocity, and J and H α are the zeroth and first neutrino moments in the fluid frame (see [36] for detailed expressions of J and H α in terms of E and F r and the closure relation).
(iii)[S α scatter ] Inelastic scattering interactions, where neutrinos scatter on electrons and appreciably change their energy and direction. This interaction necessitates a coupling of neutrino energy bins within a neutrino species. For inelastic neutrino-electron scattering, we use the source terms described in Shibata et al. [32]. In this study, we ignore inelastic scattering on nucleons.
Finally, (iv) [S α pair ], pair-production interactions where a neutrino-antineutrino pair is emitted. In GR1D, we only consider pair-production interactions involving heavy-lepton neutrinos (ν µ ,ν µ , ν τ , andν τ ) since the interactions involving electron type neutrino-antineutrino pairs are dwarfed by the charged-current rates for these neutrinos. With GR1D, there are two ways of including S α pair into the evolution equations. The first is a simplified method where we generate simplified emissivities (η νν eff ) and absorption coefficients (κ a,eff νν ) for each neutrino energy group and treat these terms like the emission and absorption interactions in (i) above. The precise form of these coefficients depends on the particular pair-production process and are described in the following section. This method is computational efficient as it does not require coupling neutrinos of different energies together when performing the implicit solution of the evolution equations and the interaction rates depend only on the temperature, electron fraction, density and neutrino energy. However, in general, these neutrino pairproduction processes do depend on the occupation density of more than one neutrino and therefore this method is an approximation. The second method is more complete, but also more computationally expensive. It uses kernels to describe the interaction between two neutrinos of different energies (and species) and takes into account the final state neutrino occupation (for emission) and initial state neutrino occupation (for annihilation) hence, coupling different energy groups. The source term for this method is based on [32] and follows from taking the appropriate angular moments of the full Boltzmann collision integral for neutrino-antineutrino annihilation, where ν is the neutrino energy, u α is the fluid fourvelocity, α is a unit vector perpendicular to u α , and B(ν, Ω) is, where for clarity we have suppressed the ν, and ν as well as Ω and Ω dependence in each of the occupation probabilities, f and f , respectively. µ, which is a function of both the prime and unprimed angular variables, is the cosine of the angle between the neutrino and antineutrino. As is typically done, we assume an angular expansion form of the production and annihilation kernels, R pro/ann ∼ R pro/ann 0 + µR pro/ann 1 , where R pro/ann 0/1 only depends on the energies of the two neutrinos involved and the underlying interaction (see the following section). Following [32], Eqs. 7 and 8 are reduce to a single integral over ν where the integrand depends only on the primed and unprimed, zeroth, first, and second neutrino moments and the R pro/ann 0/1 kernels, where h αβ = g αβ + u α u β is the projection operator and L αβ is the traceless L αβ , the second-moment tensor in the fluid frame (analogous to P αβ above, which is the coordinate frame second moment).

B. Implementation in NuLib
NuLib (http://www.nulib.org) is an open-source neutrino interaction library [36] that we use to produce tables of the neutrino-matter interaction coefficients for interpolation during our simulations. For this work, we utilized the interactions described in Table I, which are divided into the four main interaction types described above. In this work, we focused on the heavy-lepton pair-production processes and the accuracy of the prescriptions used in the transport for these interactions. For this reason, we describe these in detail below.
The two main neutrino pair-production processes in a CCSN environment are electron-positron pair annihilation and nucleon-nucleon bremsstrahlung. As discussed in § II A, we consider both an simplified prescription for these interactions and a kernel treatment. For the electron-positron pair annihilation the underlying interaction is the same in these two methods, described in [39,40]. We used NuLib to compute R ann/pro 0 and R ann/pro 1 for use in Eq. 9, which gives the neutrino pair annihilation and production rates as a function of the two neutrino energies, ν and ν , for a given value of the matter temperature and electron chemical potential. For the simplified version of neutrino emission from electron-positron annihilation (see [36] for more details), we compute η e − e + eff (ν) by assuming R pro/ann 1 = 0 (i.e. isotropic emission), no final state neutrino blocking, and integrating over all possible ν . We construct where BB is the black body intensity for heavy-lepton neutrinos with energy ν in a medium with temperature T . This ensures there is no net emission in regions where the neutrino field is the same as the equilibrium neutrino field and no absorption in regions where the neutrino field is negligible. This is an approximation.
The other main neutrino pair-production process of importance in CCNSe is nucleon-nucleon bremsstrahlung. Before this work, this interaction was included in NuLib only via an simplified way taken from Burrows et al. (2006) [39]. The simplified single neutrino emissivity (with units of erg where is the total energy emission rate for a pair of neutrinos, ζ is a correction factor (taken to be 0.5 [39]), x n/p is the mass fraction of neutron and protons, ρ 14 is the density scaled to 10 14 g cm −3 , and T is the matter temperature.
As is the case for electron-positron annihilation, we construct a simplified absorption by invoking Kirchhoff's law, κ N N eff (ν) = η N N eff (ν)/BB(ν, T ). This simplified emissivity was made in the non-degenerate-medium limit assuming an OPE potential. It is only dependant on the nucleon number densities and the temperature of the medium. In the early phases of a CCSN explosion, the nucleons at the densities of interest are rarely degenerate, however at latter stages, during the cooling of the PNS for example, the densities where the nucleon-nucleon bremsstrahlung rates can impact the evolution and emission may be in the degenerate regime, therefore this method may need to be reconsidered.
In this work, we extend NuLib to include kernels for the nucleon-nucleon bremsstrahlung process in addition to the electron-positron annihilation process. The nucleonnucleon bremsstrahlung kernels follow the form of, where n B is the baryon density, G F ∼ 1.166 × 10 −11 MeV −2 is the weak coupling constant, C a = g A /2 with g A ∼ −1.26 is the axial vector coupling constant, ω = ν + ν is the sum of the two neutrino energies, and S σ (ω) is the structure function. As for the electronpositron annihilation kernel, we decompose R pro into Legendre moments, R pro 0/1 . Given the dependence on µ in Eq. 12, this is a trivial decomposition and R pro 1 = −R pro 0 /3. In order to obtain R ann 0/1 in accordance with detailed balance, we use R pro 0/1 = e −ω/T R ann 0/1 . The exact definition of S σ (ω) depends on the underlying interaction and in this work we consider two different models. First, we include the classic nucleon-nucleon bremsstrahlung rates described in Hannestad and Raffelt (1998) [37]. Similar to the parametrization above, this interaction is derived from the OPE potential, but also includes in the structure function effects such as a nonvanishing pion mass, effects from multiple-scatterings, and is valid for both the degenerate and non-degenerate limits with an interpolation for semi-degenerate regions. The structure function is [37], This structure function is for an arbitrary nucleon interacting with a like nucleon with a nucleon density (n N ), temperature (T ), and the degeneracy factor, η = p 2 F /(2m N T ) (where p 2 F = (3π 2 n N ) 1/3 is the Fermi momentum of the nucleons with mass m N ). The spin-fluctuation rate (Γ), gives the strength of the bremsstrahlung. Also present in the structure function are dimensionless functions g(y, η) and s(ω/T, y) which are a function representing the multi-scattering effect and the interpolation of the nucleon structure function between degenerate and non-degenerate medium, respec-tively. For completeness, where α π and m π are the pion fine-structure constant and the pion mass, respectively. For detailed expressions for g(y, η) and s(ω/T, y), see [37] or the bremsstrahlung routines in NuLib (http://www.nulib.org). In NuLib, we compute a table of R pro/ann 0/1 as a function of an arbitrary nucleon density n N , the temperature T , for the pair of neutrino energies ν and ν . During our simulation we interpolate this table for three values of the nucleon density, n n , n p , √ n n n p , and combine the rates with weights of 1, 1, and 28/3, respectively [39]. The second nucleon-nucleon bremsstrahlung interaction we consider is the recent formalism from Guo and Martinez-Pinedo [38]. They calculate the structure function (S σ (ω)) used in Eq. 12 by using the T-matrix element based on the χEFT potential presented in Entem et al. [45]. A similar method was previously explored in Ref. [46]. A followup in [47] where the T-matrix formalism was shown to give modestly different results in supernova simulations from the OPE prescription above.
The T-matrix formalism used in [38] is an improvement over [46] with the inclusion of off-shell T-matrix elements in addition to on-shell elements. In NuLib, we utilize the table of S σ (ω) values provided by the authors. We interpolate this four dimensional table (ρ, T , Y e , and ν + ν ) for use in Eq. 12 in order to construct our tables.
We conclude this section by comparing each of the pair-production processes and prescriptions utilized in this work at different CCSN-like conditions. The results are shown in Fig.1, where we compare the single neutrino number isotropic emissivities, ignoring any final state neutrino blocking, as a function of energy at four densities. Following [46], we use the following relationship between density and temperature typically found in CCSN environments, and adopt an electron fraction of Y e = 0.2.
The emissivities themselves, as well as the difference between the emissivities, are strongly dependent on the density and temperature. The increase in the bremsstrahlung rates with increasing density is due to both the ρ 2 dependence and the roughly T 4.5 dependence of the number emission rate where the electronpositron annihilation number emission rates increase only to due to the increase in the temperature, scaling roughly as T 8 . Therefore we expect the importance of bremsstrahlung over electron-positron pair annihilation to scale with the density. Indeed, when the density reaches the typical values of the PNS interior the nucleon-nucleon bremsstrahlung emission dominates. In practice, the core temperatures at densities larger then a few times 10 13 g cm −3 do not reach the values predicted from Eq. 15, and therefore bremsstrahlung rates dominate over the electron-positron annihilation even more at the highest densities. For the electron-positron pair annihilation, the emissivity derived from the parametrization and the one from the kernel treatment are the same, as expected since the underlying interaction is the same.
We briefly comment on the differences between the bremsstrahlung treatments. The difference between the T-matrix and OPE treatment is very obvious for densities over 10 14 g cm −3 . There, the prescriptions derived from the OPE potential give emissivities more than 10 times greater than the T-matrix prescription. This suppression of the rates at high densities, and also the more modest enhancement of the rates at low density when compared to the OPE interaction is a consequence of the T-matrix treatment [38,46]. The parametrization, which is based on the non-degenerate limit of the OPE generally produces comparable rates for the conditions used here. However, we note that the high temperature at nuclear densities resulting from Eq. 15 are higher then expected during the cooling phase and therefore under those conditions we would expect a larger deviation of the simplified rate from the OPE results. The rates that are expected to be important during the CCSN evolution are the ones near and around the neutrinospheres where the neutrinos are decoupling from the matter. At high densities, the neutrinos are in equilibrium and the precise rate does not matter, and at low densities the rate is so low that it does not contribute appreciable to the overall neutrino emission. As pointed out in [47], the key densities are around ρ 10 12 g cm −3 during the early core-collapse phase and upward of ρ ∼ 10 14 g cm −3 for the cooling phase. Over and above this, it is important to note that the many competing neutrino rates, and their strong temperature dependence, like electron-positron annihilation, often reduce the impact of changes in any one rate.
In addition to the differences that arise from the different interactions (in the case of bremsstrahlung), differences in the actual dynamical evolution can stem from the differences in the transport treatment. As discussed above in § II A, for the simplified methods, the final state neutrino blocking is not taken in account properly for the emission, nor is the precise form of the annihilation interaction used, rather simplified emission and absorption coefficients are used. With our systematic exploration of these interactions we aim to decipher these differences.

C. Setup
We performed a set of six simulations with two different progenitors for a total of 12 simulations. A 20-M , solar-metallicity, iron-core progenitor [48] and a 9.6-M zero-metallicity iron-core progenitor [49] are used. We utilize the 20-M progenitor as it is the same as the one For the electron-positron annihilation we show the emissivity based on the kernels (solid blue) and our parametrization of them (dashed blue), both from Bruenn (1985) [40]. We note that for the two electron-positron interactions we expect the same emissivities as the underlying interaction is the same.
studied in [6] where the evolution was computed using a variety of state-of-the-art evolution codes. The variation done in our study is on the transport treatment of the neutrino pair processes, the remaining physics is held constant. This is an interesting first step to gauge the influence of the different treatments and allows us to quantify the variations against the variations seen between different codes. For this progenitor, we used a grid containing 600 zones with the inner grid spacing being fixed at 300 m for the inner 20 km and increasing logarithmically outwards until ∼ 1.3×10 10 cm. This progenitor has been explored in many studies, but in particular, Ref. [50] also consider variations on the neutrino pairproduction processes. The other progenitor we consider has a ZAMS mass of 9.6M . Unlike most iron-core pro-genitors, this one has the peculiarity to explode in 1D. Although multidimensional effects can and do impact the development of the explosion in this model [51], these spherically symmetric simulations give us general insight on the behaviour of the explosion energy development over time and on the neutrino-interaction dependence of the early cooling phase. For this progenitor, we used a spherically symmetric grid of 800 zones with a constant grid spacing of 300 m in the inner 20 km and then a logarithmically increasing zone size until ∼ 1.3 × 10 9 cm. This progenitor has been used in multidimensional studies [51][52][53]. For all of the simulations, we used the SFHo equation of state from Steiner et al. [54] with the same neutrino physics (other than the pair-production treatments) as [6]. The simulation time step is set by the radiation and is equal to the light crossing time of the smallest zone and a CFL condition of 0.4 before bounce, 0.1 near bounce and 0.5 from 20 ms after bounce for all the simulations. We used a logarithmically spaced energy grid for the neutrinos from 1 MeV to 250 MeV with 18 energy groups.

III. RESULTS
In this section, we explore the impact of the different treatments of heavy-lepton neutrino pair-production described in Sec.II A and Sec.II B on the supernova evolution. For this, we apply the six different combinations of the pair processes treatments described in Tab. II. We will first explore the impact on the 20-M progenitor evolution and follow with the exploding 9.6-M progenitor evolution.

A. s20 progenitor
The 20-M progenitor does not lead to an explosion. The shock radius evolutions are plotted in the top panel of Fig. 2. The different colors correspond to the different models in Tab.II. The blue, green, and red solid lines refer to the three simulations using the electron-positron annihilation kernels with the bremsstrahlung fit, OPE kernel, and T-matrix kernel, respectively, while the three dashed lines refer to the electron-positron annihilation simplified emissivity for the three different bremsstrahlung treatments. All of the models give qualitatively similar results. Bounce occurs at ∼298 ms after the onset of collapse. The shock then expands for ∼90 ms after bounce and reaches a radius of ∼150 km where it stalls for ∼10 ms and starts to recede. The shock radius shows a short expansion phase again at ∼230 ms after bounce, which is due to the silicon-oxygen shell interface accreting through the shock front. The shock radius then continues to recede to attain ∼50 km at 500 ms after bounce. For reference, we show with grey lines the shock radius evolution from simulations with various codes for the same progenitor and setup taken from the comparison study of [6]. The shock evolution of all our models generally agree with these simulations and the level of variation between our simulations is slightly less than that observed between the simulation codes. The different neutrino-pair production treatments only modestly impact the shock radius evolution. For the sim-ulations using the full kernel treatment for the electronpositron-annihilation to neutrino-pair process (models 4, 5, and 6) there is a consistently lower shock radius (∼5 km) compared to the models with the simplified emissivity for this process (models 1, 2, and 3). This hierarchy is correlated with the properties of the heavylepton neutrino emission (bottom panels of Fig.2). As we discuss below, during the first ∼150 ms after bounce, the largest heavy-lepton neutrino luminosities and mean energies arise from the simulations using the full kernel treatment of the electron-positron annihilation pairproduction process. These simulations give enhanced cooling, smaller PNS radii, and a smaller shock radii. This cause-and-effect is commonly seen in this model (and for which the explosion properties in multidimensional simulations of this model are particularly sensitive too), for example with modifications on the neutralcurrent scattering opacities [55,56].
In the bottom two panels of Fig.2, we show the heavylepton neutrino mean energy (middle panel) and the heavy-lepton neutrino luminosity (bottom panel) as measured in the coordinate frame at 500 km. For the luminosities, after a small peak at bounce, there is a short rise to a plateau around ∼ 35 × 10 51 erg s −1 at the time of the peak shock radius. The heavy-lepton luminosities then decrease as the PNS contracts reaching values ∼ 10 52 erg s −1 at 500 ms after bounce. For the heavylepton neutrino mean energies, after a short peak at bounce, the mean neutrino energy rises from ∼30 ms after bounce from ∼ 14.5−15 MeV to a peak of ∼16-16.5 MeV. With the accretion of the silicon/oxygen interface the heavy-lepton neutrino mean energy drops ∼ 1 MeV and generally plateaus at ∼15-15.5 MeV until the end of the simulation at 500 ms. As we have shown for the shock radius, we show the neutrino luminosities and mean energies from [6] in grey. We can see that the different neutrino pair-production formalisms create differences which are comparable to the variability seen across different transport methods and hydrodynamics. It is worth noting that in [6], the prescriptions of the treatment of heavy-lepton neutrinos also varied among the codes. The impact of the different pair-production treatments on the electron-type neutrino luminosities and mean energies (not shown) is small.
During all stages of the evolution the quantities in Fig. 2 are within ∼10% of each other for the luminosities and within ∼3% for mean energies. However, the differences seen do correlate with the different pair-production treatments. Models 4, 5, and 6, where we use the full kernel-based treatment for the electron-positron annihilation process, have the largest neutrino luminosities and mean energies during the first ∼150 ms after bounce, while the simplified electron-positron annihilation treatment (models 1, 2, and 3) shows consistently lower luminosities and energies during this time. As we mentioned above, this causes increased PNS contraction and lower shock radii for the former models. However, also as a consequence of the increased contraction, there is increased electron neutrino mean energies, and an increased specific neutrino heating (although less overall heating due to the smaller gain region). For all models, the luminosity differences mostly disappear starting at ∼200 ms after bounce, although some differences in the mean energy remain as we will discuss below.
The impact of the nucleon-nucleon bremsstrahlung treatment on the evolution is less obvious. We do observe that among the different bremsstrahlung treatment, the use of T-matrix formalism (models 3 and 6) systematically creates a higher neutrino mean energy throughout the entire simulation, but especially after ∼250 ms after bounce. This is due to the lower emissivity of this interaction at higher densities (see Fig. 1) which gives an earlier decoupling radius and therefore a harder spectrum, since the matter temperatures are higher. The luminosities also tend to be the lower soon after bounce when using this formalism. The differences between the use of the OPE potential kernel-based formalism and the simplified emissivity based on this same potential only appear in the luminosities, and even there it is minimal. It has the effect of reducing the luminosity for ∼160 ms following bounce, analogous to the use of the simplified emissivity for the electron-positron annihilation, but smaller in magnitude.
From these observations we conclude that the differences created by the use of the simplified emissivities mainly lie in simplistic treatment of the neutrino transport (i.e. ignoring the functional form of the neutrino and antineutrino distributions and their angular dependence as well as any final state blocking, as explained in § II A) rather than differences in the underlying neutrino interaction model. A previous study, [50], explored the impact of a simplified heavy-lepton neutrino pairproduction treatments as well. They find similar changes on the luminosity, mean energy and shock radius evolution as the ones we find comparing models 1 and 5. They suggest that the differences seen are a result of the implicit assumption of the angular dependence (i.e. that it is isotropic) of the neutrino annihilation partner, rather than the in situ distribution, which is forward peaked (we note we from Eq. 12 that the annihilation strength is minimal for co-travelling neutrinos). This over-predicts neutrino-antineutrino annihilations within the simplified emissivity assumption. While this is certainly true, we note that since the neutrino annihilations are occurring well below the scattering surface, the distribution function is very isotropic. We therefore suggest it is rather the overall magnitude of the occupation density of the annihilation partner (which is implicitly assumed to be the black body distribution) that causes the simplified emissivity to over-predict annihilation and thus lead to smaller emergent heavy-lepton neutrino luminosities. We show this in the following. In both cases, we use the T-matrix kernel treatment for the bremsstrahlung interaction. In red we show the radial evolution of the total outgoing neutrino luminosity, normalized to the value to the kernel treatment at 500 km. In blue we show the relative difference between the actual and equilibrium neutrino distribution function, while in green we show the flux factor. For the former radial profile we select an energy bin with ∼15 MeV. Vertical lines denote the peak of the neutrino luminosity for the simplified treatment. electron-positron treatment; solid line). In red, we show the growing outward going heavy-lepton neutrino luminosity, i.e. 4πr 2 F r , normalized so that the luminosity from the full kernel treatment is 100% at 500 km. In blue, again for model 3 and 6 with dashed and solid line, respectively, we show the difference between the equilibrium heavy-lepton neutrino distribution and the actual heavy-lepton neutrino distribution relative to the heavylepton equilibrium distribution, i.e. (f eq − f )/f eq = 1 − f /f eq for an energy bin corresponding to ∼15 MeV.
Here we take f = J/(4πν 3 ) and f eq = 1/(exp(ν/T ) + 1). In green, we show the energy averaged heavy-lepton flux factor (= F r /E). This allows us to highlight the proportional importance of the anisotropy of the neutrino field (the green lines) and the deviation of the actual neutrino distribution from equilibrium (blue lines). At the early time (left panel), at the radii where the neutrinos luminosity is rising (∼20-40 km) the distributions are almost isotropic (F r /E ∼ 15%). However at these same radii, the occupation density of the neutrinos significantly deviates from the black body, falling short of the equilibrium occupation density by ∼ 60% at ∼40 km. Since the simplified treatment implicitly assume a black body distribution for the annihilation partners, this leads to an over prediction of the annihilation rate in this regime. The result is a further out decoupling radius, and ultimately a lower heavy-lepton neutrino luminosity and mean energy. This implies that a significant factor in the difference between a full and simplified treatment is linked to the assumed distribution of the pair neutrino rather than the intrinsic anisotropy of the radiation field. For late times (right panel), at the radii where the luminosity is rising (∼16-22 km), the distribution is almost completely isotropic (rising to F r /E ∼ 5% at 22 km) and the deviation of the distribution function from the equilibrium distribution is negligible at all radii except for the last 20% of the emission, even there, the deviation is at most ∼30%. As a result we do not see any excess annihilation at these late times. It is worth noting that the simplified treatment does actually predict a larger luminosity at these late times, by a few percent. We suspect this is due to lack of final state blocking in the emission of the neutrinos within the simplified treatment. At late times the value of the heavy-lepton neutrino distribution function near the peak of the emission is several times the value seen at earlier times, raising the impact of the final state blocking. Some oscillations appear in the luminosity for late times. This is a consequence of the explicit-in-time, first-order, forward-Euler scheme used for the evaluation of the spatial fluxes. It tends to occur for energy groups that are close to free streaming (with the characteristic speeds entering the Riemann problem approaching c) in zones with a CFL condition c∆t/∆x ∼ 0.5. It causes the small oscillations seen near ∼30 km and the instabilities in the outward going neutrino flux under 15 km, which at this point is constrained to the lowest energy groups as the rest are fully trapped.

B. z9.6 progenitor
In order to show the impact of the different heavylepton pair process neutrino treatments on an exploding model, we evolved one of the few progenitors known to explode in 1D [51,53], a 9.6M , zero-metallicity star evolved with the KEPLER stellar evolution code [49]. We simulate, using GR1D, the six combinations of thermalinteraction models listed in Table II. All models successfully explode. For all models, the shock expansion does not stall as in traditional iron-core collapse progenitors like it does in Fig. 2 for the 20M progenitor. The shock radius evolution can be seen in the top panel of Fig. 4. Here we plot the location of the maximum velocity gradient which flags accurately the position of the, sometimes multiple, different shocks. The explosion times vary across the models, showing the sensitivity of the explosion to these neutrino interactions. The earliest explosion time, here arbitrarily defined as when the shock passes 1000 km for the last time, is ∼310 ms after bounce while the latest explosion is ∼470 ms after bounce, ∼50% longer. It is also notable that the initial shock formed at bounce is not the one that ultimately leads to the final explosion. In all cases, the first shock expands but is not energetic enough to runaway. The accretion of this material onto the PNS creates a burst of neutrino heating (bottom panel of Fig. 4) and a secondary shock which will ultimately lead to an explosion. The formation time of this second shock corresponds to the heating peaks in the lower panel as the in-falling matter is compressed and the neutrino heating increases. In general, the models using the simplified treatment for electron-positron annihilation show a strong initial shock expansion phase, but a later ultimate explosion time, while the models with the full kernel treatment of electron-positron an- nihilation show a lower initial shock expansion and an earlier explosion. We discuss the difference seen between the different interaction models in more detail below.
In Fig. 5, we show the heavy-lepton neutrino quantities for this progenitor model. In the top and bottom panel we show the mean energy and luminosity, respectively for each of the six neutrino pair-production treatments shown in Table. II. After a sharp and short peak at bounce, the mean energies plateau for the first ∼100 ms after bounce around 15 MeV with a spread ∼0.5 MeV. The mean energies decrease over the remaining 400 ms of the simulations by ∼1-1.5 MeV, depending on the treatment used, with the T-matrix treatment maintaining the highest mean energy, as was the case in the 20M progenitor above. The luminosities present a peak just after bounce and then slowly decrease until the end of the simulation. The heavy-lepton neutrino luminosities and energies present a similar dependence on the explored interactions as observed for the 20M progenitor with the models using the full kernel treatment of the electron-positron pair annihilation having higher luminosities soon after bounce and then reaching similar, but slightly lower, values to the simplified formalism at later times. As for the case of the 20-M progenitor, we attribute these differences to the treatment of the transport in the simplified models. Particularly, the differences at early times are attributed to the form of the effective absorption coefficient which incorrectly treats the distribution of the annihilation partner and the smaller difference at late times due to the assumption of no final state blocking for the neutrinos during the emission process. At late times, the largest difference in the neutrino quantities arises from the use of the T-matrix bremsstrahlung kernels where the lower opacity at higher densities gives rise to higher neutrino energies, by ∼3%, this is seen in both the simplified and the full kernel treatment. We conclude our discussion of the 9.6-M progenitor by examining the impact of the different neutrino interaction models in Tab. II on the development of the explosion and the shock propagation. The largest systematic difference we observe between the models using different interactions is the shock evolution for the three simulations that use the simplified treatment of electronpositron annihilation to neutrino pairs versus the three models that use the kernel-based treatment. This is seen in Fig. 4, but we show further evidence for this in Fig. 6. In the six panels of Fig. 6, we show the accretion history for each simulation. Blue colors denote negative accretion rates (infalling material) while red colors show positive accretion rates (expanding material). The three models for which we use the simplified treatment of electron-positron annihilation (models 1, 2, 3) are shown on the top row while the three models using the full kernel treatment (models 4, 5, 6) are shown on the bottom row. The left, middle, and right columns include the simplified bremsstrahlung treatment, the OPE kernel treatment, and the T-matrix kernel treatment, respectively. The three models with the simplified electron-positron annihilation pair production treatment have a slightly faster initial shock expansion from higher neutrino heating (see Fig. 4), with the T-matrix bremsstrahlung treatment having the largest such expansion. This causes these three models to undergo a strong initial shock acceleration phase starting around ∼130 ms after bounce. Although, except for a small region directly behind the shock during this expansion phase, matter is mainly accreting in the postshock region. These shocks eventually fail. It is worth noting that in multidimensional simulations of this progenitor [51,53] the explosions typically set in during this period as the added role that multidimensional effects like convection and turbulence play is enough to initiate the explosion. However, in our spherically symmetric simulations this is not the case. Eventually, secondary shocks form at the surface of the PNS between ∼200-240 ms after bounce concomitant with the increased accretion rate from the failing shock (see the dark blue regions around ∼100-150 km at this time) which eventually give the ultimate explosion and the beginnings of a neutrino driven wind. The three models with the full kernel treatment for electron-positron annihilation do not undergo this accelerated expansion and continue to mildly expand until ∼160 ms at which point the shock fails and the secondary shock forms at ∼180 ms after bounce. These models are the first to ultimately explode.
We compute the diagnostic explosion energy for each of the six models simulated, where φ is the gravitational potential, v is the fluid velocity, is the specific internal energy, and 0 is a reference zero-point taken for simplicity to be the value of the internal energy of the EOS for the same density and Y e but for T = 0. We only consider contributions to the diagnostic energy from outflowing matter and where the integrand of Eq. 16 is positive, this is a proxy for unbound material. The results are shown Fig.7. In general, for the 9.6-M progenitor the explosion energy depends on the time of the explosion. For early explosions in multiple dimensions, for example see [51,53], the explosion energies can reach several to 10 times 10 49 erg. In our relatively late spherically-symmetric explosions, we FIG. 6: Mass accretion rates for the six different interaction sets explored vs. time. Blue colors denote accreting material, red denotes expanding material. Overlaid in grey are the location of the (sometimes multiple) shock fronts. see explosion energies ∼ 10 49 erg for the five models that explode first. The last model to explode achieves only 0.4 × 10 49 erg owing to the much lower neutrino heating rates present at the later times. We see that while the rise of the explosion energy is strongly correlated with the onset of the explosion, the ultimate explosion energy is not strictly dependent on the explosion time. We see that differences can arise due to the complex interplay of the neutrino heating of the fallback material from the initial failed shock (see Fig. 6).
Finally, we note that − 0 in Eq. 16 is an estimate of the sum of the thermal energies and the recombination energy that will be converted to thermal energy as the matter recombines. We tested this against an explicit calculation similar to that of [51]. If we assume that all the free nucleons and alpha particles present in the matter would recombine to iron this sets the specific recombination energy. We use this, along with the specific internal energy from the Helmholtz EOS [57] at a matter temperature of T to replace − 0 in Eq. 16. We find comparable, to within 5%, diagnostic explosion energies as shown in Fig. 7.

IV. CONCLUSION
In this paper, we highlighted and explored the importance of heavy-lepton (ν µ , ν τ , and their antiparticles) neutrinos and the interactions that produce them in the scope of spherically-symmetry, fully general-relativistic neutrino-driven CCSNe simulations. We performed systematic simulations on two different progenitors, a solar metallicity progenitor star with a ZAMS mass of 20M which has been extensively used in the literature for CCSN simulations, including a recent comparison study and a zero metallicity progenitor star with a ZAMS mass of 9.6M that has the peculiarity to explode even in spherically symmetric simulations. We simulate the core collapse and the early post-bounce phase using the open-source software GR1D and Nulib, which we update accordingly with the work presented in this paper. In particular, we test the importance of the two main heavy-lepton neutrino pair-production processes in CCSNe, electron-positron annihilation to a neutrino-antineutrino pair and nucleon-nucleon scattering that radiates a neutrino-antineutrino pair (i.e. nucleon-nucleon bremsstrahlung). We explore two main effects. First, we study the neutrino transport implementation of both of these interactions by utilizing a simplified approached with effective emission and absorption coefficients and a complete treatment utilizing complete scattering kernels and in situ neutrino distribution functions. The aim for this part of the study is to assess the impact and quantify systematic effects of the simplified treatment on CCSN simulations with the goal of providing a robust prescription for use in multidimensional simulations. Second, we explored two independent nuclear-physics based prescriptions for the nucleonnucleon bremsstrahlung interaction (in addition to a simplified approach as mentioned above). One of these interactions is commonly used throughout the literature for CCSNe simulations and is based on the one-pion exchange formalism [37]. The other interaction [38], formulates the nucleon-nucleon bremsstrahlung scattering kernel based on a T-matrix formalism where the underlying interaction is constrained by experimental phase shifts, see also [46]. From these variations we arrive at six combinations of interactions to explore with our CCSN simulations. We find that overall the simplified neutrino interactions do a fair job at reproducing the neutrino quanti-ties (within ∼10%) and the dynamics of the core collapse event, although potentially important differences are present that need to be considered when employing the simplified treatment for precision simulations. We find the simplified method under predicts the heavylepton neutrino luminosity in the early stages, by ∼10%. We show that this is the result of the simplistic treatment of the distribution of the annihilation partner, i.e. the assumption that it follows the blackbody distribution. This leads to an overestimate of the annihilation rate as the neutrinos begin to free stream away from the CCSN core. In the region where this excess annihilation is occurring the distribution function is quite isotropic and therefore, as opposed to that suggested in previous works, this difference is unlikely due to the assumption of isotropy in the simplified treatment. While the largest impact is seen with different treatments for the electronpositron annihilation interaction, we see the same (but much smaller in impact) trend with the simplified treatment of nucleon-nucleon bremsstrahlung. At later times the simplified treatment agrees much better with the full kernel treatment, owing to the fact that the distribution of the annihilation partner is much closer to the blackbody distribution in the regions where the neutrinos are decoupling. We do see a slight overestimate of the luminosity for the simplified treatment, which we attribute to the assumption of no final-state neutrino blocking in the simplified emission treatment. Finally, we comment on the impact of using a different microphysical interaction for the nucleon-nucleon bremsstrahlung. The use of the T-matrix formalism over the standard one-pion exchange treatment systematically increases the heavylepton neutrino energy by ∼5% which we infer is due to the reduced interaction strength of the T-matrix kernel at larger densities compared the one-pion exchange kernel. This causes the heavy-lepton neutrinos to begin to decouple deeper into the PNS core, where the temperature is higher.
For cases where the dynamics can sensitively depend on the neutrino physics, for example with the 9.6M progenitor studied here, we find that the different interactions explored can impact the heating enough at the earlier stages to quantitatively effect the development of the explosion, including the ultimate explosion time and the explosion energy. To what extent this carries over to multidimensional simulations remains to be seen, although it is important to note that even in multidimensional simulations we know there is sensitivity to the neutrino physics at the ∼10% level demonstrated here [55,56].