# Neutrino transport with Monte Carlo method: II. Quantum Kinetic Equations

Chinami Kato  Department of Science and Technology, Tokyo University of
Science, 2641 Yamazaki, Noda-shi, Chiba Prefecture 278-8510, Japan  Department
of Aerospace Engineering, Tohoku University, 6-6-01 Aramaki-Aza-Aoba, Aoba-ku,
Sendai 980-8579, Japan  Hiroki Nagakura  National Astronomical Observatory of
Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Department of
Astrophysical Sciences, Princeton University, Princeton, NJ 08544  Taiki
Morinaga  Graduate School of Advanced Science and Engineering, Waseda
University, 3-4-1 Okubo, Shinjuku, Tokyo 169-8555, Japan

###### Abstract

Neutrinos have an unique quantum feature as flavor conversions. Recent studies
suggested that collective neutrino oscillations play important roles in high-
energy astrophysical phenomena. Quantum kinetic equation (QKE) is capable of
describing the neutrino flavor conversion, transport and matter collision
self-consistently. However, we have experienced many technical difficulties in
their numerical implementation. In this paper, we present a new QKE solver
based on Monte Carlo (MC) approach. This is an upgraded version of our
classical MC neutrino transport solver; in essence, a flavor degree of freedom
including mixing state is added into each MC particle. This extension requires
updating numerical treatments of collision terms, in particular for scattering
processes. We deal with the technical problem by generating a new MC particle
at each scattering event. To reduce statistical noise inherent in MC methods,
we develop the effective mean free path method. This suppresses a sudden
change of flavor state due to collisions without increasing the number of MC
particles. We present a suite of code tests to validate these new modules with
comparing to the results reported in previous studies. Our QKE-MC solver is
developed with fundamentally different philosophy and design from other
deterministic- and mesh methods, suggesting that it will be complementary to
others, and potentially provide new insights into physical processes of
neutrino dynamics.

supernova:general — neutrinos —

††journal: ApJS

##  1 Introduction

In hot and dense media, neutrinos have important roles as effective carriers
of energy and lepton-number. During the advanced stages of massive star
evolution ($M\gtrsim 10M_{\sun}$), neutrino cooling dictates the thermal
evolution of the stellar core (Farag et al., 2020; Odrzywolek & Heger, 2010).
In the subsequent core-collapse supernova (CCSN), neutrinos are produced
abundantly via various weak interaction processes. The neutrino transport and
their matter interactions may be central to the explosion mechanism (Janka,
2012; Kotake et al., 2012; Foglizzo et al., 2015; Mirizzi et al., 2016;
Müller, 2019, 2020; Burrows & Vartanyan, 2021), and of nucleosynthesis of
heavy elements (Martínez-Pinedo et al., 2014; Harris et al., 2017; Wanajo et
al., 2018; Sieverding et al., 2020; Witt et al., 2021). Regardless of the
detailed explosion mechanism, the emitted neutrinos from Galactic events are
observable signals by current and future-planned terrestrial detectors. The
direct detection of neutrinos from these nearby CCSNe will provide high
statistics neutrino data, that will shed light on the inner dynamics of CCSN
(Horiuchi & Kneller, 2018; Gallo Rosso et al., 2018; Suwa et al., 2019;
Nagakura et al., 2021a, b; Li et al., 2021; Mori et al., 2021). The remnant of
binary neutron star merger (BNSM) is another exciting production site of
neutrinos, that has received huge interest with the dawn of multi-messenger
astrophysics (Abbott et al., 2017). Although neutrino dynamics is a key
ingredient to account for the time evolution of the system, the detailed
flavor-dependent features of energy spectrum and angular distributions of
neutrinos are still uncertain. Accurate treatments of neutrino transport and
neutrino-matter interactions are indispensable towards understanding of the
global dynamics, ejecta composition and observable signals including kilonova
transients (Sekiguchi et al., 2015; Thielemann et al., 2017; Siegel & Metzger,
2017; Radice et al., 2018; Miller et al., 2019; Kawaguchi et al., 2021).

There is growing experimental evidence that neutrinos’ flavor states oscillate
from one to another during the propagation, known as neutrino oscillation. The
flavor mixing is attributed to the fact that neutrinos are massive particles,
and that their mass- and flavor eigenstates are not identical each other. In a
vacuum, the flavor conversion is dictated by the difference of eigenvalues of
the mass eigenstates (i.e., the mass difference of neutrinos), and non-zero
mixing angles. Accurate determination of these oscillation parameters is the
major challenge in experimental physics, that has made remarkable progress
over the past decades (Abe et al., 2014; Adamson et al., 2014; T2K
Collaboration & Abe, 2014; Abe et al., 2018a, b; Aartsen et al., 2018; Acero
et al., 2019; Esteban et al., 2020).

Neutrinos experience refractive effects when they propagate in matter. This
changes the dispersion relation of neutrinos, implying that flavor conversions
are altered. The matter effect dominates over the vacuum oscillation in CCSNe
and remnants of BNSM. This may induce strong flavor mixing at the resonance
MSW matter density (Wolfenstein, 1978). Another important consequence of the
matter effect is that the mixing angles between flavor- and effective mass
eigenstates are substantially reduced at a high density region, indicating
that the flavor conversion is strongly suppressed. This is a rationale behind
the assumption that no flavor mixing occur in a CCSN/BNSM core (e.g., Fuller
et al., 1987; Dighe & Smirnov, 2000); indeed, almost every numerical
simulation has been carried out with the assumption.

However, refractive effects by neutrino-self interactions potentially induce
strong flavor conversion, also known as collective neutrino oscillations
(Samuel, 1993; Sigl & Raffelt, 1993; Sigl, 1995; Sawyer, 2005). They are
genuinely a non-linear phenomenon with large off-diagonal components of the
Hamiltonian potential in the flavor-basis, suggesting that the flavor
conversion is qualitatively different from that induced by the vacuum and
matter potentials. The collective oscillation involves two distinct modes:
fast- and slow flavor conversions. The time scales of these flavor conversions
are very short, in particular for the fast mode; the oscillation frequency
increases proportionally to the neutrino number density. Although the growth
rate of the slow mode is relatively mild, various interesting phenomena have
been found, for instances, spectral swaps/splits, bipolar and synchronized
oscillations, although the occurrence of these phenomena in CCSN and BNSM
environments are still being debated (Duan et al., 2006a, b; Dasgupta et al.,
2009; Tian et al., 2017; Martin et al., 2020a).

In the last several years, there has been significant progress of our
theoretical understanding of collective neutrino oscillations. Some of the
progress are owed to the linear stability analysis that provides not only
insights into characteristics of collective neutrino oscillations but also
criteria to assess the onset of flavor conversions (Dasgupta et al., 2017;
Nagakura et al., 2019; Delfan Azari et al., 2019; Abbar et al., 2020; Morinaga
et al., 2020; Delfan Azari et al., 2020). Regarding the fast mode, electron
neutrino lepton number crossings (ELN crossings) in angular distributions are
associated with the instability (Morinaga, 2021). The recent studies suggested
that ELN crossings have been frequently observed in the neutrino data of CCSN
and BNSM simulations (Tamborra et al., 2017; Wu & Tamborra, 2017; Dasgupta et
al., 2018; Glas et al., 2020; Abbar, 2020; Capozzi et al., 2020; Abbar, 2021;
Abbar et al., 2021), that increases the possibility of occurrence of fast
flavor conversion in real CCSN/BNSM environments. Concerning the slow mode,
some recent studies have demonstrated that this type of flavor oscillation
affects observational consequences of CCSNe such as nucleosynthesis and
neutrino signals (Chakraborty et al., 2010; Zaizen et al., 2018; Sasaki et
al., 2017; Wu et al., 2017; George et al., 2020). These results motivate more
detailed studies of collective neutrino oscillations.

Aside from a flavor mixing, neutrinos experience absorption and momentum-
exchanged scattering111We note that the term ”momentum-exchanged scattering”
is used to distinguish from the forward scattering generating neutrino
refraction. during the propagation in medium. The macroscopic properties of
neutrino radiation field in CCSN and BNSM strongly hinge on these collision
term. They also affect the dynamics of collective neutrino oscillations; for
instances, both fast and slow modes may be triggered by neutrino scatterings
with nucleons or nuclei (Cherry et al., 2013, 2020; Morinaga et al., 2020;
Zaizen et al., 2020); the characteristics of flavor conversions may be altered
by these scatterings (Capozzi et al., 2019; Johns et al., 2020; Shalgar &
Tamborra, 2021a, b; Martin et al., 2021); Johns (2021) suggested very recently
that a new type of instability, collisional flavor instability, appears due to
collision terms. Hence, the coupling between collective neutrino oscillations
and collision terms is currently of great interest. On the other hand, time
scales of neutrino-matter collisions are usually longer than that of neutrino
oscillations, indicating that the investigation of the interplay to neutrino-
matter collisions should inevitably be in non-linear phase of collective
oscillations. The direct numerical simulation is, hence, the most useful
approach for the investigation.

The consistent treatment of neutrino transport, flavor conversion and
neutrino-matter collisions requires the quantum kinetic description of the
neutrino-radiation field. The quantum kinetic equation (QKE) provides the
basic equation, that can be concisely written in the mean-field approximation
(see also Eq. 1). However, the direct numerical simulation is not an easy
task, but rather a grand challenge in computational physics. This is mainly
attributed to the large disparity of time- and spatial scales between flavor
conversions and global (macroscopic) system. Recent studies in Johns et al.
(2020); Richers et al. (2021) also found that fast flavor conversions
generally involve cascade phenomena, in which fast flavor conversions trigger
the energy exchange between large- to small scales. This suggests that small
scale structures of both spatial- and momentum spaces need to be resolved,
posing another challenge in numerical simulations. For these reasons, the
currently feasible approach is either local simulations (Abbar & Volpe, 2019;
Richers et al., 2019; Martin et al., 2020b; Shalgar & Tamborra, 2021a; Richers
et al., 2021; Martin et al., 2021; Zaizen & Morinaga, 2021) or large-scale
simulations with employing crude approximations and simplifications (see,
e.g., Duan et al., 2006a; Dasgupta & Dighe, 2008; Cherry et al., 2012; Zaizen
et al., 2018; Stapleford et al., 2020; Li et al., 2021).

Regarding the methodologies, there are mainly two options: deterministic- and
probabilistic approaches. The former has been popular in QKE neutrino
transport. This may be because the numerical technique can be simply extended
from classical radiation transport; indeed, it is a good compatibility with
finite volume methods (but see a PIC method in Richers et al. (2021)). On the
other hand, there are no previous work that employ the probabilistic way in
QKE neutrino transport. The main reason of the absence of demonstration is
that probabilistic treatments are not necessary in QKE neutrino transport,
unless the collision term is taken into account. As already mentioned,
however, there is a growing interest in interplay between flavor conversions
and collision terms, indicating that the detailed investigation of the
coupling system will be made more vigorously in future. It should also be
mentioned that, since all numerical simulations involve a certain level of
approximation, it is necessary to assess the reliability of their results.
Independent numerical approach will help the assessment, and provide
confidence in the accuracy of each result. In fact, such collaborative code-
comparison projects can be seen in classical neutrino transport (Richers et
al., 2017; O’Connor et al., 2018).

In this paper, we present a new QKE neutrino transport solver based on a Monte
Carlo (MC) method, which is a representative probabilistic approach. This QKE-
MC code is an extension version of our classical MC neutrino transport one
(Kato et al., 2020), in which the neutrino transport with various weak
processes are well tested in the context of CCSN simulations. On the other
hand, there is a non-trivial issue: how the quantum kinetic states can be
handled in the particle method. Addressing this issue is one of the focuses in
this paper. We develop new modules to reduce the statistical noise inherent to
the MC approach, and to increase the computational efficiency. After
presenting the numerical implementation, we demonstrate a suite of code tests
including fast flavor conversions with momentum-exchanged scatterings to test
the capability of these new modules.

This paper is organized as follows. We start with providing an overview of the
QKE for neutrino transport in Sec. 2. The essential methodologies in our QKE-
MC solver are described in Sec. 3. Before entering into the detail of code
tests, we describe the basic assumption in Sec. 4. From Secs. 5 to 7, we
present code tests for vacuum, matter and collective neutrino oscillations,
respectively. Other tests regarding QKE with scatterings are demonstrated in
Sec. 8. We conclude the paper with a summary and discussions in Section 9.
Throughout this paper, we use the metric signature of $(-+++)$. Unless
otherwise stated, Greek subscripts denote the spacetime components, and we
work in units with $c=\hbar=1$, where $c$ and $\hbar$ denote the speed of
light and reduced Planck constant, respectively.

##  2 Basic equation

The QKE for neutrino transport can be written as (see also Sigl & Raffelt,
1993; Volpe, 2015),

