# Evolution of collisional neutrino flavor instabilities in spherically
symmetric supernova models

Zewei Xiong<EMAIL_ADDRESS>GSI Helmholtzzentrum für Schwerioneneforschung,
Planckstraße 1, 64291 Darmstadt, Germany  Meng-Ru Wu  Institute of Physics,
Academia Sinica, Taipei 11529, Taiwan  Institute of Astronomy and
Astrophysics, Academia Sinica, Taipei 10617, Taiwan  Physics Division,
National Center for Theoretical Sciences, Taipei 10617, Taiwan  Gabriel
Martínez-Pinedo  GSI Helmholtzzentrum für Schwerioneneforschung, Planckstraße
1, 64291 Darmstadt, Germany  Institut für Kernphysik (Theoriezentrum),
Fachbereich Physik, Technische Universität Darmstadt, Schlossgartenstraße 2,
64289 Darmstadt, Germany  Tobias Fischer  Institute of Theoretical Physics,
University of Wroclaw, Pl. M. Borna 9, 50-204 Wroclaw, Poland  Manu George
Institute of Physics, Academia Sinica, Taipei 11529, Taiwan  Chun-Yu Lin
National Center for High-performance Computing, National Applied Research
Laboratories, Hsinchu Science Park, Hsinchu City 30076, Taiwan  Lucas Johns
Departments of Astronomy and Physics, University of California, Berkeley, CA
94720, USA

###### Abstract

We implement a multigroup and discrete-ordinate neutrino transport model in
spherical symmetry which allows to simulate collective neutrino oscillations
by including realistic collisional rates in a self-consistent way. We utilize
this innovative model, based on strategic parameter rescaling, to study a
recently proposed collisional flavor instability caused by the asymmetry of
emission and absorption rates between $\nu_{e}$ and $\bar{\nu}_{e}$ for four
different static backgrounds taken from different stages in a core-collapse
supernova simulation. Our results confirm that collisional instabilities
generally exist around the neutrinosphere during the supernova accretion and
postaccretion phase, as suggested by Johns [arXiv:2104.11369.]. However, the
growth and transport of flavor instabilities can only be fully captured by
models with global simulations as done in this work. With minimal ingredient
to trigger collisional instabilities, we find that the flavor oscillations and
transport mainly affect (anti)neutrinos of heavy lepton flavors around their
decoupling sphere, which then leave imprints on their energy spectra in the
free-streaming regime. For electron (anti)neutrinos, their properties remain
nearly intact. We also explore various effects due to the decoherence from
neutrino-nucleon scattering, artificially enhanced decoherence from emission
and absorption, neutrino vacuum mixing, and inhomogeneous matter profile, and
discuss the implication of our work.

##  I Introduction

The phenomenology of neutrino flavor oscillations has been established by
experiments with solar, atmospheric, reactor, and accelerator neutrinos [1].
From these experiments, it is well understood that flavor oscillations depend
not only on the properties of neutrino mixing in vacuum, but also the coherent
forward scattering of neutrinos with electrons in medium, e.g., the earth
mantle or the solar interior.

In astrophysical environments such as core-collapse supernovae (CCSNe) and
binary neutron star mergers (BNSMs), neutrino fluxes are sufficiently intense
so that forward scattering among neutrinos is important and cannot be ignored.
The nonlinear interaction among neutrinos themselves leads to various
collective phenomena and flavor instabilities in two main categories: the so-
called “slow” mode (e.g., Refs. [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
15, 16, 17]), and “fast” flavor conversion (e.g., Refs. [18, 19, 20, 21, 22,
23, 24, 25, 26, 27, 28]); see also review papers [29, 30, 31, 32, 33] and
references therein. Both categories require a crossing from positive to
negative values in the distribution of the neutrino electron lepton numbers
($\nu$ELN) [34, 35]. The distinction is that the slow mode relies on the
crossing in energy spectrum while the fast mode requires one in the angular
distribution. For the slow mode, it typically has a flavor conversion length
scale of $\mathcal{O}(10-100)$ km and occurs far outside the neutrinosphere
after neutrinos decouple from the medium. For the fast mode, it converts
neutrino flavors in a much shorter length scale
$\sim(G_{F}n_{\nu})^{-1}\sim\mathcal{O}(1)$ cm with $G_{F}$ the Fermi constant
and $n_{\nu}$ the neutrino number density. Studies have shown that fast
instabilities generally exist in certain regions near or even inside the
neutrinosphere in CCSNe modeled by multidimensional simulations [36, 37, 38,
39, 40, 41, 42], and are even more ubiquitously present in the postmerger
environments of BNSMs [43, 44, 45, 46, 47, 48, 49, 50]. The potential
importance and the associated very short length scale of fast mode catalyzes
recent surges of local dynamical simulations in tiny boxes [51, 28, 52, 53,
54, 55, 56, 57, 58, 59, 50].

The change of neutrino flavor content due to oscillations can affect the
matter composition and the energy exchange of neutrinos with medium through
collisional processes, which can also lead to feedback effects on the flavor
evolution of the neutrino gas. A self-consistent treatment including both
coherent flavor evolution and collisions have been formulated as the neutrino
quantum kinetic equation (QKE) [60, 61, 62, 63]. Studies that solve simplified
QKE by including reduced set of collisions have been done in recent years [64,
65, 66, 67, 68, 69, 70]. Among those, Ref. [66, 67, 69, 70] examined the
impact of neutrino-nucleon scattering (NNS) on the evolution of fast flavor
oscillations under the condition that fast instabilities exist a priori.

Recently, a novel kind of flavor instability caused by the asymmetry of
emission and absorption (EA) rates between electron neutrinos $\nu_{e}$ and
electron antineutrinos $\bar{\nu}_{e}$ has been proposed [71]. Together with
an additional assumption that heavy lepton neutrinos may have a smaller number
density than that of $\bar{\nu}_{e}$ in regions where neutrinos decouple, it
was found that this collisional flavor instability can trigger flavor
conversion within a length scale in the order of inverse mean free path from
EA processes, i.e., $\sim\mathcal{O}$(1) km near the neutrino sphere. Although
this length scale is much longer than that of the fast mode, the required
condition for the collisional instability to occur completely differs from
that of fast instabilities such that both mechanisms can work independently.
In particular, it can potentially trigger flavor conversions in regions with
low electron number fraction ($Y_{e}$) where neutrinos decouple in CCSNe [71].

The pioneer works in Refs. [71, 72] on the collisional instability were based
on the assumption of spatial homogeneity. One general question, however, can
be raised: How can a simulation including collisions near the neutrino sphere
be performed self-consistently without including the global advection of
neutrinos? This question is not restricted to the collisional instability but
applicable to all long-term simulations including both collisions and
anisotropic angular $\nu$ELN distribution. The fundamental reasons that
trigger the flavor instabilities such as the diluted heavy lepton neutrino
number flux due to diffusion or the anisotropic angular distributions are
consequences of neutrino advection in the presence of collisional neutrino
processes and inhomogeneity in the length scale of astrophysical environment
[73]. Thus, when the advection of neutrinos is neglected, those conditions for
flavor instabilities may not be maintained self-consistently such that the
outcome of the long-term neutrino flavor evolution simulations can be
inaccurate. To address this issue, several simulations that aimed to examine
flavor oscillations triggered by fast and/or slow modes with the inclusion of
global advection in spherically symmetric geometry under given hydrostatic
radial profiles have been reported lately [74, 75, 76, 77]111Ref. [78]
implemented flavor oscillations in spherically symmetric CCSN simulations, but
taking the assumptions that flavor oscillations only occur much above the
region where neutrinos decouple. Moreover, the flavor evolution history for
all radially outgoing neutrinos were assumed to be independent of their
propagation angles..

In this work, we extend previous analyses and consider for the first time the
collisional instability and its effects in spherically symmetric hydrostatic
CCSN backgrounds. Therefore, we solve the quantum kinetic transport equations
using a multi-energy and multi-angle collective neutrino oscillation
simulator, an extended version of COSE$\nu$ [79]. As spherically symmetric
neutrino transport with realistic collisional rates do not lead to angular
spectrum crossing and hence no fast instability [80], it provides a clean
background for us to probe the consequence of collisional flavor instability,
which should exist in regions where the fast instabilities do not exist. We
consider several radial profiles of thermal quantities at selected post-bounce
times, obtained from a CCSN simulation. For each snapshot, we then simulate
the evolution of neutrinos up to $\sim 1$ ms in our simulator including
advection, collisions, and flavor oscillations based on state-of-the-art weak
rates determined from the background CCSN profiles, however, without feedback
on the medium. We start with fiducial models by including minimal but
essential ingredients that are able to trigger the collisional instability.
Additional ingredients are added case by case to explore their possible
impact.

This paper is organized as follows. In Sec. II, we describe our models and the
adopted parameters. We present the linear stability analyses for collisional
instabilities in Sec. III. We discuss our simulation results of fiducial
models in Sec. IV and those with additional parameters in Sec. V. Further
discussions and conclusions are given in Sec. VI. We adopt natural units and
$\hbar=c=k_{B}=1$ throughout the paper.

##  II Models

###  II.1 CCSN simulation and neutrino collisional processes

We use AGILE-BOLTZTRAN to simulate a CCSN, launched from a 18 $M_{\odot}$
progenitor star, based on general relativistic neutrino radiation-
hydrodynamics in spherical symmetry [81, 82, 83]. A comoving baryon mass mesh
is used for 208 radial grid points [84, 85], which features an adaptive mesh
refinement method [86] and implements as coordinate the enclosed baryon mass
instead of the actual mass mesh location [87].

For the current study, we use the nuclear equation of state of Ref. [88],
which is based on the nuclear statistical equilibrium approach for the
composition of nuclei based on several 1000 species, in combination with the
density dependent DD2 relativistic mean-field model of Ref. [89]. For low
temperatures, $T<0.45$ MeV, the ideal silicon- and sulfur gas approximation is
applied. Electron, positron, photon, and Coulomb contributions are added
following the equation of state of [90]. Note that in the current study muons
are not included in the SN simulation.

The Boltzmann neutrino-transport module of AGILE-BOLTZTRAN employs the
discrete ordinate method, and evolves four species of neutrinos, $\nu_{e}$,
$\bar{\nu}_{e}$, $\nu_{x}$ (for $\nu_{\mu}$ or $\nu_{\tau}$), and
$\bar{\nu}_{x}$ (for $\bar{\nu}_{\mu}$ or $\bar{\nu}_{\tau}$) without coherent
flavor oscillations. Although several weak collisional processes associated
with muons were employed recently [91, 92], we only include the processes from
(1a) to (4b) listed in Table 1: EA [93], NNS [94, 81], neutrino-electron
scattering (NES) [83], and neutrino-pair reactions (NPR), including both
leptonic process [94] and nucleon-nucleon bremsstrahlung [95, 96].

Table 1:  Set of weak processes considered in BOLTZTRAN (short as B) or COSE$\nu$ (C), where $\nu$ and $\bar{\nu}$ are for all neutrino flavors and $N=n,p$. Label | Weak process | Abbreviation | Adoption
---|---|---|---
(1a) | $\nu_{e}+n\leftrightarrows p+e^{-}$ | EA |  B,C
(1b) | $\bar{\nu}_{e}+p\leftrightarrows n+e^{+}$ | EA |  B,C
(1c) | $\bar{\nu}_{e}+p+e^{-}\leftrightarrows n$ | EA |  B,C
(2a) | $\nu+N\leftrightarrows\nu+N$ | NNS |  B,C
(2b) | $\bar{\nu}+N\leftrightarrows\bar{\nu}+N$ | NNS |  B,C
(3) | $\nu+e^{\pm}\leftrightarrows\nu+e^{\pm}$ | NES | B
(4a) | $\nu+\bar{\nu}\leftrightarrows e^{-}+e^{+}$ | NPR | B
(4b) | $\nu+\bar{\nu}+N+N\leftrightarrows N+N$ | NPR | B
(5a) | $\nu_{\mu}+n\leftrightarrows p+\mu^{-}$ | EA (muonic) | C
(5b) | $\bar{\nu}_{\mu}+p\leftrightarrows n+\mu^{+}$ | EA (muonic) | C

The roles of those collisional processes on the spectra of various neutrino
species are different. The energy spectra of $\nu_{e}$ and $\bar{\nu}_{e}$
near the neutrino sphere of last inelastic collision are predominantly
determined by the EA processes 222Neutrinos with different energy decouple
from matter at different radii; for details, c.f. Ref. [97].. For $\nu_{x}$
and $\bar{\nu}_{x}$, they primarily interact with the medium through the
neutral-current NNS. The NNS are considered here in the elastic approximation,
known as iso-energetic NNS. In particular near the $\nu_{x}$ and
$\bar{\nu}_{x}$ sphere of last elastic scattering, where the nucleon mass is
much larger than neutrino energy, little energy exchange between neutrinos and
medium happens. Thus, the NNS only play important role in trapping $\nu_{x}$
and $\bar{\nu}_{x}$ and define their “transport neutrino sphere”, which
locates outside the “energy sphere”, within which the NES and NPR processes
thermalize $\nu_{x}$ and $\bar{\nu}_{x}$ [94, 98, 99]. For instance, when NES
are included, the mean energy of $\nu_{x}$ and $\bar{\nu}_{x}$ outside their
neutrino spheres can be as low as $\sim 60\%$ of the values obtained without
including NES [99].

The numerical implementation of EA and NNS processes are very different from
NES and NPR, because the later two are inelastic in nature such that the
scattering kernels connect both the energies and angles of incoming and
outgoing neutrinos, which significantly increase the computational cost. When
more angular grids are used in, e.g., the simulations by COSE$\nu$ (see next
subsection), the computational complexity increases dramatically. This
prevents us from including NES and NPR unless approximated prescriptions for
them are employed. However, since these rates are subleading for $\nu_{e}$ and
$\bar{\nu}_{e}$ whose rate difference is the key to trigger collisional
instability, we expect that the results would remain qualitatively true when
the NES and NPR rates are omitted. We note that for $\nu_{\mu}$ and
$\bar{\nu}_{\mu}$, if we only include the elastic NNS, there is no other
process responsible for thermalizing them. To partly compensating for this, we
additionally include the muonic EA processes of (5a)-(5b) in Table 1 in
extended COSE$\nu$ simulations described in the next section, which can
thermalize the spectrum of muonic neutrinos with energy $\gtrsim
65\,\mathrm{MeV}$ [92].

