Neutrinos play a critical role of transporting energy and changing the lepton density within corecollapse supernovae and neutron star mergers. The quantum kinetic equations QKEs combine the effects of neutrinomatter interactions treated in classical Boltzmann transport with the neutrino flavorchanging effects treated in neutrino oscillation calculations. We present a method for extending existing neutrino interaction rates to full QKE source terms for use in numerical calculations. We demonstrate the effects of absorption and emission by nucleons and nuclei electron scattering electronpositron pair annihilation nucleonnucleon bremsstrahlung neutrinoneutrino scattering. For the first time we include all these collision terms selfconsistently in a simulation of the full isotropic QKEs in conditions relevant to corecollapse supernovae and neutron star mergers. For our choice of parameters the longterm evolution of the neutrino distribution function proceeds similarly with and without the oscillation term though with measurable differences. We demonstrate that electron scattering nucleonnucleon bremsstrahlung processes and fourneutrino processes dominate flavor decoherence in the protoneutron star PNS absorption dominates near the shock and all of the considered processes except elastic nucleon scattering and neutrinoneutrino processes are relevant in the decoupling region. Finally we propose an effective decoherence opacity that at most energies predicts decoherence rates to within a factor of 10 in our model PNS and within 20 outside of the PNS.
['Astrophysics - High Energy Astrophysical Phenomena', 'High Energy Physics - Phenomenology', 'Nuclear Theory']
Neutrino quantum kinetics in compact objects
{'Neutrino Quantum Kinetics in Compact Objects': 'Sherwood A. Richers, ∗ Gail C. McLaughlin, James P. Kneller, and Alexey Vlasenko North Carolina State University  \nNeutrinos play a critical role of transporting energy and changing the lepton density within corecollapse supernovae and neutron star mergers. The quantum kinetic equations (QKEs) combine the effects of neutrino-matter interactions treated in classical Boltzmann transport with the neutrino flavor-changing effects treated in neutrino oscillation calculations. We present a method for extending existing neutrino interaction rates to full QKE source terms for use in numerical calculations. We demonstrate the effects of absorption and emission by nucleons and nuclei, electron scattering, electron-positron pair annihilation, nucleon-nucleon bremsstrahlung, neutrino-neutrino scattering. For the first time, we include all these collision terms self-consistently in a simulation of the full isotropic QKEs in conditions relevant to core-collapse supernovae and neutron star mergers. For our choice of parameters, the long-term evolution of the neutrino distribution function proceeds similarly with and without the oscillation term, though with measurable differences. We demonstrate that electron scattering, nucleon-nucleon bremsstrahlung processes, and four-neutrino processes dominate flavor decoherence in the protoneutron star (PNS), absorption dominates near the shock, and all of the considered processes except elastic nucleon scattering and neutrino-neutrino processes are relevant in the decoupling region. Finally, we propose an effective decoherence opacity that at most energies predicts decoherence rates to within a factor of 10 in our model PNS and within 20% outside of the PNS.', 'I. INTRODUCTION': "Neutrinos interacting with matter are central to the mechanism by which massive ( > ∼ 10 M ⊙ ) stars explode as core-collapse supernovae (CCSNe, see, e.g., [1-3]). When the iron core of such a star exceeds its effective Chadrasekhar mass, it becomes unstable to collapse and the center reaches nuclear densities within a couple hundred milliseconds. At this point, nuclear forces cause the equation of state to quickly stiffen, halting collapse and launching a shock wave into the still supersonicallyinfalling iron core. This initial shock is not strong enough to lead to an explosion. Neutrinos emitted as matter accretes onto the protoneutron star (PNS) have the potential to transport enough energy to the fluid under the shock front to lead to an explosion. However, a detailed accounting of precisely how this happens remains elusive.  \nNumerical simulations have become the primary tool for studying the nonlinear and multi-physics dynamics in CCSNe. Current state of the art simulations account for neutrino transport by treating neutrinos as classical radiation and employing a variety of implementations of the relativistic Boltzmann equation or approximations thereof (e.g., [4-10]). Over the years since the first numerical model of a neutrino-driven supernova by Colgate and White [1], many different neutrino-matter interaction processes have been found to be important, including absorption and emission by nucleons and nuclei, scattering by nucleons and electrons, pair production and annihilation, nucleon-nucleon bremsstrahlung radiation, and neutrino-neutrino pair annihilation and scattering (see [11, 12] for overviews).  \nHowever, it has long been known that neutrinos are able to change flavor in-flight[13-16], though many of the relevant parameters are not yet well constrained [17]. Understanding how the neutrino flavor evolves from the neutrino's emission in the CCSN to its detection on earth is necessary in order to use neutrino signals to gain valuable insight into the depths of the CCSN explosion mechanism (e.g., [18]). Historically, calculations of flavor conversion of CCSN neutrinos have been performed by evolving neutrinos moving outward from the neutrinosphere (see [19, 20] for recent reviews). In contrast to neutrino transport calculations, simulations of neutrino flavor conversion neglect collisions but include some aspects of neutrino-matter interactions in the form of a potential through which neutrinos propagate. Neutrino-neutrino interactions also contribute to this potential, making the evolution equations very nonlinear and difficult to simulate over timescales relevant to the CCSN.  \nIf neutrinos are experiencing flavor oscillations and collisions in widely separated regions of the supernovae then the no-oscillation approximation in neutrino transport calculations and the no-collision approximation in neutrino flavor conversion calculations at first glance appear very reasonable. Neutrinos are emitted in flavor states inside of the supernova shock wave and the large matterinduced contribution to the neutrino potential suppresses flavor oscillations in regions where collisions are important [21]. However, including the contribution of neutrinos themselves to this potential (namely, the 'selfinteraction' part of the potential) makes the problem much more rich. The first calculations including the full potential but only one neutrino energy and propagation direction suggested that flavor oscillations could occur within the shock [22]. Subsequent calculations that include multiple energies and/or propagation directions indicate that, though still occurring inside the shock, oscil-  \nlations begin at too large of a radius to significantly affect the CCSN explosion mechanism [23-25]. Once the neutrinos are beyond the strong matter potential, the matter density is so low that scattering events are exceedingly rare, making collision terms unimportant where flavor oscillations take place.  \nOne could take this as evidence that a full quantum kinetic treatment including both collisions and oscillations is not needed, but recent developments have suggested that there may be scenarios where neutrino oscillations and scattering may not be so separable. So-called fast flavor conversions may occur deep within the supernova shock wave if the electron neutrino and electron antineutrino angular distributions differ enough that they intersect in momentum space. In this way, in some directions the contribution to the potential from antineutrinos overwhelms that from neutrinos, allowing neutrino flavor states to mix [26-32]. Though the lepton-emission selfsustained asymmetry [5, 33] may be conducive to fast flavor conversions[28, 34], searches for conditions conducive to fast flavor conversions in simulations of CCSNe have had mixed results [33, 35, 36].  \nThe neutrino halo effect is another such scenario where collision and oscillation physics are simultaneously relevant. Even though neutrinos are decoupled outside of the shock, the small changes to the neutrino distribution from rare scattering events can enhance the neutrino selfinteraction potential and significantly change the location and strength of the flavor conversions [37-41]. Though the halo effect is not expected to have an impact on the CCSN explosion mechanism, it could significantly change the neutrino signal we might expect to see from the next nearby CCSN.  \nNeutron star mergers also create an environment with complex neutrino radiation fields, and the details of this radiation field have profound effects on the amount of mass ejected from the merger, the eventual composition of this ejecta, the fate of the central remnant, and potential launching of a relativistic jet (e.g., [42-50]). It has been shown that the neutrinos are also likely subject to a number of interesting flavor transformation effects that are expected to occur near the decoupling region (e.g., [51-60]). However, the effects of quantum kinetics in these systems is yet to be investigated in detail.  \nTo account for both neutrino flavor conversions and collisions consistently, one must evolve the neutrino quantum kinetic equations (QKEs) [61-65]. Though the collision physics for neutrinos in flavor eigenstates has been extensively explored (see, e.g., [11, 12] for overviews) and the general form of the collision integral was worked out in [61], how neutrinos in general flavor states interact with matter through the specific collisional processes relevant to CCSNe and neutron star mergers was only recently fleshed out by [65]. However, these detailed interaction physics have yet to be explored thoroughly in numerical simulations.  \nIn this paper, we begin in Sec. II with a brief introduction of the QKEs and the concept of distribution flavor  \nvectors that are critical in visualizing the results. We discuss our numerical method for solving the QKEs and our choice of initial conditions in Sec. III. Then in Sec. IV we describe a method by which the existing wealth of understanding of interaction rates for neutrinos in flavor eigenstates can be extended in a straightforward way to the full collision kernel for the QKEs. We also show results from nonoscillating calculations in this section to build intuition about the impact of each process on the flavor dynamics. We present our main results in Sec. V. Sec. V A shows the evolution of flavor decoherence resulting from a combination of all of the interactions discussed in this paper without the oscillation term. In Sec. V B we describe the same with the oscillation term turned on and compare it to the nonoscillating calculation. Finally, we perform nonoscillating calculations using input parameters from many points along a 1D core-collapse supernova simulation snapshot and suggest and approximate decoherence length scale formula in Sec. V C. To facilitate future simulations of the QKEs in core-collapse supernovae and neutron star mergers, we write the general relativistic moment form of the QKEs and the corresponding oscillation and collision source terms in Appendix A.", 'II. QUANTUM KINETIC EQUATIONS': "In this section, we present the QKEs in general relativistic form to allow future implementations in codes that account for spacetime curvature and velocity dependence. The QKEs describe the evolution of the (dimensionless) Wigner transform of the neutrino two-point correlation function [62]  \nF ABab ( ν, Ω , x µ ) := [ f LLab f LRab f ∗ LRab f RRab ] , (1)  \nwhere f LL and f RR are N f × N f matrices representing left- and right-handed neutrinos, respectively and N f is the number of neutrino flavors. Throughout this paper, uppercase latin indices represent handedness indices, lowercase latin indices are flavor indices, greek indices are spacetime indices, and we use the ( -, + , + , +) sign convention for the metric. The flavor-diagonal ( a = b ) components of f LL and f RR are the occupation probabilities at particular spacetime coordinates x µ of neutrinos with frequency ν and direction 3-vector Ω in a comoving orthonormal tetrad. That is, f AAaa are the ordinary neutrino distribution functions and the number density n Aa of neutrino flavor a and handedness A in such a frame is  \nn Aa ( x µ ) = ∫ d 3 ν c 3 f AAaa , (2)  \nwhere d 3 ν = ν 2 dνd Ω . The off-diagonal elements of f LL and f RR describe quantum flavor coherence. f LR is also a N f × N f matrix that represents quantum coherence between left- and right-handed neutrinos. In the ultrarelativistic limit, the QKEs can be written to order ϵ 2 ,  \nwhere ϵ ≪ 1 is the ratio of the neutrino mass, mass splitting, forward-scattering potentials, or gradients to the neutrino frequency, as  \nd F dλ +force + drift = -p µ u µ ( C i ¯ hc [ H , F ] ) , (3)  \nwhere u α is the dimensionless fluid four-velocity. λ is an affine parameter such that the neutrino four-momentum is p µ = dx µ /dλ . In the comoving orthonormal tetrad, the momentum is also p µ = (1 , Ω ) hν/c , where the first component is the time component. The derivative expands in a general curved spacetime to d/dλ = p µ ∂/∂x µ -Γ µ αβ p α p β ∂/∂p µ , where Γ µ αβ ( x µ ) are connection coefficients (units of cm -1 ). The 2 N f × 2 N f collision integral C ABab ( ν, Ω , x µ ) (units of cm -1 ) can be decomposed as [62]  \nC = C + -C -= { 1 -F , ˜ Π + } -{ F , ˜ Π -} , (4)  \nwhere the calculation and use of the self-energies ˜ Π ± ABab ( ν, Ω , x µ ) will be described in more detail in Sec. IV. The oscillation potential H ABab ( ν, Ω , x µ ) (units of ergs) is described below. If neutrinos are Majorana particles, only F needs to be evolved, but if neutrinos are Dirac particles, an additional antineutrino field ¯ F must be evolved with an analogous equation. Throughout the rest of this paper, we will refer to the N f × N f matrices f ab and ¯ f ab as simply the neutrino and antineutrino distribution functions, respectively.  \nThe terms labeled 'force' and 'drift' in Eq. 3 are O ( ϵ 2 ) corrections to the d/dλ term that account for refractive effects due to the finite masses of neutrinos (see, e.g., [62, 64]). We neglect them in this work for simplicity, and more work needs to be done to assess the importance of these terms in the context of core-collapse supernovae. We will also assume for simplicity that there is no spin coherence ( f LR = ¯ f LR = 0) and, in the case of Dirac neutrinos, there are no right-handed neutrinos and no left-handed antineutrinos ( f RR = ¯ f LL = 0). Spin coherence effects are not expected to be important in CCSNe, since the neutrino-antineutrino-mixing contributions to the potential are suppressed by an additional factor of ϵ and neutrinos in typical CCSN profiles likely pass through the associated resonance too rapidly [64, 66, 67]. However, spin coherence effects may not be negligible in other environments such as neutron star mergers [54]. Under the assumption of no spin coherence, we can separate left-handed neutrinos from righthanded (anti)neutrinos, and Dirac and Majorana neutrinos evolve identically. The QKEs for f = f LL and ¯ f = f RR (Majorana) or ¯ f = ¯ f RR (Dirac) then become  \np µ ∂f ∂x µ -Γ µ αβ p α p β ∂f ∂p µ = -p µ u µ ( C -i ¯ hc [ H,f ] ) p µ ∂ ¯ f ∂x µ -Γ µ αβ p α p β ∂ ¯ f ∂p µ = -p µ u µ ( ¯ C -i ¯ hc [ ¯ H, ¯ f ] ) , (5)  \nwhere C ab and H ab are the N f × N f collision integral and Hamiltonian, respectively, from the LL quadrants  \nof C and H . Similarly, ¯ C ab and ¯ H ab come from the RR quadrants of C and H (Majorana) or ¯ C and ¯ H (Dirac).  \nThe Hamiltonian operator is often decomposed as  \nH = H vacuum + H matter + H neutrino . (6)  \nThe Hamiltonian for antineutrinos is related to that for neutrinos by ¯ H vacuum = H ∗ vacuum , ¯ H matter = -H ∗ matter , and ¯ H neutrino = -H ∗ neutrino . The vacuum Hamiltonian is  \nH vacuum = UH ( m ) vacuum U † , (7)  \nwhere H ( m ) vacuum = diag( √ h 2 ν 2 + m 2 l c 4 ) is the vacuum Hamiltonian in the neutrino mass basis, m l is the mass of the neutrino corresponding to lepton flavor l . The unitary matrix U describes the mixing between the flavor and mass bases [68, 69]. The matter potential in the local comoving frame is  \nH matter = √ 2 G F ¯ h 3 c 3 diag( n l -n ¯ l ) , (8)  \nwhere n l and n ¯ l are the number density of charged lepton and antilepton of flavor l , though in the astrophysical systems of interest electrons are the only lepton with a significant abundance. Neutral current interactions with nucleons also technically contribute to the potential, but since they affect all flavors equally, the potential offset does nothing to modify oscillations and can be ignored. Finally, the neutrino self-interaction potential is  \nH neutrino = √ 2 G F ¯ h 3 ∫ d 3 ν ' (1 -cos θ )( f ' -¯ f '∗ ) , (9)  \nwhere cos θ = Ω · Ω ' . This is analogous to the matter potential, except that neutrinos are not in general isotropic and the cos θ term is needed to account for the angular dependence of the anisotropic neutrino distributions.", 'A. Distribution flavor vector': "It is useful to visualize a neutrino quantum state as a vector in flavor isospin space. We can also visualize the distribution functions f ab and ¯ f ab with a similar flavor isospin vector, but unlike that for an individual neutrino spin 1/2 quantum state, the distribution flavor vector can represent states with a trace different from unity. When working in a two-flavor system, the vector components are written as  \nf ab = f ( t ) δ ab + f ( x ) σ ( x ) ab + f ( y ) σ ( y ) ab + f ( z ) σ ( z ) ab f ( t ) = 1 2 ( f ee + f µµ ) f ( x ) = 1 2 ( f eµ + f µe ) f ( y ) = -i 2 ( f eµ -f µe ) f ( z ) = 1 2 ( f ee -f µµ ) , (10)  \nwhere σ ( α ) are Pauli matrices. We can define the length of the spatial part of the distribution flavor vector as  \nL := √ ⃗ f · ⃗ f , (11)  \nwhere in components ⃗ f = ( f ( x ) , f ( y ) , f ( z ) ). We recall that for two arbitrary matrices A and B , one can derive the following identities:  \n-i 2 [ A,B ] = ( 0 , ⃗ A × ⃗ B ) 1 2 { A,B } = A ( t ) B ( t ) ( 1 + ⃗ A · ⃗ B A ( t ) B ( t ) , ⃗ A A ( t ) + ⃗ B B ( t ) ) , (12)  \nwhere the first term in the parentheses is the 'time' component of the distribution flavor vector. These identities are useful in gaining intuition for how the oscillation and collision terms affect the neutrino distribution. Applying the first identity to Eq. 5, we see that the Hamiltonian term is incapable of changing the number of neutrinos (the t component of the anticommutator is zero) and that it is incapable of changing the length of the flavor vector ( d ⃗ f/dλ is perpendicular to ⃗ f ). The second identity applied to the same equation together with Eq. 4 shows that the collision terms can change both.  \nFinally, given values for f ee and f µµ , one can determine the maximum magnitude of the off-diagonal component of the distribution function (i.e., √ f 2 ( x ) + f 2 ( y ) ) by requiring that no diagonal component of the distribution function be smaller than 0 (i.e., L f ≤ f ( t ) ) or larger than 1 (i.e., L f ≤ 1 -f ( t ) ) in any basis. Since f ( z ) is determined by f ee and f µµ , the off-diagonal components must satisfy  \nf 2 ( x ) + f 2 ( y ) ≤ min( f ( t ) , 1 -f ( t ) ) 2 -f 2 ( z ) . (13)  \nThis will be useful for setting the initial conditions for our test calculations.", 'III. ISOTROPICSQA': 'We evolve the QKEs assuming isotropy, homogeneity, and flat spacetime in a stochastic, operator-split manner with the new, open-source code IsotropicSQA 12 . Under these assumptions, the neutrino QKEs [Eqs. 3] take the form of  \n1 c ∂f ∂t = C -i ¯ hc [ H,f ] . (14)  \nIn the absence of collisions, this equation can be rewritten as  \n1 c ∂S ∂t = -i ¯ hc HS, (15)  \nwhere the (unitary) evolution operator S determines the distribution function at a later time as  \nf ( t ) = Sf (0) S † . (16)  \nThe (Hermitian) Hamiltonian operator H was described in Sec. II. The corresponding equations for antineutrinos are exactly analogous. We assume a two-flavor system under the normal hierarchy with a mass splitting of ∆ m 21 = 2 . 43 × 10 -3 eV 2 and mass eigenstates rotated from flavor eigenstates by 9 · . The assumptions of isotropy and homogeneity preclude important spatial transport and multiangle effects (e.g., [41, 56]), but allow us to get a handle on the effects of the various contributions to the collision term on the flavor evolution of a neutrino distribution.', 'A. Maximally mixed initial conditions': 'To demonstrate the effects of each of the contributions to the collision term, we will evolve a flavor-mixed isotropic neutrino distribution function in time without oscillations. We configure the initial distribution function to be  \nf ( ν ) = [ FD( T, µ ν e , ν ) f ( x ) , max ( ν ) f ( x ) , max ( ν ) FD( T, 0 , ν ) ] ¯ f ( ν ) = [ FD( T, -µ ν e , ν ) ¯ f ( x ) , max ( ν ) ¯ f ( x ) , max ( ν ) FD( T, 0 , ν ) ] , (17)  \nwhere the Fermi-Dirac distribution is given by  \nFD( T, µ, ν ) = 1 e ( hν -µ ) /k B T +1 . (18)  \nThe electron neutrino chemical potential is µ ν e = µ e + µ p -µ n , where µ n , µ p , and µ e are the neutron, proton, and electron chemical potentials determined by an input equation of state (EOS). We use the HShen equation of state [70] to determine the chemical potentials used in the initial conditions and collision reactions because it is consistent with the Furusawa EOS [71, 72] used to produce the core-collapse supernova simulation snapshots used in Sec. V. We set the imaginary part of the off-diagonal components to zero and give the real parts positive values that maximally mix the distribution function at each energy according to Eq. 13. Our fiducial fluid parameters for the test calculations are ρ = 10 12 g cm -3 , T = 10MeV, and Y e = 0 . 3, which yields µ n = 1 . 39 MeV, µ p = -8 . 56 MeV, µ e = 10 . 1 MeV, and µ ν e = 0 . 0977 MeV (including contributions from masses).', 'B. Stochastic integration': 'In our implementation of the QKEs depicted in Fig. 1, we evolve neutrino oscillations without collisions using Eq. 15 punctuated by interaction steps at random intervals, as described below. We discretize the distribution  \nFIG. 1. IsotropicSQA Method Summary. Neutrino oscillations are evolved explicitly via the evolution matrix S (bottom horizontal lines) separately from collisions (top horizontal lines) and the effects are combined at time intervals of dt block . This allows us to take very large time steps for the expensive and slowly-acting interactions, and very short steps for the inexpensive and rapidly varying oscillations.  \n<!-- image -->  \nfunctions and evolution matrices into 50 energy bins centered on integer multiples of 2 MeV, each with a width of 2 MeV. This choice allows a direct implementation of the neutrino-neutrino collision terms presented in [65]. We initialize the evolution matrices S and ¯ S at each energy to the identity and calculate the Hamiltonian based on f and ¯ f [Eq. 6, depicted as vertical arrows in Fig. 1]. We then evolve S and ¯ S for a block in time dt block without collisions using an adaptive sixth-order Cash-Karp Runge Kutta integrator (depicted as the bottom horizontal arrows in Fig. 1). At the end of this block, we map the evolution matrix back onto the distribution function using Eq. 16 (depicted as diagonal dashed arrows in Fig. 1). We then evolve f starting from this updated value for an identical time dt block using the same sixth-order integrator (depicted as the upper line in Fig. 1). Then the process repeats for subsequent blocks until the end of the simulation.  \nThe oscillation term time step dt oscillate , the collision term time step dt interact , and the block time step dt block are all independent. dt oscillate and dt interact are dynamically adjusted to keep the maximum relative difference between fifth- and sixth-order solutions for each step within a critical value of 10 -12 . Such high accuracy is required to maintain an accurate solution over the large number of time steps in our simulations. The time step dt oscillate is typically much shorter than dt block , but dt interact is typically very large and is restricted by dt block . This separation of oscillation and collision integration allows us to take much larger time steps for the collision term than for the oscillation term, both for computational efficiency and for preventing the conversion between f and S , which suffers from truncation errors, from ruining the accuracy of the solution (see Appendix B).  \nIf dt block is left constant or is changed deterministically, there will be artificial correlations between successive interaction steps. For example, if dt block is a multiple of the oscillation timescale, the collision step will always act on the same point in the oscillation cycle, and this would not represent the physical solution where collisions affect the distribution continuously throughout the oscillations. To  \nprevent such aliasing errors, we randomize dt block based on a target time step dt block , target as  \ndt block = min( -log( U ) , 5) dt block , target , (19)  \nwhere U is a uniform random number between 0 and 1. The resulting exponential random time step makes the interaction sampling a Poisson process, which causes the times of the collision events to be uncorrelated and uniformly sampled in time. The limiter of 5 prevents an excessively large time step. The first dt block , target is an input parameter, and subsequent values of dt block are determined in order to keep the impact of the collision term within a block relatively small. We define the impact over the previous block as  \nI = max ∣ ∣ ∣ ∣ f ( t ) -f ( t -dt block ) L ( t ) ∣ ∣ ∣ ∣ , (20)  \nwhere L ( t ) is the length of the distribution flavor vector [Eq. 11] at the same energy and helicity as f , and the maximum is taken over all matrix components, helicities, and energies. Comparing to L rather than the trace allows us to adapt dt block to quantities relevant to oscillations, since the oscillation term rotates this distribution flavor vector. After the first block, subsequent values of dt block , target are set to  \ndt block , target ← dt block , target × min ( 1 . 1 , I target dt block I dt block , target ) (21)  \nin order to drive the impact of the collisions over dt block toward the target impact. The factor of 1 . 1 prevents dt block from growing quickly following a step with serendipitously low impact. We choose a target impact of I target = 10 -4 . See Appendix B for convergence tests for the integrator accuracy, the target impact, and the number of neutrino energy bins.  \nWe use detailed flavor-incoherent neutrino interaction rates from the open-source neutrino rate library NuLib [73]. However, we must generalize these interaction rates to flavor-coherent neutrinos. This is the subject of Sec. IV.  \nFIG. 2. Two-point diagram for electron neutrino absorption and emission processes. The equivalent diagram for electron antineutrinos is the same diagram in reverse. The lack of internal neutrino lines makes this a particularly straightforward process to transform into QKE collision terms.  \n<!-- image -->', 'IV. COLLISION TERMS': 'Calculating the neutrino self-energies in Eq. 4 requires evaluating two-point irreducible Feynman diagrams starting at two loops as in [62, 65]. As in Sec. II, we neglect spin coherence and calculate Π ± ab = ˜ Π ± LLab and ¯ Π ± ab = ˜ Π ± RRab (Majorana) or ¯ Π ± ab = ˜ ¯ Π ± RRab (Dirac), though without spin coherence the Dirac and Majorana terms are identical. In the following sections, we break down the various contributions to the collision term and evaluate their effects individually. In many cases, these diagrams have already been evaluated for applications in core-collapse supernovae assuming neutrinos remain in flavor eigenstates (see, e.g., [11, 12]). We will demonstrate how to extend them to full QKE source terms without reevaluating the diagrams.', 'A. Absorption and emission': 'Let us first consider charged-current neutrino absorption by free neutrons ( ν + n → e -+ p ) and the reverse emission process (Fig. 2). Typical CCSN temperatures of around 10 MeV are not high enough to produce significant amounts of µ or τ leptons. Thus, we are able to ignore all charged-current interactions involving leptons other than electrons. This is accounted for by only considering the ee component of Π ± , such that  \nΠ + ∼ Π -∼ I ( e ) ab , (22)  \nwhere I ( c ) ab = δ ca δ cb . Expanding out Eq. 4, we find that  \nC + ∼   (1 -f ee ) -f eµ / 2 -f eτ / 2 -f µe / 2 0 0 -f τe / 2 0 0   C -∼   f ee f eµ / 2 f eτ / 2 f µe / 2 0 0 f τe / 2 0 0   . (23)  \nRecall that the contributions to the standard collision integral for the ee component in the absence of phase  \nFIG. 3. Absorption - Evolution of the flavor-off-diagonal components of a maximally mixed Fermi-Dirac neutrino distribution [Eq. 17] due to absorption and emission by nucleons and nuclei. On-diagonal components are not shown since they simply remain at their initial values. Interaction rates are based on a background described by ρ = 10 12 g cm -3 , T = 10MeV, and Y e = 0 . 3. We will show in Sec. V C that absorption dominates decoherence of higher-energy neutrinos outside of the PNS.  \n<!-- image -->  \ncoherence are  \nC + ee = j ( ν e ) (1 -f ee ) C -ee = κ ( ν e ) f ee , (24)  \nwhere j ( ν e ) ( ν, x µ ) is the electron neutrino emissivity and κ ( ν e ) ( ν, x µ ) is the electron neutrino absorption opacity, both in units of cm -1 . After matching terms and applying the same process to other lepton species, we arrive at  \nC ab = j ( ν a ) δ ab -( ⟨ j ⟩ ab + ⟨ κ ⟩ ab ) f ab , (25)  \nwhere ⟨ j ⟩ ab = ( j ( ν a ) + j ( ν b ) ) / 2, and similarly for ⟨ κ ⟩ ab . In line with the above assertion that there are very few heavy leptons present in conditions relevant to corecollapse supernovae and neutron star mergers, it is often assumed that j ( ν µ ) = j ( ν τ ) = κ ( ν µ ) = κ ( ν τ ) = 0. The collision integral for antineutrinos is exactly analogous. While for this process it is particularly straightforward to write down the collision terms, it serves to illustrate the matching procedure that we will apply to the more complex processes in the next sections.  \nThe on-diagonal collision terms are exactly identical to the standard transport collision terms, while the off-diagonal ones contain only the absorption component with interaction rates averaged between two flavors. Thus, absorption and emission will always cause flavor coherence to decay. Fig. 3 shows the evolution of the off-diagonal components of our fiducial maximally mixed Fermi-Dirac distribution [Eq. 17] relative to their initial values due only to absorption onto and emission from free nucleons and nuclei. All curves demonstrate that the off-diagonal components decay exponentially with a timescale determined by the absorption opacities according to Eq. 25. As one would expect, neutrinos at higher  \nFIG. 4. Two-point diagrams for neutrino-electron scattering processes. The charged-current processes lead to flavor decoherence and all processes allow redistribution in energy.  \n<!-- image -->  \nenergies (yellow curves) interact more strongly than those at low energies (black cuves), resulting in more rapid flavor decoherence. Similarly, neutrinos (solid lines) interact more strongly than antineutrinos (dashed lines), resulting in more rapid flavor decoherence. The ondiagonal components (not plotted) remain constant at their initial Fermi-Dirac values, since the collision terms for these components do not depend on the off-diagonal components of f .', 'B. Electron scattering': "Charged-current reactions are once again kinematically suppressed for heavy-lepton neutrinos, so we only consider scattering by electrons. The two-point diagrams contributing to electron processes are shown in Fig. 4. These give rise to four terms, each with a different flavor structure. They are [65]  \nΠ -= ∫ d 3 ν ' c 4 { A -Y L (1 -f ' ) Y L + B -Y R (1 -f ' ) Y R + D -[ Y L (1 -f ' ) Y R + Y R (1 -f ' ) Y L ] } Π + = ∫ d 3 ν ' c 4 { A + Y L f ' Y L + B + Y R f ' Y R + D + [ Y L f ' Y R + Y R f ' Y L ] } (26)  \nwhere Y R = sin 2 θ W and Y L = sin 2 θ W -1 / 2+ I ( e ) , and we use sin 2 θ W ≈ 0 . 22343 [69]. We will now use our matching procedure to find the values of the constants A ± , B ± , and D ± , These constants can depend on both the ingoing and outgoing neutrino momenta. In the absence of flavor coherence, the collision rates have long been used for neutrino transport simulations ([74] and corrections by [11]). They are usually written in the form of  \nC ( ν a ) = ∫ d 3 ν ' c 4 [ f ' ( ν a ) R + ( ν a ) (1 -f ( ν a ) ) -f ( ν a ) R -( ν a ) (1 -f ' ( ν a ) ) ] (27)  \nwhere primed quantities refer to the final state of the neutrino after scattering. R ± ( ν a ) ( ν, ν ' , cos θ, x µ ) (units of cm 3 s -1 ) are the inscattering (+) and outscattering (-) kernels for neutrino species a that depend on both the ingoing and outgoing neutrino frequencies and the cosine of the angle between the ingoing and outgoing neutrinos in a comoving orthonormal tetrad. We indicate dependence on the background fluid quantities generally with a dependence on the spacetime coordinates x µ . Detailed balance requires that R + ( ν a ) ( ν, ν ' , cos θ, x µ ) = R -( ν a ) ( ν ' , ν, cos θ, x µ ) = exp[ h ( ν ' -ν ) /kT ] R -( ν a ) ( ν, ν ' , cos θ, x µ ). When arguments are suppressed, we assume that ν ' is the second argument  \nFollowing [11], we note that the coefficients D ± are smaller than A ± and B ± by approximately a factor of the ratio of the electron mass to the electron energy, so we neglect them here. When f and ¯ f are flavor diagonal, the standard neutrino rates and the QKE collision terms should be identical, so we can simply match terms to arrive at  \nR ± ( ν e ) = ( A ± Y L Y L + B ± Y R Y R ) ee R ± ( ν µ ) = ( A ± Y L Y L + B ± Y R Y R ) µµ . (28)  \nNote that the scattering kernels for ν τ are the same as for ν µ in the absence of µ and τ leptons. We can solve for A ± and B ± . The result is  \nA ± = R ± ( ν e ) -R ± ( ν µ ) 2 sin 2 θ W B ± = (2 sin 2 θ W +1) 2 R ± ( ν µ ) -(2 sin 2 θ W -1) 2 R ± ( ν e ) 8 sin 6 θ W . (29)  \nWe can then write the self-energies very compactly as  \nΠ + ab = ∫ d 3 ν ' c 4 R + ab f ' ab Π -ab = ∫ d 3 ν ' c 4 R -ab ( δ ab -f ' ab ) , (30)  \nwhere  \nR ± ab ( ν, ν ' , cos θ ) := ⟨ R ⟩ ± ab -˜ R ± ab , ˜ R ± ab ( ν, ν ' , cos θ ) := ϵ ab R ± ( ν a ) -R ± ( ν b ) 4 sin 2 θ W , ⟨ R ⟩ ± ab ( ν, ν ' , cos θ ) := R ± ( ν a ) + R ± ( ν b ) 2 , (31)  \nand ϵ ab is the rank-two Levi-Civita symbol that makes the kernel symmetric on the flavor indices.  \nPlugging this into the collision term [Eq. 4], we get  \nC + ab = ∫ d 3 ν ' c 4 [ R + ab f ' ab -ς + ab ] C -ab = ∫ d 3 ν ' c 4 [ ⟨ R ⟩ -ab f ab -ς -ab ] , (32)  \nwhere  \nς ± ab ( ν, ν ' , cos θ, x µ ) := 1 2 ∑ c ( R ± cb f ac f ' cb + R ± ac f ' ac f cb ) (33)  \n(units of cm 3 s -1 ). The first term in each line in Eq. 32 accounts for unblocked in- and outscattering, respectively, and the second term accounts for Pauli blocking. As with absorption processes, the outscattering rate in these off-diagonal elements is based on the average of the outscattering rates of two neutrino flavors. For the flavordiagonal elements ( a = b ), the c = a terms in Eq. 33 sum to the ordinary flavor-incoherent blocking terms, but we see that there is an additional blocking or enhancement [depending on the sign of Re( f ' ac f ca )] coming from flavor coherence.  \nFor numerical applications, we expand the interaction kernels in Legendre polynomials as in [11] and consider only the first two terms. That is,  \nR ± ab ≈ 1 2 Φ ± 0 ab ( ν, ν ' , x µ ) + 3 2 Φ ± 1 ab ( ν, ν ' , x µ ) cos θ . (34)  \nInelastic scattering rates in this form are available for flavor-incoherent neutrinos as part of the opensource neutrino rate library NuLib [73]. Fig. 5 shows Φ 0 eµ / ⟨ Φ ⟩ 0 eµ over a range of incoming and outgoing neutrino frequencies at a particular temperature of 2 . 74 MeV and electron degeneracy µ e = 3 . 16 MeV chosen to clearly show the features. A value of 1 (red) means that a scattered neutrino carries over its quantum state to its new direction and energy. A value of 0 (white) means the offdiagonal flavor component is erased and that a neutrino collapses to a flavor eigenstate when scattering. Negative values (blue) mean the scattered quantum state receives a phase inversion.  \nIn the case of neutrinos scattering on electrons, the scattering kernel for a neutrino of ingoing momentum p ν and outgoing momentum p ' ν involves an integral over the ingoing electron momentum p e and outgoing electron momentum p ' e . To interpret Fig. 5, we first note that in  \nFIG. 5. Flavor structure of the isotropic part Φ ± 0 of electron scattering collision kernel R ± [Eq. 34] as used in Eq. 32. The reaction rates are calculated at an electron chemical potential of µ e = 3 . 16 MeV and a temperature of kT = 2 . 74 MeV as well as using θ W = 0 . 2223. The top plot is for neutrinos and the bottom for antineutrinos. Red regions indicate that scattering events mostly preserve the sign and magnitude of the scattered neutrino's off-diagonal elements, white regions collapse them to zero, and blue regions flip their sign.  \n<!-- image -->  \ncomputing R ( ν ) , there are three terms in this integral that come from combinations of the momentum structures of the different diagrams in Fig. 4. One term is proportional to ( p e · p ν )( p ' e · p ' ν ), another is proportional to ( p ' e · p ν )( p e · p ' ν ), and a third that we neglect (again, because it is suppressed by a factor of m e /E e ) is proportional to p ν · p ' ν (see, e.g., the integrand in Eq. C49 in [11]). If the reaction rate is dominated by the first integral, one can use the coefficients of this term for electron and muon neutrinos to show that Φ ± 0 eµ / ⟨ Φ ⟩ ± 0 eµ ≈ (4 sin 4 θ W -1) / (4 sin 4 θ W ) ≈ -0 . 67. Similarly, if the reaction is dominated by the second term, the appropriate combination of coefficients yields Φ ± 0 eµ / ⟨ Φ ⟩ ± 0 eµ = 1. As can clearly be seen in the top panel of Fig. 5, the first integral dominates for large energy transfers, causing the top left and bottom right parts of the plot to be blue. The reaction rates for antineutrinos can be  \nρ  \nFIG. 6. Diagram demonstrating how the first term in Eq. 32 (electron scattering collision integral) affects the phase of 2flavor mixed-state neutrinos after scattering when we follow the neutrino from its initial to its final energy with each scatter. We show a maximal-length single particle flavor isospin vector ⃗ ρ for an initial neutrino flavor state described by the density matrix ρ ab (unrelated to the fluid density ρ ) that is mostly electron neutrino (Initial 1) and one that is mostly muon neutrino (Initial 2) for illustration purposes ( ρ ( x ) = √ 7 / 4, ρ ( y ) = 0 and ρ ( z ) = ± 3 / 4). We set the oscillation Hamiltonian to zero and only consider electron scattering. Red regions in Fig. 5 drive the distribution function closer to the flavor axis with each scatter, as shown by the series of red arrows. Blue regions in Fig. 5 follow the same pattern, except that the neutrino phase is negated on each scatter, as shown by the series of red arrows. In both cases, the density matrix approaches the flavor states more quickly as Φ ± 0 eµ / ⟨ Φ ⟩ ± 0 eµ approaches 0. To see this, compare the Initial 1 case to the Initial 2 case, where the magnitude of the off-diagonal component of the scattering kernel is larger for the latter. The dashed circle represents the sphere along which oscillations can rotate the vector, described by ρ 2 ( x ) + ρ 2 ( y ) + ρ 2 ( z ) = 1.  \n<!-- image -->  \ndetermined by replacing p ν ↔ -p ' ν , effectively swapping the coefficients of the two energy integrals. Thus, the top left and bottom right parts of the bottom panel are red. When there are equal contributions from both, Φ 0 eµ / ⟨ Φ ⟩ 0 eµ = (8sin 4 θ W -1) / (8 sin 4 θ W +1) = -0 . 43. This is visible along the diagonal in both plots.  \nThe phase inversion process is depicted in Fig. 6 using flavor vectors for an initial single particle flavor state described by a density matrix ρ ab (unrelated to the fluid density ρ ). We define vector components of ⃗ ρ based on ρ ab in the same way that we define the distribution flavor vector components in Sec. II A. Focus first on the series of lines in the top half starting at initial state 'Initial 1'.  \nOne can interpret this initial flavor vector as representing a neutrino at a single frequency that is more likely to be measured as a neutrino than an antineutrino ( ρ ( z ) > 0) and has some flavor coherence ( | ρ ( x ) | > 0). The series of red lines depict the neutrino flavor structure after a series of scattering events where 0 < Φ ± 0 eµ / ⟨ Φ ⟩ ± 0 eµ < 1. This could be realized, for instance, by antineutrinos scattering back and forth between high and low energies (red regions in the bottom panel of Fig. 5) without blocking. During each scatter, the neutrino is not destroyed and the relative probability of finding the neutrino in an electron or muon flavor remains unchanged ( ρ ( t ) and ρ ( z ) are constant). However, since the ratio is positive but not 1, some of the flavor coherence decays away with each scattering event ( ρ ( x ) decreases). Similarly, the series of blue lines depict a scattering event where -1 < Φ ± 0 eµ / ⟨ Φ ⟩ ± 0 eµ < 0. Once again, the neutrino is not destroyed and the relative probability of being measured as each flavor remains unchanged, but the flavor coherence ( ρ ( x ) in this diagram) is negated and loses magnitude with each scattering event. In both cases, the distribution decays down to a flavor-diagonal state at the same rate. The series of lines in the bottom half start at a similar initial condition 'Initial 2' that is mostly muon neutrino for illustration purposes. The evolution of the flavor vectors proceed very similarly to those in the top half, but decay to flavor diagonal states more slowly since the magnitude of Φ ± 0 eµ / ⟨ Φ ⟩ ± 0 ab is closer to unity. In realistic QKE calculations, the impact of phase inversions on the distribution function will depend on the relative directions of the distribution flavor vectors at the initial and final momenta.  \nThe feature in the bottom left corner of Fig. 5 (below hν ∼ hν ' ∼ 10 MeV) is due to the degeneracy of electrons. Along the diagonal ( ν ≈ ν ' ), where scattering events exchange no energy, electron blocking for inscattering exactly cancels that for outscattering. However, in scattering from high energy to low energy (bottom right corner), electron blocking strongly reduces the reaction rate, whereas in scattering from low to high energies the terms that are traditionally associated with blocking effectively amplify the rate. This is to say that the feature is located precisely where one expects degeneracy effects to be important. Increasing the temperature or electron degeneracy increases the size of the feature.  \nThe leftmost and rightmost plots in Fig. 7 show the standard anisotropic component of the scattering kernel for the ee (left panels) and µµ (right panels) components of the neutrino (top panels) and antineutrino (bottom panels) distributions. As we would expect, the large positive values along the diagonal indicate that for small energy transfer, electron scattering is largely forward peaked, and the negative values away from the diagonal mean that scattering events with large energy transfer are backward peaked. However, the center plot describing the anisotropy of the scattering kernel for offdiagonal components is strikingly different because it contains both angular information similar to the plots  \nFIG. 7. Flavor structure of the linearly anisotropic part Φ ± 1 of electron scattering collision kernel R ± as used in Eq. 32. The reaction rates are calculated at an electron chemical potential of µ e = 3 . 16 MeV and a temperature of kT = 2 . 74 MeV using θ W = 0 . 2223. The top plot is for neutrinos and the bottom for antineutrinos. The anisotropic collision term for off-diagonal elements (middle panels) reflect both the anisotropy of the flavor-diagonal parts (left and right panels) and the flavor structure of Fig. 5.  \n<!-- image -->  \non the left and right and flavor information similar to Fig. 5. Thus, in the collision term for f eµ (top center panel), large energy transfers are both backward peaked and induce a phase flip, resulting in a net positive (red) value. Similarly for antineutrinos (bottom center panel), the collision term is backward peaked and does not induce a phase flip, resulting in a net negative (blue) value. For small energy transfers (along the diagonal), the scattering is forward peaked and induces a phase flip, resulting in a net negative value for both neutrinos and antineutrinos. The anisotropic (Φ 1 ) terms do not enter into the isotropic evolution equations in this paper, but will be important for quantum kinetics calculations with spatial transport.  \nGiven this intuition, we can now understand the action of the inelastic scattering kernel on our fiducial maximally mixed Fermi-Dirac distribution [Eq. 17]. The top and middle panels of Fig. 8 depict evolution of the flavordiagonal components for both neutrinos and antineutrinos. The first thing to notice is that even though the flavor-diagonal components are in Fermi-Dirac distributions, the flavor coherence temporarily drives them out  \n̸  \nof equilibrium. Inspecting Eq. 32, it becomes apparent that the only terms that modify the flavor-diagonal components are the a = b = c part of the blocking terms ς ± ab . This can be seen after noting that the flavor-diagonal components of the scattering kernel are R aa = ⟨ R aa ⟩ = R ( ν a ) and that the first terms in Eq. 32 combined with the a = b = c part of the ς ± ab terms constitute the standard flavor-diagonal collision term for incoherent neutrinos. The a = b = c part of ς ± ab is the blocking term for incoherent neutrinos, and ς ± ab itself is just the extension of the blocking term to non-flavordiagonal neutrino distributions. Eq. 12 shows that the blocking terms, which are quadratic in f , have an effect smaller than the standard flavor-incoherent terms when f and f ' are antialigned and a larger effect when aligned. Thus, what we see in Fig. 8 is that the Fermi-Dirac distribution, though an equilibrium for incoherent neutrinos, is not an equilibrium in this mixed state. The blocking terms redistribute neutrinos in energy toward a new mixed equilibrium that has more low-energy neutrinos and fewer high-energy neutrinos. However, after a short period of time the terms linear in f (i.e., the regular  \nFIG. 8. Inelastic electron scattering - Evolution of maximally mixed Fermi-Dirac neutrino distribution [Eq. 17] due to inelastic electron scattering ( ν + e -↔ ν + e -) interactions. The top panel contains the flavor-diagonal neutrino distribution components, the middle panel constains the flavordiagonal antineutrino distribution components, and the bottom panel contains the flavor-off-diagonal neutrino and antineutrino components. Interaction rates are based on a background described by ρ = 10 12 g cm -3 , T = 10MeV, and Y e = 0 . 3. The charged-current parts of the electron scattering interaction drive the off-diagonal elements to zero, and the diagonal elements redistribute in energy to non-FermiDirac values as long as the off-diagonal elements are nonzero.  \n<!-- image -->  \nscattering terms without blocking) drive the distribution to be flavor-diagonal, since R + ab is present in the inscattering terms but ⟨ R ⟩ -ab is present in the outscattering terms. This is demonstrated by the sharp decline of the off-diagonal components of both neutrino and antineutrino distributions in the bottom panel of Fig. 8. As the flavor-off-diagonal components decay, the blocking terms increasingly resemble the classic flavor-incoherent terms and the distribution settles back down to a flavordiagonal Fermi-Dirac distribution. As expected, neutrinos interact more strongly than antineutrinos and highenergy neutrinos interact more strongly than low-energy neutrinos, and stronger interaction rates lead to quicker returns to the flavor-diagonal Fermi-Dirac distribution.  \nMany modern neutrino transport codes treat electron scattering as elastic for computational efficiency. If the collisions are treated as elastic,the scattering kernels simplify as R + ab = R -ab . The collision integral greatly reduces to  \nC ab = 1 4 π ∫ d Ω ' [( κ 0 ab + κ 1 ab cos θ ) f ' ab -( ⟨ κ ⟩ 0 ab + ⟨ κ ⟩ 1 ab cos θ ) f ab ] , (35)  \nFIG. 9. Elastic electron scattering - Evolution of maximally mixed Fermi-Dirac neutrino distribution [Eq. 17] due to elastic electron scattering ( ν + e -↔ ν + e -) interactions. The top panel contains the flavor-off-diagonal neutrino and antineutrino components. The bottom panel shows the difference between the elastic ( f ) and inelastic ( f inelastic from Fig. 8) results. On-diagonal components are not shown since they simply remain at their initial values. Interaction rates are based on a background described by ρ = 10 12 g cm -3 , T = 10MeV, and Y e = 0 . 3. Though these approximate results are qualitatively similar to the bottom panel of Fig. 8, there are significant quantitative differences.  \n<!-- image -->  \nwhere the elastic but anisotropic scattering opacity is  \ndκ ( ν a ) d Ω ' = 1 4 π ( κ 0( ν a ) + κ 1( ν a ) cos θ ) (36)  \nand the matrix form of the opacities are constructed in the same way as the scattering kernels in Eq. 31.  \nThe top panel of Fig. 9 shows the evolution of the flavor off-diagonal parts of a maximally mixed FermiDirac distribution under this assumption of elastic electron scattering. The on-diagonal components f aa are not shown, since the assumption of elastic scattering causes ς ± ab (blocking) terms to cancel and the on-diagonal components do not change. We see that the elastic treatment produces results rather consistent with the full inelastic treatment. The bottom panel shows the difference between the elastic and inelastic treatments, and in this particular choice of fluid parameters elastic scattering causes the distribution to evolve more slowly, leading to differences of at most ∼ 15%. Of course, greater differences could occur for other fluid and neutrino distributions.  \nFIG. 10. Two-point diagram for neutrino-nucleon scattering.This is very similar to Fig. 4, but the lack of chargedcurrent processes makes the flavor structure less complex.  \n<!-- image -->", 'C. Nucleon scattering': "Next, we consider scattering by nucleons and nuclei (Fig. 10), which at the two-loop level only undergo scattering reactions with neutrinos through neutral-current processes. The outscattering and inscattering selfenergies are identical in form to Eq. 30, though we note that due to the neutral current nature of nucelon scattering that the flavor structure of the scattering kernels is simpler ( R ± ab = ⟨ R ⟩ ± ab and ˜ R ± ab = 0). Thus, the scattering kernels do not themselves impart any flavor information on the collision term and could be represented by a scalar instead of a flavor matrix (i.e., R ± ab = R ± ). However, in order to minimize the number of representations used in this paper we leave in the flavor indices. Note that weak magnetism corrections cause the interactions to be weaker for antineutrinos than for neutrinos [75], but leave the flavor structure is unaffected.  \nThe resulting self-energies are then identical in form to Eq. 30. The contribution to the collision term is then also given by Eq. 32. Matching the diagonal elements with f ex = f ' ex = 0 to the known scattering interaction rates for incoherent neutrinos, we see that R ± aa are the standard scattering kernels. Though inelastic nucleon scattering is not yet implemented in NuLib , we note that in the electron scattering case the finite value of ˜ R ± ab is what drives the distribution to a flavor equilibrium and leads to interesting structure in Figs. 5 and the center panels of Fig. 7. When ˜ R ± ab = 0 as it is for nucleon scattering the ratio Φ ± ab / ⟨ Φ ⟩ ± ab is everywhere identically 1. Thus, inelastic nucleon scattering will redistribute neutrinos and flavor coherence in energy and direction, but will not cause the flavor coherence to decay.  \nIf we consider the collisions to be elastic, which is a much better approximation for nucleon or nucleus scattering than electron scattering because the masses are much larger, then R + = R -. Once again, the corresponding collision term is identical to that for electron scattering given by Eq. 35, though we note that for neutral-current scattering κ ab = ⟨ κ ⟩ ab and ˜ κ ab = 0. For an isotropic radiation field and elastic nucleon scattering, f = f ' and C = 0. Thus, we do not plot the uninteresting case of evolving an isotropic mixed distribution function in the presence of elastic nucleon scattering because all quantities remain constant.", 'D. Pair annihilation': "Electron-positron pair annihilation ( e -+ e + ↔ ν + ¯ ν ) is a cross diagram of electron scattering and the collision terms have a similar structure. That is [65],  \nΠ + = ∫ d 3 ν ' c 4 { A + Y L (1 -¯ f ' ) Y L + B + Y R (1 -¯ f ' ) Y R + C + [ Y L (1 -¯ f ' ) Y R + Y R (1 -¯ f ' ) Y L ] } , Π -= ∫ d 3 ν ' c 4 { A -Y L ¯ f ' Y L + B -Y R ¯ f ' Y R + C -[ Y L ¯ f ' Y R + Y R ¯ f ' Y L ] } (37)  \nwhere ¯ f is the antineutrino distribution function. Once again, we neglect the third term in each line. One can go through the term-matching procedure in Sec. IV B and show that A ± and B ± are the same as in the electron scattering case [Eq. 29], but with values of R ± ( ν a ) from annihilation processes. This leads us to  \nΠ + ab = ∫ d 3 ν ' c 4 R + ab ( δ ab -¯ f ' ab ) Π -ab = ∫ d 3 ν ' c 4 R -ab ¯ f ' ab (38)  \nPlugging this into the collision term [Eq. 4], we arrive at  \nC + ab = ∫ d 3 ¯ ν ' c 4 [ R + ab δ ab -⟨ R ⟩ + ab f ab -R + ab ¯ f ' ab + ς + ab ] C -ab = ∫ d 3 ¯ ν ' c 4 ς -ab , (39)  \nwhere R ± ab and ⟨ R ⟩ ± ab are defined with respect to values of R ± ( ν a ) from annihilation processes in the same way as in Eq. 31. The ς ± ab is also similar to the electron scattering case [Eq. 33], but replacing f ' → ¯ f ' since we are integrating over the antineutrino distribution instead of an outgoing neutrino distribution. The first term on the first line of Eq. 39 is the emission term in the absence of Fermi blocking, which creates flavor-diagonal pairs of neutrinos. Similarly, the second line describes the neutrino pair annihilation rate, where Fermi blocking of final state electrons and positrons is already accounted for in the calculation of R -ab . We can then understand the final three terms on the first line as accounting for neutrino Fermi blocking. The second term in C + ab behaves like outscattering and will always decrease the magnitude of the flavor-off-diagonal components. The third term has flavor structure like the inscattering term in Eq. 32, but acts with the opposite sign.  \nFig. 11 shows the isotropic part [see Eq. 34] of the pair production kernel R ± eµ , similar to Fig. 5. This term does not affect the unblocked pair production or annihilation rates, but does describe how the presence of neutrinos and antineutrinos contribute to the blocking term differently and how they annihilate. When the ratio Φ ± eµ / ⟨ Φ ⟩ ± eµ is negative (blue in the top panel of the  \nFIG. 11. Flavor structure of the isotropic part Φ ± 0 of the pair annihilation kernel R ± as used in Eq. 39. The reaction rates are calculated at an electron chemical potential of µ e = 17 . 8 MeV and a temperature of kT = 2 . 74 MeV using θ W = 0 . 2223. The top plot is for neutrinos and the bottom for antineutrinos. Red and dark blue regions show where electron degeneracy has a significant impact on the flavor structure of the scattering kernel.  \n<!-- image -->  \nfigure), the flavor off-diagonal component of the blocking term will increase the neutrino distribution in the direction of the antineutrino distribution's off-diagonal component. Similarly, when the ratio is positive (red in the top panel of the figure), the neutrino's off-diagonal components will be increased in the opposite direction of the antineutrino's off-diagonal components. The same argument can be applied to antineutrinos being blocked by neutrinos (bottom panel of the figure). Notice that the antineutrino collision term is simply the transpose in ingoing/outgoing energy of the neutrino collision term. This is due to the fact that the interaction diagrams are closely related to the interaction diagrams for antineutrinos, so R ± ab ( ν, ¯ ν, cos θ, x µ ) = ¯ R ± ab (¯ ν, ν, cos θ, x µ ).  \nSimilarly, Fig. 12 shows the linearly anisotropic part of the electron-positron pair annihilation kernels. The left and right panels show the classic incoherent kernels for electron neutrinos (left) and heavy lepton neutrinos (right) and for antineutrinos (bottom left and right). The  \npredominantly blue color means that neutrinos annihilate more strongly when colliding at large angles, except for a small range of energies near the electron degeneracy energy where dependence on the collision angle is much weaker. This feature is only present when the electron degeneracy is larger than the fluid temperature. The middle plots for neutrinos (top) and antineutrinos (bottom) show the off-diagonal part of the kernel, which is a convolution of angular information in the plots to the left and right and the flavor information in Fig. 11. In the blue regions, the blocking term from neutrino-antineutrino pair creation will increase the neutrino's off-diagonal components in the same (flavor) direction as the antineutrino's if they are born moving in the same (spatial) direction, and will do the opposite if born moving in opposite (spatial) directions. In red regions, the opposite is true.  \nNote that the features in Figs. 11 and 12 only appear when the electrons are very degenerate. We chose very degenerate parameters to show these effects, but when kT > ∼ µ e there is almost no dependence of the flavor structure on either neutrino or antineutrino energy (the plots are a solid color). There are two relevant energy scales. The first is the inner edge of the features at around 10 MeV in the figures. This is related to the temperature, since in the range of 0 < h ( ν + ν ' ) < ∼ 10 MeV a significant fraction of the annihilating electrons will have energies below our example fluid temperature of 2 . 74 MeV. The second is the outer edge of the features at around 100 MeV in the figures. This energy scale is related to the electron degeneracy, since in the range of 10 < ∼ h ( ν + ν ' ) < ∼ 100 MeV a significant fraction of the reacting electrons will be degenerate.  \nFig. 13 shows the evolution of our isotropic maximally mixed Fermi-Dirac distribution of neutrinos [Eq. 17] due to electron-positron pair production and annihilation only. Since the distribution is isotropic, only the isotropic part of the kernel Φ 0 enters into the calculation. Similar to the electron scattering case, we see that the pair production kernels modify both the flavor-diagonal (top two panels of Fig. 13) and the flavor-off-diagonal (bottom panel of Fig. 13) components of the distribution. Once again, only the nonlinear ς ± ab terms lead to the initial change in the on-diagonal components, since the other terms cancel out by construction of the initial FermiDirac distribution. Unlike electron scattering, electronpositron pair production and annihilation can change the total number of neutrinos present, which occurs as all flavor-diagonal components grow in the first ∼ 0 . 1 ms. The ee component of both the neutrino and antineutrino distributions (solid lines in the top two panels) grow on a shorter timescale than the µµ components (dashed lines in the top two panels), since the collision rates for those components are larger due to the contribution of chargedcurrent processes. While this rapid change occurs, the off-diagonal components (bottom panel) of both neutrino and antineutrino distributions decline, and decline more rapidly for high-energy anti/neutrinos.  \nFor all components, the antineutrino distribution  \nFIG. 12. Flavor structure of the linearly anisotropic part Φ ± 1 of the pair annihilation kernel R ± as used in Eqs. 39. The reaction rates are calculated at an electron chemical potential of µ e = 17 . 8 MeV and a temperature of kT = 2 . 74 MeV using θ W = 0 . 2223. The top plot is for neutrinos and the bottom for antineutrinos. The anisotropic annihilation term for off-diagonal elements (middle panels) reflects both the anisotropy of the flavor-diagonal parts (left and right panels) and the flavor structure of Fig. 11.  \n<!-- image -->  \nevolves more quickly than the neutrino distribution. The initial neutrino distribution has a higher average energy and a higher overall number density than the initial antineutrino distribution due to the large degeneracy of the electrons. This causes both R ± ab and ¯ R ± ab to be larger when ν > ¯ ν than when ν < ¯ ν . If we look at the third term in Eq. 39, since the average energy of ¯ f is lower than that of f and R + ab is smaller when ¯ ν < ν , the contribution of that term to the integral is relatively small. Similarly, the corresponding term for antineutrinos is an integral over f , which has a larger average energy than ¯ f , and ¯ R + , which is larger when ν ' > ¯ ν . Thus, it is natural to expect a larger collision term for antineutrinos than neutrinos and correspondingly more rapid evolution of the antineutrino distribution.  \nIt is common for neutrino transport codes to treat pair processes as an effective absorption and emission process by applying Kirchoff's law to the pair emissivity under the assumption of no Fermi blocking. This allows the neutrino distributions to reach the correct equilibrium and is less computationally expensive. To extend this idea to coherent flavor transport, we take the flavor struc-  \nture of the term linear in f .  \nC ab = j ( ν a ) δ ab -( ⟨ j ⟩ ab + ⟨ κ ⟩ ab ) f ab (40)  \nwhere  \nj ( ν a ) = ∫ d 3 ν ' Φ + 0( ν a ) κ ( ν a ) = j ( ν a ) 1 -FD ( T, µ ν a , ν ) . (41)  \nFig. 14 shows the evolution of the same initial condition under this treatment of pair processes. The on-diagonal components are in equilibrium and do not change due to the lack of the nonlinear ς ± ab terms in Eq. 40, and so are not plotted. The top panel shows the evolution of the off-diagonal components of the neutrino (solid lines) and antineutrino (dashed lines) distributions relative to their initial values. We see that they evolve on a timescale that is significantly longer than that seen in the bottom panel of Fig. 13. When we subtract the full kernel solution from this approximated solution (bottom panel of Fig. 14), we see that the solutions differ by as  \nFIG. 13. Pair Annihilation - Evolution of maximally mixed Fermi-Dirac neutrino distribution [Eq. 17] due to electronpositron pair ( e + + e -↔ ν + ¯ ν ) interactions. The top panel contains the flavor-diagonal neutrino distribution components, the middle panel contains the flavor-diagonal antineutrino distribution components, and the bottom panel contains the flavor-off-diagonal neutrino and antineutrino components. Interaction rates are based on a background described by ρ = 10 12 g cm -3 , T = 10MeV, and Y e = 0 . 3. The nonlinear terms cause the flavor-diagonal elements to redistribute away from Fermi-Dirac values as long as the offdiagonal elements are nonzero. The long evolution timescale makes this process a subdominant driver of flavor decoherence.  \n<!-- image -->  \nmuch as ∼ 80%. The effective absorption approach is problematic, but we will see in future sections that pair annihilation is not a dominant source of flavor decoherence for these background fluid parameters.", 'E. Nucleon-nucleon bremsstrahlung': "Comparing nucleon-nucleon bremsstrahlung to electron-positron pair annihilation is very similar to comparing nucleon scattering to electron scattering. Nucleon-nucleon bremsstrahlung is a neutral-current process that requires too many diagrams to show them all here. The general structure of the diagrams is depicted in Fig. 15, but without the gauge bosons connecting the lines in many permutations (see, e.g., [76, 77] and references therein). However, in all diagrams there is a single unbroken neutrino line, so the self-energy is linear in f ' . The resulting expressions for the self-energies are identical in form to Eq. 38, except that the purely neutral-current nature of the interactions leads to the simplifications R ± ab = ⟨ R ⟩ ± ab and ˜ R ± ab = 0. The resulting collision term is then identical to Eq. 39.  \nHowever, many codes treat nucleon-nucleon  \nFIG. 14. Effective absorption pair annihilation - Evolution of maximally mixed Fermi-Dirac neutrino distribution [Eq. 17] due to electron-positron pair ( e + + e -↔ ν + ¯ ν ) interactions treated as an effective absorption/emission process. The top panel contains the flavor-off-diagonal neutrino and antineutrino components. The bottom panel shows the difference between the elastic ( f ) and inelastic ( f inelastic from Fig. 8) results. On-diagonal components are not shown since they simply remain at their initial values. Interaction rates are based on a background described by ρ = 10 12 g cm -3 , T = 10MeV, and Y e = 0 . 3. There are significant differences between how an effective absorption treatment leads to flavor decoherence compared to the full treatment.  \n<!-- image -->  \nFIG. 15. Two-point diagram for nucleon-nucleon bremsstrahlung radiation. This diagram represents many diagrams with various permutations of Z bosons connecting the lines. A single internal neutrino line makes the structure similar to electron-positron pair annihilation (Fig. 4).  \n<!-- image -->  \nbremsstrahlung as an effective absorption and emission process by applying Kirchoff's law to the pair emissivity under the assumption of no Fermi blocking. This leads to reaction rates that are approximately correct, settle to the correct equilibrium, and are much more computationally efficient. To extend this idea to coherent flavor transport, we take the flavor structure of the term linear in f , which results in a collision term as in Eq. 40. With this construction, all species have the same emissivity but different absorption opacities such that the collision term drives the distribution to Fermi-Dirac values.  \nFIG. 16. Effective absorption bremsstrahlung -Evolution of the flavor-off-diagonal components of a maximally mixed Fermi-Dirac neutrino distribution [Eq. 17] due to nucleonnucleon bremsstrahlung ( N + N ↔ N + N + ν +¯ ν ) interactions. On-diagonal components are not shown since they simply remain at their initial values. Interaction rates are based on a background described by ρ = 10 12 g cm -3 , T = 10MeV, and Y e = 0 . 3. Bremsstrahlung processes decohere low-energy neutrinos more rapidly than high-energy neutrinos and become a dominant driver of decoherence at nuclear densities. The significant errors associated with the effective absorption treatment of pair processes (Fig. 14) calls for a more realistic treatment of bremsstrahlung processes.  \n<!-- image -->  \nFig. 16 shows the evolution of our maximally mixed Fermi-Dirac distribution [Eq. 17] using an approximate absorption-emission treatment, since full kernels are not readily available. As with any simple absorption process (e.g., Sec. IV A), the off-diagonal component decays exponentially. However, we see that high-energy neutrinos evolve more slowly than low-energy neutrinos, since the bremsstrahlung emissivity decreases with neutrino energy (e.g., [12]). The neutrino distribution at a given energy also evolves more slowly than the antineutrino distribution because the values of FD(¯ ν, -µ ν e , T ) that go into Eq. 41 are smaller than FD( ν, µ ν e , T ) at the same neutrino energy. In any case, the effect of bremsstrahlung collision processes on neutrino flavor coherence is much weaker than other processes for this choice of fluid parameters. However, the effects can be much more conspicuous at the nuclear densities in a protoneutron star.", 'F. Four-neutrino processes': "Finally, we describe the contribution of four-neutrino processes to the collision term, the diagrams for which are shown in Fig. 17. The full integrals [65] are rather complicated and we do not reproduce them here. Since this reaction requires integrating over all of the neutrino distributions, we do not extend existing reaction rates and instead directly integrate Eq. 96 in [65].  \nFig. 18 shows the evolution of our fiducial maximally mixed Fermi-Dirac distribution [Eq. 17] due to neutrinoneutrino scattering and pair processes. It is apparent  \nν  \nFIG. 17. Two-point diagrams for neutrino-neutrino scattering and pair annihilation. The multiple internal neutrino lines make this a difficult process to treat in the QKEs. The lack of charged-current interactions precludes phase decoherence.  \n<!-- image -->  \nFIG. 18. Four-neutrino processes - Evolution of maximally mixed Fermi-Dirac neutrino distribution [Eq. 17] due to neutrino-neutrino scattering ( ν + ν ↔ ν + ν ) and pair ( ν +¯ ν ↔ ν +¯ ν ) interactions. The top panel contains the flavor-diagonal neutrino distribution components, the middle panel contains the flavor-diagonal antineutrino distribution components, and the bottom panel contains the flavor-off-diagonal neutrino and antineutrino components. Interaction rates are based on a background described by ρ = 10 12 g cm -3 , T = 10MeV, and Y e = 0 . 3. The neutrinos never decay to flavor-diagonal FermiDirac distributions, but instead redistribute such that the distribution flavor vector at all neutrino energies are aligned, and the neutrino vectors are antialigned with the antineutrino vectors.  \n<!-- image -->  \nthat the initial distribution is not an equilibrium distribution, and the distribution redistributes itself. Initially, the distribution flavor vector length is 1.6-2.7 (depending on the neutrino energy) times the length of the flavor vector of a flavor-diagonal Fermi-Dirac distribution due to our imposed off-diagonal components [Eq. 17]. However, over time the flavor vector length decays down to  \na nearly constant factor of 1.1-1.4 times the Fermi-Dirac values (again, depending on the neutrino energy). This is done by redistributing flavor coherence between neutrinos and helicities. Indeed, the initial flavor phase angle (the angle between the distribution flavor vector and the positive flavor axis) is initially energy and helicity dependent, but at the end of the calculation the distribution flavor vector is constant in energy at 50 . 6 · and the antineutrino distribution flavor vector at 230 . 6 · . Intuitively, neutrino-neutrino scattering will redistribute flavor phase within a helicity to make the flavor phase angle constant in energy and neutrino-neutrino pair processes will redistribute flavor phase across helicity to make the neutrino and antineutrino flavor phase angles antiparallel. The direction of this equilibrium distribution flavor vector can be in any direction and depends on the initial conditions.  \nIf we assume that distribution functions 2 and 3 are flavor-diagonal Fermi-Dirac distributions (similar to what is done by [78] for standard neutrino transport), the self-energies reduce to  \nΠ + ab = ∫ d 3 ν ' 1 c 4 ⟨ R ⟩ + ab f ' 1 ab Π -ab = ∫ d 3 ν ' 1 c 4 ⟨ R ⟩ -ab ( δ ab -f ' 1 ab ) (42)  \nwhere  \nR + ( ν a ) = ∑ c (1 + δ ac ) ∫ d 3 ν ' 2 c 3 d 3 ν ' 3 c 3 × r ( p 1 + p 3 → p + p 2 ) (1 -f ' 2 cc ) f ' 3 cc R -( ν a ) = ∑ c (1 + δ ac ) ∫ d 3 ν ' 2 c 3 d 3 ν ' 3 c 3 × r ( p + p 2 → p 1 + p 3 ) f ' 2 cc (1 -f ' 3 cc ) (43)  \nare the ordinary Boltzmann neutrino-neutrino scattering rates. This would lead to similar evolution timescales, but assuming that one of the reacting distributions is flavor diagonal breaks the flavor symmetry that allows the full scattering and annihilation kernels to relax the distributions to an arbitrary flavor phase angle. Instead, the neutrino and antineutrino distributions would always be driven to flavor-diagonal Fermi-Dirac distributions.", 'V. RESULTS': 'Armed with an understanding of the effects of individual processes, we can make sense of more complete QKE simulations. We will first outline the major points before taking a deeper dive into the results. We have published the most important data from this section on Zenodo 3 , and the rest is available on request.  \nIn Sec. VA we will combine all of the interactions in Sec. IV while still neglecting oscillations. We demonstrate for our fiducial fluid parameters of ρ = 10 10 g cm -3 , T = 10MeV, and Y e = 0 . 3 that neutrino flavor coherence decays on timescales near a microsecond depending on energy, as seen in the bottom panel of Fig. 19. Bremsstrahlung processes drive the decoherence of the lowest energy neutrinos, absorption and electron scattering drive decoherence of the highest energies, and inelastic scattering and pair processes drive changes in the flavor-diagonal components.  \nThe keystone to our discussion of the QKEs is the oscillation term, which we discuss in Sec. V B. We show in Fig. 20 that for our fiducial fluid parameters there is a nutation due to the total neutrino density (top two panels) and a precession due to the electron density, both of which occur on timescales of around a picosecond. Extending this calculation by several microseconds, we see the combined action of oscillations and collisions in Fig. 21. The neutrino distributions (thick shaded regions covering the oscillation amplitude) decohere to flavordiagonal Fermi-Dirac distributions on timescales similar to but measurably different from the calculations without oscillations (green curves).  \nGiven this result, we feel justified to do a wider parameter sweep in Sec. V C using the less-expensive nonoscillating QKE calculations to understand flavor decoherence in core-collapse supernovae. We take background matter parameters from a one-dimensional (1D) CCSN simulation (Fig. 22) and determine the decoherence time using an isotropic QKE calculation at each radial point (Fig. 23). We show the results using different interaction sets in different panels to demonstrate that in the PNS decoherence is dominated by electron scattering, nucleonnucleon bremsstrahlung, and neutrino-neutrino scattering, while near the shock decoherence is dominated by absorption, and in the decoupling region all of these processes except elastic nucleon scattering and neutrinoneutrino processes have a significant impact. Finally we are able to demonstrate that an empirically determined effective decoherence opacity κ effective [Eq. 45] is a much better predictor of flavor decoherence rates than the mean free path. In Fig. 24, one can see that the decoherence timescale from κ effective matches the simulated one generally to within a factor of 10 in the PNS and to within 20% outside of the PNS.', 'A. Quantum kinetics without oscillations': 'First, we will discuss how neutrinos evolve with a full suite of collision processes but without oscillations to provide a basis for understanding what happens when evolving the full QKEs. We begin with the same initial maximally mixed neutrino distribution [Eq. 17] and background matter parameters ( ρ = 10 12 g cm -3 , T = 10 MeV, Y e = 0 . 3) used in Sec. IV. The resulting evolution of the neutrino distributions is shown in Fig. 19, and  \nFIG. 19. All processes -Evolution of maximally mixed Fermi-Dirac neutrino distribution [Eq. 17] due to absorption/emission, inelastic electron scattering, electron-positron pair annihilation, nucleon-nucleon bremsstrahlung, neutrinoneutrino scattering, and neutrino-antineutrino pair interactions. The top panel contains the flavor-diagonal neutrino distribution components, the middle panel constrains the flavordiagonal antineutrino distribution components, and the bottom panel contains the flavor-off diagonal neutrino and antineutrino components. Interaction rates are based on a background described by ρ = 10 12 g cm -3 , T = 10MeV, and Y e = 0 . 3. This quantitatively provides an estimate of decoherence timescales on the order of a millisecond for these parameters depending on energy.  \n<!-- image -->  \nit is immediately clear that it behaves unlike that due to any single collision process.  \nThe initial reactions of the flavor-diagonal elements (top two panels of Fig. 19) are very reminiscent of the evolution due to neutrino-neutrino reactions (Fig. 18), though with different amplitudes and on a shorter timescale. In all but the lowest energy electron neutrinos (solid lines, top panel of Fig. 19), the impact of electron scattering (Fig. 8) is apparent. Similarly, at intermediate energies, the electron and muon antineutrino curves (center panel of Fig. 19 evolve similar to those caused by pair annihilation in Fig. 13. The offdiagonal elements (bottom panel) show that high-energy neutrinos lose flavor coherence more quickly than lowenergy neutrinos and that neutrinos lose flavor coherence more quickly than antineutrinos. This reflects the actions of absorption of neutrinos on nucleons (Fig. 3) and electron scattering (Fig. 8). Even though neutrinoneutrino interactions themselves support long-lived flavor coherence, when combined with coherence-destroying reactions they simply accelerate the demise of coherence. The rapid initial adjustment from neutrino-neutrino interactions quickly shuffles neutrinos between energy bins, and energy is also shuffled between neutrinos and an-  \nFIG. 20. Oscillation terms - Evolution of all components of a two-flavor neutrino distribution function due to the action of the oscillation Hamiltonian in Eq. 3 in a background fluid described by ρ = 10 12 g cm -3 , T = 10MeV, and Y e = 0 . 3. There is a short precession timescale associated with the number density of electrons visible in the flavor-off-diagonal components (bottom two panels) and a longer nutation timescale associated with the total neutrino density most obvious in the flavor-diagonal components (top two panels). Both are much shorter than the collision timescale.  \n<!-- image -->  \ntineutrinos, allowing flavor-destructive processes to operate more efficiently.', 'B. Quantum kinetics with oscillations': 'Now we can finally add the oscillation term back in. Before discussing the combined effects of oscillations and collisions, we note that the oscillations occur on an extremely short timescale of ∼ 10 -16 s. Fig. 20 shows this in detail. The top two panels show the evolution of the flavor-diagonal components of the neutrino (solid) and antineutrino (dashed) distributions. When one flavor becomes more abundant, the other flavor becomes less abundant, conserving the total number of neutrinos. Lower energy neutrinos (dark colors) are more degenerate and have larger occupation probabilities. The lowestenergy electron neutrinos (black solid curves at the top of the top panel) actually oscillate into a completely degenerate state with occupation numbers of 1. This empha-  \nzes the importance of taking care to choose only physically motivated initial conditions since a larger amount of initial flavor mixing in the initial conditions than used here would allow the oscillation term to push the occupation numbers in these degenerate states higher than 1. With our choice of initial conditions, the blocking terms in our calculations naturally prevent the neutrino distribution from developing into a state with f > 1 in any basis. For the realization of the QKEs [64] and our adaption of them, we have not, however, proven that the collision terms are incapable of developing such an unphysical distribution under all circumstances. This is trivially true for absorption and elastic scattering, for which evolution of the diagonal components is unchanged by flavor coherence, but the nonlinear processes are less straightforward.  \nThe bottom two panels of Fig. 20 show the real and imaginary parts of the distribution, respectively. First, we notice that the real and imaginary parts are out of phase with each other, indicating a circular oscillation of the distribution flavor vector around the flavor axis. Also, the direction of the antineutrino flavor vector oscillation is opposite to the the neutrino vector, indicated by the opposite sign of the imaginary components. This is expected based on the relationship between the neutrino and antineutrino matter Hamiltonians (Sec. II). Finally, we see that when the difference between the values of the on-diagonal components (the z component of the distribution flavor vector) grows toward t = 0 . 3 ps, the magnitude of the off-diagonal components ( x and y components of the distribution flavor vector) shrink. This is an expression of the fact that the oscillation term is unable to change the length of the distribution flavor vector. Effects from interactions are technically present in these data, but they do not have have significant impact over such a short time.  \nThe rapid oscillations in the bottom two panels show the effect of the matter potential. At the relevant density of 10 12 g cm -3 and electron fraction of 0.3, the matter potential [Eq. 6] is 0 . 023 eV, a factor of 1 . 8 × 10 7 larger than the vacuum potential for the 1 MeV bin. Because of the large matter potential, the oscillation axis is nearly aligned with the flavor axis. One can also see that the period of the oscillations is approximately h/ √ 2 G F n e = 1 . 8 × 10 -13 s. Deviations from this are due to the additional contribution of the neutrino selfinteraction potential. For the diagonal components in the top panel, the period of oscillation, also due to the neutrino potential, is h/ √ 2 G F n ν = 6 . 9 × 10 -13 s, where n ν is the total neutrino number density.  \nOnly on timescales much longer than the oscillation time do collisions have a chance to impact the neutrino distribution, so we continue the simulation for an additional 100 µ s. The top left panel of Fig. 21 makes it clear that the amplitude of the nutations at neutrino energies of 8 MeV are damped on a timescale of around a microsecond and that the distributions are driven to flavor states on a timescale of several microseconds. In particular, this panel shows the z component of the neutrino  \n(dark) and antineutrino (light) distribution flavor vectors relative to their equilibrium Fermi-Dirac values. In flavor-diagonal equilibrium, the distribution flavor vector is entirely in the z direction with a magnitude of L eq . As a reminder, the initial conditions were constructed as distribution flavor vectors described by f ( z ) (0) = L eq and an imposed large x component, resulting in a value of L (0) > L eq . On the left side of the plot, the oscillations from the top two panels of Fig. 20 are visible as a vertical spread in the colored regions. The neutrinos oscillate from their initial latitude (bottom edge of the dark region) toward the electron flavor axis (top edge of the dark region), indicated by a z component that grows to values of f ( z ) /L eq > 1. Antineutrinos start with a negative latitude (more muon antineutrinos that electron antineutrinos, bottom edge of the light region) and oscillate across the equator a bit beyond the opposite latitude (top edge of the light region). The oscillation term drives the flavor-diagonal components much farther from their initial values than the collision terms themselves were able to do in the absence of oscillations (plotted in solid green for neutrinos and dashed green for antineutrinos). As time progresses, two effects are evident. First, the amplitude of the oscillations in the z direction decreases as the shaded areas compress to lines. Second, the location of the average value of f ( z ) decays to the equilibrium value. Stated another way, the distribution flavor vectors nutate within a band between two latitudes, and with time the width of the band (amplitude of the nutations) shrinks and the average latitude approaches the (muon) electron flavor axis for (anti)neutrinos.  \nThe top left plot only reflects the behavior of the flavordiagonal components of the distribution. However, in the leftmost center plot in Fig. 21 we see that oscillations in the flavor off-diagonal components are damped similarly to the flavor-diagonal components and that they damp on roughly the timescale predicted by the nonoscillating calculations. This can be seen as follows. To avoid overloading the plot with the precession motion, we only plot the part of the distribution flavor vector perpendicular to the flavor axis  \nf ( ⊥ ) := √ f 2 ( x ) + f 2 ( y ) , (44)  \nstill at 8 MeV. Though this quantity also oscillates at the left-hand side of the plot, the neutrinos (dark) exhibit a much broader band than antineutrinos (light). This corroborates the interpretation of the top panel, namely that the neutrino distribution flavor vector is oscillating from its initial latitude (large f ( ⊥ ) ) to very near the electron flavor axis (small f ( ⊥ ) ) and back. Antineutrinos, on the other hand, oscillate from their initial latitudes across the flavor equator without getting very close to either axis (opposite sign f ( z ) and small variation in f ( ⊥ ) ). As with the z component, the collisions simultaneously decrease the average value of f ( ⊥ ) and the range of the oscillation. At t ≈ 0 . 4 µ s, the lower edge of the f ( ⊥ ) band is able to reach a value of zero, meaning the neutrino distribution flavor vector is able to nutate directly in line with  \nFIG. 21. Top panels : Evolution of the z component of the neutrino distribution flavor vectors relative to the equilibrium vector length at several neutrino energies. Dark (light) shades show (anti)neutrino values. The equivalent values from the nonoscillating calculations are plotted in green for neutrinos (solid lines) and antineutrinos (dashed lines). The spread is due to violent oscillations on much shorter timescales, and one can see collisions decreasing the amplitude of the oscillations (width of the spread) and driving the distributions toward the flavor axis. Center panels : Component of the distribution flavor vector perpendicular to the z direction [Eq. 44] relative to the equilibrium vector length using the same color/line conventions as above. The spread is due to the same oscillations that cause the spread in the top panels. The results roughly agree with those of the nonoscillating calculations, plotted as solid green for neutrinos and dashed green for antineutrinos. Bottom panels : Length of the neutrino distribution flavor vector relative to the equilibrium length using the same color and line conventions as above. There is not spread because the distribution flavor vector length is not directly sensitive to oscillations.  \n<!-- image -->  \nthe electron flavor axis. However, the minimum f ( ⊥ ) increases once again before slowly decaying back to zero. Compared to the decay of the corresponding quantities in the absence of oscillations (plotted in green), the overall evolution is somewhat more rapid for neutrinos and slower for antineutrinos.  \nAs time progresses, the collisions work to remove offdiagonal components of the neutrino distribution, or equivalently to reduce f ( ⊥ ) . In the case of antineutrinos, which oscillate back and forth over the flavor equator, this rapidly shortens the total length of the distribution flavor vector. The bottom left panel of Fig. 21 shows the length of the distribution flavor vector relative to the equilibrium Fermi-Dirac length at the same neutrino energy of 8 MeV. The light curves show that the collisions decrease the antineutrino distribution flavor vector length to less than the equilibrium value, since the flavordiagonal parts of the collision term are unable to replenish the on-diagonal distribution components as quickly as the off-diagonal parts of the collision term decrease the off-diagonal distribution components. Once f ( ⊥ ) is depleted, the diagonal components and hence the vector length slowly decay back to equilibrium values. Since the  \noscillations cause the neutrino distribution flavor vector to spend more time than the antineutrino vector spends near the flavor axis, the collisions are less effective at shortening the length of the distribution flavor vector and so it decreases more slowly than the case when oscillations are switched off (solid green lines). In contrast, we see that since oscillations cause the antineutrino vector to spend more time near the equator, the vector length decreases more rapidly than when oscillations are not present (dashed green lines).  \nThis story remains true at higher energies as depicted in farther right plots in each row, but at higher energies the rate that at which diagonal components of the collision term equilibrate f ( z ) to L eq become comparable to the rate at which the off-diagonal components equilibrate f ( ⊥ ) to zero. This can be seen in the top panels as the centers of the bands for both neutrinos and antineutrinos approach equilibrium more quickly, though the bands in the bottom panel approach 0 at a similar rate. In the case of low-energy neutrinos, the distribution flavor vector decays to the flavor axis, and when f ( ⊥ ) is depleted f ( z ) > L eq and the vector decays down the flavor axis until it reaches f ( z ) = L eq . On the other hand, at high  \nenergies f ( z ) is driven to L eq before f ( ⊥ ) reaches 0. As a result, the range of latitudes covered by the oscillations extends below the starting latitude, and the off-diagonal components are able to push the vector toward the flavor axis with both f ( z ) (top panels) and L (bottom panels) less than L eq . Switching now to the antineutrinos, at low energies the distribution flavor vector decays to the flavor axis with f ( z ) and L less than L eq , as we described earlier. Just as the neutrinos, at higher energies the antineutrinos are driven to f ( z ) = L eq before f ( ⊥ ) decays to zero, resulting in values of L that do not dip as far below L eq as at low energies. For both neutrinos and antineutrinos, the continued presence of off-diagonal components when f ( z ) first reaches L eq allows continued evolution of f ( z ) and L away from L eq until f ( ⊥ ) decays to 0, and the interesting fluctuations in the top and bottom panels result. These behaviors are apparent in the evolution of the nonoscillating distributions (green curves) as well, but without oscillations driving the distribution flavor vectors through a wide range of latitudes, they tend not to be as pronounced.  \nFinally, the thickness of the bands in the top and center plots are much larger at low energies than at high energies. This is a result of our choice of initial conditions, since the differences between the distributions for different neutrino flavors at low energies permit a larger initial f ( ⊥ ) relative to f ( z ) [Eq. 17]. We remind the reader that our initial conditions also place all neutrinos and antineutrinos at the same phase angle (i.e., f ( x ) > 0 and f ( y ) = 0) and with Fermi-Dirac values for diagonal components. Understanding how the evolution of the neutrino distributions proceeds under different initial conditions is certainly worth further study.', 'C. Flavor decoherence in core-collapse supernovae': "In order to understand flavor decoherence in conditions throughout a CCSN explosion, we perform a series of isotropic homogeneous QKE calculations using a range of input parameters relevant to CCSNe. We see in Fig. 23 that a detailed accounting of all of the processes described in Sec. IV (except for elastic nucleon scattering) are required to describe flavor decoherence everywhere under the shock. However, before describing the results in detail, we will outline the details of the calculations.  \nThe range of parameters come from a snapshot of a one-dimensional neutrino radiation hydrodynamics CCSN simulation [79] shown in Fig. 22. This snapshot was taken at 100 ms after core bounce, by which time the shock has stalled at a radius of around 168 km. The spatial grid of the original simulation has 384 radial points extending out to 5000 km, but we only perform QKE calculations on each of the 241 radial points inside of the shock. Inside of the PNS (dashed green line) the temperatures are only ∼ 8 MeV (black curve), but electrons and electron neutrinos are trapped and very degenerate  \nFIG. 22. Background fluid snapshot from a one-dimensional neutrino radiation hydrodynamics core-collapse supernova simulation at 100 ms after core bounce [79]. We perform separate isotropic, homogeneous quantum kinetics calculations without the oscillation term at each radial point using the matter density ρ , electron fraction Y e , and temperature T at that radius as input. The bottom panel also shows for reference the chemical potentials of electrons µ e and electron neutrinos µ ν e as determined by the HShen equation of state [70]. The location of the shock is at the right edge of the plot at 168 km, and the approximate location of the outer edge of the protoneutron star is shown at 10 km with a vertical dashed green line.  \n<!-- image -->  \n(yellow and red curves). Just outside of the PNS the temperature reaches a maximum of 18 . 5 MeV and then continues to decrease with radius. The equilibrium electron neutrino degeneracy drops more quickly with radius than does the temperature, causing thermal electron neutrinos to become nondegenerate at r ≈ 20 km.  \nWe perform an isotropic calculation for the conditions associated with each location in the CCSN shown in Fig. 22, i.e. between the center of the PNS and the position of the shock at 168 km. In these calculations, we use a large energy grid spanning a domain of 200 MeV with energy grid spacing of 1 MeV in order to resolve the neutrino distribution anywhere under the shock. The large degeneracies within and near the PNS require an energy domain extending to 200 MeV to contain the neutrino distributions, and when the temperature and electron neutrino degeneracy drop at larger radii we need to adequately resolve the more compact distribution. We elect to initialize the distribution function in each calculation with maximally mixed Fermi-Dirac values [Eq. 17] based on the background fluid temperature and chemical potentials for the initial neutrino distributions, even though outside of the neutrinosphere the neutrino distribution is more sparse and anisotropic than a thermal distribu-  \ntion. Without spatial transport, even flavor-diagonal distributions will evolve toward Fermi-Dirac distributions if they begin out of equilibrium with the fluid. It would be difficult to disentangle this known effect from the new off-diagonal parts of the collision terms if we did not start with Fermi-Dirac initial conditions on the diagonals. In addition, we do not see a way to consistently map anisotropic distributions into isotropic calculations in a meaningful way, so we construct our initial conditions using Eq. 17 rather than from output of the CCSN simulation. Effects due to anisotropy should be addressed with anisotropic calculations with spatial transport, which will be the subject of future work.  \nUsing the results of each of these simulations, we can define a flavor decoherence timescale τ decohere as the amount of time it takes for the magnitude of the flavor off-diagonal component of the neutrino distribution to decrease to a factor of e -1 of its initial value. We show τ decohere throughout the CCSN profile in Fig. 23, color coded by neutrino energy. The left panels show the results when we only consider absorption onto nucleons and nuclei 'a la Sec. IV A. The interpretation of decay time is especially accurate in this case, since the off-diagonal components truly decay exponentially according to the absorption opacities. As expected, we see once again that neutrinos (top left panel) lose flavor coherence more quickly than antineutrinos (bottom left panel), and high-energy neutrinos decohere faster than low-energy neutrinos. Going inward from the shock (right side of the panel) to the protoneutron star, the decoherence timescales decrease as increasing matter densities and temperatures increase the interaction rates. Progressing into the PNS, the temperature decreases as the electron and neutron degeneracies increase, blocking neutrino and antineutrino absorption processes, respectively. In the PNS the electron neutrinos are also very degenerate and final-state neutrino blocking increases the effective electron neutrino absorption opacity (corrected for stimulated absorption), but the electrons more than compensate for this with an even larger degeneracy that greatly reduces the opacity.  \nIn the second panels from the left of Fig. 23 we perform the same calculations, but include inelastic electron scattering and electron-positron pair annihilation processes. Outside of the PNS the results are rather similar to those in the first panel, but electron scattering totally dominates absorption and lepton pair processes as a driver of flavor decoherence in the PNS. Just as with the absorption processes, the high electron degeneracy strongly blocks electron-positron pair production, and the lack of positrons due to the same degeneracy makes pair production inefficient in the PNS. Note that the decoherence times for the highest energy neutrinos (gold curves) discontinuously jump within the PNS. This is an artifact of how we measure decoherence times. At the lower branch, the blocking terms in the off-diagonal component of the scattering term redistribute neutrinos in a way that causes the off-diagonal component | f eµ ( t ) | to  \nquickly decrease below f eµ (0) /e . On the upper branch, f eµ ( t ) decreases past this point, continues through zero, and becomes large and negative; we must then wait for | f eµ ( t ) | to decrease below f eµ (0) /e a second time. With larger radii and smaller degeneracies, this overshoot effect becomes less severe. Eventually we reach a radius where the size of the overshoot is smaller than 1 /e of the initial value and our measured decoherence time drops to the lower branch.  \nIn the third panels from the left, we once again repeat the calculations now with the addition of nucleonnucleon bremsstrahlung interactions, implemented as effective absorption and emission of heavy lepton neutrinos. The bremsstrahlung reaction predominantly decreases the decoherence timescale of low-energy neutrinos and antineutrinos, leaving the high energies for the most part unchanged. In the antineutrino case, the effect is strong enough to make decoherence times increase with energy up to 24 MeV, above which decoherence times again decrease with increasing energy. The importance of the bremsstrahlung reaction increases with density and is especially important within the PNS. More work needs to be done to develop detailed bremsstrahlung kernels to replace the effective absorption treatment, but once developed for flavor-diagonal neutrinos, they can be extended and developed into the QKE collision term in the same way as the electron-positron pair annihilation kernels.  \nFinally, in the fourth panels from the left we repeat the calculations one more time and include the full suite of reactions described in this paper. Within the PNS, neutrino-neutrino scattering increases the decoherence times of low-energy neutrinos (black curves, top right panel), decreases those of neutrinos above 100 MeV (yellow curves, top right panel), and decreases those of all antineutrinos (bottom-right panel). Outside of the PNS, there is some overal decrease in decoherence times due to neutrino-neutrino interactions, though primarily at low energies. However, we must once again consider that an isotropic thermal distribution is not realistic outside of the decoupling region, and we have not tested the effects of different initial conditions. More work is required to understand the quantum kinetics of non-equilibrium distributions.  \nThe high computational cost of QKE calculations incites us to attempt to predict the flavor decoherence rates using only the flavor-diagonal interaction rates, i.e., without needing to run QKE calculations. Based on our QKE calculations, we find that the neutrino mean free path is not an accurate predictor of decoherence rates. However, we can empirically construct an effective decoherence opacity to reflect flavor decoherence rates as  \nκ effective := κ abs ,eµ + 1 2 ˜ κ scat ,eµ . (45)  \nThe value of κ abs we use here contains bremsstrahlung radiation and electron-positron pair annihilation as effective absorption, and κ scat is the elastic electron scattering  \nFIG. 23. Decoherence times of flavor off-diagonal components of the neutrino distribution function in conditions relevant to CCSNe. We perform an isotropic, homogeneous nonoscillating quantum kinetic calculation using fluid values from each radial point in Fig. 22. Top panels show the decoherence times for neutrinos and bottom panels show those for antineutrinos. Neutrino energies from 1 MeV to 100 MeV are colored according to the color bar on the right, while energies from 101 to 200 MeV are all colored gold. Each panel shows the results using a different set of interactions as indicated by the text in the panel. The approximate outer edge of the protoneutron star is depicted with a dashed green line at 10 km. For reference, we also show the location where the radial coordinate is equal to the decoherence length scale with a dashed blue line. Electron scattering, bremsstrahlung processes, and four-neutrino processes dominate decoherence in the PNS, while absorption dominates decoherence just under the shock. All processes except neutrino-nucleon elastic scattering and neutrino-neutrino processes contribute in the decoupling region.  \n<!-- image -->  \nopacity described in Sec. IV B. In Fig. 24 we compare the decoherence times from the third panels in Fig. 23 (i.e., including absorption, inelastic electron scattering, electron-positron pair annihilation, and nucleon-nucelon bremsstrahlung) to those predicted by κ effective . Within the PNS, the effective opacity predicts flavor decoherence rates for the 1 -100 MeV within a factor of 10, while the results from higher energies are thrown off by the overshoot discussed in reference to Fig. 23. Outside of 12 km κ effective predicts decoherence rates to within 20%, except for the four lowest antineutrino energies (1 -4 MeV) at ∼ 90 km. One might expect the charged-current component of the total opacity ( κ abs ,eµ + ˜ κ scat ,eµ ) to describe flavor decoherence, but this combination is only accurate to within a factor of 2 outside of 12 km and is significantly worse than κ effective inside the PNS. For completeness, we also examined the absorption opacity by itself and found that it consistently underpredicts the decoherence rates. Another combination, the total off-diagonal opacity ( κ abs ,eµ + κ scat ,eµ ), is also less successful than κ effective , since it is only good to within a factor of 10 outside of 12 km and much worse inside the PNS.", 'VI. DISCUSSION AND CONCLUSIONS': 'CCSN simulations have long included neutrino transport with a detailed set of collision rates, but until now the technology to simultaneously and self-consistently also treat neutrino flavor oscillations did not exist. In this work, we take the first steps toward developing the technology to simulate neutrino quantum kinetics in regions relevant to the CCSN explosion mechanism. In Sec. IV we demonstrate a means of converting existing neutrino interaction rates used in CCSN simulations into full collision terms for the quantum kinetic equations. For reference, we also describe the moment-integrated form of these source terms in Appendix A for application to moment-based QKE calculations. To demonstrate the use of these collision terms, we developed a novel method and open-source code IsotropicSQA for explicitly evolving both oscillations and collisions using stochastic integration to randomly sample the impact of collisions during the evolution (Sec. III). In our code and derivations, we include absorption of neutrinos by nucleons and nuclei, inelastic scattering by electrons, elastic scattering by nucleons, kernel-based electron-positron pair annihilation, kernel-based neutrino-neutrino scattering and an-  \nFIG. 24. Ratio of the computed flavor decoherence length scale to the effective decoherence opacity [Eq. 45] in the nonoscillating calculations including absorption, inelastic electron scattering, electron-positron pair annihilation, and effective-absorption nucleon-nucleon bremsstrahlung. This effective opacity is a much better predictor of flavor decoherence rates than the mean free path.  \n<!-- image -->  \nnihilation, and nucleon bremsstrahlung as an effective absorptive process.  \nIn Sec. V, we demonstrated the use of this new method and source terms to perform the first, albeit isotropic and homogeneous, direct evolution of the QKEs in conditions relevant to the CCSN explosion mechanism. In particular, we chose conditions of ρ = 10 12 g cm -3 , T = 10MeV, and Y e = 0 . 3, and evolved an initially strongly flavormixed neutrino distribution for 20 µ s. This was long enough to demonstrate that the distribution relaxes to a flavor-diagonal Fermi-Dirac distribution on timescales similar to but measurably different from those observed in calculations where the oscillation term is neglected.  \nGiven this insight, we performed a parameter sweep using nonoscillating calculations in conditions relevant to CCSNe. We find that electron scattering, nucleonnucleon bremsstrahlung, and four-neutrino processes are the dominant drivers of decoherence within the PNS, while just inside the shock front absorption is dominant. In the decoupling region, all of the processes discussed in this paper except for elastic nucleon scattering and neutrino-neutrino processes had a significant impact on flavor decoherence times. This demonstrates the need for a more sophisticated estimate of decoherence rates than the neutrino mean free path. To address this, we defined an effective decoherence opacity in Eq. 45 that predicts decoherence rates within ∼ 20% outside of 12 km and within a factor of 10 everywhere in our parameter set. In the future, we will test this estimate for different stages of a CCSN and under more realistic treatments of neutrino angular distributions.  \nNumerical neutrino quantum kinetics is a nascent field and there is a great deal of work to do. Our calculations are limited to isotropic and homogeneous neutrino and matter distributions, but we are hopeful that this work will encourage development of QKE codes with spa-  \ntial transport. This will be required to gain a quantitative understanding of how flavor coherence is transported throughout a CCSN, especially in the decoupling region where neutrino densities and collision rates are high. The importance of nucleon-nucleon bremsstrahlung in and near the PNS begs for an improved kernel-based treatment of these reactions rather than our effective absorption method. Fortunately, the same framework for generating QKE source terms from existing electronpositron pair production rates can also do the same for bremsstrahlung interaction rates. Inelastic nucleon scattering kernels could also be implemented in the framework of inelastic electron scattering. Finally, our treatment of neutrino-neutrino scattering and annihilation could become computationally prohibitive for more largescale calculations. Since this term requires integrating over the phase space of four neutrino distributions, the cost of calculating this term increases as N 4 E , where N E is the number of energy bins. In addition, our particular discretization of neutrino energy (bins having equal widths centered on integer multiples of the first bin center), though allowing for a straightforward implementation of the four-neutrino processes, does not allow a single energy grid to cover neutrino distributions in the range of temperatures and chemical potentials seen in CCSNe unless an intractable number of energy bins are used. Developing a more efficient way to accurately treat fourneutrino processes will be necessary for QKE calculations covering the wide range of conditions seen in a CCSN.  \nThe full QKE calculation presented in Sec. V required collision term time steps on the order of 2 × 10 -11 s and oscillation time steps on the order of 3 × 10 -15 s, and required about four days on as many cores. The cost of the calculation is in the large number of time steps needed to follow the oscillations and the high accuracy required of each step to prevent numerical artifacts from appearing. Since these isotropic and homogeneous calculations evolve a relatively small number of variables, the parallelizability is limited. However, future calculations including spatial transport will have many more quantities to evolve, and thus will be much more suited to taking advantage of larger computing resources.  \nWhile simulating the QKEs in CCSNe presents new challenges, doing so is essential to understanding the neutrino signal from and potentially also to the explosion mechanism of CCSNe. In addition to CCSNe, mergers of two neutron stars or of a neutron star and a black hole are astrophysical environment that are a particularly interesting home to potential quantum kinetic effects, as neutrino flavor transformations have already been suggested to be significant very close to the decoupling region [51, 52, 54]. There is much to be explored on the frontier of numerical quantum kinetics.', 'VII. ACKNOWLEDGMENTS': "SR is supported by the N3AS Fellowship under National Science Foundation Grant No. PHY-1630782 and Heising-Simons Foundation Grant No. 2017-228. G.M.,  \n- [1] S. A. Colgate and R. H. White, 'The hydrodynamic behaviour of supernoave explosions,' Astrophys. J. 143 , 626 (1966).\n- [2] H. A. Bethe, 'Supernova mechanisms,' Rev. Mod. Phys 62 , 801-866 (1990).\n- [3] H.-T. Janka, 'Explosion Mechanisms of Core-Collapse Supernovae,' Ann. Rev. Nucl. Part. Sci. 62 , 407-451 (2012).\n- [4] E. O'Connor, R. Bollig, A. Burrows, S. Couch, T. Fischer, H. T. Janka, K. Kotake, E. J. Lentz, M. Liebendorfer, O. E. Messer, A. Mezzacappa, T. Takiwaki, and D. Vartanyan, 'Global comparison of corecollapse supernova simulations in spherical symmetry,' J. Phys. G 45 , 104001 (2018).\n- [5] E. O'Connor and S. Couch, 'Exploring Fundamentally Three-dimensional Phenomena in High-fidelity Simulations of Core-collapse Supernovae,' Astrophys. J. 865 , 81 (2018).\n- [6] A. Summa, H. T. Janka, T. Melson, and A. Marek, 'Rotation-supported Neutrino-driven Supernova Explosions in Three Dimensions and the Critical Luminosity Condition,' Astrophys. J. 852 , 28 (2017).\n- [7] S. Richers, H. Nagakura, C. D. Ott, J. Dolence, K. Sumiyoshi, and S. Yamada, 'A Detailed Comparison of Multi-Dimensional Boltzmann Neutrino Transport Methods in Core-Collapse Supernovae,' Astrophys. J. 847 , 133 (2017).\n- [8] D. Vartanyan, A. Burrows, D. Radice, M. Aaron Skinner, and J. Dolence, 'A successful 3D core-collapse supernova explosion model,' Mon. Not. R. Astron. Soc. 482 , 351369 (2019).\n- [9] K. Kotake, T. Takiwaki, T. Fischer, K. Nakamura, and G. Mart'ınez-Pinedo, 'Impact of Neutrino Opacities on Core-Collapse Supernova Simulations,' Astrophys. J. 853 , 170 (2018).\n- [10] K.-C. Pan, A. Perego, T. Kuroda, R. M. Cabez'on, O. Heinimann, K. Ebinger, M. Liebendorfer, and F.K. Thielemann, 'Core-collapse supernovae in the hall of mirrors,' Astron. Astrophys. 619 , A118 (2018).\n- [11] S. W. Bruenn, 'Stellar Core Collapse: Numerical Model and Infall Epoch,' Astrophys J. Suppl. Ser. 58 , 771-841 (1985).\n- [12] A. Burrows, S. Reddy, and T. A. Thompson, 'Neutrino opacities in nuclear matter,' Nucl. Phys. A 777 , 356-394 (2006).\n- [13] B. Pontecorvo, 'Neutrino Experiments and the Problems of Conservation of Lepton Charge,' Zh. Eksp. Teor. Fiz 26 , 1717-1725 (1968).\n- [14] L. Wolfenstein, 'Neutrino Oscillations in Matter,' Phys. Rev. D 17 , 2369-2374 (1978).\n- [15] S. P. Mikheyev and A. Y. Smirnov, 'Resonant neutrino oscillations in matter,' Progress in Particle and Nuclear  \nA.V. and J.K. are supported at NC State by DOE Grant No. DE-FG02-02ER41216. We thank Hiroki Nagakura for access to core-collapse supernova simulation results. Some computations were performed on the Wheeler cluster at Caltech, which is supported by the Sherman Fairchild Foundation and by Caltech. Feynman diagrams are generated with the TIKZ-FEYNMAN package [80].  \nPhysics 23 , 41-136 (1989).  \n- [16] J. Pantaleone, 'Neutrino oscillations at high densities,' Phys. Lett. B 287 , 128-132 (1992).\n- [17] F. Capozzi, E. Lisi, A. Marrone, and A. Palazzo, 'Current unknowns in the three-neutrino framework,' Progr. Part. Nucl. Phys. 102 , 48-72 (2018).\n- [18] S. Horiuchi and J. P. Kneller, 'What can be learned from a future supernova neutrino detection?' J. Phys. G 45 , 043002 (2018).\n- [19] H. Duan, G. M. Fuller, and Y.-Z. Qian, 'Collective Neutrino Oscillations,' Ann. Rev. Nucl. Part. Sci. 60 , 569594 (2010).\n- [20] G. Bellini, L. Ludhova, G. Ranucci, and F. L. Villante, 'Neutrino oscillations,' Advances in High Energy Physics 2014 , 191960 (2014).\n- [21] L. Wolfenstein, 'Neutrino oscillations and stellar collapse,' Phys. Rev. D 20 , 2634-2635 (1979).\n- [22] H. Duan, G. M. Fuller, J. Carlson, and Y. Z. Qian, 'Coherent development of neutrino flavor in the supernova environment,' Phys. Rev. Lett. 97 , 241101 (2006).\n- [23] S. Chakraborty, T. Fischer, A. Mirizzi, N. Saviano, and R. Tom'as, 'No collective neutrino flavor conversions during the supernova accretion phase,' Phys. Rev. Lett. 107 , 151101 (2011).\n- [24] H. Duan and A. Friedland, 'Self-induced suppression of collective neutrino oscillations in a supernova,' Phys. Rev. Lett. 106 , 091101 (2011).\n- [25] B. Dasgupta, E. P. O'Connor, and C. D. Ott, 'Role of collective neutrino flavor oscillations in core-collapse supernova shock revival,' Phys. Rev. D 85 , 065008 (2012).\n- [26] R. F. Sawyer, 'Speed-up of neutrino transformations in a supernova environment,' Phys. Rev. D 72 , 045003 (2005).\n- [27] R. F. Sawyer, 'Neutrino Cloud Instabilities Just above the Neutrino Sphere of a Supernova,' Phys. Rev. Lett. 116 , 081101 (2016).\n- [28] S. Chakraborty, R. Hansen, I. Izaguirre, and G. Raffelt, 'Collective neutrino flavor conversion: Recent developments,' Nucl. Phys. B 908 , 366-381 (2016).\n- [29] I. Izaguirre, G. Raffelt, and I. Tamborra, 'Fast Pairwise Conversion of Supernova Neutrinos: A Dispersion Relation Approach,' Phys. Rev. Lett. 118 , 021101 (2017).\n- [30] F. Capozzi, B. Dasgupta, E. Lisi, A. Marrone, and A. Mirizzi, 'Fast flavor conversions of supernova neutrinos: Classifying instabilities via dispersion relations,' Phys. Rev. D 96 , 043016 (2017).\n- [31] B. Dasgupta, A. Mirizzi, and M. Sen, 'Simple method of diagnosing fast flavor conversions of supernova neutrinos,' Phys. Rev. D 98 , 103001 (2018).\n- [32] S. Abbar and H. Duan, 'Fast neutrino flavor conversion: roles of dense matter and spectrum crossing,' Phys. Rev. D 98 , 043014 (2018).\n- [33] I. Tamborra, L. Hudepohl, G. G. Raffelt, and H.-T. Janka, 'Flavor-dependent Neutrino Angular Distribution in Core-collapse Supernovae,' Astrophys J. 839 , 132 (2017).\n- [34] B. Dasgupta, A. Mirizzi, and M. Sen, 'Fast neutrino flavor conversions near the supernova core with realistic flavor-dependent angular distributions,' J. Cosmol. Astropart. Phys. 2 , 019 (2017).\n- [35] S. Abbar, H. Duan, K. Sumiyoshi, T. Takiwaki, and M. C. Volpe, 'On the occurrence of fast neutrino flavor conversions in multidimensional supernova models,' arXiv:1812.06883 .\n- [36] M. D. Azari, S. Yamada, T. Morinaga, W. Iwakami, H. Nagakura, and K. Sumiyoshi, 'Linear Analysis of Fast-Pairwise Collective Neutrino Oscillations in CoreCollapse Supernovae based on the Results of Boltzmann Simulations,' Phys. Rev. D 99 , 103011 (2019).\n- [37] J. F. Cherry, J. Carlson, A. Friedland, G. M. Fuller, and A. Vlasenko, 'Neutrino scattering and flavor transformation in supernovae,' Phys. Rev. Lett. 108 , 261104 (2012).\n- [38] S. Sarikas, I. Tamborra, G. Raffelt, L. Hudepohl, and H. T. Janka, 'Supernova neutrino halo and the suppression of self-induced flavor conversion,' Phys. Rev. D 85 , 113007 (2012).\n- [39] A. Mirizzi and P. D. Serpico, 'Flavor stability analysis of dense supernova neutrinos with flavor-dependent angular distributions,' Phys. Rev. D 86 , 085010 (2012).\n- [40] J. F. Cherry, J. Carlson, A. Friedland, G. M. Fuller, and A. Vlasenko, 'Halo modification of a supernova neutronization neutrino burst,' Phys. Rev. D 87 , 085037 (2013).\n- [41] V. Cirigliano, M. Paris, and S. Shalgar, 'Collective neutrino oscillations with the halo effect in single-angle approximation,' J. Cosmol. Astropart. Phys. 11 , 019 (2018).\n- [42] S. Rosswog and M. Liebendoerfer, 'High resolution calculations of merging neutron stars. 2: Neutrino emission,' Mon. Not. R. Astron. Soc. 342 , 673 (2003).\n- [43] L. Dessart, C. D. Ott, A. Burrows, S. Rosswog, and E. Livne, 'Neutrino signatures and the neutrino-driven wind in binary neutron star mergers,' Astrophys. J. 690 , 1681-1705 (2009).\n- [44] B. D. Metzger and R. Fernandez, 'Red or blue? A potential kilonova imprint of the delay until black hole formation following a neutron star merger,' Mon. Not. R. Astron. Soc. 441 , 3444-3453 (2014).\n- [45] S. Richers, D. Kasen, E. O'Connor, R. Fernandez, and C. D. Ott, 'Monte Carlo Neutrino Transport Through Remnant Disks From Neutron Star Mergers,' Astrophys. J. 813 , 38 (2015).\n- [46] S. Fujibayashi, Y. Sekiguchi, K. Kiuchi, and M. Shibata, 'Properties of Neutrino-driven Ejecta from the Remnant of Binary Neutron Star Merger : Purely Radiation Hydrodynamics Case,' Astrophys. J. 846 , 114 (2017).\n- [47] A. Perego, H. Yasin, and A. Arcones, 'Neutrino pair annihilation above merger remnants: Implications of a long-lived massive neutron star,' J. Phys. G 44 , 084007 (2017).\n- [48] M. R. Wu, I. Tamborra, O. Just, and H. T. Janka, 'Imprints of neutrino-pair flavor conversions on nucleosynthesis in ejecta from neutron-star merger remnants,' Phys. Rev. D 96 , 123015 (2017).\n- [49] D. Radice, A. Perego, S. Bernuzzi, and B. Zhang, 'Longlived remnants from binary neutron star mergers,' Mon. Not. R. Astron. Soc. 481 , 3670-3682 (2018).  \n- [50] F. Foucart, M. D. Duez, L. E. Kidder, R. Nguyen, H. P. Pfeiffer, and M. A. Scheel, 'Evaluating radiation transport errors in merger simulations using a Monte Carlo algorithm,' Phys. Rev. D 98 , 063007 (2018).\n- [51] Y. L. Zhu, A. Perego, and G. C. McLaughlin, 'Matterneutrino resonance transitions above a neutron star merger remnant,' Phys. Rev. D 94 , 105006 (2016).\n- [52] A. Malkus, G. C. McLaughlin, and R. Surman, 'Symmetric and standard matter neutrino resonances above merging compact objects,' Phys. Rev. D 93 , 045021 (2016).\n- [53] D. Vaananen and G. C. McLaughlin, 'Uncovering the matter-neutrino resonance,' Phys. Rev. D 93 , 105044 (2016).\n- [54] A. Chatelain and C. Volpe, 'Helicity coherence in binary neutron star mergers and nonlinear feedback,' Phys. Rev. D 95 , 043005 (2017).\n- [55] A. Chatelain and M. C. Volpe, 'Neutrino propagation in binary neutron star mergers in presence of nonstandard interactions,' Phys. Rev. D 97 , 023014 (2018).\n- [56] A. Vlasenko and G. C. McLaughlin, 'Matter-neutrino resonance in a multiangle neutrino bulb model,' Phys. Rev. D 97 , 083011 (2018).\n- [57] M. B. Deaton, E. O'Connor, Y. L. Zhu, A. Bohn, J. Jesse, F. Foucart, M. D. Duez, and G. C. McLaughlin, 'Elastic scattering in general relativistic ray tracing for neutrinos,' Phys. Rev. D 98 , 103014 (2018).\n- [58] V. Cirigliano, M. W. Paris, and S. Shalgar, 'Effect of collisions on neutrino flavor inhomogeneity in a dense neutrino gas,' Phys. Lett. B 774 , 258-267 (2017).\n- [59] J. Y. Tian, A. V. Patwardhan, and G. M. Fuller, 'Neutrino flavor evolution in neutron star mergers,' Physical Review D 96 , 043001 (2017).\n- [60] M. R. Wu and I. Tamborra, 'Fast neutrino conversions: Ubiquitous in compact binary merger remnants,' Phys. Rev. D 95 , 103007 (2017).\n- [61] G. Sigl and G. Raffelt, 'General kinetic description of relativistic mixed neutrinos,' Nucl. Phys. B 406 , 423451 (1993).\n- [62] A. Vlasenko, G. M. Fuller, and V. Cirigliano, 'Neutrino quantum kinetics,' Phys. Rev. D 89 , 105004 (2014).\n- [63] C. Volpe, 'Neutrino Quantum Kinetic Equations,' Int. J. Mod. Phys. E 24 , 1541009 (2015).\n- [64] V. Cirigliano, G. M. Fuller, and A. Vlasenko, 'A new spin on neutrino quantum kinetics,' Phys. Lett. B 747 , 27-35 (2015).\n- [65] D. N. Blaschke and V. Cirigliano, 'Neutrino quantum kinetic equations: The collision term,' Phys. Rev. D 94 , 033009 (2016).\n- [66] J. Serreau and C. Volpe, 'Neutrino-antineutrino correlations in dense anisotropic media,' Phys. Rev. D 90 , 125040 (2014).\n- [67] J. Y. Tian, A. V. Patwardhan, and G. M. Fuller, 'Prospects for neutrino spin coherence in supernovae,' Phys. Rev. D 95 , 063004 (2017).\n- [68] Z. Maki, M. Nakagawa, and S. Sakata, 'Remarks on the Unified Model · of Elementary Particles,' Progr. Theor. Phys. 28 , 870 (1962).\n- [69] Particle Data Group, 'Review of Particle Physics,' Phys. Rev. D 98 , 030001 (2018).\n- [70] H. Shen, H. Toki, K. Oyamatsu, and K. Sumiyoshi, 'Relativistic Equation of State for Core-Collapse Supernova Simulations,' Astrophys J. Suppl. Ser. 197 , 20 (2011).\n- [71] S. Furusawa, S. Yamada, K. Sumiyoshi, and H. Suzuki, 'A new baryonic equation of state at sub-nuclear densities for core-collapse simulations,' Astrophys J. 738 , 178 (2011).\n- [72] S. Furusawa, K. Sumiyoshi, S. Yamada, and H. Suzuki, 'New equations of state based on the liquid drop model of heavy nuclei and quantum approach to light nuclei for core-collapse supernova simulations,' Astrophys. J. 772 , 95 (2013).\n- [73] E. O'Connor, 'An Open-Source Neutrino Radiation Hydrodynamics Code for Core-Collapse Supernovae,' Astrophys J. 219 , 24 (2015).\n- [74] W. Yueh and J. Buchler, 'Neutrino transport in supernova models-SN method,' Astrophys. J. 217 , 565 (1977).\n- [75] C. J. Horowitz, 'Weak magnetism for antineutrinos in supernovae,' Phys. Rev. D 65 , 043001 (2002).\n- [76] S. Hannestad and G. Raffelt, 'Supernova Neutrino Opacity from Nucleon-Nucleon Bremsstrahlung and Related Processes,' Astrophys J. 507 , 339-352 (1998).\n- [77] Y. Li, M. K. Liou, W. M. Schreiber, and B. F. Gibson, 'Neutrino-pair bremsstrahlung from nucleonnucleon scattering,' Phys. Rev. C 92 , 015504 (2015).\n- [78] R. Buras, H.-T. Janka, M. T. Keil, G. G. Raffelt, and M. Rampp, 'Electron Neutrino Pair Annihilation: A New Source for Muon and Tau Neutrinos in Supernovae,' Astrophys. J. 587 , 320-326 (2003).\n- [79] H. Nagakura, W. Iwakami, S. Furusawa, K. Sumiyoshi, S. Yamada, H. Matsufuru, and A. Imakura, 'Threedimensional Boltzmann-Hydro Code for Core-collapse in Massive Stars . II . The Implementation of Moving-mesh for Neutron Star Kicks,' Astrophys J. Suppl. Ser. 229 , 42 (2017).\n- [80] J. P. Ellis, 'TikZ-Feynman: Feynman diagrams with TikZ,' Comput. Phys. Commun. 210 , 103-123 (2017).\n- [81] M. Shibata, K. Kiuchi, Y. I. Sekiguchi, and Y. Suwa, 'Truncated Moment Formalism for Radiation Hydrodynamics,' Progr. Theor. Phys. 125 , 1255-1287 (2011).\n- [82] C. Y. Cardall, E. Endeve, and A. Mezzacappa, 'Conservative 3+1 general relativistic variable Eddington tensor radiation transport equations,' Phys. Rev. D 87 , 103004 (2013).\n- [83] L. F. Roberts, C. D. Ott, R. Haas, E. O'Connor, P. Diener, and E. Schnetter, 'General-Relativistic Three-Dimensional Multi-Group Neutrino RadiationHydrodynamics Simulations of Core-Collapse Supernovae,' Astrophys. J. 831 , 98 (2016).\n- [84] F. Foucart, E. O'Connor, L. Roberts, M. D. Duez, R. Haas, L. E. Kidder, C. D. Ott, H. P. Pfeiffer, M. A. Scheel, and B. Szilagyi, 'Post-merger evolution of a neutron star-black hole binary with neutrino transport,' Phys. Rev. D 91 , 124021 (2015).\n- [85] S. Hannestad, G. G. Raffelt, G. Sigl, and Y. Y. Wong, 'Self-induced conversion in dense neutrino gases: Pendulum in flavor space,' Phys. Rev. D 74 , 105010 (2006).", 'Appendix A: Moment form': "Following [81, 82], we can take moments of the QKEs in order to generate evolution equations for a small number of angular moments of the distribution function. In order to do this, we decompose the neutrino four-momentum  \ninto  \np α = hν c ( u α + l α ) , (A1)  \nwhere l α is a unit normal four-vector ( l α l α = 1) orthogonal to u α ( u α l α = 0). We can define the first few comoving-frame angular moments (units of Hz 3 ) of the neutrino distribution function as  \nJ ab ( ν, x µ ) := ν 3 ∫ d Ω f ab ( ν, Ω , x µ ) H α ab ( ν, x µ ) := ν 3 ∫ d Ω f ab ( ν, Ω , x µ ) l α L αβ ab ( ν, x µ ) := ν 3 ∫ d Ω f ab ( ν, Ω , x µ ) l α l β N αβγ ab ( ν, x µ ) := ν 3 ∫ d Ω f ab ( ν, Ω , x µ ) l α l β l γ . (A2)  \nNote that the oscillation Hamiltonian H ab can be distinguished from the first moment of the radiation field H α ab by the presence of a spacetime index. The two-moment evolution equations from [81] Eq. 3.19 are  \n∇ β M αβ ab -∂ ∂ν ( νM αβγ ab ∇ γ u β ) = S α ab , (A3)  \nwhere  \nM αβ ab = J ab u α u β + H α ab u β + H β ab u α + L α β ab M αβγ ab = J ab u α u β u γ + H α ab u β u γ H β ab u α u γ + H γ ab u α u β + L αβ ab u γ + L αγ ab u β + L βγ ab u α + N αβγ ab . (A4)  \nIn this work, we will focus on the source terms, which can be expressed in moment form as  \nS α ab ( ν, x µ ) := ν 3 ∫ d Ω ( C ab -i ¯ hc [ H,f ] ab ) ( u α + l α ) (A5)  \n(units of Hz 3 cm -1 ). We leave the full implementation of relativistic QKEs to future work, but note that, when neglecting the O ( ϵ 2 ) force and drift term corrections, all of the flavor coherence effects lie in the source term. These evolution equations need to be massaged into a numerically implementable form (e.g., [73, 81 ? -84]).  \nOscillation Term - Performing the angular integrals on the contribution to the source term from the vacuum and matter potentials yields  \nS α ab = -i ¯ hc [ H vacuum + H matter , Ψ α ] ab . (A6)  \nSimilarly, performing the angular integral on the contribution from the neutrino self-interaction term, we get  \nS α ab = -i √ 2¯ h 2 G F ν 3 c ∫ dν ' ν ' { [( J ' -¯ J '∗ ) , Ψ α ] ab -[( H ' β -¯ H ' β ∗ ) , Ξ α β ] ab } , (A7)  \nwhere we have defined the following Lorentz-invariant objects (units of Hz 3 ):  \nΨ α ab ( ν, x µ ) := J ab u α + H α ab Ξ α βab ( ν, x µ ) := g βµ ( H µ ab u α + L µα ab ) . (A8)  \nAbsorption and Emission - To get the source term for the two-moment form of the QKEs [Eq. A3], we perform the integral in Eq. A5 for the contribution to the source term in Eq. 25 and arrive at  \nS α ab = 4 πν 3 j ( ν a ) δ ab u α -( ⟨ j ⟩ ab + ⟨ κ ⟩ ab )Ψ α ab . (A9)  \nThe source term for antineutrinos is exactly analogous.  \nScattering - Angular moments of the inelastic scattering kernel [Eq. 31] can be integrated as  \nν ' 3 ν 3 ∫ d Ω d Ω ' f ( u α + l α ) R ± = 4 πν ' 3 2 Φ ± 0 Ψ α ν ' 3 ν 3 ∫ d Ω d Ω ' f ' ( u α + l α ) R ± = 4 πν 3 2 (Φ ± 0 J ' u α +Φ ± 1 H ' α ) ν ' 3 ν 3 ∫ d Ω d Ω ' ff ' ( u α + l α ) R ± = 1 2 Φ ± 0 J ' Ψ α + 3 2 Φ ± 1 H ' β Ξ α β . (A10)  \nThe resulting source terms from performing the angular integral in Eq. A5 on Eq. 32 are  \nS α ab = ∫ dν ' c 4 ν ' { 4 π 2 [ ν 3 (Φ + 0 ab J ' ab u α +Φ + 1 ab H ' α ab ) -ν ' 3 ⟨ Φ ⟩ -0 ab Ψ α ab ] -ζ α ab } , (A11)  \nwhere  \nζ α ab := ν ' 3 ν 3 ∫ d Ω d Ω ' ( ς + ab -ς -ab )( u α + l α ) = 1 2 ∑ c [ 1 2 (∆Φ 0 cb Ψ α ac J ' cb +∆Φ 0 ac J ' ac Ψ α cb ) + 3 2 ( ∆Φ 1 cb Ξ α βac H ' β cb +∆Φ 1 ac H ' β ac Ξ α βcb ) ] , (A12)  \n(units Hz 7 cm 3 ) and ∆Φ ab := Φ + ab -Φ -ab . The contribution to S α aa from the terms inside the square brackets in Eq. A11 and the c = a terms in ζ α aa make the well-known source term for noncoherent neutrinos as in [81]. If the scattering process is treated as elastic [Eq. 35], the source terms simplify to  \nS α ab = [ -H α ab ( ⟨ κ ⟩ 0 ab -κ 1 ab / 3) -J ab ˜ κ 0 ab u α ] (A13)  \nPair Processes - The pair process (electron-positron and nucleon-nucleon bremsstrahlung) source terms can be integrated from Eq. 39 as  \nS α ab = ∫ d ¯ ν ' ν ' { 4 π 2 [4 πν 3 ¯ ν ' 3 Φ + 0 ab δ ab -ν 3 (Φ + 0 ab ¯ J ' ab u α /c +Φ + 1 ab ¯ H ' α ab ) -¯ ν ' 3 ⟨ Φ ⟩ + 0 ab Ψ α ab ] + ζ α ab } . (A14)  \nζ α ab is defined in Eq. A12, though for Eq. A14 one should use values of R ± ab for this process, and primed quantities should be primed and barred. The source term from integrating the effective absorption collision term [Eq. 40] looks identical to Eq. A9, but here the emissivities j ( ν a ) and opacities κ ( ν a ) for heavy lepton neutrinos are nonzero.", 'Appendix B: Code tests': 'To ensure that IsotropicSQA produces realistic results, we perform a few basic tests.', '1. Equilibrium test': 'In the limit of no flavor mixing, the interactions should reduce to the well-known flavor-diagonal interactions, which drive the neutrino distributions to Fermi-Dirac distributions described by the fluid temperature and chemical potentials of µ ν e = µ p + µ e -µ n , µ ¯ ν e = -µ ν e , and µ ν x = 0. Though we see the collision kernels in the main text driving the distribution function to this equilibrium, we quantify how well the collision terms maintain this equilibrium. To test this, we initialize the diagonal components of the distribution function to their thermal equilibrium values, set off-diagonal elements to zero, and evolve including all interactions in Sec. IV and without oscillations. All tests are performed under the conditions of ρ = 10 12 g cm -3 , T = 10MeV, Y e = 0 . 3, and a neutrino energy domain of 1 -101 MeV as in Sec.s IV and V. We run three tests: (a) 50 energy bins and an integration accuracy of 10 -12 , (b) 25 energy bins an an integration accuracy of 10 -12 , and (c) 50 bins and an integration accuracy of 10 -11 . In all cases, within 25 µ s the maximum deviation of any component of the distribution function settles to within one part in 2 × 10 -15 of its initial value. The 25 bin case briefly has maximum errors that grow to 5 × 10 -12 before settling down to the previously mentioned error. We also tested individual processes with similar results. These tests show that the equilibrium Fermi-Dirac distribution is reproduced to within numerical error by our collision processes.', '2. Vacuum oscillations': 'To ensure that the oscillations evolve correctly, we first simulate vacuum oscillations by setting the matter and interaction potentials to zero. The probability of a neutrino transitioning from one flavor to another is given by the well-known formula (e.g., [20])  \nP T ( t ) = sin 2 (2 θ 12 ) sin 2 ( c 4 ∆ m 2 12 t 4 E ¯ h ) (B1)  \nFIG. 25. Vacuum oscillations test . - The top panel shows the evolution of the ν e and ν µ distribution functions for the 10 MeV and 49MeV energy bins, along with the analytic solution. The bottom panel shows the maximum error in the solution for simulations with three different time step sizes. The error increases with decreasing step size. See text for details.  \n<!-- image -->  \nThus, we expect the distribution function values to follow  \nf ee ( t ) = [1 -P T ( t )] f ee (0) + P T ( t ) f µµ (0) f µµ ( t ) = [1 -P T ( t )] f µµ (0) + P T ( t ) f ee (0) (B2)  \nThe good agreement between the computed (colored lines) and analytic (dotted lines) results is shown in the top panel of Fig. 25. The bottom panel shows the maximum relative error among all energy groups and neutrino species between the computed and analytic solution as a function of the time step size. The error actually increases when the time step decreases. This is due to the hybrid representation of the distribution function matrices. After every oscillation step, the unitary evolution operator is applied to the distribution functions. Though the operator itself is represented very accurately, some accuracy is lost in applying the evolution operator to the distribution function. This emphasizes the importance of applying the operator as sparsely as possible, and it is the reason we allow the evolution operator to evolve for many steps between applications to the distribution function.', '3. MSW resonance': 'Oscillations in a constant matter background are also straightforward to compute. Eq. B2 still applies, but the  \nFIG. 26. MSW oscillations test . - The top panel shows the evolution of the ν e and ν µ distribution functions for the 10 MeV and 49MeV energy bins, along with the analytic solution. The bottom panel shows the maximum error in the solution for simulations with three different time step sizes. The error increases with decreasing step size. See text for details.  \n<!-- image -->  \nmixing angle and mass squared difference in Eq. B1 are replaced by effective values of  \nsin 2 (2 ˜ θ 12 ) = sin 2 (2 θ 12 ) sin 2 (2 θ 12 ) + C 2 ∆ ˜ m 2 12 = ∆ m 2 12 √ sin 2 (2 θ 12 ) + C 2 C = cos(2 θ 12 ) -2 V E ∆ m 2 12 c 4 (B3)  \nwhere V = ± √ 2 G F ¯ h 3 c 3 n e is the matter potential (positive for neutrinos, negative for antineutrinos). Fig. 26 shows the error associated with evolution in a matter potential that puts the 20 MeV neutrinos on the Mikheyev-Smirnov-Wolfenstein (MSW) resonance ( C = 0).', '4. Bipolar oscillations': 'Fig. 27 shows a calculation of bipolar oscillations and associated errors. Bipolar oscillations occur for pure electron anti/neutrinos with the inverted mass hierarchy and for pure anti/muon neutrinos with the normal mass hierarchy. In this test, we set up an initial distribution of pure anti/muon neutrinos with f µµ = ¯ f µµ = 1 in a single energy bin at an energy of hν . The background matter density is set to 0, the angle between the mass and flavor  \nFIG. 27. Bipolar oscillation test . -First panel: evolution of a neutrino distribution function f initially consisting of pure 50 MeV ν µ and ¯ ν µ in equal amounts ignoring collision terms. Second panel: Maximum error between the fiducial calculation (integration accuracy of 10 -13 and dt block = 10 -3 s) and that with lower integration accuracy (red line) and longer dt block (blue line). Third panel: Relative change in the length of the distribution flavor vector, which should be 0 since collision terms are set to zero. Fourth panel: How non-Hermitian the distribution function becomes at each step before enforcing Hermitivity.  \n<!-- image -->  \nbases to 0 . 01, and the self-interaction Hamiltonian to be  \nH neutrino = 10( m 2 2 -m 2 1 ) c 4 2 hν ( f -¯ f ) . (B4)  \nWe also set the mass difference to m 2 2 -m 2 1 = (2 hν )¯ h/c 4 , such that the approximate frequency for bipolar oscillations is then κ ≈ 0 . 995 s -1 , corresponding to the µ = 10 , ω = 1 case in [85]. The top panel of Fig. 27 shows the evolution of the neutrinos over a period of 8 s. The oscillations do indeed occur on approximately this timescale and match Fig. 1 in [85]. However, the period of the oscillations in the numerical solution decreases over time due to numerical errors caused predominantly by the mapping of the evolution matrix onto the distribution function; applying the mapping less often causes the period to decrease more slowly. The third panel shows  \nthe relative change in the length of the distribution flavor vector with time, which should always be zero in this test without collision terms. It is once again clear that a longer dt block prevents accumulation of error, and that the integration accuracy parameter does not significantly affect the solution on these timescales. Similarly, the distribution function should always be Hermitian. The nonHermitivity of f is shown in the bottom panel after each step prior to making f Hermitian again. Once again, the largest errors occur during the flavor transitions.', '5. Resolution': 'We demonstrate here that errors associated with our numerical treatment of the QKEs are within acceptable bounds. To do this, we rerun the combined oscillation and collision QKE calculation in Sec. V under varied parameter choices. The black curves in the top two panels of Fig. 28 show the length of the distribution flavor vectors for neutrinos and antineutrinos, run with 25 energy bins and an energy grid spacing of 4 MeV. Each curve represents one neutrino energy. The distribution flavor vector length is a particularly useful quantity to compare to, since it changes on collision timescales and does not oscillate violently as individual components of the neutrino distribution do. We aim to ensure that this macroscopic behavior converges and do not expect the instantaneous phases of the underlying oscillations to match between calculations.  \nUnderneath the black curves in Fig. 28 are the corresponding curves from each of the simulations with numerical tweaks. The relative errors between the fiducial and tweaked calculations are shown in the same color in the bottom panel. First we look at the green curves, which result from decreasing the accuracy of the integrator by a factor of 10 from the fiducial case. This is clearly not the dominant source of error. The solid green curve in the bottom panel shows that the relative error maximized over all neutrino energies and both helicities is always less than 1 . 4%. The dashed and dotted lines show the neutrino and antineutrino errors at 28 MeV, respectively, which is near the average energy of the distribution. The errors there are negligibly small.  \nNext, we increase the target impact by a factor of 10 and show the results as blue curves. These results sampled values of dt block that are on average also a factor of 10 larger, causing the collisions to more sparsely sample the rapidly varying distributions. This also leads to at most a 1.3% error. The sparse sampling seems only to make the evolution of L noisier without significantly changing its long-term evolution.  \nOf course, we also need to check our energy grid resolution. The red curves (data from Sec.V) show that doubling the energy resolution leads to a much more systematic error of at most ∼ 1 . 4%, and an even lower error at 28 MeV. Here, we only calculate and plot the errors at neutrino energies at bin centers present in both grids  \nFIG. 28. Fidelity tests . - Calculations evolving the QKEs including oscillations and collisions. Top panel: Neutrino distribution flavor vector length relative to the flavor-diagonal Fermi-Dirac length. Center panel: antineutrino distribution flavor vector length. Bottom panel: Relative error compared to the fiducial calculation. Black curves show the results of a fiducial calculation using 25 energy bins and an energy grid spacing of 4 MeV, while red curves show the results from Sec. V with 50 energy bins and an energy spacing of 2 MeV. Green curves result from worsening the accuracy of the integrator to 10 -11 . Blue curves result from increasing the target impact to 10 -3 . Gold curves result from increasing the upper bound of the neutrino energy domain and the number of energy bins by 20%, keeping the energy grid spacing constant. In the bottom panel, solid lines take the maximum over all energies and helicities, dashed lines show only errors for neutrinos at 28 MeV, and dotted lines show only errors for antineutrinos at 28 MeV.  \n<!-- image -->  \n(i.e., integer multiples of 4 MeV). Interestingly, the results seem to be much more sensitive to the size of the energy domain. The gold curves demonstrate the effect of extending the grid out to 120 MeV with an additional five energy zones in order to keep the grid spacing constant. This leads to significant errors of up to 3.7%, again where errors are only calculated for neutrino energies contained in both domains. However, we note that the largest errors are at the highest energies close to the energy boundary. If we instead look at the errors at 28 MeV (near the average neutrino energy of 31 . 5 MeV), the error is always under 0.27% for antineutrinos and 0.12% for neutrinos. The fiducial grid spans a domain of 2 -102 MeV, corresponding to 99 . 7% of the neutrino number and 99 . 0% of the neutrino energy. The extended grid spans a domain of 2 -122 MeV, corresponding to 99 . 9% of the neutrino number and 99 . 8%of the neutrino energy. Given the small number of neutrinos in these high-energy bins, it is unclear why the energy domain size is so important.'}