$\displaystyle i\left(p^{\alpha}\left.\frac{\partial\rho}{\partial
x^{\alpha}}\right|_{p^{i}}-\Gamma^{i}_{\alpha\beta}p^{\alpha}p^{\beta}\left.\frac{\partial\rho}{\partial
p^{i}}\right|_{x^{\mu}}\right)$ $\displaystyle\ \ \ \ \
=(-p^{\mu}\zeta_{\mu})\left[H_{\rm{vac}}+H_{\rm{mat}}+H_{\nu\nu},\rho\right]+C,$
$\displaystyle i\left(p^{\alpha}\left.\frac{\partial\bar{\rho}}{\partial
x^{\alpha}}\right|_{p^{i}}-\Gamma^{i}_{\alpha\beta}p^{\alpha}p^{\beta}\left.\frac{\partial\bar{\rho}}{\partial
p^{i}}\right|_{x^{\mu}}\right)$ $\displaystyle\ \ \ \ \
=(-p^{\mu}\zeta_{\mu})\left[H_{\rm{vac}}^{*}-H_{\rm{mat}}^{*}-H_{\nu\nu}^{*},\bar{\rho}\right]+\bar{C},$
(1)

where $\rho$ $(\bar{\rho})$ is the density matrix for neutrinos (anti-
neutrinos); $x^{\mu}$ denotes the space-time coordinate; $p^{\mu}$ denotes the
four-momentum of neutrinos; $C$ $(\bar{C})$ is the collision term for
neutrinos (anti-neutrinos); $\Gamma$ is the Christoffel symbols; $\zeta^{\mu}$
denotes the unit vector normal to the spatial hypersurface of $t={\rm const}$,
where $t$ is the laboratory time. In the expression, Hamiltonian operators are
measured in the laboratory frame222We note that it can be measured with
arbitrary time. In this case, we replace $(-p^{\mu}\zeta_{\mu})$ by
$(-p^{\mu}\zeta^{\prime}_{\mu})$ in Eq. 1, where $\zeta^{\prime}_{\mu}$
denotes the arbitrary timelike vector.. Each potential can be given as

$\displaystyle H_{\rm{vac}}$ $\displaystyle=$ $\displaystyle
U\frac{M^{2}}{2E_{\nu}}U^{\dagger},$ $\displaystyle H_{\rm{mat}}$
$\displaystyle=$ $\displaystyle\sqrt{2}G_{F}\frac{p_{\mu}}{E_{\nu}}{\rm
diag}\left[j^{\mu}_{l}-\bar{j}^{\mu}_{l}\right],$ $\displaystyle H_{\nu\nu}$
$\displaystyle=$ $\displaystyle\sqrt{2}G_{F}\int
dV_{p^{\prime}}\frac{p^{\mu}p_{\mu}^{\prime}}{E_{\nu}E_{\nu}^{\prime}}\left(\rho^{\prime}-\bar{\rho}^{\prime\ast}\right).$
(2)

In the expression, $M^{2}$ can be written as $M^{2}_{ij}=m^{2}_{a}\delta_{ij}$
(the index $a$ specifies the mass eigenstate of neutrinos) in the relativistic
limit, where $m_{a}$ denote the rest-mass of $\nu_{a}$.
$E_{\nu}=-p^{\mu}\zeta_{\mu}$ represents the energy of neutrinos measured by
an observer of $\zeta_{\mu}$. $U$ is the unitary matrix describing the mixing
between flavor and mass bases. $G_{F}$ is the Fermi coupling constant.
$j^{\mu}_{l}$ represents the number current of charged-leptons, which can be
written as $j^{\mu}_{l}=n_{l}u^{\mu}_{l}$, where $n_{l}$ and $u^{\mu}_{l}$
denotes the the number density (measured at each lepton’s rest frame) and
four-velocity of each charged-leptons, respectively. $dV_{p}$ denotes the
momentum volume element for neutrino momentum space, which can be written as

$\displaystyle dV_{p}$ $\displaystyle=$
$\displaystyle\frac{1}{(2\pi)^{3}}\hskip 1.42262ptd\Omega\hskip
1.42262ptd\left(\frac{E_{\nu}^{3}}{3}\right),$ (3)

where $d\Omega$ denotes the solid angle of neutrinos measured in $\zeta_{\mu}$
observer333We note that $\zeta_{\mu}$ can be replaced any time-like unit
vector. This depends on the choice of tetrad basis to describe the momentum
space of neutrinos (see e.g., Shibata et al., 2014, for more detail)..

##  3 Numerical Implementations

Many modules in our QKE-MC solver are imported from classical one (Kato et
al., 2020). However, the flavor states in MC particles are no longer constant
during the flight due to flavor conversions. A possible solution is that a
flavor degree of freedom is added in each MC particle. This change also
requires an update of treatments for collision terms. In this section, we
describe our new ideas of the formulation and numerical techniques for them.

###  3.1 Transport and oscillation

Each MC particle of position and momentum state follows the geodesic equation
during propagation,

$\displaystyle\frac{dx^{\mu}}{d\lambda}=p^{\mu},$
$\displaystyle\frac{dp^{\mu}}{d\lambda}-\Gamma^{\mu}_{\hphantom{\mu}\alpha\beta}p^{\alpha}p^{\beta}=0,$
(4)

where $\lambda$ denotes the affine parameter for a trajectory of each
particle. We assume neutrinos as massless particles, i.e., Eq. 4 is null
geodesics. It is a reasonable assumption for CCSN and BNSM system, since the
typical energy of neutrinos are above MeV and it is much larger than the rest
mass of neutrinos ($\lesssim 0.1$ eV) (Aker et al., 2019). We note that
neutrinos are treated as massive particles in the evolution of flavor state
(see below).

Different from conventional MC neutrino transport, the flavor degree of
freedom is embedded in each MC particle, that allows us to treat the neutrino
flavor conversion self-consistently. More specifically, we define the
particle-matrix, $q$, on each MC particle. The matrix structure of $q$ is the
same as that of density matrix ($\rho$), i.e., $N_{f}\times N_{f}$ matrix
($N_{f}$ denotes the number of flavor). The matrix element represents the
number of particles in the flavor basis, i.e., this is a natural extension
from ”sample particles” used in classical MC method. The density matrix
($\rho$) at each time can be computed from $q$ as

$\displaystyle\rho(\mathbb{x},\mathbb{p})$ $\displaystyle=$
$\displaystyle\lim_{\Delta V_{x}\to 0}\lim_{\Delta V_{p}\to 0}\frac{1}{\Delta
V_{x}\hskip 1.42262pt\Delta V_{p}}$ (5)
$\displaystyle\times\sum_{s}q_{s}\int_{\mathbb{x}-0.5\mathbb{\Delta{x}}}^{\mathbb{x}+0.5\mathbb{\Delta{x}}}\int_{\mathbb{p}-0.5\mathbb{\Delta{p}}}^{\mathbb{p}+0.5\mathbb{\Delta{p}}}dV_{x}^{\prime}dV_{p}^{\prime}\hskip
1.42262pt$ $\displaystyle\times\
\delta^{3}\left(\mathbb{x}^{\prime}-\mathbb{x_{s}}(t)\right)\hskip 1.42262pt$
$\displaystyle\times\left(2\pi\right)^{3}\delta^{3}\left(\mathbb{p}^{\prime}-\mathbb{p_{s}}(t)\right),$

where $s$ denotes the index of a MC particle, and $\bar{\rho}$ (anti-
neutrinos) can be computed by replacing $q\rightarrow\bar{q}$. In the
expression, $x_{s}$ and $p_{s}$ are the spatial position and momentum of
$s$-th QKE-MC particle, respectively; $dV_{x}$ denotes a spatial three-
dimensional volume element,

$\displaystyle dV_{x}$ $\displaystyle=$ $\displaystyle\sqrt{-g}\hskip
1.42262ptdx_{1}\hskip 1.42262ptdx_{2}\hskip 1.42262ptdx_{3},$ (6)

where $g$ denotes the determinant of the space-time metric.

In practice, $\rho$ is computed by setting finite but small values of $\Delta
V_{x}$ and $\Delta V_{p}$ in Eq. 5. Here we take the on-shell condition for
neutrinos. The number density of neutrinos (measured in the laboratory frame)
can be computed through $\rho$,

$\displaystyle n_{ij}=\int\rho_{ij}dV_{p}^{\prime},$ (7)

where $i$ and $j$ specify the flavor state. The anti-neutrino counterpart can
be given by replacing $\rho\rightarrow\bar{\rho}$.

We evolve $q$ and $\bar{q}$ as

$\displaystyle i\frac{dq_{s}}{d\lambda}$ $\displaystyle=$
$\displaystyle(-p^{\mu}n_{\mu})[H_{\rm{vac}}+H_{\rm{mat}}+H_{\rm{\nu\nu}},q_{s}],$
$\displaystyle i\frac{d\bar{q}_{s}}{d\lambda}$ $\displaystyle=$
$\displaystyle(-p^{\mu}n_{\mu})[H_{\rm{vac}}^{\ast}-H_{\rm{mat}}^{\ast}-H_{\rm{\nu\nu}}^{\ast},\bar{q}_{s}],$
(8)

where $d/d\lambda$ in Eq. 8 is an ordinary differentiation, which is solved
with forth-order explicit Runge-Kutta method. We note that the self-
interaction potential, $H_{\rm{\nu\nu}}$, can be computed from $q$ and
$\bar{q}$ by using Eqs. 2 and 5.

###  3.2 Collision term

Figure 1: The schematic picture illustrating our numerical treatments of
scattering process. As an example, we display the case that $\nu_{e}$ ($ee$
component of $q$) experiences scatterings. We note that $a$ corresponds to a
control parameter relevant to effective mean free path (EMFP) method. See text
for more detail.

The implementation of collision terms is well established in classical MC
method, and its extension to QKE is straightforward. Below, we first describe
the essence of the extension.

Let us start with emission processes. Although the neutrino matter interaction
needs to be changed with the extension of QKE (Vlasenko et al., 2014; Blaschke
& Cirigliano, 2016), the implementation is the same as that used in classical
MC transport. The blocking factor can also be taken into account by the same
prescription used in our classical MC method. We generate new MC particles
with specifying the neutrino energy and flight direction. Following the
neutrino emissivity, we insert the value of neutrino number of production into
each matrix element of $q_{s}$ (and $\bar{q}_{s}$) of newly generated MC
particles.

For absorption process, we probabilistically determine the absorption length,
the distance to a position where the neutrinos experience absorption; the
procedure is exactly the same as that used in our classical MC transport. We
note, however, that the length can not be assigned into each MC particle,
since the reaction rate in general varies with the neutrino flavor. Therefore,
it is assigned to each flavor state of MC particles. MC particles are evolved
following the geodesic equation until the travel distance reaches the
absorption length. At the absorption point, the corresponding matrix element
of $q$ set to be zero. If $q$ becomes completely null due to the absorption,
we exclude the MC particle from our simulations.

The scattering process can be handled as the similar procedure as that used in
neutrino absorptions. The scattering length, a distance to a point where the
MC particle scatters with matter, is determined probabilistically. Since it is
varied with neutrino flavor state, the scattering length is assigned to each
matrix element of MC particle. At the scattering point, we generate a new MC
particle, in which the neutrino energy and flight direction is chosen so as to
coincide with the scattered neutrinos444These momentum space values, i.e.,
$p^{\mu}$, can be determined through the same probabilistic way used in our
classical MC method (see Kato et al. (2020)).; i.e., it plays a role of the
particle after the scattering (see Fig. 1). The matrix element of $q$ in the
new MC particle is copied from that in the original particle, meanwhile we set
zero in the matrix component of $q$ of the original MC particle (but see below
for a corrected treatment to reduce the statistical noise). It should be
emphasized that this prescription guarantees the exact conservation of
neutrino numbers.

Although these treatments are capable of providing accurate solutions in
principle, the statistical noise inherent in MC methods is a concern. It would
affect self-interaction potentials, and then unphysical flavor oscillations
may be induced. Indeed, it potentially generate artificial ELN crossings, that
would kick in artificial fast neutrino flavor conversions. It should be
stressed that the problem is unavoidable as long as the number of MC particle
is finite555In any particle-based methods, we can always take an infinitesimal
solid angle so that the number of particle is zero, indicating that artificial
ELN crossings are ubiquitous at the small scale., indicating that the special
treatment is mandatory.

To overcome the difficulty, we adopt a smoothing prescription in evaluating
the self-interaction potential. Similar as mesh methods, we discretize the
6-dimensional phase space, and compute $\rho$ on each mesh by using Eq. 5. The
momentum integral of neutrinos in self-interaction potential is calculated by
using $\rho$. More specifically, we adopt the values of energy and flight
directions at the cell center value of the momentum space grid in computations
of the self-interaction potential. Although the prescription smears out a
small scale structure of neutrino momentum space, the statistical noise can be
controlled by changing the mesh resolution. This should be set appropriately
depending on the problem.

We develop another prescription of collision terms to further reduce the
statistical noise, namely effective mean free path method (hereafter we refer
to it as EMFP method). In general, large statistical noise potentially appears
when the MC particle experiences an absorption or scattering with matter. This
is attributed to the fact that a MC particle represents a bundle of neutrinos,
and that the physical states of all neutrinos change suddenly at each
absorption/scattering point. Although the most straightforward way to reduce
the noise is to increase the number of MC particles or to run many simulations
for the increase of statistics, we propose another way to address the issue;
the essence is as follows. We introduce a control parameter, $a$ ($0\leq a\leq
1$), that determines the rate of change in each MC particle by neutrino matter
interactions. Let us provide an example how we use $a$ for scattering
processes. At the scattering point, $a\times q_{ij}$ (where $i,j$ specify the
flavor of scattered neutrinos) is transferred from the original MC particle to
the new one (see Fig. 1), i.e., the change of $q$ is softened by $a$. To
compensate this effect, we reduce the scattering length by the same factor of
$a$. As we show in Sec. 8.1, this does not affect the physical solution of QKE
system. The merit of this prescription is that the statistical noise can be
reduced without pressing memory capacity in simulations.