###  II.2 Neutrino transport including flavor oscillations and collisions

####  II.2.1 Quantum kinetic equation

Although a full description of oscillation phenomena in CCSNe may require to
consider three active neutrino flavors [100, 101, 102, 103], we take the
approximation by considering only electron and muon flavors. The general
equations governing the spatial and temporal evolution of Wigner transformed
density matrices $\varrho$ (for neutrinos) and $\bar{\varrho}$ (for anti-
neutrinos), are given by

$\displaystyle(\partial_{t}+v_{r}\partial_{r}+\frac{1-v_{r}^{2}}{r}\partial_{v_{r}})\varrho(E,v_{r},r,t)=$
$\displaystyle-i[\mathbf{H}(E,v_{r},r,t),\varrho(E,v_{r},r,t)]+\mathbf{C}(E,v_{r},r,t),$
(1)
$\displaystyle(\partial_{t}+v_{r}\partial_{r}+\frac{1-v_{r}^{2}}{r}\partial_{v_{r}})\bar{\varrho}(E,v_{r},r,t)=$
$\displaystyle-i[\bar{\mathbf{H}}(E,v_{r},r,t),\bar{\varrho}(E,v_{r},r,t)]+\bar{\mathbf{C}}(E,v_{r},r,t),$
(2)

with

$\displaystyle\varrho(E,v_{r},r,t)$
$\displaystyle=\begin{bmatrix}\varrho_{ee}(E,v_{r},r,t)&\varrho_{e\mu}(E,v_{r},r,t)\\\
\varrho_{e\mu}^{*}(E,v_{r},r,t)&\varrho_{\mu\mu}(E,v_{r},r,t)\end{bmatrix},$
$\displaystyle\bar{\varrho}(E,v_{r},r,t)$
$\displaystyle=\begin{bmatrix}\bar{\varrho}_{ee}(E,v_{r},r,t)&\bar{\varrho}_{e\mu}(E,v_{r},r,t)\\\
\bar{\varrho}_{e\mu}^{*}(E,v_{r},r,t)&\bar{\varrho}_{\mu\mu}(E,v_{r},r,t)\end{bmatrix},$
(3)

in the flavor basis and normalized by the neutrino number density so that
$n_{\nu_{i}}(r,t)=\int dE\,dv_{r}\,\varrho_{ii}(E,v_{r},r,t)$ where
$i=e,\,\mu$333The dependence on spacial and temporal indices will not be
explicitly written unless emphasized in the following discussions.. Clearly,
the diagonal elements $\varrho_{ii}$ relate to the neutrino phase-space
distribution function $f_{\nu_{i}}$ by $\varrho_{ii}=f_{\nu_{i}}\times
E^{2}/(2\pi^{2})$ while the off diagonal elements characterize the degree of
flavor mixing. In Eqs. (1) and (2), the advection term
$[v_{r}\partial_{r}-(1/r)(1-v_{r}^{2})\partial_{v_{r}}]$ takes a
simplification from the full QKE [63, 104], and the second part accounts for
the aberration as the radial projection of velocity $v_{r}$ (the cosine of the
angle between the traveling direction of neutrino and the radial direction)
changes for neutrinos propagating non-radially.

On the right-hand side, $\mathbf{H}$ and $\mathbf{C}$ are the coherent
propagation Hamiltonian and collisional term for neutrinos respectively. Three
different contributions to $\mathbf{H}$ include the vacuum mixing term

$\mathbf{H}_{\mathrm{vac}}(E)=\frac{\delta m^{2}}{4E}\begin{bmatrix}-\cos
2\theta_{V}&\sin 2\theta_{V}\\\ \sin 2\theta_{V}&\cos
2\theta_{V}\end{bmatrix},$ (4)

where $\delta m^{2}$ is the vacuum mass-squared difference and $\theta_{V}$ is
the mixing angle, the diagonal matrix of the matter term

$\mathbf{H}_{\mathrm{mat}}=\mathrm{diag}[V_{\rm mat},0],$ (5)

with $V_{\rm mat}=\sqrt{2}G_{F}\rho Y_{e}/m_{u}$ corresponding to neutrino
forward scattering on $e^{\pm}$, neutrons, and protons where $m_{u}$ is the
atomic mass, and the neutrino self-induced term

$\displaystyle\mathbf{H}_{\nu\nu}(v_{r})=$ $\displaystyle\sqrt{2}G_{F}\int
dE^{\prime}\,dv_{r}^{\prime}\,\times$
$\displaystyle(1-v_{r}v_{r}^{\prime})[\varrho(E^{\prime},v_{r}^{\prime})-\bar{\varrho}^{*}(E^{\prime},v_{r}^{\prime})],$
(6)

corresponding to neutrino forward scattering on other neutrinos.

For the collisional term $\mathbf{C}$, we include the processes of EA and NNS
in our calculations. The contribution from EA is given by

$\displaystyle\mathbf{C}_{\mathrm{EA}}(E)=$
$\displaystyle\frac{1}{2}\left\\{\mathrm{diag}[j_{e}(E),j_{\mu}(E)],\varrho_{\rm
FO}(E)-\varrho(E)\right\\}$
$\displaystyle-\frac{1}{2}\left\\{\mathrm{diag}[\chi_{e}(E),\chi_{\mu}(E)],\varrho(E)\right\\},$
(7)

where $j_{e}(E)$ and $j_{\mu}(E)$ are the emissivities for the reactions for
processes of (1a)-(1c) and (5a)-(5b) in Table 1, $\chi_{e}(E)$ and
$\chi_{\mu}(E)$ are the opacities for their inverse reactions, $\varrho_{\rm
FO}(E)=\mathcal{I}\times E^{2}/(2\pi^{2})$, with $\mathcal{I}$ being the
identity matrix, is the fully-occupied differential number density for a
specific $E$, and the curly bracket is the anti-commutator. For the process of
NNS, we have

$\displaystyle\mathbf{C}_{\mathrm{NNS}}(E,v_{r})=\int dv_{r}^{\prime}\,$
$\displaystyle R_{\rm NNS}(E,v_{r},v_{r}^{\prime})\times$
$\displaystyle[\varrho(E,v_{r}^{\prime})-\varrho(E,v_{r})],$ (8)

where $R_{\rm NNS}(E,v_{r},v_{r}^{\prime})$ is the scattering kernel
transferring neutrino of energy $E$ with a radial velocity $v_{r}$ to the same
energy but a different velocity $v_{r}^{\prime}$. The special feature of the
scattering kernel in NNS, $R_{\rm NNS}(E,v_{r},v_{r}^{\prime})=\chi_{\rm
NNS}(E)/2+v_{r}v_{r}^{\prime}\tilde{\chi}_{\rm NNS}(E)/2$, allows further
simplification

$\displaystyle\mathbf{C}_{\mathrm{NNS}}(E,v_{r})=$ $\displaystyle-\chi_{\rm
NNS}(E)\left[\varrho(E,v_{r})-\int\frac{dv_{r}^{\prime}}{2}\,\varrho(E,v_{r}^{\prime})\right]$
$\displaystyle+v_{r}\tilde{\chi}_{\rm
NNS}(E)\int\frac{dv_{r}^{\prime}}{2}\,v_{r}^{\prime}\varrho(E,v_{r}^{\prime}),$
(9)

where $\chi_{\rm NNS}(E)$ and $\tilde{\chi}_{\rm NNS}(E)$ are the opacity in
NNS. Notice that the Pauli blocking in iso-energetic NNS has no impact because
$R_{\rm NNS}(E,v_{r},v_{r}^{\prime})=R_{\rm NNS}(E,v_{r}^{\prime},v_{r})$ as
shown in Appendix A.

For anti-neutrinos, $\bar{\mathbf{H}}$ and $\bar{\mathbf{C}}$ on the right-
hand side are defined similarly:
$\bar{\mathbf{H}}_{\mathrm{vac}}(E)=\mathbf{H}_{\mathrm{vac}}(E)$,
$\bar{\mathbf{H}}_{\mathrm{mat}}=-\mathbf{H}_{\mathrm{mat}}$,
$\bar{\mathbf{H}}_{\nu\nu}(v_{r})=-\mathbf{H}^{*}_{\nu\nu}(v_{r})$,

$\displaystyle\bar{\mathbf{C}}_{\mathrm{EA}}(E)=$
$\displaystyle\frac{1}{2}\left\\{\mathrm{diag}[\bar{j}_{e}(E),\bar{j}_{\mu}(E)],\bar{\varrho}_{\rm
FO}(E)-\bar{\varrho}(E)\right\\}$
$\displaystyle-\frac{1}{2}\left\\{\mathrm{diag}[\bar{\chi}_{e}(E),\bar{\chi}_{\mu}(E)],\bar{\varrho}(E)\right\\},$
(10)

and

$\displaystyle\bar{\mathbf{C}}_{\mathrm{NNS}}(E,v_{r})=\int dv_{r}^{\prime}\,$
$\displaystyle\bar{R}_{NNS}(E,v_{r},v_{r}^{\prime})\times$
$\displaystyle[\bar{\varrho}(E,v_{r}^{\prime})-\bar{\varrho}(E,v_{r})],$ (11)

where $\bar{j}_{e}(E)$, $\bar{j}_{\mu}(E)$, $\bar{\chi}_{e}(E)$,
$\bar{\chi}_{\mu}(E)$, and $\bar{R}_{\rm NNS}(E,v_{r},v_{r}^{\prime})$ are the
emissivities, opacities, and scattering kernel for antineutrinos.

We use emissivities and opacities from BOLTZTRAN and calculate the kernel of
NNS based on the method in Refs. [94, 105] for the simulations with COSE$\nu$
(see Appendix A for details).

####  II.2.2 Attenuation factors used in simulations

The rates involved in the flavor evolution equation span a wide range in
magnitudes. Due to the nature of weak interactions, the collisional rates are
much lower than those of coherent forward scatterings by a factor of $\sim
G_{F}E^{2}$ that can be $\sim 10^{-9}$–$10^{-5}$ for relevant neutrino
energies. For example, Fig. 1 compares the scales of all rates for neutrino
energy $E=21\,\mathrm{MeV}$ computed at a post-bounce time
$t_{\mathrm{pb}}\approx 247$ ms of the simulated CCSN. The magnitude of
neutrino self-induced term $\mathbf{H}_{\nu\nu}$ can be evaluated by taking
the difference of its diagonal elements
$V_{\nu\nu}=H_{\nu\nu,ee}-H_{\nu\nu,\mu\mu}$. It varies from $\sim
10^{9}~{}\mathrm{km}^{-1}$ at 10 km to $\sim 10^{4}~{}\mathrm{km}^{-1}$ at 85
km, while the opacities of $\nu_{e}$ or $\bar{\nu}_{e}$ are $\sim
10~{}\mathrm{km}^{-1}$ at 10 km and as low as $\sim
10^{-3}~{}\mathrm{km}^{-1}$ at 85 km. This huge gap poses a computational
challenge to properly handle all those rates. One possible approach is to
focus on a local simulation within $\mathcal{O}(10~{}\mathrm{m})$ or even
smaller range [66], but it loses the capability to account for global
advection, and thus cannot model the long-term evolution of the system self-
consistently. Another constraint comes from the minimal radial interval
$\Delta r$. Simulating in a larger radial range
$\sim\mathcal{O}(10~{}\mathrm{km})$ implies that the number of uniformly
distributed radial grids $N_{r}=L/\Delta r$ in a radial range $L=75$ km must
be $\gtrsim 10^{8}$ to properly resolve the length scale of oscillating flavor
wave imposed by $\mathbf{H}_{\nu\nu}$ developed at $\sim 30$ km (see Sec. IV).
Taking such a large number of radial grids together with other required
resolutions in angular and energy distributions is not feasible for our
extended version of COSE$\nu$ at present. Therefore, for the trilemma among
self-consistency, advection, and exact rates, at least one component must be
compromised. For this work, we take the approach of introducing artificial
attenuations to some of the rates as follows.

Figure 1:  Comparison of the scales of all the included collision rates and
different potentials in the Hamiltonian for neutrino energy $E=21$ MeV in the
transport equation [Eqs. (1) and (2)] for the background snapshot at
$t_{\mathrm{pb}}\approx 247$ ms (our Model II) from the CCSN simulation. Also
shown is $(\Delta r)^{-1}$ with the radial grid size $\Delta r=3$ m adopted in
our simulations. The black dotted curve shows the attenuated potential
$V_{\nu\nu}$ (see Sec. II.2.2 and Table 2) while the blue dotted curve shows
the parametrized matter potential $V_{\rm mat}$ (see Sec. V.4). We take
$\delta m^{2}_{\rm atm}=2.3\times 10^{-3}~{}{\rm eV}^{2}$ and $\delta
m^{2}_{\odot}=8\times 10^{-5}~{}{\rm eV}^{2}$ corresponding to the measured
values in atmospheric and solar neutrino experiments respectively.

First, we scale down all elements in neutrino self-induced term
$\mathbf{H}_{\nu\nu}$ by multiplying them with an attenuating factor [76],

$a_{\nu\nu}(r)=\frac{a_{1}}{1+e^{(a_{2}-r)/a_{3}}},$ (12)

where $a_{1}$, $a_{2}$, and $a_{3}$ are parameters, for which we take
different values for different snapshot models (see Table 2). For instance,
adopting $a_{1}=10^{-3}$, $a_{2}=35~{}{\rm km}$, and $a_{3}=3.0~{}{\rm km}$
for the snapshot shown in Fig. 1 results in attenuated $V_{\nu\nu}$ (black
dotted curve) smaller than the value of $(\Delta r)^{-1}$, when taking
$N_{r}=25000$ uniformly distributed radial grids. This form attenuates
$V_{\nu\nu}$ more at smaller radii as the neutrino number densities are
significantly higher. Notice that we also make sure our choice of
$a_{\nu\nu}(r)$ always maintains the hierarchy between $V_{\nu\nu}$ and other
collisional terms for most relevant neutrino energies. For example, the
attenuated $V_{\nu\nu}$ shown in Fig. 1 remains larger than the collisional
rates by at least a factor of 10–100 for neutrino energy of 21 MeV.

Second, for the matter potential, the high density around the neutrino sphere
also leads to too large value of $V_{\rm mat}$ to be directly included in our
simulation without modifications (see Fig. 1). Thus, for models that we
include $\mathbf{H}_{\rm mat}$ and $\mathbf{\bar{H}}_{\rm mat}$, we instead
adopt a parametrized function to explore the effect due to the presence of an
inhomogeneous matter term (see Sec. V.4).

Third, the collisional rates for neutrinos with high energy $E\gtrsim 90$ MeV
can exceed $(\Delta r)^{-1}$. To avoid that, we also apply attenuation
functions

$a_{\rm EA}(E,v_{r})=\left[1+\frac{\chi_{\rm
EA,max}^{2}(E)}{{a^{2}_{4}}}\right]^{-1/2},$ (13)

to all emissivities and opacities in EA and

$a_{\rm NNS}(E,v_{r})=\left[1+\frac{\chi_{\rm
NNS}^{2}(E)}{{a^{2}_{5}}}\right]^{-1/2},$ (14)

to the scattering kernels in NNS, where $\chi_{\rm EA,max}(E)$ is the maximal
between $j_{e}(E)+\chi_{e}(E)$ and $j_{\mu}(E)+\chi_{\mu}(E)$, $a_{4}$ and
$a_{5}$ are parametric saturation rates. We adopt values of $a_{4}$ and
$a_{5}$ close to $(\Delta r)^{-1}$ to make sure all rates are smaller than
$(\Delta r)^{-1}$ for all energies. Note that for neutrinos with energy
$\lesssim 90$ MeV, their rates are practically not attenuated since $\chi_{\rm
EA,max}/a_{4}$ and $\chi_{\rm NNS}/a_{5}$ are much smaller than 1. For those
with $E\gtrsim 90$ MeV, they remain strongly trapped during our simulation
duration.

####  II.2.3 Boundary and initial conditions

We numerically solve the QKE under spherical symmetry within a radial range
between an inner boundary $r_{i.b.}=10$ km and an outer boundary $r_{o.b.}=85$
km. Right above the inner boundary, we set up a region with a length $\sim
l_{\rm i.b.}$ and decaying width $w_{i.b.}$ wherein the EA rates for all
energies are artificially increased by the amounts $\Delta j_{i}(E)$ and
$\Delta\chi_{i}(E)$ as follows ($i=e,\mu$),

$\displaystyle\Delta j_{i}(E)$
$\displaystyle=\frac{N_{r}/L}{1+e^{(r-r_{i.b.}-l_{i.b.})/w_{i.b.}}}\frac{\varrho_{{\rm
eq},ii}(E)}{\varrho_{\rm FO}(E)+\varrho_{{\rm eq},ii}(E)},$
$\displaystyle\Delta\chi_{i}(E)$
$\displaystyle=\frac{N_{r}/L}{1+e^{(r-r_{i.b.}-l_{i.b.})/w_{i.b.}}}\frac{\varrho_{\rm
FO}(E)}{\varrho_{\rm FO}(E)+\varrho_{{\rm eq},ii}(E)},$ (15)

before we apply the attenuation factor described in the above subsection, such
that neutrinos in this zone reach equilibrium state $\varrho_{{\rm eq},ii}$
determined by the local temperature and the equilibrium neutrino chemical
potential $\mu_{\nu}^{\rm(eq)}=\mu_{e}+\mu_{p}-\mu_{n}$ within $\sim 0.5$ ms.
The purpose of having this equilibrium zone is mainly to prevent the
artificially fast leakage of $\bar{\nu}_{e}$ below $\sim 10$ MeV and
$\nu_{\mu}$ below $\sim 20$ MeV due to the lack of their main production and
opacity source of (inverse) bremsstrahlung process that are not included here
[106, 99]. At the outer boundary, we employ the free-streaming boundary
condition for forward propagating neutrinos with $v_{r}\geq 0$, and force
neutrino density matrices to be zero for back propagating ones. A different
outer boundary condition was used in Ref. [76] where neutrinos with
equilibrium number density multiplied by a dilution factor are injected
inward. We also tested this boundary condition and found negligible impact to
our results when the dilution factor is $<10^{-4}$ given that we focus on the
transport near the neutrino sphere rather than the possible halo effect [107]
at large radii.

Since we evolve the neutrino density matrices over static matter backgrounds
provided by snapshots from AGILE-BOLTZTRAN, but use a slightly different set
of neutrino reactions, we do not directly take the detailed neutrino
distributions from BOLTZTRAN and evolve them with Eqs. (1) and (2). Instead,
we adopt a two-step approach as follows. First, we take the local thermal
equilibrium state $\varrho_{\rm eq}$ and $\bar{\varrho}_{\rm eq}$ at each
radius (with nonzero values in diagonal components only) as the initial
condition, and evolve them with extended COSE$\nu$ by neglecting $\mathbf{H}$
and $\mathbf{\bar{H}}$. Typically, $\varrho$ and $\bar{\varrho}$ reach a
stationary state after a simulation time $t\sim 1$ ms. The stationary state
derived in this way automatically contains anisotropic angular distribution
near and outside the neutrino sphere due to the advection and collisions. Note
that the EA and NNS rates are all evaluated based on a given background
profile with the exception that now we assume the condition of
$\mu_{\nu_{\mu}}=0$ in the beginning instead of $Y_{\mu}=0$ (i.e.,
$\mu_{\mu}=0$). It is a good approximation to have
$\mu_{\nu_{\mu}}\approx\mu_{\mu}\approx 0$ because initially $\nu_{\mu}$ and
$\bar{\nu}_{\mu}$ are produced by pair processes and the same for $\mu^{+}$
and $\mu^{-}$.

After obtaining the stationary state distribution of $\varrho$ and
$\bar{\varrho}$ without flavor oscillations, we take this as the initial
condition for simulations fully taking into account all terms in Eqs. (1) and
(2) for each snapshot.

When the vacuum term $\mathbf{H}_{\rm vac}$ is ignored, we apply an initial
seed of neutrino flavor mixing using an Gaussian perturbation for the off-
diagonal elements of $\varrho_{\nu}$ ($\bar{\varrho}_{\nu}$) to trigger the
flavor instability

$\displaystyle\frac{\varrho^{\rm
pert}_{e\mu}(E,v_{r})}{\varrho_{ee}(E,v_{r})-\varrho_{\mu\mu}(E,v_{r})}=\frac{\bar{\varrho}^{\rm
pert}_{e\mu}(E,v_{r})}{\bar{\varrho}_{ee}(E,v_{r})-\bar{\varrho}_{\mu\mu}(E,v_{r})}$
$\displaystyle=$ $\displaystyle
10^{-3}\exp\left[-\left(\frac{r~{}{\rm[km]}-47.5}{10}\right)^{2}\right].$ (16)

We also examine several cases that include $\mathbf{H}_{\rm vac}$ (see Sec.
V.3). For those cases, because $\mathbf{H}_{\rm vac}$ generates flavor mixing
automatically, no artificial perturbations are given initially.

####  II.2.4 List of evolved models

We take four hydrodynamical snapshot profiles at different post-bounce times
$t_{\rm pb}=144$ ms, 247 ms, 503 ms, and 1000 ms from the BOLTZTRAN SN
simulation, labeled by Models I–IV in Table 2. For all models, we evolve
$\varrho$ and $\bar{\varrho}$ between $r_{\rm i.b.}=10$ km and $r_{\rm
o.b.}=85$ km with the total radial range $L=75$ km. For each snapshot, we
perform three simulations using different radial resolutions. Our fiducial
models listed in Table 2 are high resolution ones with $N_{r}=25000$.
Comparison to results obtained with smaller $N_{r}=10000$ and 2500 are given
in Appendix B. For angular grids, we take $v_{r}$ uniformly between $-1$ and
$1$ with $N_{v_{r}}=50$. We have also examined cases with $N_{v_{r}}=200$ for
the lowest radial resolution runs and found that the results are nearly
identical (see Appendix B). We take $N_{E}=20$ energy grids between 2 MeV to
160 MeV spaced nearly uniformly in logarithmic scale.

Since our main goal is to examine the flavor conversions triggered by the
collisional instability without the interference of other types of flavor
instability, we first only include in our fiducial models the neutrino self-
induced term, the full form of EA, and the diagonal elements of
$\mathbf{C}_{\rm NNS}$ for the NNS collisions, which are the minimal
ingredients to induce the collisional instability. For Model II, we perform
additional simulations including $\mathbf{H}_{\rm vac}$ and $\mathbf{H}_{\rm
m}$ separately and study their effects. Moreover, we introduce two auxiliary
parameters $b_{\rm EA}$ and $b_{\rm NNS}$, which are multiplied to all off-
diagonal elements of $\mathbf{C}_{\rm EA}$ and $\mathbf{C}_{\rm NNS}$.
Different values of $b_{\rm EA}=1,2,4$ and $b_{\rm NNS}=0,1$ are taken to test
possible impact of the size of the off-diagonal elements in $\mathbf{C}_{\rm
EA}$ and $\mathbf{C}_{\rm NNS}$ on the results. Table 2 summarizes all models
presented in this paper. Each simulation takes $\sim\mathcal{O}(30000)$ CPU
hours for a duration of simulation time $t\sim 1$ ms.

Table 2:  Parameters used in each model. For all models, the radial range is from 10 km to 85 km, $l_{\rm i.b.}=6$ km, $w_{i.b.}=0.2$ km, $N_{r}=25000$, $N_{v_{r}}=50$, $N_{E}=20$, $a_{4}=1/\Delta r$ [km-1] , $a_{5}=1/(2\Delta r)$ [km-1] and $\theta_{V}=10^{-6}$. model |  $t_{\mathrm{pb}}$ [ms] | $a_{1}$ | $a_{2}~{}\mathrm{[km]}$ | $a_{3}~{}{\rm[km]}$ | $b_{\mathrm{EA}}$ | $b_{\mathrm{NNS}}$ |  $\delta m^{2}$ [eV2] |  $V_{\mathrm{mat}}(r~{}\mathrm{[km]})$ [km-1]
---|---|---|---|---|---|---|---|---
I | 144 | $10^{-3}$ | 45 | 4.5 | 1 | 0 | 0 | 0
II | 247 | $10^{-3}$ | 35 | 3 | 1 | 0 | 0 | 0
III | 503 | $10^{-3}$ | 25 | 1.5 | 1 | 0 | 0 | 0
IV | 1000 | $10^{-3}$ | 21 | 1.1 | 1 | 0 | 0 | 0
IIn | 247 | $10^{-3}$ | 35 | 3 | 1 | 1 | 0 | 0
IIe1 | 247 | $10^{-3}$ | 35 | 3 | 2 | 0 | 0 | 0
IIe2 | 247 | $10^{-3}$ | 35 | 3 | 4 | 0 | 0 | 0
IIv1 | 247 | $10^{-3}$ | 35 | 3 | 1 | 0 | $8\times 10^{-5}$ | 0
IIv2 | 247 | $10^{-3}$ | 35 | 3 | 1 | 0 | $2.3\times 10^{-3}$ | 0
IIm | 247 | $10^{-3}$ | 35 | 3 | 1 | 0 | 0 | $(500/3)\cdot\exp\left[-(r-10)^{4}/18^{4}\right]$

##  III Linear stability analysis

We perform the linear stability analysis in the regime where the off-diagonal
elements of density matrices are small and can be treated as perturbations
compared to the number densities, i.e.,
$|\varrho_{e\mu}|/|\varrho_{ee}-\varrho_{\mu\mu}|\ll 1$ and
$|\bar{\varrho}_{e\mu}|/|\varrho_{ee}-\varrho_{\mu\mu}|\ll 1$ [108, 109].
Neglecting the vacuum and matter Hamiltonian contribution as well as the
aberration term, we further assume that a collective mode of the perturbation
$\varrho_{e\mu}=Q(\Omega,K_{r},E,r,v_{r})e^{-i[\Omega t-K_{r}(r^{\prime}-r)]}$
and $\bar{\varrho}_{e\mu}^{*}=\bar{Q}(\Omega,K_{r},E,r,v_{r})e^{-i[\Omega
t-K_{r}(r^{\prime}-r)]}$ can develop locally near $r$. With these assumptions,
the off-diagonal parts of Eqs. (1) and (2) become

$\displaystyle\left[\Omega-K_{r}v_{r}-\Phi(v_{r})+iC_{\rm
e\mu,EA}(E)\right]Q(\Omega,K_{r},E,v_{r})=$
$\displaystyle-\sqrt{2}G_{F}[\varrho_{ee}(E,v_{r})-\varrho_{\mu\mu}(E,v_{r})]\int
dE^{\prime}\,dv_{r}^{\prime}\,\times$
$\displaystyle(1-v_{r}v_{r}^{\prime})[Q(\Omega,K_{r},E^{\prime},v_{r}^{\prime})-\bar{Q}(\Omega,K_{r},E^{\prime},v_{r}^{\prime})],$
(17)

and

$\displaystyle\left[\Omega-K_{r}v_{r}-\Phi(v_{r})+i\bar{C}_{\rm
e\mu,EA}(E)\right]\bar{Q}(\Omega,K_{r},E,v_{r})=$
$\displaystyle-\sqrt{2}G_{F}[\bar{\varrho}_{ee}(E,v_{r})-\bar{\varrho}_{\mu\mu}(E,v_{r})]\int
dE^{\prime}\,dv_{r}^{\prime}\,\times$
$\displaystyle(1-v_{r}v_{r}^{\prime})[Q(\Omega,K_{r},E^{\prime},v_{r}^{\prime})-\bar{Q}(\Omega,K_{r},E^{\prime},v_{r}^{\prime})],$
(18)