One may wonder if the above prescription causes increasing numerical costs and
makes the simulation unfeasible. In fact, the number of scattering events
increases proportional to $1/a$. It should be noted, however, that this
problem is under control; indeed $a=1$ recovers the original treatment.
Importantly, this prescription is suited for problems of our interest:
interplay between collision terms and collective neutrino oscillations. The
reason is as follows. The time scale of collective neutrino oscillation is
usually much shorter than that of neutrino-matter collisions, indicating that
the numerical cost is dominated by transport- and oscillation terms. Indeed,
for the case with $a=1$, almost every MC particle does not interact with
matter during each time step ($\Delta t$) of simulations. By introducing $a$,
we can manage the balance of computational time among different terms. As a
result, we successfully increase the statistics without both pressing memories
and computational cost.

We make another important remark regarding the number of MC particles. In our
treatments of collision term, the number of MC particles tends to increase
with time, since new MC particles are generated in not only emission but also
scattering processes. We address this issue by the following prescription.
When the particle number almost exceeds the limit of computational resources,
we combine adjacent MC particles so that the number, energy, and moment of
neutrinos are conserved. Although the prescription violates the causality,
this effect is expected to be minor666Regarding the causality, the computation
of self-interaction potential is also an issue, since it is computed by
collecting MC particles located at different spatial locations. However, the
error can be estimated by the convergence study; indeed, the solution
approaches to the real by increasing the number of MC particles and decreasing
the phase space volume in the computation of $\rho$..

###  3.3 Flow chart

Figure 2: The flow chart for our QKE-MC solver. We highlight the parts
associated with neutrino-matter interactions as purple color.

To facilitate readers’ understanding, we provide a flow chart of our QKE-MC
solver in Fig. 2. It summarizes the overall procedure for the time evolution
of $t=t_{n}$ to $t_{n+1}$. We first evolve the MC particles by solving Eq. 4.
If the scattering and absorption processes are included, we change the flight
direction and energy of neutrinos during the time step. The spatial transport
of generated MC particles associated with scatterings (see Sec. 3.2) are also
handled there. It should be mentioned that we do not evolve $q$ of each MC
particle at this step.

As the next step, we compute the self-interaction potential from MC particles.
As mentioned in Sec. 3.2, we employ a smoothing prescription to evaluate it if
collision terms are taken into account. We then evolve $q$ and $\bar{q}$ by
solving Eq. 8.

There are two important remarks in the above procedure. First, we treat the
time evolution of $q$ as a similar procedure of operating-splitting method; in
which neutrino oscillations and collision terms are handled separately. As we
have described in Sec. 3.2, matrix elements of scattered MC particles needs to
be corrected, and the new generated one takes over them (see also Fig. 1).
This computation is done after solving Eq. 8, i.e., completing the transport-
and oscillation processes. As long as $\Delta t$ is small enough to resolve
the system, the artifact due to the prescription is negligible. Second, the
actual procedure of our QKE-MC solver is a bit more complicated than that
described above. This is because the high-order time integration (fourth-order
Runge-Kutta method) is implemented in our QKE-MC. We update the self-
interaction potential at each sub-step. We note that the high-order precision
is necessary only for solving Eq. 8, i.e., to determine the evolution of
neutrino flavor conversions. Except for the computation of flavor conversion,
we adopt the Eulerian integration in their time integration.

##  4 Basic assumptions in code tests

Although our QKE-MC solver is capable of handling the problems with arbitrary
number of neutrino flavors, we assume a two flavor system throughout the code
tests. In the approximation, the density matrix has four independent
components that can be written as,

$\displaystyle\rho=\left(\begin{array}[]{cc}\rho_{ee}&\rho_{ex}\\\
\rho_{ex}^{*}&\rho_{xx}\end{array}\right),\hskip
5.69054pt\bar{\rho}=\left(\begin{array}[]{cc}\bar{\rho}_{ee}&\bar{\rho}_{ex}\\\
\bar{\rho}_{ex}^{*}&\bar{\rho}_{xx}\end{array}\right).$ (13)

In the two flavor system, the expression with the polarization vector
$\mathbb{P}$ and $\bar{\mathbb{P}}$ in the flavor space is useful. It is
defined as,

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\rho_{ee}+\rho_{xx}\right)I+\frac{1}{2}\mathbb{P}\cdot\mathbb{\sigma},$
$\displaystyle\bar{\rho}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\bar{\rho}_{ee}+\bar{\rho}_{xx}\right)I+\frac{1}{2}\mathbb{\bar{P}}\cdot\mathbb{\sigma^{*}},$
(14)

where $\sigma$ and $I$ denotes the Pauli matrix and Unit matrix. The component
of polarization vector can be written as,

$\displaystyle\mathbb{P}=\left(\begin{array}[]{c}2{\rm Re}{\rho_{ex}}\\\
-2{\rm Im}{\rho_{ex}}\\\ \rho_{ee}-\rho_{xx}\end{array}\right),\ \
\mathbb{\bar{P}}=\left(\begin{array}[]{c}2{\rm Re}{\bar{\rho}_{ex}}\\\ 2{\rm
Im}{\bar{\rho}_{ex}}\\\ \bar{\rho}_{ee}-\bar{\rho}_{xx}\end{array}\right).\ \
\ \ $ (21)

Using the poralization vectors, Eq. 1 can be written as Bloch-type equations:

$\displaystyle\frac{D}{d\lambda}\mathbb{P}=\left[\omega\mathbb{B}+\sqrt{2}G_{F}\left(n_{e}-\bar{n}_{e}-n_{x}+\bar{n}_{x}\right)\mathbb{u}\right.$
$\displaystyle\ \ \ \ \ \ \ \ \left.+\sqrt{2}G_{F}\int
dV_{p^{\prime}}\frac{p^{\mu}p_{\mu}^{\prime}}{E_{\nu}E_{\nu}^{\prime}}\left(\mathbb{P}^{\prime}-\bar{\mathbb{P}}^{\prime}\right)\right]\times\mathbb{P},$
$\displaystyle\frac{D}{d\lambda}\mathbb{\bar{P}}=\left[-\omega\mathbb{B}+\sqrt{2}G_{F}\left(n_{e}-\bar{n}_{e}-n_{x}+\bar{n}_{x}\right)\mathbb{u}\right.$
$\displaystyle\ \ \ \ \ \ \ \ \left.+\sqrt{2}G_{F}\int
dV_{p^{\prime}}\frac{p^{\mu}p_{\mu}^{\prime}}{E_{\nu}E_{\nu}^{\prime}}\left(\mathbb{P}^{\prime}-\bar{\mathbb{P}}^{\prime}\right)\right]\times\mathbb{\bar{P}},\
\ \ \ \ $ (22)

where $\theta$ is the angle between the propagation directions of two
neutrinos and $\omega\equiv\Delta m^{2}/2E_{\nu}$ is the vacuum oscillation
frequency with the mass difference $\Delta m^{2}\equiv m_{2}^{2}-m_{1}^{2}$.
In the expression, we use the following expressions for Hamiltonian
potentials,

$\displaystyle H_{\rm{vac}}=\frac{1}{2}\omega\mathbb{B}\cdot\mathbb{\sigma},$
$\displaystyle
H_{\rm{mat}}=\frac{\sqrt{2}}{2}G_{F}\left(n_{e}-\bar{n}_{e}-n_{x}+\bar{n}_{x}\right)\mathbb{u}\cdot\mathbb{\sigma},$
$\displaystyle H_{\nu\nu}=\frac{\sqrt{2}}{2}G_{F}\int
dV_{p^{\prime}}\frac{p^{\mu}p_{\mu}^{\prime}}{E_{\nu}E_{\nu}^{\prime}}\left(\mathbb{P}^{\prime}-\mathbb{\bar{P}}^{\prime}\right)\cdot\sigma,$
(23)

with $\mathbb{u}=(0,0,1)$ and
$\mathbb{B}=\left(\sin{2\theta_{0}},0,-\cos{2\theta_{0}}\right)$. $\theta_{0}$
denotes the mixing angle. In this expression, we subtract the trace part,
since they do not affect flavor conversions.

In the following code tests, we take the flat space-time. We also assume that
the matter distribution is static and homogeneous. It should be mentioned that
neutrino distributions are assumed to be homogeneous throughout our code
tests. More detailed tests and scientific simulations including neutrino
transport under inhomogeneous backgrounds are postponed to future work.

##  5 Vacuum oscillations

We start with tests of vacuum oscillations. By virtue of the linear system, we
can solve the time evolution of $\rho$ analytically. They can be written in
the flavor-basis as,

$\displaystyle\rho_{ee}=\left[1-\sin^{2}{2\theta_{0}}\sin^{2}{\frac{\omega
t}{2}}\right]\rho_{ee}^{0}$ $\displaystyle\ \ \ \ \
+\sin^{2}{2\theta_{0}}\sin^{2}{\frac{\omega t}{2}}\rho_{xx}^{0}$
$\displaystyle\ \ \ \ \
-2\cos{2\theta_{0}}\sin{2\theta_{0}}\sin^{2}{\frac{\omega t}{2}}{\rm
Re}\rho_{ex}^{0}$ $\displaystyle\ \ \ \ \ -\sin{2\theta_{0}}\sin{\omega t}{\rm
Im}\rho_{ex}^{0},$ (24)
$\displaystyle\rho_{xx}=\sin^{2}{2\theta_{0}}\sin^{2}{\frac{\omega
t}{2}}\rho_{ee}^{0}$ $\displaystyle\ \ \ \ \
+\left[1-\sin^{2}{2\theta_{0}}\sin^{2}{\frac{\omega
t}{2}}\right]\rho_{xx}^{0}$ $\displaystyle\ \ \ \ \
+2\cos{2\theta_{0}}\sin{2\theta_{0}}\sin^{2}{\frac{\omega t}{2}}{\rm
Re}\rho_{ex}^{0}$ $\displaystyle\ \ \ \ \ +\sin{2\theta_{0}}\sin{\omega t}{\rm
Im}\rho_{ex}^{0},$ (25) $\displaystyle{\rm
Re}\rho_{ex}=\cos{2\theta_{0}}\sin{2\theta_{0}}\sin^{2}{\frac{\omega
t}{2}}\left(\rho_{xx}^{0}-\rho_{ee}^{0}\right)$ $\displaystyle\ \ \ \ \
+\left[1-2\cos^{2}{2\theta_{0}}\sin^{2}{\frac{\omega t}{2}}\right]{\rm
Re}\rho_{ex}^{0}$ $\displaystyle\ \ \ \ \ -\cos{2\theta_{0}}\sin{\omega t}{\rm
Im}\rho_{ex}^{0},$ (26) $\displaystyle{\rm
Im}\rho_{ex}=\frac{1}{2}\sin{2\theta_{0}}\sin{\omega
t}\left(\rho_{ee}^{0}-\rho_{xx}^{0}\right)$ $\displaystyle\ \ \ \ \
+\cos{2\theta_{0}}\sin{\omega t}{\rm Re}\rho_{ex}^{0}$ $\displaystyle\ \ \ \ \
+\left[1-2\sin^{2}{\frac{\omega t}{2}}\right]{\rm Im}\rho_{ex}^{0},$ (27)

where $\rho_{ij}^{0}$ denotes the initial density matrices of neutrinos (the
$i$ and $j$ indexes take either $e$ or $x$).

In this test, we adopt the following mass and mixing parameters: $\Delta
m^{2}=(-)2.45\times 10^{-15}\ \rm{MeV}^{2}$ and
$\sin^{2}{\theta_{0}}=2.24\times 10^{-2}\ (2.26\times 10^{-2})$ for the normal
mass ordering (the inverted mass ordering) (Olive & Particle Data Group,
2014). Hereafter, we refer to each mass hierarchy as NO and IO, respectively.
The neutrino energy ($E_{\nu}$) is assumed to be $20$ MeV. As an initial
condition, we emit a single MC particle (the flight direction can be arbitrary
chosen) with setting a purely $\nu_{e}$ state, i.e., the $ee$ component of $q$
is 1 and others are set to be zero777Since it is a linear system (see also Eq.
1 or 8), the result is scale free, indicating that an arbitrary positive value
can be chosen.. During the evolution, time step ($\Delta t$) is constant with
$\omega\Delta t/2\pi=10^{-3}$.

Figures 3 show the time evolution of each matrix element of $q$ in cases with
NO (left panel) and IO (right panel). The QKE-MC results are shown in solid
lines while the analytical solutions are in black dashed lines. In the bottom
figures, the errors from analytic solutions are displayed. As shown in these
figures, QKE-MC results are good agreement with analytic solutions for both
mass hierarchies.

Figure 3: The time evolution of each matrix element of $q_{s}$ in the test of
vacuum oscillations. Left and right panels show the normal and inverted mass
ordering, respectively. Solid lines represented numerical solutions and dashed
black lines are analytical solutions. In the bottom panel, the errors from the
analytical solutions are displayed.