where $\Phi(v_{r})=\sqrt{2}G_{F}\int
dE^{\prime}\,dv_{r}^{\prime}\,(1-v_{r}v_{r}^{\prime})[\varrho_{ee}(E^{\prime},v_{r}^{\prime})-\varrho_{\mu\mu}(E^{\prime},v_{r}^{\prime})-\bar{\varrho}_{ee}(E^{\prime},v_{r}^{\prime})+\bar{\varrho}_{\mu\mu}(E^{\prime},v_{r}^{\prime})]$,
$C_{\rm e\mu,EA}(E)=[j_{e}(E)+\chi_{e}(E)+j_{\mu}(E)+\chi_{\mu}(E)]/2$, and
$\bar{C}_{\rm
e\mu,EA}(E)=[\bar{j}_{e}(E)+\bar{\chi}_{e}(E)+\bar{j}_{\mu}(E)+\bar{\chi}_{\mu}(E)]/2$.
The eigenvalues of $\Omega$ can be numerically derived for a given wave number
$K_{r}$ using the same discretization scheme for both $E$ and $v_{r}$ as in
the simulations. Solutions derived in this way will be shown in the next
section to help understand the simulation outcome.

Another useful simplification that one can take is to consider the limit where
the environment is homogeneous and all neutrinos are monochromatic. Taking
these assumptions, the linearized equations above can be further simplified
into

$\displaystyle(\Omega+iC_{\rm e\mu,EA})Q(\Omega)=$
$\displaystyle-\sqrt{2}G_{F}[\Delta n_{\bar{\nu}}Q(\Omega)-\Delta
n_{\nu}\bar{Q}(\Omega)],$ $\displaystyle(\Omega+i\bar{C}_{\rm
e\mu,EA})\bar{Q}(\Omega)=$ $\displaystyle-\sqrt{2}G_{F}[\Delta
n_{\bar{\nu}}Q(\Omega)-\Delta n_{\nu}\bar{Q}(\Omega)].$ (19)

The corresponding eigenvalue of $\Omega$ therefore satisfies a quadratic
secular equation

$\displaystyle\Omega^{2}+[i(C_{\rm e\mu,EA}+\bar{C}_{\rm
e\mu,EA})+\sqrt{2}G_{F}(\Delta n_{\bar{\nu}}-\Delta n_{\nu})]\Omega$
$\displaystyle=$ $\displaystyle i\sqrt{2}G_{F}(C_{\rm e\mu,EA}\Delta
n_{\nu}-\bar{C}_{\rm e\mu,EA}\Delta n_{\bar{\nu}})+C_{\rm e\mu,EA}\bar{C}_{\rm
e\mu,EA},$ (20)

where $\Delta n_{\nu}=n_{\nu_{e}}-n_{\nu_{\mu}}$ and $\Delta
n_{\bar{\nu}}=n_{\bar{\nu}_{e}}-n_{\bar{\nu}_{\mu}}$. The imaginary part of
analytical solutions can be obtained as

$\displaystyle{\rm Im}$ $\displaystyle(\Omega)\approx-\frac{|C_{\rm
e\mu,EA}+\bar{C}_{\rm e\mu,EA}|}{2}\pm\frac{|C_{\rm e\mu,EA}-\bar{C}_{\rm
e\mu,EA}|}{2}\times$ $\displaystyle\frac{|\Delta n_{\nu}+\Delta
n_{\bar{\nu}}|}{|\Delta n_{\nu}-\Delta
n_{\bar{\nu}}|}\left[1-\frac{1}{4G_{F}^{2}}\left(\frac{C_{\rm
e\mu,EA}-\bar{C}_{\rm e\mu,EA}}{\Delta n_{\nu}-\Delta
n_{\bar{\nu}}}\right)^{2}\right]$ (21)

to the order of $\mathcal{O}[C^{3}_{\rm e\mu,EA}/(G_{F}\Delta n_{\nu})^{2}]$.
When $\sqrt{2}G_{F}(\Delta n_{\nu}-\Delta n_{\bar{\nu}})$, i.e. $V_{\nu\nu}$,
becomes much larger than the collisional rates, the high-order term is
negligible and approaches to an asymptotic value. Although this simplified
expression cannot replace Eqs. (III) and (III) to precisely predict the growth
rates, it can help us quantitatively understand the general feature of the
neutrino system. Whether there exists a runaway solution roughly depends on
the relation between the asymmetry factor of EA rates $\alpha_{C}\equiv|C_{\rm
e\mu,EA}-\bar{C}_{\rm e\mu,EA}|/|C_{\rm e\mu,EA}+\bar{C}_{\rm e\mu,EA}|$ and
that of neutrino number densities $\alpha_{n}\equiv|\Delta n_{\nu}-\Delta
n_{\bar{\nu}}|/|\Delta n_{\nu}+\Delta n_{\bar{\nu}}|$ [71]. Clearly, when
ignoring higher order terms in Eq. (21), larger asymmetry in the EA rates and
smaller asymmetry in neutrino number densities can result in a positive
imaginary value of $\Omega$, which leads to an instability. The criteria of
having the instability is usually satisfied near the neutrino sphere. The
reason is twofold. First, around the neutrino sphere the $Y_{e}$ is usually
low ($Y_{e}\lesssim 0.15$) such that the $\nu_{e}$ EA rates can be much larger
than that of $\bar{\nu}_{e}$ by a factor of $\approx 5$–10. Second, the
$\nu_{\mu}$ and $\bar{\nu}_{\mu}$ number densities are usually higher than
that of $\bar{\nu}_{e}$ well inside the neutrino sphere but become smaller
around the sphere, due to the diffusion of $\nu_{\mu}$. At radii where
$n_{\nu_{\mu}}\simeq n_{\bar{\nu}_{\mu}}\lesssim n_{\bar{\nu}_{e}}$, the
asymmetry factor in neutrino number densities becomes the largest. Thus, a
combination of these two facts allows the instability generally exist around
the neutrino sphere.

We caution that the factors discussed above are not necessarily the only ones
determining how the collisional instability to induce flavor conversions. This
is because in the stability analysis presented above, we ignored the
aberration term and assume the neutrino number density is locally homogeneous.
To understand the role of these dynamical properties, it requires a complete
simulation as we will discuss in the next two sections. However, the linear
stability analysis still provides some insights to the underlying mechanism
that drives the flavor instability as will be discussed further below.

##  IV Results of collisional instability

For all four fiducial models (I–IV) that represent different evolution stages
of the CCSN model, we observe the occurrence of collisional instability. The
features in Model II, III, and IV are qualitatively similar while Model I
shows a unique behavior. Below, we begin by discussing results obtained in
Model II. We then address the similarities and differences of the models.

For this section, we define energy-integrated density matrices
$\langle\varrho\rangle_{E}(v_{r})=\int dE\,\varrho(E,v_{r})$, angle-integrated
density matrices $\langle\varrho\rangle_{A}(E)=\int dv_{r}\,\varrho(E,v_{r})$,
and neutrino mean energy

$\langle E_{\nu_{i}}\rangle=\frac{\int
dE\,dv_{r}\,E\varrho_{ii}(E,v_{r})}{\int dE\,dv_{r}\,\varrho_{ii}(E,v_{r})},$
(22)

where $i=e,\,\mu$ for analysis.

###  IV.1 Model II

We first show in Fig. 2 the energy-integrated diagonal and off-diagonal
elements of neutrino density matrices $\langle\varrho_{\mu\mu}\rangle_{E}$
(left panels) and $|\langle\varrho_{e\mu}\rangle_{E}|$ (middle panels), as
well as the dimensionless ratio
$s_{e\mu}\equiv|\langle\varrho_{e\mu}\rangle_{E}|/|\langle\varrho_{ee}\rangle_{E}-\langle\varrho_{\mu\mu}\rangle_{E}|$
(right panels), as functions of $r$ and $v_{r}$ taken at different times from
the COSE$\nu$ simulation for Model II. Both the left and middle panels show
the anisotropic shape of $\langle\varrho_{\mu\mu}\rangle$ and
$|\langle\varrho_{e\mu}\rangle|$ at larger radii due to the nature of neutrino
free streaming. However, the ratio $s_{e\mu}$ in right panels reveals that the
development and evolution of the flavor off-diagonal elements in fact depend
on $v_{r}$ rather weakly. Clearly, the flavor instability existing around
$r=25$-30 km in the beginning leads to the rapid growth of
$|\langle\varrho_{e\mu}\rangle_{E}|$ to the nonlinear regime within $\sim
0.008$ ms, indicated by $s_{e\mu}\gtrsim 0.1$. Afterwards, $s_{e\mu}$ stops
growing in magnitude and splits into two branches that move in two opposite
directions. The inward moving one gets damped by the interaction with denser
medium and eventually disappears at $\sim 0.04$ ms. The outward going branch
propagates with a group velocity of $\sim 0.4$ and further bifurcates into
another two sub-branches at $\sim 36$ km around $t\sim 0.08$ ms. The inner
sub-branch moves at a similar group velocity of $\sim 0.4$ and disappears at
$t\sim 0.5$ ms. On the other hand, the outer one propagates with a group
velocity of $\sim 0.9$, quick enough to reach the free-streaming region where
flavor mixing gets transported away from our simulation domain for positive
$v_{r}$ modes. Afterwards, the neutrino fields remain stable (against
oscillations) for a duration of $\sim 0.2$ ms during 0.5 ms$\lesssim t\lesssim
0.7$ ms. At $t\approx 0.7$ ms, the collisional instability appears once again
around $r\sim 28$ km after the neutrino field self-regulates its distributions
near the decoupling region. Unlike the first instability discussed above, the
flavor mixing due to the second collisional instability does not get
transported away this time and eventually freezes into a stationary state
until the end of our simulation at $t\sim 1$ ms.

Figure 2:  Evolution of energy-integrated number density
$\langle\varrho_{\mu\mu}\rangle_{E}$ (left panel), flavor mixing
$|\langle\varrho_{e\mu}\rangle_{E}|$ (middle panel), and the dimensionless
ratio
$s_{e\mu}\equiv|\langle\varrho_{e\mu}\rangle_{E}|/|\langle\varrho_{ee}\rangle_{E}-\langle\varrho_{\mu\mu}\rangle_{E}|$
(right panel) for Model II. Figure 3:  Growth rates Im$(\Omega)$ from linear
stability analysis as functions of radius $r$ for four simulation times in
Model II. We sample 31 values of $K_{r}$ from -2 km-1 to 2 km-1 and show the
maximal values for each radius. The red curve indicates $5/r$ as an empirical
criteria to determine the growth of instability against advection.  Figure 4:
Time evolution of the dimensionless ratio $s_{e\mu}$ of radial velocity
$v_{r}=1$ for four different radii in Model II. Each of them is compared with
a black dashed line determined by $s_{e\mu}(t=0)\exp[{\rm Im}(\Omega)t]$ with
the growth rate Im$(\Omega)$ in Fig. 3 at $t=0$ ms.

Figure 3 shows the maximum growth rates of flavor instabilities among all
$K_{r}$ modes, Im$(\Omega)$, as a function of radius for Model II at different
times. These Im$(\Omega)$ are derived by numerically solving the linearized
Eqs. (III) and (III). At $t=0$ ms, Im$(\Omega)$ peaks at $r\simeq 27$ km. The
growth of this instability thus dominates the evolution of the system
initially, consistent with results shown in Fig. 2. When $t=0.024$ ms after
flavor transformation occurs, Im$(\Omega)$ in 25 km$\lesssim r\lesssim 32$ km
becomes smaller than $\sim 0.1$ km-1, while its value maintain roughly the
same for $r\gtrsim 32$ km. For $t=0.42$ ms when flavor mixings around $r\sim
30$ km gets suppressed, Im$(\Omega)<0.1$ km-1 for all radii. At an even later
time $t=0.66$ ms when flavor conversion reappears (see Fig. 2), larger
Im$(\Omega)$ are found in 27 km$\lesssim r\lesssim 31$ km again.

We compare the maximum growth rates obtained from the stability analysis at 0
ms with the numerical evolution for Model II at four different radii. Figure 4
shows that the time evolution of $s_{e\mu}$ of radial velocity $v_{r}=1$ in
the linear regime perfectly agree with the prediction determined by
$s_{e\mu}(t=0)\exp[{\rm Im}(\Omega)t]$ at 28, 30, and 32 km respectively,
which is expected since the growth of collisional instability dominates over
the disturbance from advection.

Although positive Im$(\Omega)$ are found for nearly all radii larger than 20
km at all times, not all of them lead to the growth of
$|\langle\varrho_{e\mu}\rangle|$ in the simulation. This is because the
stability analyses can only tell how a perturbation evolves around where the
local condition can be maintained. However, in realistic simulations where
advection occurs in the presence of inhomogenous neutrino number density, the
instability growth rate needs to compete with advection for a perturbation to
grow before it being transported away. To illustrate this, we show in Fig. 2 a
characteristic value of advection rate as $5/r$ for Model II by the red solid
line. We take this function in an empirical way motivated by the advection
term in Eqs. (1) and (2) being generally proportional to $1/r$. Comparing
Figs. 2 and 3, it seems to suggest that when Im$(\Omega)$ is roughly less than
$5/r$, the growth rate of the instability is too small against the advection,
such that no significant flavor conversion can develop. For example, although
the stability analysis yields a positive Im$(\Omega)$ at $r=40$ km and $t=0$
ms, $s_{e\mu}$ decreases in the simulation as shown in Fig. 4. In addition,
Im$(\Omega)$ are positive from 0.024 ms to 0.42 ms at $r=27$ km in Fig. 3,
while $s_{e\mu}$ decreases in Fig. 2.

Next, we examine the impact of flavor conversion due to collisional
instability on the property of neutrinos of all flavors. Figure 5 shows the
angular-integrated neutrino energy spectra $\langle\varrho_{\nu}\rangle_{A}$
for two radii $r=32$ and $55$ km at three different simulation times $t=0.0$,
0.02, and $0.16$ ms. For $r=32$ km close to where the initial flavor
instability occurs, the maximum amount of flavor conversion happens when
$t\sim 0.02$ ms. At this time, the $\nu_{\mu}$ spectrum is significantly
shifted towards low energy compared to the one at $t=0$, because high-energy
$\nu_{\mu}$ are converted to $\nu_{e}$ while low-energy $\nu_{\mu}$ receive
contribution from flavor conversion of $\nu_{e}$ in return. However, the
spectra of both $\nu_{e}$ and $\bar{\nu}_{e}$ are only affected marginally
(reduced by $\sim 5\%$–$8\%$ near 15 MeV). This is mainly because the
relatively large EA processes [(1a)-(1c) in Table 1] that act on $\nu_{e}$ and
$\bar{\nu}_{e}$ quickly bring them to a state close to thermal equilibrium. At
later times after the flavor mixing gets transported to outer radii, e.g.,
$t=0.16$ ms shown in Fig. 5, the EA processes restore the $\nu_{e}$ and
$\bar{\nu}_{e}$ spectra fully back to the equilibrium state, while the
advection of $\nu_{\mu}$ from denser region gradually shifts the peak of its
spectrum to higher energy again. For spectra at a larger radius of $r=55$ km,
the $\nu_{\mu}$ spectrum remains identical as the initial one in the beginning
(e.g., at $t=0.02$ ms) since there is no flavor instability at this radius.
Differences only show up at later times (e.g., at $t=0.16$ ms) when the
$\nu_{\mu}$ affected by collisional instability propagate here from inner
radii, leading to a clearly visible shift of $\nu_{\mu}$ spectra to low
energy. For $\nu_{e}$ and $\bar{\nu}_{e}$, EA collisions wipe out all the
effects due to flavor conversions such that their spectra remain unaffected by
collisional instability at all.

(a)

(b)

Figure 5:  Angle-integrated energy spectrum for $\nu_{e}$(solid),
$\bar{\nu}_{e}$ (dashed), and $\nu_{\mu}$ (dash-dotted) at two radii [32 km in
panel (a) and 55 km in panel (b)]. Blue, orange, and green colors correspond
to three different simulation times 0.0 ms, 0.02 ms, and 0.16 ms respectively.
Note that in panel (b), the blue curve of overlaps with the orange curve for
$\nu_{\mu}$, while the blue and orange curves overlap with the green curve for
$\nu_{e}$ and $\bar{\nu}_{e}$.

Our results here suggest that there exists a fundamental difference in models
assuming spatial homogeneity (e.g., Ref. [71]) and more realistic ones
including advection. For a model where spatial homogeneity is assumed and
advection is neglected, flavor conversion due to collisional instability can
lead to $\nu_{\mu}$ spectrum identical to that of $\nu_{e}$ or $\bar{\nu}_{e}$
given enough time. This is because in homogeneous models, $\nu_{\mu}$
(decoupled from medium collisionally) continuously acquires supply from
electron neutrinos due to flavor conversion while electron flavors can be
repopulated by EA processes, so that their spectra becomes identical given
long enough time. However, similar to what argued above regarding the growth
of instability, in a realistic environment, the advection introduces another
timescale so that the transport of $\nu_{\mu}$ away from where the instability
occurs can prevent their spectra become identical to that of $\nu_{e}$.

In addition to the energy spectra, we examine two energy-integrated
quantities, the neutrino number density and mean energy, to characterize the
transport of collision-induced neutrino flavor conversion. Figure 6 shows the
radial profile of these two quantities at different times. The $\nu_{e}$ and
$\bar{\nu}_{e}$ number densities are barely affected except for the tiny dip
in $\nu_{e}$ around 31 km at $t=0.02$ ms. For $\nu_{\mu}$, slight changes of
$n_{\nu_{\mu}}$ due to flavor conversion around 32 km are present. More
clearly visible is the larger enhancement of $n_{\nu_{\mu}}$ at $t=0.16$ ms
for $r\gtrsim 35$ km due to the extra contribution from the converted ones at
the low energy part of the spectrum discussed above. Although the impact on
the number densities seem to be small, flavor conversion and transport due to
collisional instability does lead to significant changes to the mean neutrino
energies, particularly for $\nu_{\mu}$. Figure 6(b) shows that a large
decrease of $\langle E_{\nu_{\mu}}\rangle$ from $\simeq 40$ MeV locally around
$r\simeq 30$ km to as low as $\lesssim 30$ MeV at $0.02$ ms due to flavor
conversion. Once again, this effect gets transported outward and leaves behind
a lowered $\langle E_{\nu_{\mu}}\rangle$ for $r\gtrsim 25$ km at e.g.,
$t=0.16$ ms.

(a)

(b)

Figure 6:  Radial profiles of neutrino number densities [panel (a)] and mean
energies [panel (b)] for $\nu_{e}$(solid), $\bar{\nu}_{e}$ (dashed), and
$\nu_{\mu}$ (dash-dotted) at two radii (32 km and 55 km). Blue, orange, and
green colors correspond to three different simulation times 0.0 ms, 0.02 ms,
and 0.16 ms respectively.

###  IV.2 Models I, III, and IV

After examining the simulation results in Model II, we now discuss the
similarity and differences obtained in Model I, III, and IV, based on
different background SN profiles taken from the BOLTZTRAN simulation. Figure 7
shows the evolution of $s_{e\mu}$ for these three models. For Model I, the
collisional instabilities are not strong enough at the beginning so the
Gaussian perturbation simply gets transported outward. At $t\sim 0.2$ ms, an
unstable mode grows around $r\sim 38$ km. The growth rate is significantly
lower than that in Model II since the dimensionless ratio reaches at maximum
$\sim 5\times 10^{-3}$ at $t=0.5$ ms and $\sim 5\times 10^{-2}$ at $t=0.74$
ms.

Figure 7:  Evolution of the dimensionless ratio $\log_{10}(s_{e\mu})$ in
models I (left panel), III (middle panel) and IV (right panel) at different
simulation times. The lower-right corners in the middle and right panels have
too small values of
$|\langle\varrho_{ee}\rangle_{E}-\langle\varrho_{\mu\mu}\rangle_{E}|<10^{29}~{}{\rm
cm}^{-3}$ and hence not shown.

For Models III and IV, the flavor conversion and transport behave similarly to
Model II. Collisional instabilities develop immediately around 20–25 km and
18–22 km in these two models, respectively. Afterwards, they grow to nonlinear
regime and bifurcate. Once again, the inward going mode gets suppressed and
the outward going mode sustains and propagate to larger radii. Also, second
onsets of flavor instability happen in both models at $t\simeq 0.36$ ms and
$t\simeq 0.2$ ms and seem to reach stationary states without propagating
outwards at the end of our simulations.

As in Model II, the effect of flavor conversions on radial profile of number
densities in all three models are limited and we do not present them
explicitly. For the neutrino mean energies, we show results from Model III and
IV in Fig. 8. Similar to Fig. 6, the collisional instabilities lead to
significant reductions of $\langle E_{\nu_{\mu}}\rangle$ by
$\sim\mathcal{O}(10)$ MeV and slightly increases $\langle E_{\nu_{e}}\rangle$
and $\langle E_{\bar{\nu}_{e}}\rangle$ at where the instabilities occur. The
changes in $\nu_{\mu}$ are carried to large radii, while strong EA collisions
of $\nu_{e}$ and $\bar{\nu}_{e}$ refrains their local changes from propagating
outwards.

(a)

(b)

Figure 8:  Radial profiles of neutrino mean energies for $\nu_{e}$(solid),
$\bar{\nu}_{e}$ (dashed), and $\nu_{\mu}$ (dash-dotted) in model III [panel
(a)] and model IV [panel (b)]. Blue, orange, and green colors correspond to
three different simulation times 0.0 ms, 0.02 ms, and 0.08 ms respectively.

The general behaviors in those four models can be understood with the simple
criteria of asymmetries in neutrino number densities and the collision rates
in Eq. (21). Given that $C_{\rm e\mu,EA}/\bar{C}_{\rm
e\mu,EA}\approx\chi_{e}/\bar{\chi}_{e}\approx(1-Y_{e})/Y_{e}$, the asymmetry
factor of EA rate is $\alpha_{C}\simeq|1-2Y_{e}|$. Neglecting the higher order
correction, Eq. (21) predicts a local instability when $\alpha_{C}>\alpha_{n}$
with a growth rate proportional to $|C_{\rm e\mu,EA}+\bar{C}_{\rm
e\mu,EA}|\times(\alpha_{C}/\alpha_{n}-1)$.

Figure 9 compares $|1-2Y_{e}|\simeq\alpha_{C}$ (dashed curves) to $\alpha_{n}$
(solid curves) and shows the corresponding $Y_{e}$ profiles in the bottom
panel. The deleptonization in CCSN leads to more neutron-rich condition near
the neutrinosphere at the later time so that the maximal value of
$|1-2Y_{e}|\simeq\alpha_{C}$ increases with time. For $\alpha_{n}$, it is
larger than 1 in denser region where
$n_{\nu_{e}}>n_{\nu_{\mu}}>n_{\bar{\nu}_{e}}$ but decreases with radius and
asymptotic to a value smaller than 1 for all four snapshots. More importantly,
the asymptotic $\alpha_{n}<1$ is lower at a later SN snapshot considered. This
is because the SN evolves from the initial stage of neutronization burst
dominated by $\nu_{e}$ emission toward the accretion phase where the
difference between $\nu_{e}$ and $\bar{\nu}_{e}$ becomes smaller, thus leading
to smaller values of the asymptotic $\alpha_{n}$. Consequently,
$\alpha_{C}/\alpha_{n}>1$ become larger around the neutrinosphere in a later
SN profile. Moreover, the collision rates also increases as the matter density
around the neutrinosphere increases over time. Combining all these factors
described above, it is clear that later SN stages contain more favorable
conditions for the collisional instabilities to grow against advection,
consistent with our simulation results.

Figure 9:  Profiles for the asymmetry factor $\alpha_{n}$ (solid),
$|1-2Y_{e}|\simeq\alpha_{C}$ (dashed) in the upper panel, and $Y_{e}$ in the
lower panel at four CCSN snapshots. An approximate condition for collisional
instability is $\alpha_{C}>\alpha_{n}$.

Our results obtained above in Models I–IV tend to suggest that the collisional
instability in spherically symmetric SN models can affect the properties of
heavy lepton flavor neutrinos at and beyond their decoupling region, but does
not lead to major impact on $\nu_{e}$ and $\bar{\nu}_{e}$ spectra above their
neutrinospheres. Thus, the effect of collisional instability are more likely
to manifest in the CCSN dynamics and the emitted neutrino signals, and may
moderately affect neutrino (induced) nucleosynthetic processes that are
sensitive to the spectra of heavy-lepton neutrino flavors. For the condition
of nucleosynthesis in the neutrino driven wind, the impact of collisional
instability may be minor since the outcome depend more on the properties of
$\nu_{e}$ and $\bar{\nu}_{e}$.

##  V Effects of other terms on collisional instability

As mentioned earlier, the typical timescale for collisional instability is in
similar order as the inverse of collisional rates. This makes this problem
different from e.g., the fast flavor conversion for which effects from the
vacuum mixing, matter inhomogeneity, and collisions may be ignored to a good
approximation for local simulations [51, 28, 52, 54, 56, 57]. In the previous
section, we ignored the vacuum term $\mathbf{H}_{\rm vac}$, the matter term
$\mathbf{H}_{\rm mat}$, and the impact of the NNS collision on the off-
diagonal elements of the neutrino density matrices, for the purpose of purely
investigating the outcome of collision instability. In this section, we
explore consequences of including these terms each by each, as well as effects
due to the artificially enhanced decoherence in EA based on the SN background
used for Model II. These six additional models are listed in the second part
of Table 2. The simulation results for all of them are summarized in Figs. 10
and 11: The former shows the radial profile of
$|\langle\varrho_{e\mu}\rangle_{E}|$ for the angular mode $v_{r}=1$ at
different times, while the latter shows the mean energies of $\nu_{\mu}$ as a
function of radius at the same simulation times. Below, we discuss the results
of the impact due to the NNS induced decoherence, the artificially enhanced EA
decoherence, the matter term, and the vacuum term in each of the subsections.

(a) IIn
(b) IIe1
(c) IIe2
(d) IIv1
(e) IIv2
(f) IIm

Figure 10:  Radial profiles of dimensionless ratio $s_{e\mu}$ of radial
velocity $v_{r}=1$ in exploratory models IIn, IIe1, IIe2, IIv1, IIv2, and IIm
[from panel (a) to (f) respectively]. The simulation times are indicated by
different colors as listed in the legend. Results from model II for each
simulation time are presented by transparent and dashed curves for comparison.

(a) IIn
(b) IIe1
(c) IIe2
(d) IIv1
(e) IIv2
(f) IIm

Figure 11:  Radial profiles of neutrino mean energies for $\nu_{e}$(solid),
$\bar{\nu}_{e}$ (dashed), and $\nu_{\mu}$ (dash-dotted) in exploratory models
IIn, IIe1, IIe2, IIv1, IIv2, and IIm [from panel (a) to (f) respectively]. The
simulation times are indicated by different colors as listed in the legend.
Results from model II are presented in more transparent colors for each
simulation time for comparison.

###  V.1 NNS decoherence