Figure 4 portrays the convergence of our code with respect to $\Delta t$. In
this check, we compare $q_{ee}$ at the time of $\omega t/2\pi=10$ between the
numerical (QKE-MC) and analytical solutions. We confirm that our code achieved
a 4th-order convergence. We note that, the round-off-error is seen in the case
with $\omega\Delta t/2\pi=10^{-4}$, which results in degrading the
convergence.

Figure 4: Numerical error as a function of $\Delta t$ for tests of vacuum
oscillation in the case of the normal mass ordering (NO).

The norm of polarization vector, $|P|$, is a conserved quantity in this test;
hence, it is also useful to measure the accuracy of our code. In Fig. 5, we
show the time evolution of $|1-|P||$ (we note that $|P|$ is 1 at the initial
condition) for two different time step ($\omega\Delta t/2\pi=10^{-2}$ and
$10^{-3}$). We display the results with NO oscillation in the figure. This
result suggests that $|P|$ is well conserved in the test, and that the
accuracy is improved with decreasing $\Delta t$.

Figure 5: The time evolution of $|1-|P||$ for the test of vacuum oscillation.
We show the results for simulations with two different time resolutions:
$\omega\Delta t/2\pi=10^{-3}$ (solid line) and $10^{-2}$ (dashed line). In
this figure, we assume NO for mass hierarchy.

##  6 Neutrino oscillation in matter

As long as the matter potential, $H_{\rm mat}$, is static and homogeneous, we
can obtain exact solutions of QKE (with respect to the time evolution of the
system) with vacuum and matter potentials. If only electrons surround
neutrinos, the exact solution can be obtained by replacing
$\theta_{0}\rightarrow\theta_{m}$ and $\omega\rightarrow\omega_{m}\equiv\Delta
m_{m}^{2}/2E_{\nu}$ in Eqs. 24-27, where $\theta_{m}$ and $\Delta m_{m}^{2}$
can be written as,

$\displaystyle\sin{2\theta_{m}}$ $\displaystyle=$
$\displaystyle\frac{\sin{2\theta_{0}}}{\sqrt{\left(\cos{2\theta_{0}}\pm
A\right)^{2}+\sin^{2}{2\theta_{0}}}},$ (28) $\displaystyle\Delta m_{m}^{2}$
$\displaystyle=$ $\displaystyle\Delta m^{2}\sqrt{\left(\cos{2\theta_{0}}\pm
A\right)^{2}+\sin^{2}{2\theta_{0}}},$ (29)

with $A=2\sqrt{2}G_{F}n_{e}E_{\nu}/\Delta m^{2}$, where $n_{e}$ denotes the
number density of electrons. Regarding the sign in front of $A$ in Eqs. 28 and
29, positive and negative signs are taken in the cases of neutrinos and anti-
neutrinos, respectively. As shown in Eq. 28, the large mixing (resonance)
occurs when the denominator of right hand side of Eq. 28 is zero. This is
equivalent to the condition of $n_{e}=(-)n_{e0}$ for (anti-)neutrinos, where

$\displaystyle n_{e0}=\frac{\Delta
m^{2}\cos{2\theta_{0}}}{2\sqrt{2}G_{F}E_{\nu}}.$ (30)

In Fig. 6, we show the time evolution of $q_{ee}$ for various cases of
$n_{e}$. In these tests, the initial condition of neutrinos is the same as
that used in the test of vacuum oscillation (see Sec. 5). The result of
analytic solutions are plotted with black dashed lines. This figure suggests
that our code accurately demonstrates the neutrino oscillation in matter.

Figure 6: We show the time evolution of $\nu_{e}$ for the test of MSW neutrino
oscillation. Color lines show the results for different $n_{e}$, and the black
dashed lines represent the analytic solutions

. In this figure, we assume NO for mass hierarchy.

It may be interest to see the resonance case ($n_{e}=n_{e0}$) in more detail.
Figs. 7 display the same quantities as displayed in Figs. 3. Fig. 8 shows the
result of convergence check. The error is measured at $t=2\pi/\omega\times
10$, which is the same as that in the vacuum test (see also Fig. 4). These
tests suggest that the MSW resonance can be handled accurately in our QKE-MC
solver.

Figure 7: The same as Fig. 3 but for MSW effects. We adopt $n_{e}=n_{e0}$
(resonant oscillation). Left and right panels show the results with normal and
inverted mass ordering, respectively. Figure 8: The same as Fig. 4 but for a
MSW test. We adopt $n_{e}=n_{e0}$ (resonant oscillation) in these simulations.

##  7 Collective oscillation

We perform three representative tests for collective neutrino oscillations.
One of them is a collective neutrino oscillation in isotropic angular
distributions. This is the same problem as done in Hannestad et al. (2006).
The result is summarized in Sec. 7.1. In Sec. 7.2, we test a two-beam
collective oscillation following Chakraborty et al. (2016). We compare the
results with the linear analysis quantitatively. In Sec. 7.3, we carry out a
test of fast flavor conversion. We evolve the system, in which the anisotropic
neutrino distribution with a single ELN-crossing is set as an initial
condition, by using QKE-MC solver. It should be mentioned that, although the
physical setup is exactly the same as one of the simulations in Shalgar &
Tamborra (2021a), we do not reproduce their results. To address the issue,
another simulation is carried out for the same problem but using our
independent code (finite-difference method), that strengthened the reliability
of QKE-MC solver. We note that, thanks to the cooperation by the authors of
Shalgar & Tamborra (2021a), we identify the source of discrepancy. We finally
reached an agreement with them.

###  7.1 isotropic distribution

In this test, we assume that neutrinos and anti-neutrinos have the same
isotropic distributions. Regarding the matter potential, we assume
$\bar{n}_{e}=n_{x}=\bar{n}_{x}=0$. Following the convention of Hannestad et
al. (2006), we describe the QKE by using polarization vectors (see also Eq.
22),

$\displaystyle\frac{\partial}{\partial
t}\mathbb{P^{\prime}}=\left[\omega\mathbb{B}+\lambda\mathbb{L}+\mu\left(\mathbb{P^{\prime}}-\bar{\mathbb{P}}^{\prime}\right)\right]\times\mathbb{P^{\prime}},$
$\displaystyle\frac{\partial}{\partial
t}\mathbb{\bar{P}}^{\prime}=\left[-\omega\mathbb{B}+\lambda\mathbb{L}+\mu\left(\mathbb{P^{\prime}}-\bar{\mathbb{P}}^{\prime}\right)\right]\times\mathbb{\bar{P}}^{\prime},$
(31)

with $\lambda=\sqrt{2}G_{F}n_{e}$, $\mu=\sqrt{2}G_{F}n_{\nu}$
($n_{\nu}=n_{\nu_{e}}=n_{\bar{\nu}_{e}}$) and $\mathbb{L}={\rm diag}[1,0]$. In
the expression, we normalize the polarization vectors as
$\mathbb{P}^{\prime}=\mathbb{P}/n_{\nu}$ and
$\mathbb{\bar{P}}^{\prime}=\mathbb{\bar{P}}/n_{\nu}$.

As the first test, we ignore the matter potential, i.e., $\lambda$ is set to
be zero. In this case, Eq. 31 can be rewritten as,

$\displaystyle\mathbb{\dot{Q}}$ $\displaystyle=$
$\displaystyle\mu\mathbb{D}\times\mathbb{Q},$ $\displaystyle\mathbb{\dot{D}}$
$\displaystyle=$ $\displaystyle\omega\mathbb{B}\times\mathbb{Q},$ (32)

with $\mathbb{Q}=\mathbb{P}+\mathbb{\bar{P}}-\omega\mathbb{B}/\mu$ and
$\mathbb{D}=\mathbb{P}-\mathbb{\bar{P}}$. As shown in Hannestad et al. (2006),
we have four conserved quantities in the system: $|\mathbb{Q}|$,
$\mathbb{D}\cdot\mathbb{B}$, $\mathbb{D}\cdot\mathbb{Q}$,
$H\equiv\omega\mathbb{B}\cdot{Q}+\mu|\mathbb{D}|^{2}/2$. They are useful to
check the accuracy of our QKE-MC solver. Setting the initial condition with
$\mathbb{P}^{\prime}(0)=\mathbb{\bar{P}^{\prime}}(0)=(0,0,1)$, we obtain

$\displaystyle\mathbb{D}=\frac{1}{\mu}\frac{\mathbb{Q}\times\mathbb{\dot{Q}}}{|\mathbb{Q}|^{2}}.$
(33)

Eq. 33 is an analogous equation to that describing the dynamics of charged
particles propagating in electric field. It can be further cast into the
following form,

$\displaystyle\ddot{\phi}=-\omega\mu Q\sin{\left(\phi+2\theta_{0}\right)},$
(34)

with $\mathbb{Q}=Q(\sin{\phi},0,\cos{\phi})$. This suggests that the dynamics
is analogous to a pendulum. We refer readers to Hannestad et al. (2006) for
more detail.

In the numerical setup of QKE-MC solver, we inject one MC particle with
$E_{\nu}=20$ MeV initially with purely $\nu_{e}$ and $\bar{\nu}_{e}$ states.
We note that, their flight direction can be arbitrary chosen, meanwhile we
ignore the angular dependence of neutrinos in the evaluation of the neutrino
self-interaction potential during the time evolution. By virtue of the
prescription, only the isotropic component of the potential remains. It should
be noted that it is possible to set isotropic angular distributions of
neutrinos by generating many MC particles, and compute the self-interaction
potential without the above prescription. On the other hand, neutrino flight
directions of each MC particle are constant with time, and the neutrino
oscillation is identical among all MC particles. This indicates that
computations of many MC particles are wasteful. Our prescription does not
compromise the analysis and saves the computational resource. The number of
neutrinos, i.e., the $ee$ component of $q$ (and $\bar{q}$) is chosen so that
$\mu$ becomes $10\omega$ ($\Delta m^{2}=-2.45\times 10^{-15}\ {\rm MeV^{2}}$).
The mixing angle is set as $\theta_{0}=0.01$ to reproduce the result of
Hannestad et al. (2006). In this test, we adopt a constant time step with
$\mu\Delta t=10^{-3}$.

Figure 9 portrays the time evolution of the z-components for $\mathbb{P}$ and
$\mathbb{\bar{P}}$ computed by our QKE-MC solver. This is in good agreement
with the result of Hannestad et al. (2006) (see Fig.1 in their paper). To see
the capability of our code more quantitatively, we compare to the solutions
obtained by solving Eq. 34; the results are summarized in Fig. 10. We confirm
that our QKE-MC solver provides a consistent result with the solutions of Eq.
34; indeed, the deviation ($|\Delta Q_{z}|$) is $\lesssim 10^{-3}$ up to the
time of $\omega t=80$ (see the bottom panel of Fig. 10).

Figure 9: The time evolution of the z-components in polarization vectors in a
test of isotropic collective neutrino oscillation without matter potential. We
show the results for $\mu\Delta t=10^{-3}$ in this figure. See also Fig.1 in
Hannestad et al. (2006) that corresponds to a counterpart of this figure.
Figure 10: Time evolution of $Q_{z}$ for a test of isotropic collective
neutrino oscillation without matter potential. In the top panel, we compare
the results between QKE-MC (orange-solid line) and those obtained by solving
Eq. 34 (black-dashed line). The difference, $|\Delta Q_{z}|$, is displayed in
the bottom panel. We adopt a constant time step ($\mu\Delta t=10^{-3}$) in our
QKE-MC simulation.

Figure. 11 displays the time evolution of conserved quantities. In this
figure, we add the result with $\mu\Delta t=10^{-2}$ to see the dependence of
time step. All the quantities, or $|\mathbb{P}|$, $|\mathbb{Q}|$, $H$,
$\mathbb{D}\cdot\mathbb{B}$ and $\mathbb{D}\cdot\mathbb{Q}$, are well
conserved in our results, and we confirm that the smaller $\Delta t$ tends to
improve the accuracy of the conserved quantities.

Figure 11: The time evolution of errors of conserved quantities for a test of
isotropic collective neutrino oscillation without matter potential. From top
to bottom, we display the results of $1-|\mathbb{P}|$,
$1-|\mathbb{Q}|/|\mathbb{Q_{0}}|$, $1-H/H_{0}$, $\mathbb{D}\cdot\mathbb{B}$
and $\mathbb{D}\cdot\mathbb{Q}$, respectively. We show the results with two
different time steps in our QKE-MC simulation: $\mu\Delta t=10^{-3}$ (solid-
line) and $10^{-2}$ (dashed-line).

As another test of isotropic collective neutrino oscillations, we include the
matter term in addition to the self-interaction potential (see Eq. 31). Same
as Hannestad et al. (2006), we vary the value of $\lambda$:
$\lambda=10^{2}\omega$, $10^{3}\omega$ and $10^{4}\omega$. Figure 12 shows the
time evolution of $P_{z}$ and $\bar{P}_{z}$ obtained by our QKE-MC (see also
Fig. 9 for comparison to the case without matter potential). We confirm that
our results reproduce the same results reported in Hannestad et al. (2006)
(see Fig.2 in the paper). This suggests that our code is capable of handling
the collective neutrino oscillations in isotropic distributions.