We investigate the decoherence effect due to the NNS by setting the parameter
$b_{\rm NNS}=1$ in Model IIn, which allows nonzero off-diagonal elements in
$\mathbf{C}_{\rm NNS}$ and $\bar{\mathbf{C}}_{\rm NNS}$ that were artificially
suppressed in Model I–IV. Figure 10(a) shows that the evolution during the
initial phase of $t\lesssim 0.16$ ms is extremely insensitive to the
decoherence of NNS as the curves from Model II and Model IIn overlap.
Deviations start to appear at $r\gtrsim 36$ km after $t\gtrsim 0.04$ ms when
the flavor mixing propagate outwards. The NNS decoherence effect in Model IIn
suppresses the amount of flavor mixing by a factor of $\sim 10$ for all
velocity modes than the result obtained in Model II. The main reason that the
NNS decoherence effect only shows up at later times at larger radii is related
to the angular distribution of neutrinos at locations where flavor mixing
occurs. At $r\sim 30$ km (inside the neutrinosphere) where initially the
instability grows, $\varrho_{e\mu}$ are nearly isotropic (see Fig. 2) such
that the NNS does not affect their evolution [see Eq. (9)]. However, at larger
radii where $\varrho_{e\mu}$ becomes anisotropic in $v_{r}$, the NNS
decoherence term not only isotropizes $\varrho_{e\mu}$, but also leads to
damping of $\varrho_{e\mu}$ due to the interplay with $H_{\nu\nu}$. Since the
NNS decoherence term mainly affects $\varrho_{e\mu}$, including this term does
not lead to any changes to the averaged neutrino energy shown in Fig. 11(a).

###  V.2 Enhanced EA decoherence

The next models IIe1 and IIe2 assume artificially enhanced off-diagonal
elements in $\mathbf{C}_{\rm EA}$ and $\bar{\mathbf{C}}_{\rm EA}$ to examine
how the results depend on the size of these terms, which are the major source
of the collisional instability. We take the enhancement factors $b_{\rm EA}=2$
in Model IIe1 and $b_{\rm EA}=4$ in Model IIe2. Panel (b) and (c) in Fig. 10
show that for most radii, $s_{e\mu}$ at $0.004$ ms in Models IIe1 and IIe2 are
similar to the values at $0.008$ ms and 0.016 ms in Model II, respectively.
This is consistent with the conclusion obtained with the LSA in Sec. III,
because only the off-diagonal elements of $C_{e\mu,\rm EA}$ and
$\bar{C}_{e\mu,\rm EA}$ enter the linearized equation.

After the instability grows to the linear regime and propagate outwards, the
evolution of flavor mixing are also clearly affected by the enhancement.
Overall, the mixing propagates faster with larger values of $b_{\rm EA}$ but
also gets damped earlier (see e.g., the curves at 0.08 and 0.16 ms). The
effect of faster propagation is also clearly shown in the evolution of
$\langle E_{\nu_{\mu}}\rangle$ shown in panels (b) and (c) of Fig. 11.

###  V.3 Vacuum term

Models IIv1 and IIv2 include the vacuum Hamiltonian ${\mathbf{H}}_{\rm vac}$.
In both models, we assume an effectively reduced mixing $\theta_{V}=10^{-6}$.
For $\delta m^{2}$, we take different values of $\delta m^{2}=\delta
m^{2}_{\odot}=8\times 10^{-5}~{}{\rm eV}^{2}$ for Model IIv1 and $\delta
m^{2}=\delta m^{2}_{\rm atm}=2.3\times 10^{-3}~{}{\rm eV}^{2}$ for Model IIv2,
which correspond to the measured values in solar and atmospheric neutrino
experiments, respectively. For these models, we do not apply any initial
perturbation in the off-diagonal terms of $\varrho$ and $\bar{\varrho}$ as
flavor mixing is directly generated by ${\mathbf{H}}_{\rm vac}$.

Panels (d) and (e) of Fig. 10 show that the collisional instabilities develop
in both models and reach the nonlinear regime within a similar timescale
$\sim\mathcal{O}(0.01~{}{\rm ms})$, similar to the case using artificial
initial Gaussian perturbation. Because Model IIv2 has larger $\delta m^{2}$
than Model IIv1, which generates larger flavor mixing seed, $s_{e\mu}$ grows
faster in Model IIv2 in the linear regime.

For both models, we find that the inclusion of $H_{\rm vac}$ tends to help
amplify flavor mixing when they propagate outwards, as shown by the larger
value of $s_{e\mu}$ around $r\simeq 38$ km at $t=0.048$ ms (purple line). As a
result, the corresponding $\langle E_{\nu_{\mu}}\rangle$ in these models also
become smaller; see panel (d) and (e) in Fig. 11.

In Model IIv2 with a larger $\delta m^{2}$, another instability occurs at
$r\sim 45$ km at 0.032 ms, which later grows to the nonlinear regime. This
flavor conversion is not triggered by the collisional instability, but due to
the well studied slow mode (see Sec. I). However, we note that this onset of
slow mode in our simulation domain is related to the fact that we have
attenuated $\mathbf{H}_{\nu\nu}$ and neglected $\mathbf{H}_{\rm mat}$. Unlike
the collisional instability which is not sensitive to these assumptions, the
slow mode instability can be suppressed by large values of $V_{\nu\nu}$ and
$V_{\rm mat}$ [110, 111, 112], and should only appear at much larger radii if
we had used $V_{\nu\nu}$ and $V_{\rm mat}$ from the corresponding SN
simulation snapshot.

###  V.4 Matter term

An ingredient commonly ignored in the flavor evolution equation of homogeneous
oscillation models is the matter term. In those models, the medium density is
assumed to be homogeneous such that it effectively introduces a global phase
to all considered neutrinos, which can be rotated away [3]. However, since the
length scale of collisional instability is $\sim\mathcal{O}$(1) km, comparable
to the scale length of varying matter profile, the contribution of
${\mathbf{H}}_{\rm mat}$ at different locations varies, which may modify the
dispersion relation of the flavor wave.

Since the original $V_{\rm mat}$ from SN simulation snapshot has the largest
values in our simulation domain (see Fig. 1), whose corresponding length scale
of oscillations cannot be resolved by our simulation setup, we take a
parametric function

$V_{\rm mat}=(2\Delta r)^{-1}\exp\left[-(r-10~{}{\rm km})^{4}/(18~{}{\rm
km})^{4}\right]$ (23)

shown by the blue dotted curve in Fig. 1 to probe the effect of inhomogeneous
matter profile on the evolution of collisional instability. This parametric
function is taken in a way that it mainly varies between 20 km and 50 km where
the collisional instability occurs. It becomes a constant when $r<20$ km and
vanishes to zero when $r>50$ km.

Figure 10(f) shows that although this inhomogeneous matter term does not
affect the initial evolution of flavor instability for $t\leq 0.016$ ms, it
does affect the later transport of flavor mixing by mainly reducing the group
velocity of the propagating flavor waves than that obtained in Model II
without any matter term, which is clearly demonstrated by the red and purple
curves at $t=0.08$ and $0.16$ ms shown in the same panel. Once again, this
effect is also reflected in $\langle E_{\nu_{\mu}}\rangle$ shown in Fig.
11(f), where the reduction of the average energy spreads to larger radii more
gradually compared to Model II.

##  VI Discussion and conclusions

We have implemented a multi-group and discrete-ordinate collective neutrino
oscillation simulator, and solved the neutrino QKE in a spherically symmetric
geometry, including global advection as well as realistic collisional rates in
a self-consistent way. We used this simulator to study the occurrence and
transport of collisional instability in the absence of fast flavor conversion.
Our simulations were performed within the radial range of 10 km to 85 km under
hydrostatic backgrounds taken at different stages from a CCSN simulation. We
confirmed the existence of collisional instability near the neutrinosphere,
and found that the strength of instability increases in later SN stages which
have more asymmetric EA rates due to the deleptonization of the matter and
less differences between $\nu_{e}$ and $\bar{\nu}_{e}$ number densities.

We found that the collisional instability can lead to significant flavor
conversions in three fiducial models that we examined (Models II-IV), which
included minimal ingredients to trigger the instability. Flavor mixing
developed at the location of the initial instability can be transported both
inwards and outwards. The inner branch gets damped and reset by collisions,
while the outer branch transports to the free-streaming regime. For electron
(anti)neutrinos, although their distributions are affected at the onset of the
flavor conversion, the relatively large EA rates quickly reset their spectra
close to the equilibrium ones. For heavy lepton (anti)neutrinos, flavor
conversions not only affect their distributions near their decoupling region,
but also leave imprints in their spectra at the free-streaming regime.
Overall, their neutrino number densities are slightly increased and their mean
energy are significantly reduced. In one of the four models (Model I), we
found that the growth rate of the flavor instability is too small against
advection so that flavor mixing does not reach the nonlinear regime at the end
of our simulation time $\sim 1$ ms.

The results derived in this work suggest a major difference from Ref. [71]:
Although (anti)neutrinos of heavy lepton flavors are affected, their number
densities do not converge to those of electron (anti)neutrinos as predicted in
the homogeneous model wherein global advection is absent and therefore flavor
conversion runs to completion. Our results also imply that although the
collisional instability may not directly affect the $\nu_{e}$ and
$\bar{\nu}_{e}$ heating rates behind the accretion shock, the altered emission
of heavy lepton flavors from region around neutrinosphere may still have
potential impact on supernova dynamics. In addition, the changes of the energy
spectra of heavy lepton flavors can be probed by the future detection of
nearby CCSNe with high statistics.

Beyond the fiducial models, we have also examined impacts from different terms
in the QKE, including the effects due to the NNS decoherence, artificially
enhanced EA decoherence, neutrino vacuum mixing, and inhomogeneous matter
profile for Model II. Including the off-diagonal element of the NNS
collisional term introduces little changes to inner regions where neutrino
distributions are nearly isotropic, but damps flavor mixing at larger radii.
The artificially enhanced EA decoherence leads to higher growth rates of
flavor instability. With non-zero vacuum mixing, the flavor instability
reaches the nonlinear regime earlier. It also results in further decreased
mean energy of heavy lepton flavor neutrinos. The inhomogeneous matter
potential mainly changes the group velocity of the flavor mixing wave by
alternating the dispersion relation. We note that although each of these terms
affects the quantitative behavior of the system, the qualitative features
demonstrated by our fiducial models remain robust.

Although our results confirmed and provided important insights to understand
the collisional flavor instability in SN environment, some cautions should be
noted. First, our simulations are restricted by the trilemma among self-
consistency, advection, and exact rates and hence rely on the strategical
attenuation and parametrization of coherence scattering potentials that may
give rise to inaccurate results when the vacuum mixing and matter potential
are included. Second, due to the computational limitation, we did not include
the NES contribution to collisions, which are important in determining the
exact mean energy of heavy lepton flavors despite their subdominant role in
trapping them. It will be crucial to include these rates in future to assess
the actual impact of collisional flavor instabilities on heavy lepton flavors.
Moreover, our models are based on static SN backgrounds and assume spherical
symmetry, thus neglects the potential impacts due to dynamic evolution of
background profiles and the anisotropy, as well as the feedback of flavor
instability on SN evolution. It remains to be seen how the conclusion derived
in this work holds in a more complete SN model where these simplifications and
assumptions are addressed.

###### Acknowledgements.

We thank Gang Guo and Ninoy Rahman for discussions on the neutrino collisional
rates. ZX and GMP acknowledge support by the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme (ERC
Advanced Grant KILONOVA No. 885281), the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) – Project ID No. 279384907 - SFB 1245 “Nuclei:
From Fundamental Interactions to Structure and Stars”, and the State of Hesse
within the Cluster Project ELEMENTS. MRW and MG acknowledge supports from the
National Science and Techonology Council, Taiwan under Grant No.
110-2112-M-001-050, No. 111-2628-M-001-003-MY4, and the Academia Sinica
(Project No. AS-CDA-109-M11). MRW also acknowledges supports from the Physics
Division of the National Center for Theoretical Sciences, Taiwan. TF
acknowledges support from the Polish National Science Center (NCN) under Grant
No. 2020/37/B/ST9/00691. LJ was supported by NASA through the NASA Hubble
Fellowship Grant No. HST-HF2-51461.001-A awarded by the Space Telescope
Science Institute, which is operated by the Association of Universities for
Research in Astronomy, Incorporated, under NASA Contract No. NAS5-26555. We
would like to acknowledge the use of the following software: Matplotlib [113],
Numpy [114], and SciPy [115].