Figure 12: The same as Fig. 9 but for the test including matter potential.
Different colors denote results for different values of matter potential:
$\lambda=10^{2}\omega,10^{3}\omega$, and $10^{4}\omega$. We adopt a constant
time step ($\mu\Delta t=10^{-3}$) in this test .

###  7.2 Two-beam collective oscillation

Figure 13: The time evolution of the real part of $n_{ex}$ in the test of two-
beam collective oscillation. In the left and right panels, we adopt the
different asymmetric parameter, $\alpha=0$ (left) and $0.5$ (right). Orange
and black lines denote the QKE-MC results and the those computed from eigen
frequency, respectively. We show the results with two different
$\mu^{\prime}$’s in top and bottom panels. The top and bottom panels
correspond to the stable and unstable system to flavor conversion. See text
for more detail.

We employ a two-beam model as a test of collective oscillation. In this model,
neutrinos and anti-neutrinos, propagating in the opposite direction, intersect
each other. The characteristics of flavor conversion has been well studied
(see, e.g., Chakraborty et al., 2016). We also note that the dispersion
relation of the instability can be written in a concise way (see below).
Hence, we compare the result of QKE-MC simulation to the analytic solution in
the linear regime.

Following Chakraborty et al. (2016), we give an overview of the derivation of
eigen frequency. In the system, we obtain the time evolution of the number
densities from QKE as,

$\displaystyle i\frac{\partial n_{ex,1}}{\partial t}$ $\displaystyle=$
$\displaystyle\left(-\frac{\Delta
m^{2}}{2E_{\nu}}-2\sqrt{2}G_{F}\bar{n}_{ee,2}\right)n_{ex,1}$ (35)
$\displaystyle\ \ \ \ \ \ \ +2\sqrt{2}G_{F}\bar{n}_{ex,2}^{\ast}n_{ee,1},$
$\displaystyle i\frac{\partial\bar{n}_{ex,2}}{\partial t}$ $\displaystyle=$
$\displaystyle\left(-\frac{\Delta
m^{2}}{2E_{\nu}}-2\sqrt{2}G_{F}n_{ee,1}\right)\bar{n}_{ex,2}$ (36)
$\displaystyle\ \ \ \ \ \ \ +2\sqrt{2}G_{F}n_{ex,1}^{\ast}\bar{n}_{ee,2},$

where indices 1 and 2 distinguishes the flight direction of neutrinos, i.e.,
specifying neutrinos and anti-neutrinos, respectively; we assume
$\theta_{0}=0$ and a monochromatic neutrino energy with $E_{\nu}$. It should
be mentioned that we integrate over the momentum space in these equations.
Before linearizing the equation, we define two new variables, $g$ and $S$, as

$\displaystyle n_{i}-\bar{n}_{i}^{\ast}=\frac{1}{2}{\rm
Tr}\left(n_{i}-\bar{n}_{i}^{\ast}\right)+\frac{g_{i}}{2}\left(\begin{array}[]{cc}1&S_{i}\\\
S_{i}^{\ast}&-1\end{array}\right),$ (39)

where $g_{i}=n_{ee,i}-n_{xx,i}-\bar{n}_{ee,i}+\bar{n}_{xx,i}$. In our test
case, $g_{1}=n_{ee,1}$, $g_{2}=-\bar{n}_{ee,2}$, $n_{ex,1}=g_{1}S_{1}/2$ and
$\bar{n}_{ex,2}=-g_{2}S_{2}^{\ast}/2$. By using the new variables, the
linearized equations can be written as

$\displaystyle i\frac{\partial S_{1}}{\partial t}$ $\displaystyle=$
$\displaystyle\left(-\frac{\Delta
m^{2}}{2E_{\nu}}+2\sqrt{2}G_{F}g_{2}\right)S_{1}-2\sqrt{2}G_{F}g_{2}S_{2},$
(40) $\displaystyle i\frac{\partial S_{2}^{\ast}}{\partial t}$
$\displaystyle=$ $\displaystyle-\left(\frac{\Delta
m^{2}}{2E_{\nu}}+2\sqrt{2}G_{F}g_{1}\right)S_{2}^{\ast}+2\sqrt{2}G_{F}g_{1}S_{1}^{\ast}.\
\ \ \ \ \ \ $ (41)

Adopting plane wave expansion with $S_{i}=Q_{i}\exp{(-i\Omega t)}$, the eigen
frequency can be written as,

$\displaystyle\Omega$ $\displaystyle=$
$\displaystyle\sqrt{2}G_{F}\left(g_{1}+g_{2}\right)$ (42)
$\displaystyle\pm\sqrt{2G_{F}^{2}\left(g_{1}+g_{2}\right)^{2}+\omega^{\prime}\left(\omega^{\prime}+2\sqrt{2}G_{F}\left(g_{1}-g_{2}\right)\right)},\
\ \ \ \ \ \ \ $

with $\omega^{\prime}=-\omega$. In this test, neutrino distributions are
assumed to be homogeneous (see Chakraborty et al., 2016, for more general
cases). For convenience, we introduce an asymmetry parameter between neutrinos
and anti-neutrinos, $\alpha$. We rewrite the number densities of neutrinos and
anti-neutrinos by using $\alpha$ as $n_{ee}=(1+\alpha)n_{t}$ and
$\bar{n}_{ee}=(1-\alpha)n_{t}$, respectively, where $n_{t}$ represents their
average, i.e., $n_{t}=0.5(n_{ee}+\bar{n}_{ee})$. The dispersion relation can
be rewritten as,

$\displaystyle\Omega_{\pm}=2\alpha\mu^{\prime}\pm\sqrt{\left(2\alpha\mu^{\prime}\right)^{2}+\omega^{\prime}\left(\omega^{\prime}-4\mu^{\prime}\right)},$
(43)

with $\mu^{\prime}=\sqrt{2}G_{F}n_{t}$. The imaginary part of $\Omega$ has a
non-zero value in cases with
$\left(2\alpha\mu^{\prime}\right)^{2}+\omega^{\prime}\left(\omega^{\prime}-4\mu^{\prime}\right)<0$.
In this condition, the slow flavor instability kicks in. In other cases, the
flavor conversion occurs periodically with constant frequencies. The
corresponding eigenvector of $\Omega_{-}$, for example, can be written as

$\displaystyle\left(\begin{array}[]{c}Q_{1}\\\
Q_{2}\end{array}\right)=\left(\begin{array}[]{c}2\mu^{\prime}\left(\alpha-1\right)\\\
\omega^{\prime}+2\mu^{\prime}\left(\alpha-1\right)-\Omega_{-}\end{array}\right)b.$
(48)

This is used for an initial condition in our QKE-MC solver.

To reproduce the two-beam environment in our QKE-MC simulation, we prepare two
MC particles. One of them is for neutrinos, and the other is for anti-
neutrinos. The flight direction of each particle is opposite each other. The
neutrino energy is set to be $E_{\nu}=50$ MeV. To determine $\Delta t$, we
compare the vacuum- and collective oscillation timescale at each time step,
and then adopt the smaller one. We then multiply it by a factor of $10^{-4}$.
We demonstrate four simulations by changing $\alpha$ ($\alpha=0$ and $0.5$)
and $\mu^{\prime}$ ($\mu^{\prime}=\omega^{\prime}/8$ and $\omega^{\prime}/2$).
In this test, we adopt the value of mass difference: $\Delta m^{2}=-2.45\times
10^{-3}{\rm eV}^{2}$. The corresponding time scale of $1/\omega$ is
2.69$\times 10^{-5}$ s, and the number density of $n_{t}$ with
$\mu^{\prime}=\omega^{\prime}$ is $1.93\times 10^{26}{\rm cm^{-3}}$. The
initial condition of $q$ and $\bar{q}$ is set so as to coincide with an
eigenvector with $\Omega_{-}$ with ${\rm Re}b={\rm
Im}b=10^{-6}/2\mu^{\prime}(\alpha-1)$ in Eq. 48.

Figure 13 shows the time evolution of the real part of $n_{ex}$ in the linear
regime ($|S|\ll 1$). Orange and black lines represent the QKE-MC results and
those computed from the linearlized QKE. Top and bottom panels show the
results for $\mu^{\prime}=\omega^{\prime}/8$ and $\omega^{\prime}/2$, where
the former and latter systems are the stable and unstable with respect to slow
flavor conversion. Regardless of the asymmetric parameter and stability, our
QKE-MC code can accurately reproduce the analytic solutions in the linear
phase.

###  7.3 fast flavor conversion in axisymmetric momentum space

As a test of fast flavor conversion, we carry out simulations with the same
physical condition of Shalgar & Tamborra (2021a), in which the angular
distributions of $\nu_{e}$ and $\bar{\nu}_{e}$ has initially a single crossing
in momentum space, indicating that it is unstable to fast flavor conversion.
We apply our QKE-MC solver to evolve the system under homogeneous
approximation, and then we compare the results to those in Shalgar & Tamborra
(2021a).

In this test, we assume a monochromatic neutrino energy with $E_{\nu}$ and
axial symmetry in momentum space. For convenience, we define the azimuthal-
integrated density matrix as,

$\displaystyle\rho_{a}(\theta_{\nu})=\int\rho(\theta_{\nu},\phi_{\nu})d\phi_{\nu}.$
(49)

The governing equation in the system can be written as

$\displaystyle i\frac{\partial\rho_{a}}{\partial t}$ $\displaystyle=$
$\displaystyle[H_{\rm vac}^{\prime}+H_{\rm
mat}+H_{\nu\nu}^{\prime},\rho_{a}],$ $\displaystyle
i\frac{\partial\bar{\rho}_{a}}{\partial t}$ $\displaystyle=$
$\displaystyle[H_{\rm vac}^{\prime}-H_{\rm
mat}^{\ast}-H_{\nu\nu}^{\prime\ast},\bar{\rho}_{a}],$ (50)

with

$\displaystyle H_{\rm vac}^{\prime}$ $\displaystyle=$ $\displaystyle
U\frac{1}{2E_{\nu}}\left(\begin{array}[]{cc}m_{1}^{2}&0\\\
0&m_{2}^{2}\end{array}\right)U\dagger,$ (53) $\displaystyle
H_{\nu\nu}^{\prime}$ $\displaystyle=$
$\displaystyle\sqrt{2}G_{F}\int\frac{d\cos{\theta_{\nu}^{\prime}}}{\left(2\pi\right)^{3}}\left(1-\cos{\theta_{\nu}}\cos{\theta_{\nu}^{\prime}}\right)$
(54)
$\displaystyle\times\left(\rho_{a}^{\prime}-\bar{\rho}_{a}^{\prime\ast}\right).$

We set an initial condition,

$\displaystyle\rho_{ee,a0}$ $\displaystyle=$ $\displaystyle 0.5n_{\nu},$
$\displaystyle\bar{\rho}_{ee,a0}$ $\displaystyle=$
$\displaystyle\left[0.47+0.05\exp{\left(-\left(\cos{\theta_{\nu}}-1\right)^{2}\right)}\right]n_{\nu},$
(55)

where $n_{\nu}$ denotes the total number density of $\nu_{e}$. We assume
$\mu=10^{5}{\rm km^{-1}}$, $\theta_{0}=10^{-6}$, $\Delta m^{2}=2.50\times
10^{-6}{\rm eV^{2}}$, and $E_{\nu}=50$ MeV, which leads to $\omega=1.27\times
10^{-4}{\rm km}^{-1}$. We note that this setup is the same as the case A in
Shalgar & Tamborra (2021a).

For convenience, we provide the time scale under the choice of these
parameters. They are $3.33\times 10^{-11}$s and $2.67\times 10^{-2}$s for
neutrino self-interaction and vacuum oscillation, respectively. Since
$\mu\gg\omega$ is satisfied, the neutrino oscillation is dominated by fast
flavor conversion888However, the vacuum term is important to generate a seed
perturbation in the system. Since there are no neutrinos in the mixing flavor-
states at initial time, the flavor conversion does not occur without the
vacuum contribution..

In QKE-MC simulations, we inject MC particles isotropically. To see the
convergence, we change the total number of MC particle as 64, 128, 256 and
512. Following Eq. 55, we determine the values of $q$ and $\bar{q}$ in each MC
particle. During the time evolution, $\Delta t$ is not constant. It is
essentially determined by the maximum difference of $\rho_{a}$ and
$\bar{\rho}_{a}$ at each time step. We compute the effective ELN number
density by assuming that the maximum difference is distributed isotropically.
This gives a conservative time scale of the fast flavor conversion. We further
divide it by a factor of $dt_{\rm fac}$. We adopt $dt_{\rm fac}=50$ as a
representative case. It should be mentioned that the MC particles do not
change their momentum state in this test due to no collision terms, indicating
that the statistical noise is not an issue in this problem. Hence, we do not
need a special treatment to reduce the noise (see also Sec. 3.2).

Figure 14 shows the time evolution of a transition probability, which is
defined by

$\displaystyle<P_{ex}>=1-\frac{\int\rho_{ee,a}d\cos{\theta_{\nu}}}{\int\rho_{ee,a0}d\cos{\theta_{\nu}}}.$
(56)

Solid lines represent our results and color distinguishes the total number of
MC particles. In the linear phase ($\lesssim 10^{-6}$ s), the resolution
dependence (with respect to the total number of neutrinos) is minor; indeed
the first peak profile in Fig. 14 is almost identical among all simulations.
However, the deviation is obvious at the latter phase; the peak profile in the
low resolution tends to appear earlier. On the other hand, we find no clear
differences between the simulation of 256 and 512 MC particles up to the third
peak.

In Shalgar & Tamborra (2021a), there is a counterpart figure, the blue line in
the top left panel of Fig. 2 in their paper, to our Fig. 14. By comparing
between the two results, the flavor conversion in their simulation tends to
evolve slower than ours; for instance, the first peak of $<P_{ex}>$ appears at
$t\sim 0.6\times 10^{-5}$ s. We address the issue by the following way. At
first, we perform another simulation, in which we solve QKE with poralization
vectors (see Eq. 22) by using finite-difference method. This is a
deterministic- and mesh-based solver, indicating that the numerical artifacts
inherent in our MC solver are not involved. We carry out two simulations with
two different resolutions; the high resolution simulation adopts $\Delta
t=0.5{\rm cm}/c$ and 256 angular grids, while the low resolution one adopts
$\Delta t=1{\rm cm}/c$ and 128 angular grids. If the MC code has some issues,
the detailed comparison between the two independent code illustrates the
origin of the discrepancy from Shalgar & Tamborra (2021a).

We found, however, that the finite difference method shows a consistent result
with QKE-MC, which can be seen in Fig. 14 (black dashed- and dashed lines). It
should also be mentioned that the difference between both codes becomes
smaller with increasing resolutions that lends confidence to our QKE-MC code.
In the wake of the result, we have contacted the authors of Shalgar & Tamborra
(2021a), and they mentioned that there are typos in the paper999In their
simulation, the neutrino self-interaction potential was four time smaller than
that described in Shalgar & Tamborra (2021a) (private communication)., and
they reproduced our results. Hence, the issue of discrepancy is completely
resolved.

Figure 14: The time evolution of transition probability $<P_{ex}>$ for a test
of fast flavor conversion (Sec. 7.3). Solid lines represent QKE-MC results and
the color distinguishes the total number of MC particles. Black lines denote
results of finite-difference solver. The line type distinguishes the
resolution. The high resolution case adopts 256 angular grids and $\Delta
t=0.5{\rm cm}/c$, while the low resolution one uses 128 angular grids and
$\Delta t=1{\rm cm}/c$. $c$ denotes the speed of light in cgs unit.

To see the flavor conversion in more detail, we display the time evolution of
neutrino angular distributions computed by our QKE-MC solver up to the time of
the first peak in Fig. 15. In this figure, we show the result with 256 MC
particles. As shown in the plot, the flavor conversion occurs around the ELN
crossing point; the qualitative trend is the same as that found in Shalgar &
Tamborra (2021a).

Figure 15: The time evolution of neutrino angular distributions for a test of
fast flavor conversion (Sec. 7.3). We display the result up to the time of the
first peak. Color distinguishes the time. Solid and dashed lines denote
$\rho_{a,ee}$ and $\bar{\rho}_{a,ee}$, respectively. In this figure, we show
the result of QKE-MC simulation with 256 MC particles.

Finally, we check the dependence of $\Delta t$, that can be controlled by
$dt_{\rm fac}$. The results are displayed in Figure 16. We vary $dt_{\rm fac}$
from 1 to 50. We find that $dt_{\rm fac}=5$ provides the sufficient time
resolution to this model. This also provides the consistent result with that
obtained by the high resolution simulation of finite-difference method.

Figure 16: The $\Delta t$ dependence of $<P_{ex}>$ evolution. Different colors
denote the results with different time resolution: $dt_{\rm fac}=1,5,10$ and
50 (see text for the definition of $dt_{\rm fac}$). For comparison, the result
of the high resolution simulation of finite-difference solver is also shown in
black dashed line.

##  8 QKE with scattering process

In this section, we apply our QKE-MC solver to problems including collision
term. We pay special attention to scattering processes in this paper; there
are mainly two reasons. One of them is that the scattering process is
essentially the combined process of emission and absorption, indicating that
this test provides a good validation for these processes. The other reason is
that roles of scattering in collective neutrino oscillations is gaining
increased attention from the community (Shalgar & Tamborra, 2021a; Martin et
al., 2021; Johns, 2021). Following the previous studies, we focus on an
isoenergetic scattering in QKE. We note that the blocking factor in
isoenergetic scatterings is exactly canceled between in- and out-
scatterings101010The blocking factor needs to be evaluated appropriately when
we apply our code to emission-, absorption- and non-isoenergetic scattering
processes. To do this, the energy volume in momentum space needs to be taken
into account (to compute the density matrix), indicating that the
monochromatic approximation can not be applied in the system. In fact, the
delta function in the energy spectrum of neutrinos (corresponding to the
monochromatic approximation) violates the characteristics of Fermions..
Although we postpone more detailed investigations on QKE with collision terms,
this study has a role of the pilot study for them.

###  8.1 Test for scattering module

Before presenting a code test for a coupling system between neutrino
oscillations and scatterings, we carry out another test to see the capability
of scattering module implemented in our QKE-MC solver. As we described in Sec.
3.2, we change numerical treatments of scattering from our classical MC
solver, indicating that the capability needs to be tested. We perform a cross-
validation study by using our classical MC solver.

The detail of this test is as follows. We set anisotropic neutrino
distributions on the monochromatic neutrino energy as an initial condition.
The angular distribution of neutrinos is the same as that used for
$\bar{\nu}_{e}$ in another test of fast flavor conversion (see Eq. 55 in
Section 7.3). We follow the evolution of the system by using two codes: QKE-MC
solver and classical MC solver. We consider the homogeneous system and
incorporate isoenergetic scattering (see below for the reaction rate). We note
that the oscillation module is turned off in the simulation by the QKE-MC
solver.

In QKE-MC solver, we adopt 2,000 MC particles. In each particle, $ee$
component of $q$ is finite but others are zero. The inverse mean free path
($C$) is set as $0.05\ {\rm km^{-1}}$. In this simulation, we apply EMFP
method to reduce the statistical noise. We vary the parameter $a$ as
$10^{-2}$, $10^{-3}$ and $10^{-4}$. We note that simulations with the smaller
value of $a$ reduce more noise but become more computationally expensive
(since the frequency of scattering for each MC particle increases with
decreasing $a$, see Sec. 3.2 for the detail.). On the other hand, the memory
capacity is not pressed, that enables us to run the simulation on standard
workstations. To reduce the computational cost, the newly generated MC
particles due to scatterings on each time step are combined promptly to the
adjacent ones. Although this prescription smears out the small angular
structure of neutrinos, such a detailed structure is out of the scope of this
test. In the classical MC transport (Kato et al., 2020), on the other hand, we
use $10^{5}$ MC particles and run 200 simulations with the identical initial
condition. The resultant statistics is high enough to quantitatively capture
the time evolution of angular distributions of neutrinos.

Figure 17 shows the time evolution of the neutrino angular distribution,
illustrating that scatterings isotropize the angular distribution of
neutrinos. We find that the statistics is improved with reducing $a$. The
result with $a=10^{-4}$ is almost identical to that obtained from classical MC
solver. This is a strong evidence that the EMFP method works well without
degrading physical fidelity.

Figure 17: The time evolution of neutrino angular distributions in a test of
scattering module (Sec. 8.1). Color distinguishes the model: purple
($a=10^{-2}$ with QKE-MC), orange ($a=10^{-3}$ with QKE-MC), light blue
($a=10^{-4}$ with QKE-MC), and black (classical MC). Line-type distinguishes
the time snapshot: $t/t_{m}$=0, 0.3, 0.75, 1.5, where $t_{m}$ denotes the mean
free time ($t_{m}=6.67\times 10^{-5}$ s).

###  8.2 Fast flavor conversions with scatterings

Figure 18: The time evolution of the angular distributions for $\nu_{e}$
(solid) and $\bar{\nu}_{e}$ (dashed) for a test of fast flavor conversion with
scatterings. Different colors denote different time steps.

As a final test, we demonstrate simulations for fast flavor conversion with
scatterings. As shown in Sec. 7.3, however, our simulations did not reproduce
the published result of Shalgar & Tamborra (2021a) in cases without
scatterings111111As mentioned in Sec. 7.3, the inconsistency to the results of
Shalgar & Tamborra (2021a) was resolved by making detailed comparison between
the two results; see text for more detail.. We, hence, take another way to
validate our results. In this test, we focus on the early phase (the flavor
conversion growth exponentially), and provide clear physical interpretations
to the dynamics. For a quantitative assessment, we compare the growth rate of
flavor conversion in our QKE-MC simulation to an empirical (analytic) formula
given by Morinaga et al. (2020).

We start with describing the governing equation of the system. Following the
convention used in Sec. 7.3, we rewrite the QKE as,

$\displaystyle i\frac{\partial\rho_{a}}{\partial t}$ $\displaystyle=$
$\displaystyle\left[H^{\prime}_{\rm vac}+H_{\rm
mat}+H^{\prime}_{\nu\nu},\rho_{a}\right]$ (57)
$\displaystyle+i\int_{-1}^{1}\tilde{C}\rho_{a}^{\prime}d\cos{\theta_{\nu}^{\prime}}$
$\displaystyle-i\int_{-1}^{1}\tilde{C}\rho_{a}d\cos{\theta_{\nu}^{\prime}},$
$\displaystyle i\frac{\partial\bar{\rho}_{a}}{\partial t}$ $\displaystyle=$
$\displaystyle\left[H^{\prime}_{\rm vac}-H^{\ast}_{\rm
mat}-H^{\prime\ast}_{\nu\nu},\bar{\rho}_{a}\right]$ (58)
$\displaystyle+i\int_{-1}^{1}\tilde{C}\bar{\rho}_{a}^{\prime}d\cos{\theta_{\nu}^{\prime}}$
$\displaystyle-i\int_{-1}^{1}\tilde{C}\bar{\rho}_{a}d\cos{\theta_{\nu}^{\prime}},\
\ \ \ $

where $\tilde{C}$ represents a scattering rate. In this test, we assume that
$\tilde{C}$ is constant (no angular dependence). We note that the inverse mean
free path $C$ is 2 times larger than $\tilde{C}$ in this expression (see also
Shalgar & Tamborra, 2021a) and we adopt $C=10\ {\rm km}^{-1}$ in this
study121212This parameter corresponds to the highest rate among models in
Shalgar & Tamborra (2021a).. It should be mentioned that a large number of MC
particles is necessary compared to the case without scatterings (see Sec.
7.3). This is due to the fact that the statistical noise makes the simulation
troublesome in this test. Following the prescriptions presented in Sec. 3.2,
we reduce the noise by a smoothing method (in computations of self-interaction
potential) and EMFP method. As a reference, we adopt the number of MC
particles, the number of angular grids ($N_{\theta}$) for the smoothing, and
EMFP parameter ($a$) are 32,000, 64, and $10^{-5}$, respectively. We show that
this provides a sufficient resolution up to the time of end of our simulation
($3\times 10^{-7}$ s). The resolution dependence is also discussed in this
section.

Figure 18 shows the time evolution of the angular distributions for $\nu_{e}$
and $\bar{\nu}_{e}$. The former has an isotropic distribution initially,
whereas the latter has an anisotropic distribution (see black lines). This
figure illustrates that $\bar{\nu}_{e}$ angular distribution is isotropized
with time, and the ELN crossing point changes with time (approaching to
$\cos\theta_{\nu}=1$). On the other hand, $\nu_{e}$ distribution does not
evolve, indicating that the flavor conversion is subtle. This suggests that
the remarkable change of angular distributions of $\bar{\nu}_{e}$ is due to
scatterings.

To see the impact of scatterings on fast flavor conversion, we show the time
evolution of off-diagonal component of $\rho_{a}$ and $\bar{\rho}_{a}$ in Fig.
19. We pick up three different angular positions in the figure. This figure
illustrates that both the oscillation amplitude and the growth rate in flavor
conversion are suppressed by scatterings regardless of neutrino and anti-
neutrinos. This may be attributed to the fact that the ELN angular
distribution evolves to an isotropic distribution by scatterings (see Fig.
18), i.e., the amplitude of ELN crossing decreases with time. As discussed in
Richers et al. (2021), the reduction of the amplitude of ELN crossing may
account for the suppression. We reckon that there is another reason. As shown
in Fig. 15, the fast flavor conversion occurs vigorously around the angular
point of ELN crossing. On the other hand, the crossing point is migrated by
scatterings (see Fig. 18), that works to disperse the growth of flavor
conversions into wide angles. This would result in suppressing flavor
conversion, unless the flavor conversion grows substantially during the
stagnation of the crossing. This suggests an interesting possibility that the
effect of scatterings hinges on the cross section of scatterings, which is in
line with results in Martin et al. (2021) and discussions in Shalgar &
Tamborra (2021a).

Figure 19: The time evolution of ${\rm Im}\rho_{a,ex}+{\rm
Im}\bar{\rho}_{a,ex}$ at three different directions: $\cos\theta_{\nu}=-0.99$
(black), $0.21$ (orange), and $0.99$ (light blue), for a test of fast flavor
conversion with scatterings. Color distinguishes angular directions. Solid and
dashed lines represent the results with and without scatterings, respectively.