## References

  * Tanabashi  _et al._ [2018] M. Tanabashi, K. Hagiwara, K. Hikasa, K. Nakamura, and et. al., Phys. Rev. D 98, 030001 (2018).
  * Pastor  _et al._ [2002] S. Pastor, G. Raffelt, and D. V. Semikoz, Phys. Rev. D 65, 053011 (2002), arXiv:hep-ph/0109035 [hep-ph] .
  * Duan  _et al._ [2006a] H. Duan, G. M. Fuller, and Y.-Z. Qian, Phys. Rev. D 74, 123004 (2006a), arXiv:astro-ph/0511275 .
  * Hannestad  _et al._ [2006] S. Hannestad, G. G. Raffelt, G. Sigl, and Y. Y. Y. Wong, Phys. Rev. D 74, 105010 (2006), [Erratum: Phys.Rev.D 76, 029901 (2007)], arXiv:astro-ph/0608695 .
  * Duan  _et al._ [2006b] H. Duan, G. M. Fuller, J. Carlson, and Y.-Z. Qian, Phys. Rev. D 74, 105014 (2006b), arXiv:astro-ph/0606616 .
  * Duan  _et al._ [2007] H. Duan, G. M. Fuller, and Y.-Z. Qian, Phys. Rev. D 76, 085013 (2007), arXiv:0706.4293 [astro-ph] .
  * Raffelt and Sigl [2007] G. G. Raffelt and G. Sigl, Phys. Rev. D 75, 083002 (2007), arXiv:hep-ph/0701182 [hep-ph] .
  * Dasgupta and Dighe [2008] B. Dasgupta and A. Dighe, Phys. Rev. D 77, 113002 (2008), arXiv:0712.3798 [hep-ph] .
  * Gava  _et al._ [2009] J. Gava, J. Kneller, C. Volpe, and G. C. McLaughlin, Phys. Rev. Lett. 103, 071101 (2009), arXiv:0902.0317 [hep-ph] .
  * Friedland [2010] A. Friedland, Phys. Rev. Lett. 104, 191102 (2010), arXiv:1001.0996 [hep-ph] .
  * Raffelt  _et al._ [2013] G. Raffelt, S. Sarikas, and D. de Sousa Seixas, Phys. Rev. Lett. 111, 091101 (2013), [Erratum: Phys.Rev.Lett. 113, 239903 (2014)], arXiv:1305.7140 [hep-ph] .
  * Duan [2013] H. Duan, Phys. Rev. D 88, 125008 (2013), arXiv:1309.7377 [hep-ph] .
  * Wu  _et al._ [2015] M.-R. Wu, Y.-Z. Qian, G. Martinez-Pinedo, T. Fischer, and L. Huther, Phys. Rev. D 91, 065016 (2015), arXiv:1412.8587 [astro-ph.HE] .
  * Mirizzi  _et al._ [2015] A. Mirizzi, G. Mangano, and N. Saviano, Phys. Rev. D 92, 021702 (2015), arXiv:1503.03485 [hep-ph] .
  * Abbar and Duan [2015] S. Abbar and H. Duan, Phys. Lett. B 751, 43 (2015), arXiv:1509.01538 [astro-ph.HE] .
  * Dasgupta and Mirizzi [2015] B. Dasgupta and A. Mirizzi, Phys. Rev. D 92, 125030 (2015), arXiv:1509.03171 [hep-ph] .
  * Martin  _et al._ [2019] J. D. Martin, S. Abbar, and H. Duan, Phys. Rev. D 100, 023016 (2019), arXiv:1904.08877 [hep-ph] .
  * Sawyer [2005] R. F. Sawyer, Phys. Rev. D 72, 045003 (2005), arXiv:hep-ph/0503013 [astro-ph] .
  * Sawyer [2016] R. F. Sawyer, Phys. Rev. Lett. 116, 081101 (2016), arXiv:1509.03323 [astro-ph.HE] .
  * Dasgupta  _et al._ [2017] B. Dasgupta, A. Mirizzi, and M. Sen, JCAP 02, 019 (2017), arXiv:1609.00528 [hep-ph] .
  * Chakraborty  _et al._ [2016] S. Chakraborty, R. S. Hansen, I. Izaguirre, and G. Raffelt, JCAP 03, 042 (2016), arXiv:1602.00698 [hep-ph] .
  * Izaguirre  _et al._ [2017] I. Izaguirre, G. Raffelt, and I. Tamborra, Phys. Rev. Lett. 118, 021101 (2017), arXiv:1610.01612 [hep-ph] .
  * Capozzi  _et al._ [2017] F. Capozzi, B. Dasgupta, E. Lisi, A. Marrone, and A. Mirizzi, Phys. Rev. D 96, 043016 (2017), arXiv:1706.03360 [hep-ph] .
  * Abbar and Duan [2018] S. Abbar and H. Duan, Phys. Rev. D 98, 043014 (2018), arXiv:1712.07013 [hep-ph] .
  * Airen  _et al._ [2018] S. Airen, F. Capozzi, S. Chakraborty, B. Dasgupta, G. Raffelt, and T. Stirner, JCAP 12, 019 (2018), arXiv:1809.09137 [hep-ph] .
  * Johns  _et al._ [2020] L. Johns, H. Nagakura, G. M. Fuller, and A. Burrows, Phys. Rev. D 101, 043009 (2020), arXiv:1910.05682 [hep-ph] .
  * Xiong and Qian [2021] Z. Xiong and Y.-Z. Qian, Phys. Lett. B 820, 136550 (2021), arXiv:2104.05618 [astro-ph.HE] .
  * Bhattacharyya and Dasgupta [2021] S. Bhattacharyya and B. Dasgupta, Phys. Rev. Lett. 126, 061302 (2021), arXiv:2009.03337 [hep-ph] .
  * Duan  _et al._ [2010] H. Duan, G. M. Fuller, and Y.-Z. Qian, Annu. Rev. Nucl. Part. 60, 569 (2010), arXiv:1001.2799 [hep-ph] .
  * Chakraborty  _et al._ [2016] S. Chakraborty, R. Hansen, I. Izaguirre, and G. Raffelt, Nucl. Phys. B. 908, 366 (2016), arXiv:1602.02766 [hep-ph] .
  * Tamborra and Shalgar [2020] I. Tamborra and S. Shalgar, Annu. Rev. Nucl. Part. 71, 165 (2020).
  * Capozzi and Saviano [2022] F. Capozzi and N. Saviano, Universe 8, 94 (2022), arXiv:2202.02494 [hep-ph] .
  * Richers and Sen [2022] S. Richers and M. Sen, arXiv e-prints (2022), arXiv:2207.03561 [astro-ph.HE] .
  * Morinaga [2022] T. Morinaga, Phys. Rev. D 105, L101301 (2022), arXiv:2103.15267 [hep-ph] .
  * Dasgupta [2022] B. Dasgupta, Phys. Rev. Lett. 128, 081102 (2022), arXiv:2110.00192 [hep-ph] .
  * Abbar  _et al._ [2019] S. Abbar, H. Duan, K. Sumiyoshi, T. Takiwaki, and M. C. Volpe, Phys. Rev. D 100, 043004 (2019), arXiv:1812.06883 [astro-ph.HE] .
  * Delfan Azari  _et al._ [2020] M. Delfan Azari, S. Yamada, T. Morinaga, H. Nagakura, S. Furusawa, A. Harada, H. Okawa, W. Iwakami, and K. Sumiyoshi, Phys. Rev. D 101, 023018 (2020), arXiv:1910.06176 [astro-ph.HE] .
  * Morinaga  _et al._ [2020] T. Morinaga, H. Nagakura, C. Kato, and S. Yamada, Phys. Rev. Res. 2, 012046(R) (2020), arXiv:1909.13131 [astro-ph.HE] .
  * Abbar  _et al._ [2020] S. Abbar, H. Duan, K. Sumiyoshi, T. Takiwaki, and M. C. Volpe, Phys. Rev. D 101, 043016 (2020), arXiv:1911.01983 [astro-ph.HE] .
  * Glas  _et al._ [2020] R. Glas, H. T. Janka, F. Capozzi, M. Sen, B. Dasgupta, A. Mirizzi, and G. Sigl, Phys. Rev. D 101, 063001 (2020), arXiv:1912.00274 [astro-ph.HE] .
  * Nagakura  _et al._ [2021] H. Nagakura, A. Burrows, L. Johns, and G. M. Fuller, Phys. Rev. D 104, 083025 (2021), arXiv:2108.07281 [astro-ph.HE] .
  * Harada and Nagakura [2022] A. Harada and H. Nagakura, Astrophys. J. 924, 109 (2022), arXiv:2110.08291 [astro-ph.HE] .
  * Wu and Tamborra [2017] M.-R. Wu and I. Tamborra, Phys. Rev. D 95, 103007 (2017), arXiv:1701.06580 [astro-ph.HE] .
  * Wu  _et al._ [2017] M.-R. Wu, I. Tamborra, O. Just, and H.-T. Janka, Phys. Rev. D 96, 123015 (2017), arXiv:1711.00477 [astro-ph.HE] .
  * George  _et al._ [2020] M. George, M.-R. Wu, I. Tamborra, R. Ardevol-Pulpillo, and H.-T. Janka, Phys. Rev. D 102, 103015 (2020), arXiv:2009.04046 [astro-ph.HE] .
  * Li and Siegel [2021] X. Li and D. M. Siegel, Phys. Rev. Lett. 126, 251101 (2021), arXiv:2103.02616 [astro-ph.HE] .
  * Just  _et al._ [2022] O. Just, S. Abbar, M.-R. Wu, I. Tamborra, H.-T. Janka, and F. Capozzi, Phys. Rev. D 105, 083024 (2022), arXiv:2203.16559 [astro-ph.HE] .
  * Richers [2022] S. Richers, arXiv e-prints (2022), arXiv:2206.08444 [astro-ph.HE] .
  * Fernández  _et al._ [2022] R. Fernández, S. Richers, N. Mulyk, and S. Fahlman, arXiv e-prints (2022), arXiv:2207.10680 [astro-ph.HE] .
  * Grohs  _et al._ [2022] E. Grohs, S. Richers, S. M. Couch, F. Foucart, J. P. Kneller, and G. C. McLaughlin, arXiv e-prints , arXiv:2207.02214 (2022), arXiv:2207.02214 [hep-ph] .
  * Martin  _et al._ [2020] J. D. Martin, C. Yi, and H. Duan, Phys. Lett. B 800, 135088 (2020), arXiv:1909.05225 [hep-ph] .
  * Bhattacharyya and Dasgupta [2020] S. Bhattacharyya and B. Dasgupta, Phys. Rev. D 102, 063018 (2020).
  * Richers  _et al._ [2021a] S. Richers, D. E. Willcox, N. M. Ford, and A. Myers, Phys. Rev. D 103, 083013 (2021a), arXiv:2101.02745 [astro-ph.HE] .
  * Wu  _et al._ [2021] M.-R. Wu, M. George, C.-Y. Lin, and Z. Xiong, Phys. Rev. D 104, 103003 (2021), arXiv:2108.09886 [hep-ph] .
  * Richers  _et al._ [2021b] S. Richers, D. Willcox, and N. Ford, Phys. Rev. D 104, 103023 (2021b), arXiv:2109.08631 [astro-ph.HE] .
  * Zaizen and Morinaga [2021] M. Zaizen and T. Morinaga, Phys. Rev. D 104, 083035 (2021), arXiv:2104.10532 [hep-ph] .
  * Abbar and Capozzi [2022] S. Abbar and F. Capozzi, J. Cosmol. Astropart. Phys. 2022, 051 (2022), arXiv:2111.14880 [astro-ph.HE] .
  * Richers  _et al._ [2022] S. Richers, H. Duan, M.-R. Wu, S. Bhattacharyya, M. Zaizen, M. George, C.-Y. Lin, and Z. Xiong, Phys. Rev. D 106, 043011 (2022), arXiv:2205.06282 [astro-ph.HE] .
  * Bhattacharyya and Dasgupta [2022] S. Bhattacharyya and B. Dasgupta,  106, 103039 (2022), arXiv:2205.05129 [hep-ph] .
  * Sigl and Raffelt [1993] G. Sigl and G. Raffelt, Nuclear Physics B 406, 423 (1993).
  * Vlasenko  _et al._ [2014] A. Vlasenko, G. M. Fuller, and V. Cirigliano, Phys. Rev. D 89, 105004 (2014), arXiv:1309.2628 [hep-ph] .
  * Volpe [2015] C. Volpe, Int. J. Mod. Phys. E 24, 1541009 (2015), arXiv:1506.06222 [astro-ph.SR] .
  * Blaschke and Cirigliano [2016] D. N. Blaschke and V. Cirigliano, Phys. Rev. D 94, 033009 (2016), arXiv:1605.09383 [hep-ph] .
  * Capozzi  _et al._ [2019] F. Capozzi, B. Dasgupta, A. Mirizzi, M. Sen, and G. Sigl, Phys. Rev. Lett. 122, 091101 (2019), arXiv:1808.06618 [hep-ph] .
  * Richers  _et al._ [2019] S. A. Richers, G. C. McLaughlin, J. P. Kneller, and A. Vlasenko, Phys. Rev. D 99, 123014 (2019), arXiv:1903.00022 [astro-ph.HE] .
  * Martin  _et al._ [2021] J. D. Martin, J. Carlson, V. Cirigliano, and H. Duan, Phys. Rev. D 103, 063001 (2021), arXiv:2101.01278 [hep-ph] .
  * Shalgar and Tamborra [2021] S. Shalgar and I. Tamborra, Phys. Rev. D 103, 063002 (2021), arXiv:2011.00004 [astro-ph.HE] .
  * Sigl [2022] G. Sigl, Phys. Rev. D 105, 043005 (2022), arXiv:2109.00091 [hep-ph] .
  * Kato and Nagakura [2022] C. Kato and H. Nagakura, arXiv e-prints , arXiv:2207.09496 (2022), arXiv:2207.09496 [astro-ph.HE] .
  * Sasaki and Takiwaki [2022] H. Sasaki and T. Takiwaki, Progress of Theoretical and Experimental Physics 2022, 073E01 (2022).
  * Johns [2021] L. Johns, arXiv e-prints , arXiv:2104.11369 (2021), arXiv:2104.11369 [hep-ph] .
  * Johns and Xiong [2022] L. Johns and Z. Xiong, arXiv e-prints , arXiv:2208.11059 (2022), arXiv:2208.11059 [hep-ph] .
  * Johns and Nagakura [2022] L. Johns and H. Nagakura, arXiv e-prints , arXiv:2206.09225 (2022), arXiv:2206.09225 [hep-ph] .
  * Shalgar and Tamborra [2022a] S. Shalgar and I. Tamborra, arXiv e-prints , arXiv:2206.00676 (2022a), arXiv:2206.00676 [astro-ph.HE] .
  * Shalgar and Tamborra [2022b] S. Shalgar and I. Tamborra, arXiv e-prints , arXiv:2207.04058 (2022b), arXiv:2207.04058 [astro-ph.HE] .
  * Nagakura [2022] H. Nagakura, arXiv e-prints , arXiv:2206.04098 (2022), arXiv:2206.04098 [astro-ph.HE] .
  * Nagakura and Zaizen [2022] H. Nagakura and M. Zaizen, arXiv e-prints , arXiv:2206.04097 (2022), arXiv:2206.04097 [astro-ph.HE] .
  * Stapleford  _et al._ [2020] C. J. Stapleford, C. Fröhlich, and J. P. Kneller, Phys. Rev. D 102, 081301 (2020), arXiv:1910.04172 [astro-ph.HE] .
  * George  _et al._ [2022] M. George, C.-Y. Lin, M.-R. Wu, T. G. Liu, and Z. Xiong, arXiv e-prints , arXiv:2203.12866 (2022), arXiv:2203.12866 [hep-ph] .
  * Tamborra  _et al._ [2017] I. Tamborra, L. Hüdepohl, G. G. Raffelt, and H.-T. Janka, Astrophys. J. 839, 132 (2017), arXiv:1702.00060 [astro-ph.HE] .
  * Mezzacappa and Bruenn [1993a] A. Mezzacappa and S. W. Bruenn, Astrophys. J. 405, 669 (1993a).
  * Mezzacappa and Bruenn [1993b] A. Mezzacappa and S. W. Bruenn, Astrophys. J. 405, 637 (1993b).
  * Mezzacappa and Bruenn [1993c] A. Mezzacappa and S. W. Bruenn, Astrophys. J. 410, 740 (1993c).
  * Liebendörfer  _et al._ [2001] M. Liebendörfer, A. Mezzacappa, and F.-K. Thielemann, Phys. Rev. D 63, 104003 (2001), arXiv:astro-ph/0012201 [astro-ph] .
  * Liebendörfer  _et al._ [2004] M. Liebendörfer, O. E. B. Messer, A. Mezzacappa, S. W. Bruenn, C. Y. Cardall, and F. K. Thielemann, Astrophys. J., Suppl. 150, 263 (2004), arXiv:astro-ph/0207036 [astro-ph] .
  * Liebendörfer  _et al._ [2002] M. Liebendörfer, S. Rosswog, and F.-K. Thielemann, Astrophys. J., Suppl. 141, 229 (2002), arXiv:astro-ph/0106539 [astro-ph] .
  * Fischer  _et al._ [2010] T. Fischer, S. C. Whitehouse, A. Mezzacappa, F. K. Thielemann, and M. Liebendörfer, Astron. Astrophys. 517, A80 (2010), arXiv:0908.1871 [astro-ph.HE] .
  * Hempel and Schaffner-Bielich [2010] M. Hempel and J. Schaffner-Bielich, Nucl. Phys. A. 837, 210 (2010), arXiv:0911.4073 [nucl-th] .
  * Typel  _et al._ [2010] S. Typel, G. Röpke, T. Klähn, D. Blaschke, and H. H. Wolter, Phys. Rev. C 81, 015803 (2010), arXiv:0908.2344 [nucl-th] .
  * Timmes and Swesty [2000] F. X. Timmes and F. D. Swesty, Astrophys. J., Suppl. 126, 501 (2000).
  * Guo  _et al._ [2020] G. Guo, G. Martínez-Pinedo, A. Lohs, and T. Fischer, Phys. Rev. D 102, 023037 (2020), arXiv:2006.12051 [hep-ph] .
  * Fischer  _et al._ [2020a] T. Fischer, G. Guo, G. Martínez-Pinedo, M. Liebendörfer, and A. Mezzacappa, Phys. Rev. D 102, 123001 (2020a), arXiv:2008.13628 [astro-ph.HE] .
  * Fischer  _et al._ [2020b] T. Fischer, G. Guo, A. A. Dzhioev, G. Martínez-Pinedo, M.-R. Wu, A. Lohs, and Y.-Z. Qian, Phys. Rev. C 101, 025804 (2020b), arXiv:1804.10890 [astro-ph.HE] .
  * Bruenn [1985] S. W. Bruenn, Astrophys. J., Suppl. 58, 771 (1985).
  * Thompson and Burrows [2001] T. A. Thompson and A. Burrows, Nucl. Phys. A. 688, 377 (2001), arXiv:astro-ph/0009449 [astro-ph] .
  * Fischer [2016] T. Fischer, Astron. Astrophys. 593, A103 (2016), arXiv:1608.05004 [astro-ph.HE] .
  * Fischer  _et al._ [2012] T. Fischer, G. Martínez-Pinedo, M. Hempel, and M. Liebendörfer, Phys. Rev. D 85, 083003 (2012), arXiv:1112.3842 [astro-ph.HE] .
  * Raffelt [2001] G. G. Raffelt, Astrophys. J. 561, 890 (2001), arXiv:astro-ph/0105250 [astro-ph] .
  * Keil  _et al._ [2003] M. T. Keil, G. G. Raffelt, and H.-T. Janka, Astrophys. J. 590, 971 (2003), arXiv:astro-ph/0208035 [astro-ph] .
  * Duan  _et al._ [2008] H. Duan, G. M. Fuller, and Y.-Z. Qian, Phys. Rev. D 77, 085016 (2008), arXiv:0801.1363 [astro-ph] .
  * Dasgupta  _et al._ [2008] B. Dasgupta, A. Dighe, A. Mirizzi, and G. G. Raffelt, Phys. Rev. D 77, 113007 (2008), arXiv:0801.1660 [hep-ph] .
  * Capozzi  _et al._ [2020] F. Capozzi, M. Chakraborty, S. Chakraborty, and M. Sen, Phys. Rev. Lett. 125, 251801 (2020), arXiv:2005.14204 [hep-ph] .
  * Capozzi  _et al._ [2022] F. Capozzi, M. Chakraborty, S. Chakraborty, and M. Sen, arXiv e-prints (2022), arXiv:2205.06272 [hep-ph] .
  * Richers  _et al._ [2019] S. A. Richers, G. C. McLaughlin, J. P. Kneller, and A. Vlasenko, Phys. Rev. D 99, 123014 (2019), arXiv:1903.00022 [astro-ph.HE] .
  * Rampp and Janka [2002] M. Rampp and H. T. Janka, Astronomy and Astrophysics 396, 361 (2002), arXiv:astro-ph/0203101 [astro-ph] .
  * Hannestad and Raffelt [1998] S. Hannestad and G. Raffelt, Astrophys. J. 507, 339 (1998), arXiv:astro-ph/9711132 [astro-ph] .
  * Cherry  _et al._ [2012] J. F. Cherry, J. Carlson, A. Friedland, G. M. Fuller, and A. Vlasenko, Phys. Rev. Lett. 108, 261104 (2012), arXiv:1203.1607 [hep-ph] .
  * Banerjee  _et al._ [2011] A. Banerjee, A. Dighe, and G. Raffelt, Phys. Rev. D 84, 053013 (2011), arXiv:1107.2308 [hep-ph] .
  * Izaguirre  _et al._ [2017] I. Izaguirre, G. Raffelt, and I. Tamborra, Phys. Rev. Lett. 118, 021101 (2017), arXiv:1610.01612 [hep-ph] .
  * Esteban-Pretel  _et al._ [2008] A. Esteban-Pretel, A. Mirizzi, S. Pastor, R. Tomas, G. G. Raffelt, P. D. Serpico, and G. Sigl, Phys. Rev. D 78, 085012 (2008), arXiv:0807.0659 [astro-ph] .
  * Duan and Friedland [2011] H. Duan and A. Friedland, Phys. Rev. Lett. 106, 091101 (2011), arXiv:1006.2359 [hep-ph] .
  * Chakraborty  _et al._ [2011] S. Chakraborty, T. Fischer, A. Mirizzi, N. Saviano, and R. Tomas, Phys. Rev. Lett. 107, 151101 (2011), arXiv:1104.4031 [hep-ph] .
  * Hunter [2007] J. D. Hunter, Computing in Science & Engineering 9, 90 (2007).
  * van der Walt  _et al._ [2011] S. van der Walt, S. C. Colbert, and G. Varoquaux, Computing in Science & Engineering 13, 22 (2011).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).