To strengthen the above argument, we compare the growth rate of fast flavor
conversion obtained from the simulation to an empirical estimation proposed in
our previous paper (Morinaga et al., 2020). We estimate the growth rate from
our simulation as follows. We take the peak amplitude of off-diagonal
component of $\rho_{a}$ and $\bar{\rho}_{a}$ on each propagation direction
(see, e.g., Fig. 19). We compute the slope (growth rate) by computing the
difference at adjacent peak points and divided by the difference of time, and
then take the angular average; the result is plotted in black line of Fig. 20.
It should be noted that the growth rate is normalized by the one without
scattering. Our result suggests that scatterings under the choice of the
parameters reduce the growth of flavor conversion by $\sim$70% at $t\sim
3\times 10^{-7}$ s. The orange line in the same figure, on the other hand,
portrays the growth rate obtained by an empirical formula,

$\displaystyle
G\sim\sqrt{\left(\int_{\Delta\rho_{a}>0}\Delta\rho_{a}d\cos{\theta_{\nu}}\right)\left(\int_{\Delta\rho_{a}<0}\Delta\rho_{a}d\cos{\theta_{\nu}}\right)},\
\ \ \ \ \ $ (59)

with $\Delta\rho_{a}\equiv\rho_{a,ee}-\bar{\rho}_{a,ee}$ (see Morinaga et al.
(2020)). We estimate the growth rate at each timestep by using
$\Delta\rho_{a}$ obtained from the simulation. Although the two lines in Fig.
20 are slightly different from each other, the qualitative trend is
essentially the same131313The deviation is mainly due to the accuracy of the
empirical relation; indeed, the exact growth rate is quantitatively different
from it. See Morinaga et al. (2020) for more detail.. This result suggests
that the growth of flavor conversion is determined mainly by the angular
distribution of $\nu_{e}$ and $\bar{\nu}_{e}$ at each instantaneous time, and
that the isotropized ELN angular distribution due to scatterings is the major
factor of the suppression (see Fig. 18). It should also be stressed that the
reasonable agreement between the two results suggests that our QKE-MC solver
works well to this problem.

Figure 20: The time evolution of the growth rate ratio of fast flavor
conversion between with and without scatterings. The ratio of $<1$ indicates
that the scattering suppresses fast flavor conversion. Orange line denotes the
result computed by the time evolution of ${\rm Im}\rho_{ex}+{\rm
Im}\bar{\rho}_{ex}$ from our QKE-MC simulation. The ratio computed from the
empirical laws by Morinaga et al. (2020) is shown in black line. See text for
more detail.

Finally, we perform resolution studies. In the top panel of Fig. 21, we show
the resolution dependence of $N_{\theta}$. The reference value is 64, but we
check the case with 128. It should be mentioned that the number of MC
particles is also increased twice from the reference in the case with
$N_{\theta}=128$. This figure suggests that the difference between the two
results is subtle. We also show the result without scatterings in the same
figure, and the difference to the cases with scatterings is remarkable. Hence,
the resolution of $N_{\theta}=64$ is sufficient to capture the impact of
scatterings qualitatively. In the middle panel of Fig. 21, we study the
influence of statistical noise. For convenience, we use another variable, $N$,
that is defined as the total number of MC particles divided by $N_{\theta}$,
i.e., it represents the total number of MC particles on each angular mesh. We
find that $N=500$ (reference value) is almost identical to the case with
$N=1000$, indicating that the statistical noise in our reference model does
not compromise our results. It should be mentioned that the statistical noise
is also reduced by the EMFP method. To see the impact more quantitatively, we
run two additional simulations with $a=10^{-4}$ and $10^{-6}$, and the results
are displayed in the bottom panel of Fig. 21. This figure suggests that the
simulation with $a=10^{-5}$ (reference value) is capable of providing a
physically accurate result. This convergence test strengthens the reliability
and capability of our new modules. It should also be mentioned that we made a
comparison to the group of Shalgar & Tamborra (2021a), and confirmed that both
simulations provide consistent results.

Figure 21: Resolution studies for the test of fast flavor conversion with
scatterings (Sec. 8.2). Top: the number of angular grids ($N_{\theta}$).
Middle: the number of sample particles in each grid ($N$). Bottom: EMFP
parameter ($a$). Color distinguishes different parameters for each resolution
test. Green line corresponds to a reference calculation ($N_{\theta}=64$,
$N=500$, $a=10^{-5}$). For comparison, the results without scatterings are
also plotted with black dashed lines.

##  9 Summary and discussions

Recent years have witnessed a number of theoretical indications that
collective neutrino oscillations commonly occur in CCSN and BNSM environments.
This suggests that quantum kinetic features as neutrino oscillations needs to
be taken into account, that potentially impact the global dynamics,
nucleosynthesis, and neutrino signals. The full QKE neutrino transport is
capable of describing both the microscopic and macroscopic dynamics. Numerical
simulations may be the most general ab-initio approach. Although a
considerable amount of efforts has been devoted to numerical simulations, the
physically consistent treatments of flavor conversion, transport, and matter
collision are technically difficult. In this paper, we propose a Monte Carlo
(MC) method, which is potentially capable of handling all elements self-
consistently. To this end, we updated our classical MC neutrino transport
solver (Kato et al., 2020) to that of QKE. It is a representative
probabilistic and particle method, indicating that our numerical models will
be complementary to those simulated by other deterministic and mesh-based
approaches.

In this QKE-MC solver, we handle neutrino flavor conversions by embedding a
flavor degree of freedom into each MC particle (see Sec. 3.1). With the
upgrades, numerical treatments of collision terms, in particular for
scattering processes, need to be changed (see Sec. 3.2). New techniques are
also developed to reduce statistical noise inherent in MC method. In this
paper, we present a suite of code tests (see Secs. 5 to 8) to see the
capability of each upgraded module by taking a cross-validation strategy, in
which we reproduce the results reported in the previous literature. We show
that our results are in good agreement with them for the tests of vacuum,
matter, and collective oscillations. On the other hand, we did not reproduce
the result of Shalgar & Tamborra (2021a) in tests of fast flavor conversions
in cases without scatterings. We, hence, validated our results by applying
another code (finite-difference scheme) to the same problem, and made a
quantitative comparison, that provides consistent results with QKE-MC solver.
In cases including scatterings, we validated our results by comparing an
empirical (analytic) formula, and we showed that the QKE-MC solver is capable
of providing physically accurate solutions. It should also be mentioned that
we and the authors of Shalgar & Tamborra (2021a) made an in-depth analysis on
each simulation. We finally reached a consensus that the QKE-MC solver is
capable of providing physically accurate solutions.

The present paper is the first report describing the essential philosophy,
design, and implementation of QKE neutrino transport by a MC approach. The
detailed investigations of non-linear dynamics entangling neutrino
oscillation, transport, and collision term are not the focus of this work but
the greatest motivation. One of our short-term targets is to reveal the
effects of collision terms, in particular for scatterings, on fast flavor
conversions. We note that the results of numerical simulations by different
groups are still at odds each other (see, e.g., Capozzi et al., 2019; Shalgar
& Tamborra, 2021a; Martin et al., 2021). Since our numerical method is
fundamentally different from theirs, the results will provide an unique
insight into the source of controversy. We have already started the long-term
simulations of such a system, and the results will be reported in a
forthcoming paper.

We are grateful to Sherwood Richers for valuable comments and discussions. We
also thank Irene Tamborra and Shashank Shalgar for carefully checking their
simulation to resolve the discrepancy between us. C.K. is supported by Grant-
in-Aid for Early-Career Scientists (No.20K14457) and Grant-in-Aid for
Scientific Research on Innovative Areas (No.20H05240). T.M. is supported by
JSPS Grant-in-Aid for JSPS Fellows (No. 19J21244) from the Ministry of
Education, Culture, Sports, Science and Technology (MEXT), Japan.

## References

  * Aartsen et al. (2018) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2018, Phys. Rev. Lett., 120, 071801, doi: 10.1103/PhysRevLett.120.071801
  * Abbar (2020) Abbar, S. 2020, J. Cosmology Astropart. Phys, 2020, 027, doi: 10.1088/1475-7516/2020/05/027
  * Abbar (2021) Abbar, S. 2021, Phys. Rev. D, 103, 045014, doi: 10.1103/PhysRevD.103.045014
  * Abbar et al. (2021) Abbar, S., Capozzi, F., Glas, R., Janka, H.-T., & Tamborra, I. 2021, Phys. Rev. D, 103, 063033, doi: 10.1103/PhysRevD.103.063033
  * Abbar et al. (2020) Abbar, S., Duan, H., Sumiyoshi, K., Takiwaki, T., & Volpe, M. C. 2020, Phys. Rev. D, 101, 043016, doi: 10.1103/PhysRevD.101.043016
  * Abbar & Volpe (2019) Abbar, S., & Volpe, M. C. 2019, Physics Letters B, 790, 545, doi: 10.1016/j.physletb.2019.02.002
  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, ApJ, 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Abe et al. (2014) Abe, K., et al. 2014, Nuclear Instruments and Methods in Physics Research A, 737, 253, doi: 10.1016/j.nima.2013.11.081
  * Abe et al. (2018a) Abe, K., Bronner, C., Haga, Y., et al. 2018a, Phys. Rev. D, 97, 072001, doi: 10.1103/PhysRevD.97.072001
  * Abe et al. (2018b) Abe, K., Akutsu, R., Ali, A., et al. 2018b, Phys. Rev. Lett., 121, 171802, doi: 10.1103/PhysRevLett.121.171802
  * Acero et al. (2019) Acero, M. A., Adamson, P., Aliaga, L., et al. 2019, Phys. Rev. Lett., 123, 151803, doi: 10.1103/PhysRevLett.123.151803
  * Adamson et al. (2014) Adamson, P., Anghel, I., Aurisano, A., et al. 2014, Phys. Rev. Lett., 112, 191801, doi: 10.1103/PhysRevLett.112.191801
  * Aker et al. (2019) Aker, M., Altenmüller, K., Arenz, M., et al. 2019, Phys. Rev. Lett., 123, 221802, doi: 10.1103/PhysRevLett.123.221802
  * Blaschke & Cirigliano (2016) Blaschke, D. N., & Cirigliano, V. 2016, Phys. Rev. D, 94, 033009, doi: 10.1103/PhysRevD.94.033009
  * Burrows & Vartanyan (2021) Burrows, A., & Vartanyan, D. 2021, Nature, 589, 29, doi: 10.1038/s41586-020-03059-w
  * Capozzi et al. (2020) Capozzi, F., Chakraborty, M., Chakraborty, S., & Sen, M. 2020, Phys. Rev. Lett., 125, 251801, doi: 10.1103/PhysRevLett.125.251801
  * Capozzi et al. (2019) Capozzi, F., Dasgupta, B., Mirizzi, A., Sen, M., & Sigl, G. 2019, Phys. Rev. Lett., 122, 091101, doi: 10.1103/PhysRevLett.122.091101
  * Chakraborty et al. (2010) Chakraborty, S., Choubey, S., Goswami, S., & Kar, K. 2010, J. Cosmology Astropart. Phys, 2010, 007, doi: 10.1088/1475-7516/2010/06/007
  * Chakraborty et al. (2016) Chakraborty, S., Hansen, R. S., Izaguirre, I., & Raffelt, G. G. 2016, J. Cosmology Astropart. Phys, 2016, 042, doi: 10.1088/1475-7516/2016/03/042
  * Cherry et al. (2012) Cherry, J. F., Carlson, J., Friedland, A., Fuller, G. M., & Vlasenko, A. 2012, Phys. Rev. Lett., 108, 261104, doi: 10.1103/PhysRevLett.108.261104
  * Cherry et al. (2013) Cherry, J. F., Carlson, J., Friedland, A., Fuller, G. M., & Vlasenko, A. 2013, Phys. Rev. D, 87, 085037, doi: 10.1103/PhysRevD.87.085037
  * Cherry et al. (2020) Cherry, J. F., Fuller, G. M., Horiuchi, S., et al. 2020, Phys. Rev. D, 102, 023022, doi: 10.1103/PhysRevD.102.023022
  * Dasgupta & Dighe (2008) Dasgupta, B., & Dighe, A. 2008, Phys. Rev. D, 77, 113002, doi: 10.1103/PhysRevD.77.113002
  * Dasgupta et al. (2009) Dasgupta, B., Dighe, A., Raffelt, G. G., & Smirnov, A. Y. 2009, Phys. Rev. Lett., 103, 051105, doi: 10.1103/PhysRevLett.103.051105
  * Dasgupta et al. (2017) Dasgupta, B., Mirizzi, A., & Sen, M. 2017, J. Cosmology Astropart. Phys, 2017, 019, doi: 10.1088/1475-7516/2017/02/019
  * Dasgupta et al. (2018) Dasgupta, B., Mirizzi, A., & Sen, M. 2018, Phys. Rev. D, 98, 103001, doi: 10.1103/PhysRevD.98.103001
  * Delfan Azari et al. (2019) Delfan Azari, M., Yamada, S., Morinaga, T., et al. 2019, Phys. Rev. D, 99, 103011, doi: 10.1103/PhysRevD.99.103011
  * Delfan Azari et al. (2020) —. 2020, Phys. Rev. D, 101, 023018, doi: 10.1103/PhysRevD.101.023018
  * Dighe & Smirnov (2000) Dighe, A. S., & Smirnov, A. Y. 2000, Phys. Rev. D, 62, 033007, doi: 10.1103/PhysRevD.62.033007
  * Duan et al. (2006a) Duan, H., Fuller, G. M., Carlson, J., & Qian, Y.-Z. 2006a, Phys. Rev. D, 74, 105014, doi: 10.1103/PhysRevD.74.105014
  * Duan et al. (2006b) —. 2006b, Phys. Rev. D, 74, 105014, doi: 10.1103/PhysRevD.74.105014
  * Esteban et al. (2020) Esteban, I., Gonzalez-Garcia, M. C., Maltoni, M., Schwetz, T., & Zhou, A. 2020, Journal of High Energy Physics, 2020, 178, doi: 10.1007/JHEP09(2020)178
  * Farag et al. (2020) Farag, E., Timmes, F. X., Taylor, M., Patton, K. M., & Farmer, R. 2020, ApJ, 893, 133, doi: 10.3847/1538-4357/ab7f2c
  * Foglizzo et al. (2015) Foglizzo, T., Kazeroni, R., Guilet, J., et al. 2015, PASA, 32, e009, doi: 10.1017/pasa.2015.9
  * Fuller et al. (1987) Fuller, G. M., Mayle, R. W., Wilson, J. R., & Schramm, D. N. 1987, ApJ, 322, 795, doi: 10.1086/165772
  * Gallo Rosso et al. (2018) Gallo Rosso, A., Vissani, F., & Volpe, M. C. 2018, J. Cosmology Astropart. Phys, 2018, 040, doi: 10.1088/1475-7516/2018/04/040
  * George et al. (2020) George, M., Wu, M.-R., Tamborra, I., Ardevol-Pulpillo, R., & Janka, H.-T. 2020, Phys. Rev. D, 102, 103015, doi: 10.1103/PhysRevD.102.103015
  * Glas et al. (2020) Glas, R., Janka, H.-T., Capozzi, F., et al. 2020, Phys. Rev. D, 101, 063001, doi: 10.1103/PhysRevD.101.063001
  * Hannestad et al. (2006) Hannestad, S., Raffelt, G. G., Sigl, G., & Wong, Y. Y. Y. 2006, Phys. Rev. D, 74, 105010, doi: 10.1103/PhysRevD.74.105010
  * Harris et al. (2017) Harris, J. A., Hix, W. R., Chertkow, M. A., et al. 2017, ApJ, 843, 2, doi: 10.3847/1538-4357/aa76de
  * Horiuchi & Kneller (2018) Horiuchi, S., & Kneller, J. P. 2018, Journal of Physics G Nuclear Physics, 45, 043002, doi: 10.1088/1361-6471/aaa90a
  * Janka (2012) Janka, H.-T. 2012, Annual Review of Nuclear and Particle Science, 62, 407, doi: 10.1146/annurev-nucl-102711-094901
  * Johns (2021) Johns, L. 2021, arXiv e-prints, arXiv:2104.11369.  https://arxiv.org/abs/2104.11369
  * Johns et al. (2020) Johns, L., Nagakura, H., Fuller, G. M., & Burrows, A. 2020, Phys. Rev. D, 101, 043009, doi: 10.1103/PhysRevD.101.043009
  * Kato et al. (2020) Kato, C., Nagakura, H., Hori, Y., & Yamada, S. 2020, ApJ, 897, 43, doi: 10.3847/1538-4357/ab97b2
  * Kawaguchi et al. (2021) Kawaguchi, K., Fujibayashi, S., Shibata, M., Tanaka, M., & Wanajo, S. 2021, ApJ, 913, 100, doi: 10.3847/1538-4357/abf3bc
  * Kotake et al. (2012) Kotake, K., Takiwaki, T., Suwa, Y., et al. 2012, Advances in Astronomy, 2012, 428757, doi: 10.1155/2012/428757
  * Li et al. (2021) Li, S. W., Roberts, L. F., & Beacom, J. F. 2021, Phys. Rev. D, 103, 023016, doi: 10.1103/PhysRevD.103.023016
  * Martin et al. (2021) Martin, J. D., Carlson, J., Cirigliano, V., & Duan, H. 2021, Phys. Rev. D, 103, 063001, doi: 10.1103/PhysRevD.103.063001
  * Martin et al. (2020a) Martin, J. D., Carlson, J., & Duan, H. 2020a, Phys. Rev. D, 101, 023007, doi: 10.1103/PhysRevD.101.023007
  * Martin et al. (2020b) Martin, J. D., Yi, C., & Duan, H. 2020b, Physics Letters B, 800, 135088, doi: 10.1016/j.physletb.2019.135088
  * Martínez-Pinedo et al. (2014) Martínez-Pinedo, G., Fischer, T., & Huther, L. 2014, Journal of Physics G Nuclear Physics, 41, 044008, doi: 10.1088/0954-3899/41/4/044008
  * Miller et al. (2019) Miller, J. M., Ryan, B. R., Dolence, J. C., et al. 2019, Phys. Rev. D, 100, 023008, doi: 10.1103/PhysRevD.100.023008
  * Mirizzi et al. (2016) Mirizzi, A., Tamborra, I., Janka, H. T., et al. 2016, Nuovo Cimento Rivista Serie, 39, 1, doi: 10.1393/ncr/i2016-10120-8
  * Mori et al. (2021) Mori, M., Suwa, Y., Nakazato, K., et al. 2021, Progress of Theoretical and Experimental Physics, 2021, 023E01, doi: 10.1093/ptep/ptaa185
  * Morinaga (2021) Morinaga, T. 2021, arXiv e-prints, arXiv:2103.15267.  https://arxiv.org/abs/2103.15267
  * Morinaga et al. (2020) Morinaga, T., Nagakura, H., Kato, C., & Yamada, S. 2020, Physical Review Research, 2, 012046, doi: 10.1103/PhysRevResearch.2.012046
  * Müller (2019) Müller, B. 2019, Annual Review of Nuclear and Particle Science, 69, 253, doi: 10.1146/annurev-nucl-101918-023434
  * Müller (2020) —. 2020, Living Reviews in Computational Astrophysics, 6, 3, doi: 10.1007/s41115-020-0008-5
  * Nagakura et al. (2021a) Nagakura, H., Burrows, A., & Vartanyan, D. 2021a, MNRAS, 506, 1462, doi: 10.1093/mnras/stab1785
  * Nagakura et al. (2021b) Nagakura, H., Burrows, A., Vartanyan, D., & Radice, D. 2021b, MNRAS, 500, 696, doi: 10.1093/mnras/staa2691
  * Nagakura et al. (2019) Nagakura, H., Morinaga, T., Kato, C., & Yamada, S. 2019, ApJ, 886, 139, doi: 10.3847/1538-4357/ab4cf2
  * O’Connor et al. (2018) O’Connor, E., Bollig, R., Burrows, A., et al. 2018, Journal of Physics G Nuclear Physics, 45, 104001, doi: 10.1088/1361-6471/aadeae
  * Odrzywolek & Heger (2010) Odrzywolek, A., & Heger, A. 2010, Acta Phys. Polon. B, 41, 1611
  * Olive & Particle Data Group (2014) Olive, K. A., & Particle Data Group. 2014, Chinese Physics C, 38, 090001, doi: 10.1088/1674-1137/38/9/090001
  * Radice et al. (2018) Radice, D., Perego, A., Hotokezaka, K., et al. 2018, ApJ, 869, 130, doi: 10.3847/1538-4357/aaf054
  * Richers et al. (2017) Richers, S., Nagakura, H., Ott, C. D., et al. 2017, ApJ, 847, 133, doi: 10.3847/1538-4357/aa8bb2
  * Richers et al. (2021) Richers, S., Willcox, D. E., Ford, N. M., & Myers, A. 2021, Phys. Rev. D, 103, 083013, doi: 10.1103/PhysRevD.103.083013
  * Richers et al. (2019) Richers, S. A., McLaughlin, G. C., Kneller, J. P., & Vlasenko, A. 2019, Phys. Rev. D, 99, 123014, doi: 10.1103/PhysRevD.99.123014
  * Samuel (1993) Samuel, S. 1993, Phys. Rev. D, 48, 1462, doi: 10.1103/PhysRevD.48.1462
  * Sasaki et al. (2017) Sasaki, H., Kajino, T., Takiwaki, T., et al. 2017, Phys. Rev. D, 96, 043013, doi: 10.1103/PhysRevD.96.043013
  * Sawyer (2005) Sawyer, R. F. 2005, Phys. Rev. D, 72, 045003, doi: 10.1103/PhysRevD.72.045003
  * Sekiguchi et al. (2015) Sekiguchi, Y., Kiuchi, K., Kyutoku, K., & Shibata, M. 2015, Phys. Rev. D, 91, 064059, doi: 10.1103/PhysRevD.91.064059
  * Shalgar & Tamborra (2021a) Shalgar, S., & Tamborra, I. 2021a, Phys. Rev. D, 103, 063002, doi: 10.1103/PhysRevD.103.063002
  * Shalgar & Tamborra (2021b) —. 2021b, arXiv e-prints, arXiv:2106.15622.  https://arxiv.org/abs/2106.15622
  * Shibata et al. (2014) Shibata, M., Nagakura, H., Sekiguchi, Y., & Yamada, S. 2014, Phys. Rev. D, 89, 084073, doi: 10.1103/PhysRevD.89.084073
  * Siegel & Metzger (2017) Siegel, D. M., & Metzger, B. D. 2017, Phys. Rev. Lett., 119, 231102, doi: 10.1103/PhysRevLett.119.231102
  * Sieverding et al. (2020) Sieverding, A., Müller, B., & Qian, Y. Z. 2020, ApJ, 904, 163, doi: 10.3847/1538-4357/abc61b
  * Sigl (1995) Sigl, G. 1995, Phys. Rev. D, 51, 4035, doi: 10.1103/PhysRevD.51.4035
  * Sigl & Raffelt (1993) Sigl, G., & Raffelt, G. 1993, Nuclear Physics B, 406, 423, doi: 10.1016/0550-3213(93)90175-O
  * Stapleford et al. (2020) Stapleford, C. J., Fröhlich, C., & Kneller, J. P. 2020, Phys. Rev. D, 102, 081301, doi: 10.1103/PhysRevD.102.081301
  * Suwa et al. (2019) Suwa, Y., Sumiyoshi, K., Nakazato, K., et al. 2019, ApJ, 881, 139, doi: 10.3847/1538-4357/ab2e05
  * T2K Collaboration & Abe (2014) T2K Collaboration, & Abe, K. 2014, arXiv e-prints, arXiv:1403.1532.  https://arxiv.org/abs/1403.1532
  * Tamborra et al. (2017) Tamborra, I., Hüdepohl, L., Raffelt, G. G., & Janka, H.-T. 2017, ApJ, 839, 132, doi: 10.3847/1538-4357/aa6a18
  * Thielemann et al. (2017) Thielemann, F. K., Eichler, M., Panov, I. V., & Wehmeyer, B. 2017, Annual Review of Nuclear and Particle Science, 67, 253, doi: 10.1146/annurev-nucl-101916-123246
  * Tian et al. (2017) Tian, J. Y., Patwardhan, A. V., & Fuller, G. M. 2017, Phys. Rev. D, 96, 043001, doi: 10.1103/PhysRevD.96.043001
  * Vlasenko et al. (2014) Vlasenko, A., Fuller, G. M., & Cirigliano, V. 2014, Phys. Rev. D, 89, 105004, doi: 10.1103/PhysRevD.89.105004
  * Volpe (2015) Volpe, C. 2015, International Journal of Modern Physics E, 24, 1541009, doi: 10.1142/S0218301315410098
  * Wanajo et al. (2018) Wanajo, S., Müller, B., Janka, H.-T., & Heger, A. 2018, ApJ, 852, 40, doi: 10.3847/1538-4357/aa9d97
  * Witt et al. (2021) Witt, M., Psaltis, A., Yasin, H., et al. 2021, arXiv e-prints, arXiv:2107.00687.  https://arxiv.org/abs/2107.00687
  * Wolfenstein (1978) Wolfenstein, L. 1978, Phys. Rev. D, 17, 2369, doi: 10.1103/PhysRevD.17.2369
  * Wu & Tamborra (2017) Wu, M.-R., & Tamborra, I. 2017, Phys. Rev. D, 95, 103007, doi: 10.1103/PhysRevD.95.103007
  * Wu et al. (2017) Wu, M.-R., Tamborra, I., Just, O., & Janka, H.-T. 2017, Phys. Rev. D, 96, 123015, doi: 10.1103/PhysRevD.96.123015
  * Zaizen et al. (2020) Zaizen, M., Cherry, J. F., Takiwaki, T., et al. 2020, J. Cosmology Astropart. Phys, 2020, 011, doi: 10.1088/1475-7516/2020/06/011
  * Zaizen & Morinaga (2021) Zaizen, M., & Morinaga, T. 2021, arXiv e-prints, arXiv:2104.10532.  https://arxiv.org/abs/2104.10532
  * Zaizen et al. (2018) Zaizen, M., Yoshida, T., Sumiyoshi, K., & Umeda, H. 2018, Phys. Rev. D, 98, 103020, doi: 10.1103/PhysRevD.98.103020