##  Appendix A Rates of collisions

The emissivity and opacity for EA processes are taken from BOLTZTRAN using
relativistic dispersion relations for nucleons in the nuclear medium and
including weak magnetism, pseudoscalar term, and form factor effects [91].

NNS is almost elastic and commonly treated as iso-energetic so that

$\displaystyle\mathbf{C}_{\rm NNS}(E,v_{r})=\frac{1}{2}\int
dv_{r}^{\prime}\,\times$ $\displaystyle\left[-R_{\rm
NNS}(E,v_{r},v_{r}^{\prime})\left\\{\varrho(E,v_{r}),\mathcal{I}-\frac{4\pi^{2}}{E^{2}}\varrho(E,v_{r}^{\prime})\right\\}\right.$
$\displaystyle\left.+R_{\rm
NNS}(E,v_{r}^{\prime},v_{r})\left\\{\varrho(E,v_{r}^{\prime}),\mathcal{I}-\frac{4\pi^{2}}{E^{2}}\varrho(E,v_{r})\right\\}\right],$
(24)

where, e.g., $4\pi^{2}\varrho_{ee}(E,v_{r})/E^{2}$ inside the anti-commutator
is the distribution function of $\nu_{e}$ accounting for the Pauli blocking.
This iso-energetic scattering kernel has the symmetry $R_{\rm
NNS}(E,v_{r},v_{r}^{\prime})=R_{\rm NNS}(E,v_{r}^{\prime},v_{r})$ so that the
blocking in the gain and loss terms cancels as
$\\{\varrho(E,v_{r}),\varrho(E,v_{r}^{\prime})\\}=\\{\varrho(E,v_{r}^{\prime}),\varrho(E,v_{r})\\}$
and effectively has no impacts.

The two moments of NNS opacities in Eq. (9) without weak-magnetism and recoil
correction are given as below

$\displaystyle\chi_{\rm NNS}(E)$
$\displaystyle=\frac{G_{F}^{2}E^{2}}{\pi}[(c_{V,n}^{2}+3c_{A,n}^{2})\eta_{n}+(c_{V,p}^{2}+3c_{A,p}^{2})\eta_{p}],$
$\displaystyle\tilde{\chi}_{\rm NNS}(E)$
$\displaystyle=\frac{G_{F}^{2}E^{2}}{\pi}[(c_{V,n}^{2}-c_{A,n}^{2})\eta_{n}+(c_{V,p}^{2}-c_{A,p}^{2})\eta_{p}],$
(25)

with weak coupling constants $c_{V,p}=1/2-2\sin^{2}\theta_{\rm W}$,
$c_{V,n}=-1/2$, $c_{A,p}=g_{A}/2$, $c_{A,n}=-g_{A}/2$, $g_{A}\approx 1.27$,
and $\sin^{2}\theta_{\rm W}\approx 0.23$. The effective number densities

$\eta_{N}=\frac{3Tn_{N}}{\sqrt{4E_{{\rm F},N}^{2}+9T^{2}}}$ (26)

account for the nucleon final state blocking with Fermi energies $E_{{\rm
F},N}=(3\pi^{2}n_{N})^{2/3}/(2m_{u})$, and number densities $n_{N}=\rho
Y_{N}/m_{u}$ for $N=n,p,$ respectively.

##  Appendix B Resolution

In this section, we discuss the dependence of our results on the adopted
simulation resolutions for Model II. In addition to the fiducial models with
highest radial resolution $N_{r}=25000$, we examine the other two different
radial resolutions with $N_{r}=10000$ and 2500. For models with lower radial
resolutions, we use the same parameters for Model II listed in Table 2, except
$a_{1}=4\times 10^{-4}$ for $N_{r}=10000$ and $a_{1}=10^{-4}$ for
$N_{r}=2500$. The different choice of $a_{1}$ here is to satisfy the
requirement that the length scale of oscillations can be resolved by $\Delta
r$. We have also examined cases with $N_{v_{r}}=200$ for the lowest radial
resolution runs. The evolution histories are shown in Fig. 12. The case with
$N_{v_{r}}=200$ undergoes almost identical evolution as the one with
$N_{v_{r}}=50$ which shows that 50 angular grids are sufficient.

Figure 12:  Evolution of the dimensionless ratio $\log_{10}(s_{e\mu})$ in
Model II at $t_{\rm pb}=247$ ms with three resolutions:
$N_{r}=10000,N_{v_{r}}=50$ (left panel), $N_{r}=2500,N_{v_{r}}=50$ (middle
panel), and $N_{r}=2500,N_{v_{r}}=200$ (right panel).

With the same number of angular grids $N_{v_{r}}=50$, the lowest radial
resolution run ($N_{r}=2500$) shows clearly visible differences from results
obtained with intermediate ($N_{r}=10000$) and high ($N_{r}=25000$)
resolutions. First of all, the grow rate in the initial linear regime is
reduced while the onset radius of instability shifts from $r\approx 29$ km to
$32$ km. Similar to runs with higher resolutions, after the flavor mixing
reaches nonlinear regime between $t=0.02$ and 0.04 ms, it bifurcates. However,
the inward-moving branch does not get damped completely, regrows around $t\sim
0.2$ ms, and eventually reaches the stationary state much earlier than the
high resolution runs by $\approx 0.5$ ms. On the other hand, results obtained
with $N_{r}=10000$ and $N_{r}=25000$ (see Fig. 2) are qualitatively more
similar, although slight differences exist.

The smaller growth rates in the initial linear regime with $N_{r}=2500$ is
related to the imposed larger attenuation on $\mathbf{H}_{\nu\nu}$. As
indicated from Eq. (21), if $V_{\nu\nu}$ and collisional rates are in similar
magnitude, the high-order term can reduces the growth rate of the instability.
On the other hand, if $V_{\nu\nu}$ is much larger than the collisional rates,
the growth rate of collisional instability converges to an asymptotic value.
The LSA analysis results obtained by solving Eqs. III and III also support
this conclusion. Figure 13 shows the growth rates as functions of radius for
all three radial resolution runs at $t=0$. The peak of growth rates with
$N_{r}=2500$ has Im$(\Omega)\approx 1.6$ km-1 at $\sim 30$ km, lower than
other two cases with peak values $\approx 5-6$ km-1 at $\sim 26-27$ km.

Figure 13:  Growth rates Im$(\Omega)$ from linear stability analysis as
functions of radius $r$ at $t=0$ ms for three choices of radial resolutions
and attenuation factors in Model II. We sample different $K_{r}$ and show the
maximal values of growth rates as in Fig. 3.

For the resolution in the energy grid, in addition to $N_{E}=20$, we have also
used $N_{E}=8$. However, we find that in this case each grid has too wide
width to accurately model the sharp dependence of the muonic EA rates at $\sim
65$ MeV. This insufficient resolution then leads to the artificial
amplification of number density $n_{\nu_{\mu}}$ near the decoupling region and
affects our simulation results.