# Symmetry breaking due to multi-angle matter-neutrino resonance in neutron
star merger remnants

Ian Padilla-Gay ID Shashank Shalgar ID and Irene Tamborra ID

###### Abstract

Neutron star merger remnants are unique sites for exploring neutrino flavor
conversion in dense media. Because of the natural excess of $\bar{\nu}_{e}$
over $\nu_{e}$, the neutrino-neutrino potential can cancel the matter
potential, giving rise to matter-neutrino resonant flavor conversion. Under
the assumption of two (anti)neutrino flavors and spatial homogeneity, we solve
the neutrino quantum kinetic equations to investigate the occurrence of the
matter-neutrino resonance within a multi-angle framework. We find that
isotropy is broken spontaneously, regardless of the mass ordering. Relying on
a hydrodynamical simulation of a binary neutron star merger remnant with a
black hole of $3\ M_{\odot}$ and an accretion torus of $0.3\ M_{\odot}$, we
find that complete flavor conversion caused by the matter-neutrino resonance
is unlikely, although the matter and neutrino potentials cancel at various
locations above the disk. Importantly, the matter-neutrino resonant flavor
conversion crucially depends on the shape of the neutrino angular
distributions. Our findings suggest that an accurate modeling of the neutrino
angular distributions is necessary to understand flavor conversion physics in
merger remnants, its implications on the disk physics and synthesis of the
elements heavier than iron.

SLAC-PUB-17761

##  1 Introduction

Core collapse supernovae and neutron star mergers are among the richest
environments in neutrinos. In neutron star merger remnants, neutrinos are
responsible for affecting the disk cooling rate [1] as well as the
nucleosynthesis of the elements heavier than iron, especially in the polar
region [2, 3]. Neutrinos change flavor as they propagate, and the rate of
flavor transformation can be modified by the medium through which neutrinos
travel; a known example of this sort is the Mikheyev-Smirnov-Wolfenstein (MSW)
effect in which the neutrino flavor transformation is modified by the coherent
forward scattering of neutrinos with the electrons in the medium [4, 5]. In
core-collapse supernovae and neutron star mergers, the neutrino number density
is so large that neutrino flavor transformation can also be affected by the
coherent forward scattering of neutrinos off other neutrinos. This causes the
flavor of neutrinos and antineutrinos with different momenta to evolve in a
correlated manner; unlike the MSW effect, this is a nonlinear phenomenon with
a rich phenomenology [6, 7, 8, 9, 10].

While neutrino self-interaction can occur in both supernovae and neutron star
mergers, there is a crucial difference between these two environments. In
neutron star mergers, there is an excess of electron antineutrinos over
neutrinos, while the opposite is true for core-collapse supernovae. The fact
that the neutrino self-interaction strength has an opposite sign compared to
the neutrino-matter interaction one in neutron star mergers can lead to a
cancellation between the two potentials; this phenomenon is called “matter-
neutrino resonance” (MNR) [11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. Due to the
MNR, electron neutrinos completely convert into non-electron flavors, while
antineutrinos return to their initial state after the MNR transition. This
resonant behavior, however, holds in the single-angle (SA) approximation,
while only a few investigations on the MNR physics have been carried out
within the multi-angle (MA) framework [20, 19]. Yet, appropriate modeling of
the flavor outcome due to the MNR can be important to predict the
nucleosynthetic production [2], as well as contribute to place robust bounds
of non-standard physics linked to neutrinos, e.g. [21].

In this paper, we demonstrate that, in an isotropic neutrino gas, the MNR
leads to spontaneous breaking of isotropy. Unlike in the case of neutrino
self-interaction without MNR [22], this spontaneous breaking of isotropy
occurs for both neutrino mass orderings. We further explore such findings
relying on the outputs of a 2D hydrodynamical simulation of a post-merger
black hole accretion disk [23].

This work is organized as follows. In Sec. 2, we introduce the neutrino
equations of motion (EOMs). In Sec. 3, we outline the concept of MNR
transitions and highlight the key differences in the MNR occurrence between
the SA and MA frameworks. In Sec. 4, we present the post-merger BH accretion
disk system and investigate the occurrence of MNR transitions. Lastly, in Sec.
5, we present our conclusions and closing remarks. In addition, basic physics
on the MNR in the SA approximation is covered in Appendix A. Details on the
numerical convergence of our simulations are provided in Appendix B. Lastly,
in Appendix C, we reproduce well-known results for the stability conditions of
symmetry-breaking solutions in the linear regime.

##  2 Neutrino equations of motion

In this section, we introduce the neutrino EOMs as well as the MNR physics.

###  2.1 Quantum kinetic equations

For the sake of simplicity, we rely on the two-flavor approximation:
$(\nu_{e},\nu_{x})$. The neutrino field can then be modeled in terms of
$2\times 2$ density matrices, $\rho(\vec{x},\vec{p},t)$ for neutrinos and
$\bar{\rho}(\vec{x},\vec{p},t)$ for antineutrinos. The neutrino and
antineutrino density matrices are respectively,

$\displaystyle\rho(\vec{x},\vec{p},t)=\begin{pmatrix}\rho_{ee}(\vec{x},\vec{p},t)&\rho_{ex}(\vec{x},\vec{p},t)\\\
\rho_{ex}^{*}(\vec{x},\vec{p},t)&\rho_{xx}(\vec{x},\vec{p},t)\end{pmatrix}\ \
\mathrm{and}\ \
\bar{\rho}(\vec{x},\vec{p},t)=\begin{pmatrix}\bar{\rho}_{ee}(\vec{x},\vec{p},t)&\bar{\rho}_{ex}(\vec{x},\vec{p},t)\\\
\bar{\rho}_{ex}^{*}(\vec{x},\vec{p},t)&\bar{\rho}_{xx}(\vec{x},\vec{p},t)\end{pmatrix}\
,$ (2.1)

where the diagonal elements represent the occupation numbers of neutrinos of
different flavors, while the off-diagonal terms contain information on flavor
coherence. Since the luminosity of the non-electron flavors is negligible with
respect to the ones of $\nu_{e}$ and $\bar{\nu}_{e}$ in the BH torus disk [24,
25, 26], we consider $\rho_{xx}(t=0)=\bar{\rho}_{xx}(t=0)=0$ hereafter.

The EOMs for neutrinos are

$\displaystyle\kern-10.00002pt(\partial_{t}+\vec{v}\cdot\vec{\nabla})\rho(\vec{x},\vec{p},t)$
$\displaystyle=$
$\displaystyle-i[H(\vec{x},\vec{p},t),\rho(\vec{x},\vec{p},t)]+\mathcal{C}[\rho(\vec{x},\vec{p},t)]\
;$ (2.2)

the term on the left-hand side is the total derivative which includes the
advective term ($\vec{v}\cdot\vec{\nabla}$) that is neglected throughout this
paper—see Refs. [27, 28] for the role of advection in the context of flavor
conversion. The last term ($\mathcal{C}$) on the right-hand side is the
collision term, including emission, absorption, and momentum-changing
scatterings of neutrinos with matter [29, 30, 31, 32]; this term is also
neglected in the rest of this work for simplicity. Note, however, that the
advective and collision terms determine the angular distribution of neutrinos;
flavor evolution further contributes to modifying the angular distributions.
Throughout this work, we use natural units ($\hbar=c=1$).

The Hamiltonian, $H(\vec{x},\vec{p},t)$, that governs the EOMs has three
terms: the vacuum, matter, and the self-interaction terms:

$\displaystyle
H(\vec{x},\vec{p},t)=H_{\mathrm{vac}}(\vec{p})+H_{\mathrm{mat}}(\vec{x},\vec{p})+H_{\nu\nu}(\vec{x},\vec{p},t)\
.$ (2.3)

A similar equation of motion holds for antineutrinos, except for
$H_{\mathrm{vac}}\rightarrow-H_{\mathrm{vac}}$. The vacuum term is

$\displaystyle
H_{\mathrm{vac}}=\frac{\omega}{2}\begin{pmatrix}-\cos{2\theta_{V}}&\sin{2\theta_{V}}\\\
\sin{2\theta_{V}}&\cos{2\theta_{V}}\end{pmatrix}\ ,$ (2.4)

where $\omega=\Delta m^{2}/2E$ is the vacuum oscillation frequency and
$\theta_{V}$ is the vacuum mixing angle. We assume that all neutrinos have the
same energy for the sake of simplicity and set $\omega\equiv{\Delta
m^{2}}/{2E}=1\ \textrm{km}^{-1}$. Note that, in order to trigger the MNR,
sufficiently large values of the mixing angle are required (see Appendix A).
Hence, we use $\theta_{\textrm{V}}=0.15\ \textrm{rad}$, unless otherwise
specified. The matter Hamiltonian due to the coherent forward scattering with
the background matter is given by

$\displaystyle H_{\mathrm{mat}}=\lambda\begin{pmatrix}1&0\\\ 0&0\\\
\end{pmatrix}\ ,$ (2.5)

where the matter potential,

$\displaystyle\lambda=\frac{\sqrt{2}G_{\rm F}\rho_{B}}{m_{N}}Y_{e}\ ,$ (2.6)

depends on the Fermi constant $G_{\rm F}$, the baryon mass density $\rho_{B}$,
the nucleon mass $m_{N}$, and the electron fraction
$Y_{e}=(n_{e^{-}}-n_{e^{+}})/n_{B}$. The neutrino-neutrino interaction term of
the Hamiltonian takes into account the coherent forward scattering of
neutrinos with their own background and is given by

$\displaystyle H_{\nu\nu}(\vec{p},t)=\mu\int
d\vec{v}^{\prime}[\rho(v^{\prime},t)-\bar{\rho}(v^{\prime},t)][1-\vec{v}\cdot\vec{v}^{\prime}]\
,$ (2.7)

where $\mu$ is the neutrino-neutrino interaction strength. The $H_{\nu\nu}$
term couples neutrinos of different momenta and is responsible for the non-
linear nature of neutrino flavor conversion. The case with
$\vec{v}\cdot\vec{v}^{\prime}=0$ corresponds to the SA approximation (Eq. 2.7
reduces to $H^{\mathrm{SA}}_{\nu\nu}(t)=2\mu[\rho(t)-\bar{\rho}(t)]$); notice
that one obtains the same reduced system from Eq. 2.7 choosing isotropic
$\rho(v,t)$ and $\bar{\rho}(v,t)$ so that the integral of the term
proportional to $\vec{v}\cdot\vec{v}^{\prime}$ vanishes [33, 34]. The term
$(1-\vec{v}\cdot\vec{v}^{\prime})$ in Eq. 2.7 depends on both polar and
azimuthal angles, $\theta$ and $\phi$, respectively. However, for the sake of
simplicity, we ignore the dependence on $\phi$, assuming axial symmetry. The
consequences of this simplification are not known for the case of MNR,
however, they might be non-trivial as found in the context of fast flavor
evolution [35].

###  2.2 Matter-neutrino resonance

The MNR is the enhancement of flavor conversion due to the cancellation in the
diagonal terms of the Hamiltonian:

$\displaystyle
H_{ee}^{\nu\nu}(\vec{x},\vec{p})-H_{xx}^{\nu\nu}(\vec{x},\vec{p})+\lambda(\vec{x})\approx
0\ .$ (2.8)

Since the self-interaction Hamiltonian ($H_{\nu\nu}$) is angle-dependent, this
condition cannot be satisfied for all angles simultaneously in the MA case,
unless the angular distributions are isotropic. When the SA approximation is
relaxed, the MNR can occur at any points in space and momentum modes for which
Eq. 2.8 is satisfied. As time evolves, the advection of neutrinos (both
forward and backward according to the direction of propagation) ensures that
the flavor-evolved neutrinos affect the flavor composition at other spatial
locations.

The full solution of the EOMs, in the presence of the MNR and as a boundary
problem involving collisions and advection, is beyond the reach of present-day
computational facilities. In order to gain an insight into the difference
between neutrino flavor evolution due to the MNR in the SA and MA cases, in
Sec. 3, we rely on a simplified system with angular distributions fixed a
priori and with the strength of the self-interaction potential and the matter
potential changing as functions of time, while ignoring the spatial dependence
of the latter.

##  3 Matter-neutrino resonance in an ideal neutrino gas

Before investigating the MNR in the context of neutron star merger remnants,
we focus on an idealized neutrino gas. Hereafter, we assume that the spatial
coordinate is directly connected to the time coordinate through the speed of
light $c$ since neutrino advection and non-forward collisions are neglected.
We consider a time-dependent self-interaction potential and constant matter
potential:

$\displaystyle\mu(t)=3\times 10^{3}e^{-qt}\ \textrm{km}^{-1}\ \ \mathrm{and}\
\ \lambda=10^{3}\ \textrm{km}^{-1}\ ,$ (3.1)

where $q=0.1$ km-1 (or equivalently $3\times 10^{4}$ s-1). We stress that
$\mu$ is proportional to the neutrino number density; however, in this
section, we choose to define the self-interaction strength as an arbitrary
function of time (this simplification is relaxed in Sec. 4).

For the MNR to occur, Eq. 2.8 should be satisfied, which is only possible if
there are more antineutrinos than neutrinos. Hence, we assume that
$\rho_{ee}(t_{0})=1$ and $\bar{\rho}_{ee}(t_{0})=4/3$ at the initial time
$t_{0}$, and the initial asymmetry parameter is

$\displaystyle\alpha\equiv\frac{\bar{\rho}_{ee}(t_{0})}{\rho_{ee}(t_{0})}=\frac{4}{3}\
.$ (3.2)

With these initial conditions, the neutrino-neutrino Hamiltonian in the SA
approximation is $H^{\mathrm{SA}}_{\nu\nu}(t_{0})=2\mu[1-\alpha]$. The same
Hamiltonian can be obtained for the MA system with initially isotropic angular
distributions. Note that, in the MA case, we solve the EOMs fixing the angular
distributions at $t_{0}$, the neutrino distributions then evolve because of
flavor conversion (but not due to advection and collisions, both neglected in
this work).

###  3.1 Single-angle and multi-angle solutions of the matter-neutrino
resonance for an isotropic neutrino gas

Most of the MNR investigations in the literature have been performed using the
SA approximation [12, 16, 14, 13, 15]. We briefly review the MNR results in
the SA approximation in order to compare them with the ones obtained in the MA
approach. The top panels of Fig. 1 show the oscillated total potential
$H_{ee}^{\nu\nu}-H_{xx}^{\nu\nu}+\lambda$, the evolution of the off-diagonal
terms of the density matrices, and the survival probabilities of $\nu_{e}$ and
$\bar{\nu}_{e}$, from left to right respectively. As we can see, a general
feature of the MNR in the SA approximation is that $\nu_{e}$’s completely
change their flavor while $\bar{\nu}_{e}$’s return to their original
configuration. Unlike the MSW effect, the MNR is independent of the sign of
$\omega$; however, $\theta_{V}$ should be large enough for the MNR to be
possible, as discussed in Appendix A (see also Fig. 14). The times $t_{1}$ and
$t_{\mathrm{end}}$ (see dotted vertical lines in the top panels of Fig. 1)
mark the times at which the cancellation of the matter and neutrino-neutrino
potentials takes place. During $\Delta t=t_{\mathrm{end}}-t_{1}\simeq
6.5\times 10^{-5}$ s, $\rho(t)$ and $\bar{\rho}(t)$ evolve holding the MNR
conditions. The analytical dependence of $\Delta t=t_{\mathrm{end}}-t_{1}$ on
the properties of the neutrino gas is provided in Ref. [12]. A direct
comparison between the analytical estimates and our numerical results is shown
in Appendix A (cf. Fig. 12), where one can see that the agreement between both
approaches is excellent.

Figure 1: Flavor evolution for an ideal neutrino gas. Single angle (top),
multi-angle with isotropic angular distributions (middle), and multi-angle
with isotropic angular distributions and a perturbation of
$10^{-22}\cos{\theta}$ in $\bar{\rho}_{ee}$ (bottom). Left: Oscillated total
potential (sum of self-interaction and matter potentials) as the MNR
transition occurs. Middle: Off-diagonal terms of the density matrices,
$\rho_{ex}$ and $\bar{\rho}_{ex}$. Right: Survival probabilities for electron
neutrinos and antineutrinos. For the SA case (top panels), the duration of the
transition $\Delta t$ can be computed analytically (see Eq. A.1). The MA
isotropic case without perturbation reproduces the SA approximation; $500$
bits per number are needed to achieve this agreement. For the MA isotropic
case with perturbation (bottom panels), the duration of the MNR transition,
$\Delta t=t_{3}-t_{1}$, is smaller than in the SA approximation. Due to the
growth of non-isotropic modes, the MNR transition is interrupted in the MA
case around $t_{3}=0.33\times 10^{-4}$ s and flavor conversion is not complete
for electron neutrinos.

The simplest comparison between the MA and SA flavor outcomes can be carried
out considering a gas of neutrinos that is initially isotropic. As in the SA
approximation, the MA Hamiltonian does not have any dependence on $\vec{v}$
initially and naively one would not expect a major difference in the flavor
outcome with respect to the SA case. However, there is a stark difference in
the flavor evolution between the SA and MA cases. This can be seen by
comparing the three cases shown in the three rows of Fig. 1. For each row, the
panels show the evolution of the oscillated potential, the evolution of the
off-diagonal components of the density matrices, and the survival
probabilities for neutrinos and antineutrinos from left to right. The plots in
the top row refer to the SA case. The same angle-averaged quantities are shown
for an isotropic system in the middle row. The evolution of all the quantities
is identical in the first and the middle rows as expected. However, we see a
drastic difference in the evolution in the bottom row compared to the middle
row even though the only difference between the two systems is that for the
system in the bottom row a perturbation of $10^{-22}\cos\theta$ is added to
the antineutrino angular distribution. We note that the size of the
perturbation is smaller than the error of a double precision number on
computers. Indeed, we use $500$ bits of precision per number and an absolute
and relative tolerance of $10^{-25}$ each. If double precision floating point
numbers are used, then the isotropy of the system is broken by the floating
point error and it is not possible to obtain the results presented in the
middle row of Fig. 1 111We rely on the support for arbitrary precision of
numbers in the DifferentialEquations.jl package of Julia [36, 37]. In the
absence of these precautions, the results depend on the solver used to solve
the differential equations.. Hence, the MA case undergoes a spontaneous
breaking of isotropy even if a tiny anisotropy is initially present.

The spontaneous breaking of isotropy has been known to exist in collective
neutrino flavor conversion [38, 22] for either mass ordering depending on
whether there is an excess of neutrinos over antineutrinos or vice versa.
However, in the MNR case, the spontaneous breaking of isotropy is present
irrespective of the mass ordering because the underlying MNR transition that
triggers the symmetry breaking is almost identical in both mass orderings, as
shown in Appendix A (Fig. 15). Moreover, this spontaneous breaking of isotropy
manifests itself only if the isotropic mode is unstable.

If the isotropy is spontaneously broken during the MNR, this is reflected in
the time evolution of the survival probability of neutrinos,

$\displaystyle P(\nu_{e}\rightarrow\nu_{e})=\frac{\int
dv[\rho_{ee}(v,t)-\rho_{xx}(v,t_{0})]}{\int
dv[\rho_{ee}(v,t_{0})-\rho_{xx}(v,t_{0})]}\ $ (3.3)

(a similar definition holds for $P(\bar{\nu}_{e}\rightarrow\bar{\nu}_{e})$
with the replacement $\rho_{ii}\rightarrow\bar{\rho}_{ii}$). In the bottom
right panel of Fig. 1, one can see that the survival probabilities are
affected by the anisotropic modes developing between $t_{1}$ and $t_{3}$. Note
that the evolution of $P(\nu_{e}\rightarrow\nu_{e})$ for the SA and MA cases
are substantially different; this is especially true for
$P(\bar{\nu}_{e}\rightarrow\bar{\nu}_{e})$, which does not describe full
flavor conversion in the MA case (cf. right bottom panel of Fig. 1).

The sudden jump in the evolution of the off-diagonal terms of the density
matrix in the middle bottom panel of Fig. 1 can be better understood by
looking at the time snapshots of the perturbed $\nu_{e}$ and $\bar{\nu}_{e}$
angular distributions in the bottom panel of Fig. 2, before and after the MNR.
The system enters the MNR at about $t_{1}$ (see Fig. 1), but the MNR is
interrupted at $t_{3}$. The neutrino and antineutrino angular distributions no
longer preserve the isotropy that was initially imposed on them. A clear
particle excess (deficit) in the backward (forward) direction dynamically
forms during the MNR in the MA solution. More precisely, while the MNR
condition holds, the $v$-dependence of the density matrices in Fig. 2 slowly
transitions from isotropic (orange and red curves) to linear in $v$ (green
curve), and later to more complex functions non-linear in $v$ (blue curve).
The $t_{3}$ and $t_{4}$ snapshots show the angular distributions in the non-
linear regime–small-scale angular structures develop within a few $10^{-6}$ s.
This diffusion from large to small angular scales is triggered by the MNR.
Animations of the temporal evolution of the potentials, survival probabilities
and (anti)neutrino angular distributions are provided as Supplemental
Material. In Appendix B we show that the breaking of isotropy, and therefore
the cascade of flavor waves to small scales, still occurs if the number of
angular bins is increased, confirming that this is a MA effect caused by tiny
anisotropies of $\mathcal{O}(10^{-22}\cos{\theta})$ and not an artifact caused
by insufficient angular resolution.

Figure 2: Snapshots of the angular distributions of $\rho_{ee}$ (left panels)
and $\bar{\rho}_{ee}$ (right panels) for the MA isotropic case without
perturbation (top panels) and MA isotropic case with perturbation (bottom
panels) shown in Fig. 1. For the MA case without perturbation (top), the
system preserves the isotropy in the neutrino angular distributions during the
MNR. For the MA case with perturbation (bottom), the angular distributions
start to evolve as the MNR condition is met. At $t_{3}$, a dipole (green
curve) forms dynamically breaking isotropy as shown in the inset; at the same
time, a “jump” in $\rho_{ex}$ occurs as visible from the bottom left panel of
Fig. 1. From $t_{1}$ to $t_{4}$, diffusion of flavor waves from large to small
angular scales takes place.  Figure 3: Temporal evolution of the Legendre
coefficients of $\rho_{ij}$ for the system in the bottom panels of Fig. 1.
Left:  Time evolution of the norm of $\rho_{ex}^{0}$, $\rho_{xx}^{0}$,
$\rho_{ex}^{1}$ and $\rho_{xx}^{1}$. During $\Delta t=t_{3}-t_{1}$, the
potentials cancel each other (see bottom panel in Fig. 1), and the system
transitions from isotropic to non-isotropic (see bottom panel in Fig. 2) as a
result of the growing dipole. Interestingly, during $\Delta t$, both
$\rho_{ex}^{1}$ and $\rho_{xx}^{1}$ grow exponentially at the same rate.
Right:  Time evolution of the first six $\ell$ modes. Modes with $\ell>1$ also
grow within the time interval $\Delta t=t_{3}-t_{1}$ during the MNR. Larger
$\ell$ modes are fully developed around $t_{3}$, when the evolution is highly
non-linear. The discontinuity around $t_{2}$ is caused by the post-process
($64$ bits) of the numerical solution ($500$ bits) which has a larger
floating-point error.

In order to understand the breaking of isotropy, it is useful to decompose the
(anti)neutrino angular distributions in terms of Legendre polynomials:

$\displaystyle\rho_{ij}=\sum_{\ell=0}^{\infty}\rho_{ij}^{\ell}P_{\ell}\ ,$
(3.4)

where $\rho_{ij}^{\ell}$ are the coefficients of the expansion that can be
determined by integrating the product of the density matrix elements and the
corresponding Legendre polynomial, i.e.

$\displaystyle\rho_{ij}^{\ell}=\frac{2\ell+1}{2}\int_{-1}^{1}dv\rho_{ij}(v)P_{\ell}(v)\
.$ (3.5)

Thus $\rho_{ij}(v,t)$ allows us to compute the coefficients
$\rho_{ij}^{\ell=0,...,5}(t)$. These coefficients encode information about the
structure of the angular distributions and help to gauge the level of
anisotropy of the neutrino gas. In the following, we refer to the $\ell=0$ and
$\ell=1$ modes as the “monopole” and the “dipole” of the angular
distributions, respectively.

Figure 3 shows the time evolution of the monopole ($\ell=0$) and dipole
($\ell=1$) of $\rho_{ex}$ in the left panel, as well as the evolution of
higher $\ell>1$ modes in the right panel. We can see that $\rho_{ex}^{1}$ and
$\rho_{xx}^{1}$ (similar behavior for antineutrinos, although not shown) start
at $\sim 10^{-18}$ and grow exponentially during $\Delta t=t_{3}-t_{1}$ (note
that it is numerically impossible to have $\rho_{ex}^{1}$ exactly equal to
zero due to the floating-point error—see also Appendix B). The discontinuity
in Fig. 3 where $|\rho_{ex}^{1}|$ jumps from $10^{-18}$ to $10^{-12}$ is
because, in the post-processing of the solution of the EOMs, we implement
double precision ($64$ bits per number) to compute the Legendre coefficients,
while the numerical solution is obtained with $500$ bits per number; the
exponential growth would also be visible from $t_{1}$ and onwards, if the
post-processing had the same level of floating-point error accuracy of the
numerical EOM solution. During $\Delta t$, a partial MNR transition occurs
(cf. the time interval between $t_{1}$ and $t_{3}$ in Fig. 1). As the dipole
modes grow during $\Delta t$, the angular distributions become non-isotropic,
as shown in Fig. 2 at $t_{1}=0.23\times 10^{-4}$ s and $t_{3}=0.33\times
10^{-4}$ s. Furthermore, larger $\ell>1$ modes also grow exponentially early
on around $t_{1}$, as displayed in the right panel of Fig. 3. Large $\ell$
modes are not visible in the angular distributions until later times
($t_{3-4}$) because they are initially negligible. In the non-linear regime,
not only the $\ell=1$ modes have fully developed but also higher ($\ell>1$)
modes. This is a crucial difference with respect to the SA approximation: the
coherence between different angular modes is not maintained in the MA case due
to the growth of high $\ell$ multipoles during the MNR.

###  3.2 Multi-angle solution of the matter-neutrino resonance for a non-
isotropic neutrino gas

The MNR cancellation found in the MA isotropic case in Sec. 3.1 is possible
due to all angle bins initially acting as one since they share the same
$H_{\nu\nu}(v)$ as the MNR begins. We now investigate whether (partial) MNR
transitions are likely to occur in a non-isotropic neutrino gas. To this
purpose, we consider initial angular distributions that are non-isotropic and
parametrize them as follows

$\displaystyle\rho_{ee}(v)=1+av\ \ \mathrm{and}\ \
\bar{\rho}_{ee}(v)=\alpha(1+\bar{a}v)\ ,$ (3.6)

with $\alpha=\int dv\bar{\rho}_{ee}/\int dv\rho_{ee}=4/3$ (as for the SA and
MA cases in the previous section, see Eq. 3.2). The parameters $a$ and
$\bar{a}$ provide the initial degree of anisotropy and are arbitrarily fixed
to $a=0.1$ and $\bar{a}=0.2$ (note that these distributions do not have
electron lepton number crossings [8]). We define $\mu(t)$, $\lambda$,
$\omega$, and $\theta_{V}$ as in Eq. 3.1 to facilitate a comparison with our
previous SA and MA results.

A crucial difference between the SA and isotropic MA calculations is that the
Hamiltonian depends on $(1-vv^{\prime})$ in the latter case and couples
different angular modes with each other even before the MNR cancellation
occurs. As a result, different angle modes have different
$H_{ee}^{\nu\nu}(v)-H_{xx}^{\nu\nu}(v)$, which result in different locations
for the MNR as a function of $v$; all possible
$H_{ee}^{\nu\nu}(v)-H_{xx}^{\nu\nu}(v)+\lambda$ for any $v$ lie in between the
corresponding values for $v=-1$ and $v=1$.

Figure 4: Multi-angle solution for an ideal neutrino gas with non-isotropic
angular distributions, cf. Fig. 1 for comparison. Left: Oscillated total
potentials as functions of time. We mark a few selected times, $t_{1-5}$, to
highlight when the MNR condition is met for different angular modes. The gray-
shaded area marks the time interval during which the MNR occurs. Middle: Off-
diagonal terms of the density matrices of neutrinos and antineutrinos. Right:
Survival probabilities for $\nu_{e}$’s and $\bar{\nu}_{e}$’s. Although the
total potential cancels out for all angular modes at different times, none of
them successfully undergoes full MNR transition. The oscillated potentials at
$v=\pm 1$ (black and blue) do not remain zero (dashed line), which is a
characteristic feature of MNR transitions. Snapshots of the (anti)neutrino
angular distributions at $t_{1-5}$ are shown in Fig. 5.

Figure 4 shows the oscillated potentials for the angular bins $v=\pm 1$ (left
panel; the envelope of all potentials obtainable for any $v$ is enclosed
between the black and blue lines), the evolution of the off-diagonal entries
of the density matrices (middle panel), and the survival probabilities for
$\nu_{e}$ and $\bar{\nu}_{e}$ (right panel). The total potential for the $v=1$
bin (black line) fulfills the MNR condition at $t_{1}$. At $t_{3}$, the angle
bin with $v=-1$ (blue line) fulfills the MNR condition. To guide the eye, the
gray-shaded area highlights the region within which any $v$ between $v=-1$ and
$v=1$ meets the MNR condition. Most of the flavor dynamics takes place within
the gray-shaded area (which highlights the time interval between $t_{1}$ and
$t_{3}$); this is clearly visible from the evolution of the survival
probabilities in the right panel of Fig. 4. For $t>t_{3}$, the flavor dynamics
seems to come to a halt, and an approximate steady state is reached. These
results are qualitatively similar to those reported in Figs. 4 and 5 of Ref.
[20].

Figure 5: Snapshots of the angular distributions of $\nu_{e}$ (left panel) and
$\bar{\nu}_{e}$ (right panel) at the selected times $t_{1-5}$ for the MA
solution for the non-isotropic neutrino gas in Fig. 4. Due to the angular
distributions of neutrinos being non-isotropic initially (dashed), different
angular modes meet the MNR condition at different times. At first, only
forward modes ($v\sim 1$) evolve as shown by the orange curve in the inset. At
later times the instability spreads across all other angle bins as seen in the
green and blue curves.

Figure 5 shows snapshots of the angular distributions of $\nu_{e}$ and
$\bar{\nu}_{e}$ at the time snapshots $t_{1-5}$. The initial non-isotropic
distributions are plotted with black dashed lines. At $t_{1}$, when the bin
$v=1$ meets the MNR condition (see left panel of Fig. 4), forward and backward
angular bins start behaving differently. Notice that, between $t_{1}$ and
$t_{2}$, only forward $v$ modes evolve (see orange and pink curves), and
later, between $t_{2}$ and $t_{3}$, both forward and backward angular modes
evolve. After all the angular bins have undergone MNR, the system reaches a
steady state at $t_{3}$. At $t\gtrsim t_{4}$ the system is highly non-linear
and fine-grained angular structure arises.

Figure 6:  Temporal evolution of the Legendre coefficients of $\rho_{ij}$ for
the non-isotropic, MA neutrino gas in Figs. 4 and 5. Left:  Time evolution of
the norm of $\rho_{ex}^{0},\rho_{xx}^{0},\rho_{ex}^{1}$ and $\rho_{xx}^{1}$.
Unlike the isotropic MA system in Fig. 3, the components $\rho_{ex}^{1}$
(orange) and $\rho_{xx}^{1}$ (dashed yellow) are of $\mathcal{O}(10^{-3})$ and
$\mathcal{O}(10^{-5})$, respectively, at $t_{0}=0$. At $t_{1}$, when the first
angular bin meets the MNR condition (see Fig. 4), both $\ell=0$ and $\ell=1$
modes are comparable to each other in magnitude and thus already coupled to
each other in the EOMs. Right:  Temporal evolution of the first six
$\ell$-modes. Modes with $\ell>1$ are of the same order of magnitude already
when the system enters the MNR at $t_{1}$. This is a crucial difference
compared to the isotropic system, whose large $\ell$ modes are zero and grow
during the MNR dynamically.

Figure 6 displays the time evolution of the monopole, dipole, and high-$\ell$
modes for the neutrino density matrices. At $t_{1}$, the total potential of
the $v=1$ bin (black curve in the left panel of Fig. 4) crosses the zero line,
i.e. meets the MNR condition. At $t_{1}$, all $\ell$ modes are comparable to
each other in magnitude, as visible from the right panel of Fig. 6. A crucial
difference with respect to the isotropic case (Fig. 3) is that larger $\ell>0$
modes are not negligible to begin with. As a result, high $\ell$ modes couple
to each other, preventing the oscillated potential
$H_{ee}^{\nu\nu}(v)-H_{xx}^{\nu\nu}(v)+\lambda$ from staying in resonance.
Between $t_{1}$ and $t_{3}$, the trend continues more or less unchanged.
During $\Delta t=t_{3}-t_{1}$, all angular bins meet the MNR condition. After
$t_{3}$, the last bin $v=-1$ (blue curve in left panel Fig. 4) undergoes MNR,
and the system reaches a steady state. For $t\gtrsim t_{4}$, the $\ell>1$
modes are equal to $\sim 10^{-2}$ until the end of the evolution at $10^{-4}$
s. This is reflected in the overall shape of the angular distributions at the
$t_{3}$, $t_{4}$ and $t_{5}$, as visible from the red, green and blue curves
in Fig. 5 which overlap with each other.

##  4 Matter-neutrino resonance in neutron star merger remnants

In order to appreciate the differences in the MNR flavor conversion physics
between the SA and MA cases, we have used a simple idealized neutrino gas in
the previous sections. In this section, we apply our findings to the BH
accretion disk resulting from the binary neutron star post-merger phase.

Figure 7: Characteristic properties of the BH torus remnant at $15$ ms in the
$x$–$z$ plane. Top:  Product of the baryon density $\rho_{B}$ and electron
fraction $Y_{e}$ (i.e., electron density, left) and the ratio of the chemical
potentials and temperature (right). Non-degenerate regions roughly correspond
to locations where $\mu_{e}/T\lesssim 1$. Bottom: Electron neutrino (left) and
antineutrino (right) flux factors as defined in Eq. 4.1. The solid lines
represent the locations where the $j_{\nu_{e}}=1/3$ (red) and
$j_{\bar{\nu}_{e}}=1/3$ (black)—inside these regions, $\nu_{e}$’s and
$\bar{\nu}_{e}$’s are trapped and outside these regions, they are
approximately in the free-streaming regime.

###  4.1 Neutron star merger remnant model

We adopt the 2D hydrodynamical simulation output carried out in Ref. [39]. In
particular, we consider the BH accretion torus model M3A8m3a5. This simulation
employs neutrino moment transport and subgrid viscosity. An idealized
equilibrium torus of $0.3\ M_{\odot}$ is simulated in the surroundings of a
central BH of $3\ M_{\odot}$, with a dimensionless BH spin parameter of $0.8$.
The newly born accretion torus starts to lose mass while accreting on the BH.
Neutrino cooling is not efficient in the first $\mathcal{O}(10)$ ms and the
environment is optically thick. Then, as the density drops, neutrino cooling
starts to balance viscous heating. The neutrino production rate decreases as
the mass and density of the torus decrease until neutrino cooling becomes
inefficient and the torus enters a phase dominated by advection, with viscous
heating driving the expansion of the torus and launching the outflows.

For simplicity and illustrative purposes, we focus on the disk properties at
$15$ ms after the merger. Figure 7 shows the following physical quantities in
the $x$–$z$ plane: the product between the baryon mass density $\rho_{B}$ and
the electron fraction $Y_{e}$ (i.e., the electron density), the ratio of the
electron potential $\mu_{e}$ and the electron temperature $T$, and the flux
factors of neutrinos and antineutrinos, $j_{\nu_{e}}$ and $j_{\bar{\nu}_{e}}$.
We can see that the region in the proximity of the polar axis is less dense
with respect to the disk one. Moreover, $\mu_{e}/T\lesssim 1$ in a large
portion of the polar region.

At a given location above the disk, the energy-integrated neutrino number flux
vector $({\vec{F}_{\nu_{\alpha}}})$ and the neutrino number density
$(n_{\nu_{\alpha}})$ are extracted from the hydrodynamical simulation (model
M3A8m3a5) to compute the location of the electron neutrino and antineutrino
decoupling surfaces. To this purpose, we require that the flux factors in the
absence of flavor conversion are [40]

$\displaystyle
j_{\nu_{e},\bar{\nu}_{e}}\equiv\frac{|{\vec{F}_{\nu_{e},\bar{\nu}_{e}}}|}{n_{\nu_{e},\bar{\nu}_{e}}}=\frac{1}{3}\
.$ (4.1)

If the flux factor is small, then $\nu_{e}$’s or $\bar{\nu}_{e}$’s are
expected to be strongly coupled to matter. On the contrary, for large values
of the flux factor, i.e., $\mathcal{O}(j_{\nu_{e},\bar{\nu}_{e}})\sim 1$,
$\nu_{e}$’s or $\bar{\nu}_{e}$’s are approximately in the free streaming
regime. The red and black solid lines in Fig. 7 are the locations where the
flux factors equal $1/3$. Note that, while (anti)neutrinos are distributed
according to an energy spectrum in model M3A8m3a5, for simplicity we consider
mono-energetic neutrinos as outlined in Sec. 2.

###  4.2 Modeling of the neutrino angular distributions

The BH torus simulation does not provide the angular distributions of
$\nu_{e}$’s and $\bar{\nu}_{e}$’s. Hence, we set them arbitrarily at the
location above the disk where we start to evolve the EOM. We adopt two
parametrizations for the angular distributions.

In the case of isotropic angular distributions, we assume

$\displaystyle\rho_{ee}(v)=\frac{1}{2}\ \mathrm{and}\
\bar{\rho}_{ee}(v)=\frac{\alpha}{2}\ ,$ (4.2)

where $\alpha=n_{\bar{\nu}_{e}}/n_{\nu_{e}}$ is extracted from our benchmark
BH torus model.

Figure 8: Left: Initial angular distributions (parametrized as in Eq. 4.3) of
$\nu_{e}$ (solid red) and $\bar{\nu}_{e}$ (solid blue) and their difference
(solid green) at $x_{0}=0.3\times 10^{6}$ cm, $z_{0}=0.5\times 10^{7}$ cm (cf.
star marker in Fig. 9). The flux factors at this location are
$j_{\nu_{e}}=0.853$ and $j_{\bar{\nu}_{e}}=0.855$, while the number densities
are such that $\alpha\approx 1.11$. For comparison, the dotted lines represent
the isotropic distributions (Eq. 4.2) for $\nu_{e}$’s (red) and
$\bar{\nu}_{e}$’s (blue). Right:  Same as the left panel, but at
$x_{0}=0.4\times 10^{6}$ cm, $z_{0}=0.2\times 10^{7}$ cm (diamond marker in
Fig. 9); since this location is closer to the central compact object, the flux
factors are smaller ($j_{\nu_{e}}=0.583$ and $j_{\bar{\nu}_{e}}=0.588$). As a
result, the angular distributions are less forward-peaked than the ones in the
left panel. The number densities are such that $\alpha\approx 1.10$. Each
marker (star, diamond and square) shows $v_{*}$ (for the cases shown in Fig.
9) which corresponds to the angular bin where (anti)neutrinos start to
plateau. In Table 1 we summarize the parameters used for each case.  Figure 9:
Isocontours of the total potential [$\sqrt{2}G_{\rm
F}(n_{\nu_{e}}-n_{\bar{\nu}_{e}})+\lambda$] in the $x$–$z$ plane in the SA (or
angle-integrated MA isotropic) case. The red and black solid lines represent
the decoupling surfaces of neutrinos and antineutrinos, respectively. The red
dashed line marks the loci where the neutrino potential $\sqrt{2}G_{\rm
F}(n_{\nu_{e}}-n_{\bar{\nu}_{e}})$ and the matter potential $\lambda$ cancel
each other (Eq. 2.8). The black dotted line represents a selected trajectory
with $z=0.5\times 10^{7}$ cm (trajectory 1), and the black dash-dotted line is
another trajectory with $z=0.2\times 10^{7}$ cm (trajectory 2) along which we
solve the EOMs. The square, star and diamond markers show the loci at which
the angular distributions are modeled.

For the scenario with anisotropic angular distributions, we parametrize the
latter with two independent variables, which can be uniquely determined given
the first two neutrino moments obtained from the simulation data [41, 42, 43].
We use a sigmoid function:

$\displaystyle\rho_{ee}(v,v_{c})=C\frac{1}{1+e^{-\beta(v-v_{c})}}\ \
\mathrm{and}\ \
\bar{\rho}_{ee}(v,\bar{v}_{c})=\bar{C}\frac{\alpha}{1+e^{-\beta(v-\bar{v}_{c})}}\
,$ (4.3)

where the constants $C$ and $\bar{C}$ ensure that $\int dv\rho_{ee}(v)=1$ and
$\int dv\bar{\rho}_{ee}(v)=\alpha$, respectively. The parameters $v_{c}$ and
$\bar{v}_{c}$ shift the distributions towards left or right in $v$ (the closer
$v_{c}$ and $\bar{v}_{c}$ are to $1$ the more forward-peaked the distributions
are), while the parameter $\beta=30$ in the exponent softens the edges of the
angular distributions; the larger $\beta$ the closer the angular distributions
resemble top-hat distributions. Varying the value of $v_{c}$ ($\bar{v}_{c}$),
we can match any neutrino (antineutrino) flux factor (Eq. 4.1) between $0$ and
$1$; notice that since for example the flux factor for $\nu_{e}$ is the ratio
between $\int dvv\rho_{ee}$ and $\int dv\rho_{ee}$ (similar for
antineutrinos), the constants $C$, $\bar{C}$ and $\alpha$ drop out and the
only tunable parameter left is $v_{c}$ ($\bar{v}_{c}$). In the limit
$v_{c}\approx 1$, the flux factor is maximum ($j_{\nu_{e}}\approx 1$), while
for $v_{c}\approx-1$ the flux factor is minimum ($j_{\nu_{e}}\approx 0$).
Intuitively, these two extremes correspond to the cases of a Delta
distribution function in the forward direction and a completely isotropic
distribution, respectively.

We define $v_{*}$ as the first angular bin at which the (anti)neutrino angular
distributions reach a plateau. Note that the value of $v_{*}$ is approximately
the same for both neutrinos and antineutrinos because the distributions peak
in the same angular region. This is a result of the (anti)neutrino flux
factors being close to each other in magnitude (see bottom panels of Fig. 7).
If the angular distributions are isotropic, we consider $v_{*}=0$.

Figure 8 shows the isotropic distributions for the square (dotted lines in the
left panel), the non-isotropic angular distributions in Eq. 4.3 for the star
(solid lines in the left panel) and diamond (solid lines in the right panel)
marker locations indicated in Fig. 9. For the non-isotropic distributions in
the left panel of Fig. 8 (trajectory 1), we use $v_{c}=0.717$ and
$\bar{v}_{c}=0.723$ given that the correspondent flux factors are
$j_{\nu_{e}}=0.853$ and $j_{\bar{\nu}_{e}}=0.855$ (see bottom panels of Fig.
7). On the other hand, in the right panel of Fig. 8, we show the angular
distributions extracted in correspondence of the diamond marker in Fig. 9
(trajectory 2), with $v_{c}=0.17$ and $\bar{v}_{c}=0.18$; the flux factors at
this location are $j_{\nu_{e}}=0.583$ and $j_{\bar{\nu}_{e}}=0.588$. Table 1
provides a summary of these parameters ordered by case. For the case of the
star marker (trajectory 1, non-isotropic), we have $v_{*}=0.9$. For the
diamond marker (trajectory 2, non-isotropic), we have $v_{*}=0.7$. (Note that
the star, diamond, and square markers are also used in Fig. 9 to indicate the
location above the merger disk whose angular distributions are used as initial
conditions to evolve the EOMs.)

Table 1: Parameters for the angular distributions shown in Fig. 8 and their location above the torus (see Fig. 9). Marker | $(x_{0},z_{0})[10^{7}\ \mathrm{cm}]$ | $(j_{\nu_{e}},j_{\bar{\nu}_{e}})$ | $(v_{c},\bar{v}_{c})$ | $v_{*}$ | $\alpha$ | Label
---|---|---|---|---|---|---
Square | $(0.20,0.50)$ | – | – | 0.0 | 1.12 | traj. 1 (iso.)
Star | $(0.03,0.50)$ | $(0.853,0.855)$ | $(0.717,0.723)$ | 0.9 | 1.11 | traj. 1 (non-iso.)
Diamond | $(0.04,0.20)$ | $(0.583,0.588)$ | $(0.170,0.180)$ | 0.7 | 1.10 | traj. 2 (non-iso.)

The angular distributions in the left panel of Fig. 8 do not cross each other,
implying the absence of an electron lepton number crossing necessary for fast
flavor evolution [44, 45, 46]. On the other hand, the angular distributions of
$\nu_{e}$ and $\bar{\nu}_{e}$ in the right panel of Fig. 8 exhibit a small
angular crossing (see solid green line); however, we have checked that the
system is stable under fast flavor instabilities. Moreover, the initial
angular distributions for trajectory 1 and trajectory 2 are such that we have
an overall excess of $\bar{\nu}_{e}$ over $\nu_{e}$, as well as an overall
excess of (anti)neutrinos along the forward ($v=1$) direction; antineutrinos
are, however, not only more abundant but also more forward-peaked, resembling
the flux factor distributions.

###  4.3 Characterizing the matter and neutrino-neutrino interaction
potentials

To investigate the neutrino flavor conversion physics, we compute the matter
potential from Eq. 2.6, while the neutrino self-interaction potential follows
from Eq. 2.7, hence for example

$\displaystyle H_{\nu\nu}^{ee}(v,t_{0})=\sqrt{2}G_{\rm F}n_{\nu_{e}}\int
dv^{\prime}\left[\rho_{ee}(v^{\prime},t_{0})-\alpha\bar{\rho}_{ee}(v^{\prime},t_{0})\right][1-vv^{\prime}]\
.$ (4.4)

We rely on the output data of the model M3A8m3a5 to compute the diagonal
elements of $H_{\nu\nu}(v,t_{0})$ as illustrated above.

Figure 9 shows isocontours of the total potential, $\sqrt{2}G_{\rm
F}(n_{\nu_{e}}-n_{\bar{\nu}_{e}})+\lambda$, in the MA isotropic (or SA)
approximation and the $x$–$z$ plane. The red (blue) regions in the $x$–$z$
plane indicate the locations where the total potential is positive (negative).
This means that the MNR condition (Eq. 2.8) is satisfied along the isocontour
which we show through the red dashed line.

In order to solve the EOMs, we pick trajectories 1 and 2, both crossing the
red dashed line in Fig. 9, i.e. the black dotted and dash-dotted lines along
$z=0.5\times 10^{7}$ cm and $z=0.2\times 10^{7}$ cm, respectively. At each
point along these trajectories, we compute the matter and neutrino-neutrino
potentials relying on the outputs of our BH disk model. We construct the
initial angular distributions at the initial $(x_{0},z_{0})$ locations above
the BH torus (see Fig. 8) as illustrated in the previous subsection and then
let the angular distributions evolve due to flavor conversion while neglecting
neutrino advection and non-forward collisions as discussed in Sec. 2.

Figure 10: Left: Matter and neutrino-neutrino potentials along trajectory 1 as
a function of $x$. The spatial domain covers a portion of the horizontal
dotted line in Fig. 9. The MNR condition is fulfilled within $\Delta x=10^{6}$
cm from positive to negative values of the total potential. The purple-shaded
area corresponds to the maximum variation of the total potential. To guide the
eye, the red dashed line shows the zero value for the total potential. Right:
Same as the left panel, but for trajectory 2 (dash-dotted line in Fig. 9).
Given the much closer location to the central BH at which trajectory 2 starts,
the MNR condition can be fulfilled at smaller values of $x$.

The “average” potential felt by neutrinos should be different for the
isotropic and anisotropic MA configurations. In order to take this into
account when exploring where the MNR condition is expected to take place, we
introduce modified total potentials in Fig. 10. For illustration, we consider
two extreme angular distributions: an isotropic angular distribution and a
Delta function [$\delta(v_{*}-v)$]. Since in the MA framework, the MNR
condition depends on $v$, these two angular distributions should meet the MNR
condition at different $x$ and $v$ on average. For the isotropic case, the MNR
condition is given by $\sqrt{2}G_{\rm
F}(n_{\nu_{e}}-n_{\bar{\nu}_{e}})+\lambda=0$; for the Delta distribution, the
MNR condition is on average given by $\sqrt{2}G_{\rm
F}(n_{\nu_{e}}-n_{\bar{\nu}_{e}})(1-v_{*}^{2})+\lambda=0$.

Figure 10 represents the sum of the matter and neutrino self-interaction
potentials as functions of $x$ for trajectory 1 (left panel) and trajectory 2
(right panel), see also Fig. 9. The purple-shaded bands in Fig. 10 represent
the variation of the total potential. The lower limit of the band corresponds
to the case of isotropic angular distributions ($v_{*}^{2}=0$). The upper
limit corresponds to the case where the first term in the modified total
potential is zero ($v_{*}^{2}=1$) and no MNR is expected because such
potential is always positive. Intermediate values of $v_{*}$ shift the
modified potential along $x$, which means that the MNR could occur close to
the central compact object depending on how forward-peaked the angular
distributions are.

For the initial conditions summarized in Table 1, the interval $\Delta
x=10^{6}$ cm is sufficient for all cases to fulfill the MNR condition for a
large number of angular bins. For all cases shown in Table 1, the total
potentials are negative at their initial locations $(x_{0},z_{0})$ but become
positive as one moves towards larger values of $x$, thereby fulfilling the MNR
conditions around $x=0.1$–$0.3\times 10^{7}$ cm depending on the specific
value of $v_{*}$.

We also investigate the MNR conditions along trajectory 2 (see dash-dotted
line in Fig. 9), which lies closer to the neutrino decoupling regions and has
less forward-peaked distributions with respect to the ones of trajectory 1,
while still maintaining a certain level of anisotropy. The initial angular
distributions of trajectory 2 are shown in the right panel of Fig. 8 with its
parameters summarized in Table 1.

###  4.4 Neutrino flavor evolution along selected trajectories

Figure 11: Neutrino flavor conversion along trajectory 1, with isotropic (top
row, Eq. 4.2) and non-isotropic (middle row, Eq. 4.3) angular distributions,
as well as for non-isotropic angular distributions and trajectory 2 (bottom
row). The flavor evolution starts at the $x_{0}$ location marked in Fig. 10
and Table 1. The left panels show the oscillated total potentials. The middle
panels show the temporal evolution of the survival probabilities for neutrinos
and antineutrinos. The right panels show snapshots of the angular
distributions at the selected time $t^{\prime}$ marked by the vertical dotted
line. For the isotropic case (top panels), our results are very similar to the
ones in Fig. 1 where a dipole mode grows during the partial MNR transition.
For this case, we show the corresponding SA system (dash-dotted lines in the
top middle panel). In the non-isotropic case (middle row), there is no MNR
transition. In the bottom panels, we show the results for trajectory 2. Here,
a large number of forward angular bins undergo MNR, and flavor conversion
begins as the first forward bin is involved in the MNR. However, the angular
distributions along trajectory 2 are less forward-peaked than the ones of
trajectory 1, and minimal flavor conversion takes place.

We now present the solution of the EOMs along the selected trajectories 1 and
2 shown in Fig. 9, whose total potentials in the absence of flavor conversion
are displayed in Fig. 10. Animations of the temporal evolution of the
characteristic quantities for the selected trajectories are provided as
Supplemental Material. Note that, because of the simplified setup adopted in
this work, the solution of the EOMs depends on where along a given trajectory
we choose to start the simulation, i.e. on $(x_{0},z_{0})$ and the initial
configuration corresponding to the square, star and diamond markers in Fig. 9
(see also Table 1). The initial location $(x_{0},z_{0})$ is chosen such that
the total potential starts from a negative value that is not too far from
zero. This means that we choose the initial location sufficiently close to the
MNR condition so that our numerical routine can cross the MNR within a
reasonable amount of CPU time. Although the exact location $(x_{0},z_{0})$ is
important for extracting the flux factors (see Table 1), small displacements
do not significantly change our angular distributions. Moreover, since there
are no ELN angular crossings, the system is stable between $(x_{0},z_{0})$ and
the location where the first angular bin ($v=1$) meets the MNR condition.

We adjust the vacuum mixing angle to $\theta_{V}=0.3$ (instead of
$\theta_{V}=0.15$; see Sec. 3). The reason is that the potentials from the BH
disk simulation are such that a slightly larger $\theta_{V}$ is needed to see
an MNR transition even in the SA scenario. Thus, if we do not find MNR for
this large $\theta_{V}$, it is not because of an inadequate choice of vacuum
mixing parameters (see Sec. 3.1 and Appendix A) but because of MA effects.

The top and middle rows of Fig. 11 show the flavor evolution of
(anti)neutrinos in the isotropic (Eq. 4.2) and non-isotropic scenarios (Eq.
4.3 as in Fig. 10) along trajectory 1–black dotted line in Fig. 9. These
results are qualitatively very similar to the ones presented for the ideal
neutrino gas in Sec. 3. The flavor outcome in the SA and MA-isotropic cases is
qualitatively different, as visible from the middle panel in the first row of
Fig. 11. The development of fine angular structures is the result of the
diffusion of flavor waves from large to small angular scales which hinders the
MNR development.

For trajectory 1, the MA isotropic case shows the least amount of flavor
conversion, while the MA non-isotropic case is prone to significant flavor
conversion. This difference in the flavor outcome is because, for forward-
peaked distributions, almost all neutrinos are clustered in a narrow angular
region and undergo the MNR roughly at the same time. In the isotropic case,
neutrinos are equally distributed across $-1<v<1$ and undergo MNR
progressively (more slowly) as $H_{ee}^{\nu\nu}(v)-H_{xx}^{\nu\nu}(v)+\lambda$
crosses zero for any $v$.

For trajectory 2, the bottom panels of Fig. 11 show the flavor outcome in the
MA non-isotropic case. In trajectory 2, (anti)neutrinos almost completely
populate the angular region with $0<v<1$, with no (anti)neutrinos in the
backward direction. Similar to the isotropic case of trajectory 1, no
significant flavor conversion is found. For both trajectory 1 (isotropic) and
trajectory 2 (non-isotropic), the survival probabilities of (anti)neutrinos
remain close to $0.9$, unlike the non-isotropic case in trajectory 1, where
survival probabilities are approximately equal to $0.7$ after $1.5\times
10^{-5}$ s.

These results highlight the crucial differences between SA and MA solutions of
the EOMs and suggest that the MNR may only be relevant far from the decoupling
regions where the angular distributions are forward-peaked. Alternatively, the
MNR is unlikely to lead to significant flavor conversion in regions where the
(anti)neutrino distributions are isotropic.

##  5 Conclusions

In neutron star merger remnants, neutrinos could undergo the MNR,
characterized by a cancellation of the matter and neutrino-neutrino
interaction potentials. In this paper, we explore the flavor conversion
outcome caused by the MNR, solving the EOMs in the SA and MA scenarios. In the
MA case, we explore the flavor conversion physics for initially isotropic and
non-isotropic angular distributions. Under the assumptions of spatial
homogeneity and mono-energetic, ideal (anti)neutrino gas, we find a
significant difference in the flavor outcome between the SA and MA cases.
Intriguingly, the MA calculation of neutrino flavor evolution in the presence
of MNR results in the spontaneous breaking of isotropy for both neutrino mass
orderings.

We further investigate the MNR physics relying on the output of a 2D
hydrodynamical simulation of a BH accretion disk, resulting from the post-
merger phase of a neutron star merger. Focusing on a representative time
snapshot, we solve the EOMs along two selected trajectories above the
accretion disk. We find that, depending on whether we employ initially
isotropic neutrino angular distributions or non-isotropic ones, the flavor
outcome is qualitatively very different. In the former case, no significant
flavor conversion is found after the MNR is partially achieved. In the latter
case, significant flavor conversion is triggered by the MNR condition being
fulfilled for some angular bins; this scenario on average tends to resemble
what would happen in the SA approximation, despite the isotropy-breaking
feature typical of the MA case. These findings suggest that the flavor outcome
resulting from the MNR is sensitive to the shape of the angular distributions,
where forward-peaked distributions lead to more substantial flavor conversion
because the MNR condition is met for the most populated angular bins around
the peak of the distribution almost simultaneously. This is not the case for
initially isotropic neutrino distributions, where the MNR condition is met
across all equally-populated bins over a wider spatial range.

We stress that our modeling of the MNR in the MA framework is simplistic as it
does not include a self-consistent evolution of the neutrino angular
distributions based on non-forward collisions with the medium and neutrino
advection. Yet, these preliminary findings highlight an intrinsic difference
between the SA and MA solutions of the flavor evolution in the presence of the
MNR. Our work also suggests a strong directional dependence of the MNR physics
in neutron star merger remnants and calls for the development of a self-
consistent modeling of the neutrino angular distributions.

## Acknowledgments

We are grateful to Oliver Just for providing access to the output data of the
neutron star merger remnant simulation adopted in this work and Alexander
Friedland for useful discussions. IPG acknowledges support from the U.S.
Department of Energy under contract number DE-AC02-76SF00515. In Copenhagen,
this project has received support from the Danmarks Frie Forskningsfond
(Project No. 8049-00038B), the Villum Foundation (Project No. 13164), and the
European Union (ERC, ANET, Project No. 101087058). Views and opinions
expressed are those of the authors only and do not necessarily reflect those
of the European Union or the European Research Council. Neither the European
Union nor the granting authority can be held responsible for them.

##  Appendix A Characteristic timescales of the matter-neutrino resonance in
the single-angle approximation

A clear feature of the MNR in the SA approximation is the complete conversion
of $\nu_{e}$, while $\bar{\nu}_{e}$ returns to its original configuration. In
Fig. 12, we show our numerical solution together with the analytical
prediction for the survival probabilities. The duration of the single-angle
MNR transition $\Delta t$ is well understood and can be estimated
analytically. Reference [12] derived the adiabaticity criterion from matching
the timescale of the transition ($\delta T_{1}$) with the timescale set by the
capacity of the system to change flavor ($\delta T_{2}$). We estimate these
timescales to validate our SA calculations. The duration of the transition is
set by the timescale

$\displaystyle\delta
T_{1}\approx\tau_{\lambda/\mu}\log{\left[{(1+\alpha)}/{(\alpha-1)}\right]}\ ,$
(A.1)

where $\tau_{\lambda/\mu}$ is the effective timescale of the ratio of the
potentials, namely

$\displaystyle\tau_{\lambda/\mu}=\left|{d\log[\lambda/\mu(t)]}/{dt}\right|^{-1}\
.$ (A.2)

Using the potentials in Eq. 3.1, it is straightforward to show that
$\tau_{\lambda/\mu}=q^{-1}$ s, which gives $\delta T_{1}\simeq 6.49\times
10^{-5}$ s. The agreement between this prediction and our numerical
simulations can be seen in the upper panel of Fig. 1, where the transition
timescale is $\Delta t\simeq 6.50\times 10^{-5}$ s and therefore $\Delta
t\approx\delta T_{1}$.

Figure 12: Left: Unoscillated self-interaction (red) and matter potentials
(green) as functions of time. Right: Survival probabilities of neutrinos
(blue) and antineutrinos (red), obtained by solving the EOMs, as functions of
time. The dashed lines are the predicted survival probabilities as shown in
Eqs. 3-4 in Ref. [12]; we can see that the numerical solution is in excellent
agreement with the theoretical prediction.  Figure 13: Duration of the
transition $\delta T_{1}$ and the timescale of the transition $\delta T_{2}$
as functions of the mixing angle $\theta_{V}$. Timescales are comparable
around our $\theta_{V}=0.15$ marked through the dashed vertical line. Choosing
$\theta_{V}\approx 0$ results in very different different timescales $\delta
T_{1}$ and $\delta T_{2}$, thus preventing MNR transitions.

For MNR transitions to occur, the timescale $\delta T_{1}$ has to be
comparable to the distance scale of the transition $\delta T_{2}$. The latter
is given by

$\displaystyle\delta
T_{2}\approx\frac{\alpha}{\omega\sin{2\theta_{V}}\langle\mathrm{Im}(\rho_{ex}+\alpha\bar{\rho}_{ex})\rangle}\
,$ (A.3)

where the term $\langle\cdot\rangle$ in the denominator is the average value
of the imaginary component of the off-diagonal entries. For $\delta T_{1}$ and
$\delta T_{2}$ to be comparable, $\langle\cdot\rangle$ must adjust to the
vacuum oscillation parameters; clearly, for $\omega\approx 0$ or
$\theta_{V}\approx 0$, this requirement cannot be fulfilled. The mass
hierarchy does not play an important role in the occurrence of MNR
transitions, as it simply shifts the total potential by a negligible amount.
However, the magnitudes of $\omega$ and $\theta_{V}$ set the timescale $\delta
T_{2}$. In Fig. 13, we show the timescales $\delta T_{1}$ (Eq. A.1) as a
function of the vacuum mixing angle $\theta_{V}$, while $\delta T_{2}$ (Eq.
A.3) is fixed by the choice of the potentials. The difference between both
timescales is largest for very small $\theta_{V}$. The dashed line marks our
choice with $\theta_{V}=0.15$, which results in comparable $\delta T_{1,2}$,
so MNR transitions are possible.

In Fig. 14, we show the flavor evolution for $\theta_{V}=0.01$ and
$\theta_{V}=0.15$ (i.e., small and large mixing angles), where the former does
not display MNR transitions while the latter does. Figure 15 confirms that the
mass ordering does not play an important role in the occurrence of MNR
transitions as it simply shifts the total potential by a small amount. The
magnitude of $\omega$, however, sets the timescale $\delta T_{2}$ together
with $\theta_{V}$.

Figure 14: Same as the upper panels of Fig. 1 but for two different values of
the vacuum mixing angle $\theta_{V}$. The survival probabilities (left panel)
do not evolve for small $\theta_{V}=0.01$ (dashed) because the timescales
$\delta T_{1}$ and $\delta T_{2}$ are too far apart (see also Fig. 13) and the
MNR transition is not possible. As a result, the total potential (right panel)
for the system with $\theta_{V}=0.01$ (dashed) continues to increase without
staying zero after $\sim 0.2\times 10^{-4}$ s. Figure 15: Same as the upper
panels of Fig. 1 but for normal ($\omega>0$) and inverted ($\omega<0$) mass
orderings, shown in solid and dotted lines, respectively. In agreement with
Eqs. 3-4 of Ref. [12], both mass orderings display qualitatively similar
dynamics and only small deviations are visible. The value of $\omega$ shifts
the total potential (right panel) by a small amount around $t=0.2\times
10^{-4}$ s and has no overall impact on the outcome.

Although not shown here, we have performed simulations where the roles of
$\mu$ and $\lambda$ are swapped [i.e. $\lambda(t)\propto\exp(-qt)$ and time
independent $\mu$]. In this case, the cancellation of the potentials is also
possible. One would naively expect that an MNR would occur in this scenario as
well. This is, however, not the case. The reason is that, if the roles of
$\mu$ and $\lambda$ are swapped, $\tau_{\lambda/\mu}=q=3\times 10^{4}$ s
(instead of $q^{-1}=0.33\times 10^{-5}$ s), which results in $\delta T_{1}$
and $\delta T_{2}$ being a few orders of magnitude apart. Our numerical
results confirm the absence of MNR transitions for this scenario in agreement
with Fig. 1(b) of Ref. [12].

Figure 16: Evolution of $\rho_{ex}$ for a similar isotropic system shown in
the lower panels of Fig. 1 for different numbers of angular bins. We compare a
low-resolution run with higher-resolution runs. Although the non-linear regime
is affected by the number of angle bins, the general behavior and key features
discussed in Sec. 3 are preserved as the angular resolution is increased. The
breaking of angular symmetry still occurs around $t\simeq 0.3\times 10^{-4}$ s
for the low and high-resolution runs.

##  Appendix B Numerical convergence of the solution of the neutrino equations
of motion

To evolve the (anti)neutrino EOMs in time, we use the Runge-Kutta-
Fehlberg(7,8) method with adaptive stepsize from the odeint library of Boost
[47]. It is well-known that MA calculations can display spurious instabilities
caused by a lack of resolution [48, 49]. To ensure that our simulations are
not prone to spurious instabilities, we adopt very high angular resolution
throughout this work, thus the number of angular bins ranges between
$N_{v}=\mathcal{O}(10^{3}$–$10^{4})$.

The breaking of symmetry in the EOMs has been recognized and well understood
in the context of slow conversion [38], but never before in the presence of
MNR. In order to cross-check that the symmetry breaking does not stem from
inadequate angular binning, we have varied the angular resolution of the
isotropic scenario in Fig. 1. The difference between results in Fig. 16 and
those in Fig. 1 is that in this Appendix we use double precision in the
numerical runs instead of the $500$ bits per number implemented before; thus,
here we show that the breaking of symmetry is not related to the lack of
angular resolution, rather to the introduction of a small but important
floating-point error from the double precision assignment.

In Fig. 16, we show the flavor evolution of $\rho_{ex}$ for the isotropic
scenario in the bottom panels of Fig. 1. Notably, the breaking of angular
symmetry still happens around $t\simeq 0.3\times 10^{-4}$ s. Although not
shown, we have repeated the calculation presented in Fig. 16 for the solvers
Bulirch-Stoer, Rosenbrock, Runge-Kutta-Dopri(5), Runge-Kutta-Cash-Karp(5,4)
from the Boost library [47] as well as using VCABM, DP8, Tsit5 solvers from
DifferentialEquations.jl package of Julia [36, 37], and found that the
breaking of symmetry occurs at the same time regardless of the choice of the
solver.

For the perturbed and non-perturbed isotropic systems in Sec 3.1, we implement
runs with 500 bits per number and $N_{v}=100$ angular bins. On the contrary,
for the non-isotropic runs of Sec. 3.2, since less precision is needed in the
floating point error, double precision is enough (64 bits per number) and one
can increase the number of angular bins to $N_{v}=2000$. For the BH disk
merger calculations in Sec. 4.4, we use double precision and $N_{v}=30000$ to
ensure a reliable final flavor configuration.

##  Appendix C Behavior of the neutrino monopole and dipole due to the matter-
neutrino resonance

In this appendix, we describe how symmetry-breaking solutions arise and under
which conditions they can be expected in different scenarios. In Appendix C.1,
we investigate the instability conditions of the monopole and dipole modes in
the simplest of cases where the matter potential is effectively zero. In
Appendix C.2, we present the linear stability analysis for the monopole and
dipole and compare the analytical predictions of the growth rates to our
numerical results.

###  C.1 Growth of the monopole and dipole

We show how the monopole and the dipole modes trigger flavor instabilities in
opposite neutrino mass orderings for isotropic systems, building on the work
of Ref. [38], we compare the analytical estimates previously known from Ref.
[38] with our results.

For simplicity, we introduce the polarization vectors defined as

$\displaystyle\bm{P}_{v}$ $\displaystyle=$
$\displaystyle\big{(}2\mathrm{Re}[\rho_{ex}(v)],-2\mathrm{Im}[\rho_{ex}(v)],\rho_{ee}(v)-\rho_{xx}(v)\big{)}\
,$ (C.1)

and similar for $\bm{\bar{P}}_{v}$ with
$\rho_{ij}(v)\rightarrow\bar{\rho}_{ij}(v)$. The EOMs in terms of the sum and
difference of the polarization vectors,
$\bm{S}_{v}=\bm{P}_{v}+\bm{\bar{P}}_{v}$ and
$\bm{D}_{v}=\bm{P}_{v}-\bm{\bar{P}}_{v}$, are given by

$\displaystyle\dot{\bm{D}}_{v}$ $\displaystyle=$
$\displaystyle-\omega\bm{B}\times\bm{S}_{v}-\mu\bm{D}_{v}\times\int(1-v\cdot
v^{\prime})\bm{D}_{v^{\prime}}d\Omega^{\prime}\ ,$ (C.2)
$\displaystyle\dot{\bm{S}}_{v}$ $\displaystyle=$
$\displaystyle-\omega\bm{B}\times\bm{D}_{v}-\mu\bm{S}_{v}\times\int(1-v\cdot
v^{\prime})\bm{D}_{v^{\prime}}d\Omega^{\prime}\ ,$ (C.3)

where $d\Omega=dvd\phi$ is the differential solid angle. Although in the main
text we have assumed axial symmetry, here we keep the $\phi$ angle to show
that both monopole and dipole modes can form regardless of the presence of
axial symmetry. Moreover, the vacuum potential is represented by the vector
$\bm{B}=(\sin{2\theta_{V}},0,\cos{2\theta_{V}})$ and we work under the
assumption that the matter term is so large that $\theta_{V}\ll 1$ [50, 51] so
that $\bm{B}\approx\bm{e}_{z}$.

In the linear regime, Eqs. C.2 and C.3 can be approximated as

$\displaystyle\dot{\bm{D}}_{v}$ $\displaystyle\approx$
$\displaystyle-\omega\bm{B}\times\bm{S}_{v}\ ,$ (C.4)
$\displaystyle\dot{\bm{S}}_{v}$ $\displaystyle\approx$
$\displaystyle\Big{[}\omega\bm{D}_{v}+\mu\int(1-v\cdot
v^{\prime})\bm{D}_{v^{\prime}}d\Omega^{\prime}\Big{]}\times\bm{B}\ .$ (C.5)

Here, we have dropped the second term in Eq. C.2 because of negligible
$\mathcal{O}(|\bm{D}_{v}|^{2})$ and used
$\bm{S}_{v}\approx\bm{B}\approx\bm{e}_{z}$. Notice that Eqs. C.4 and C.5 are
equivalent to Eqs. (6a) and (6b) of Ref. [38] under the transformation
$\bm{S}_{v}\rightarrow\bm{S}_{v}-\bm{e}_{z}\equiv\bm{q}_{v}$, which leaves the
evolution equations unchanged.

To unveil the spontaneous breaking of angular symmetry is useful to express
the polarization vectors in terms of the spherical harmonics $Y_{\ell m}(v)$:

$\displaystyle\tilde{\bm{D}}_{\ell m}$ $\displaystyle=$ $\displaystyle\int
Y_{\ell m}^{*}(v)\bm{D}_{v}d\Omega\ ,$ (C.6)
$\displaystyle\tilde{\bm{S}}_{\ell m}$ $\displaystyle=$ $\displaystyle\int
Y_{\ell m}^{*}(v)\bm{S}_{v}d\Omega\ ,$ (C.7)

where $\ell=0,1,...$ and $m=-\ell,...,\ell$. To show the decoupling of the
different $\ell,m$ modes in the linear regime, we express the factor
$(1-v\cdot v^{\prime})$ in the basis of spherical harmonics. Hence Eq. C.5
becomes

$\displaystyle\dot{\bm{S}}_{v}$ $\displaystyle\approx$
$\displaystyle\Big{[}\omega\bm{D}_{v}+4\pi\mu\int\Big{(}Y_{00}(v)Y_{00}^{*}(v^{\prime})-\frac{1}{3}\sum_{m=0,\pm
1}Y_{1m}^{*}(v)Y_{1m}^{*}(v^{\prime})\Big{)}\bm{D}_{v^{\prime}}d\Omega^{\prime}\Big{]}\times\bm{B}\
,$ (C.8)

involving the $\ell=0$ and $\ell=1$ spherical harmonics only. To turn this
into an equation for the $\ell,m$ modes we multiply both sides by $Y_{\ell
m}^{*}(v)$, integrate over $d\Omega$, and use the orthonormality of the
spherical harmonics,$\int d\Omega Y_{\ell
m}^{*}Y_{\ell^{\prime}m^{\prime}}=\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}$,
to obtain:

Figure 17: Flavor evolution in the normal mass ordering ($\omega=1\
\mathrm{km^{-1}}$, left panel) and inverted mass ordering ($\omega=-1\
\mathrm{km^{-1}}$, right panel) for a constant $\mu_{0}=3\times 10^{3}\
\mathrm{km^{-1}}$. The matter potential is set to $\lambda=0$, we also assume
a vanishing mixing angle $\theta_{V}$ and seed the system with a small off-
diagonal seed [$\mathcal{O}(10^{-8})$]. This choice of parameters is a
simplified version of the parameters used in Sec. 3 where $\mu$ does not
change in time. In the case with $\omega=1\ \mathrm{km^{-1}}$, the dipole mode
$\ell=1$ (solid orange) exponentially develops with a well-defined growth
rate, while the monopole $\ell=0$ (solid blue) is stable. In both left and
right panels, the growth rate of the monopole of $\rho_{xx}$ (blue dashed) is
exactly twice the growth rate of the unstable mode (in the case with
$\omega=-1\ \mathrm{km^{-1}}$, the unstable mode is the monopole). This is
because in the linear regime $\rho_{xx}\sim\rho_{ex}^{2}$ (see Sec. 2.2 of
Ref. [8]).

$\displaystyle\dot{\tilde{\bm{S}}}_{\ell
m}\approx\Big{[}\omega\tilde{\bm{D}}_{\ell m}+4\pi\mu\Big{(}\delta_{\ell
0}\tilde{\bm{D}}_{0m}-\frac{\delta_{\ell
1}}{3}\tilde{\bm{D}}_{1m}\Big{)}\Big{]}\times\bm{B}\ ,$ (C.9)

or more compactly

$\displaystyle\dot{\tilde{\bm{S}}}_{\ell
m}\approx\Big{[}\omega+\tilde{\mu}_{\ell}\Big{]}\tilde{\bm{D}}_{lm}\times\bm{B}\
,$ (C.10)

where we have defined the $\ell-$dependent interaction strength
$\tilde{\mu}_{\ell}$ as

$\tilde{\mu}_{\ell}=4\pi\begin{cases}\mu&\text{for $\ell=0$}\ ,\\\
-\mu/3&\text{for $\ell=1$}\ ,\\\ 0&\text{for $\ell>1$}\ .\end{cases}$

The expression above implies that $\ell=0$ and $\ell=1$ independently follow
the same EOMs with their corresponding neutrino-neutrino interaction strength
given by $\mu$ and $-\mu/3$, respectively [38]. Notice that the linearized
equations of motion cannot capture the growth of $\ell>1$ modes because they
start growing much later when the system is already in the non-linear regime.

The equivalent equation for $\dot{\tilde{\bm{D}}}_{\ell m}$ follows from
repeating the same procedure for Eq. C.4. Therefore, the complete set of
equations is then

$\displaystyle\dot{\tilde{\bm{D}}}_{\ell m}$ $\displaystyle\approx$
$\displaystyle\omega\tilde{\bm{S}}_{\ell m}\times\bm{B}\ ,$
$\displaystyle\dot{\tilde{\bm{S}}}_{\ell m}$ $\displaystyle\approx$
$\displaystyle\Big{[}\omega+\tilde{\mu}_{\ell}\Big{]}\tilde{\bm{D}}_{lm}\times\bm{B}\
.$ (C.11)

In Fig. 17, we show how the modes $\ell=0,1$ are decoupled from each other due
to the $\ell-$dependence of $\tilde{\mu}_{\ell}$ in the linear regime, which
results in a re-scaling and a flip of sign in the neutrino-neutrino
interaction strength (see Eq. C.1). In Appendix C.2, we compute the analytical
growth rates of the monopole and the dipole modes and show the agreement with
our numerical simulations of flavor evolution.

###  C.2 Linear stability analysis for the $\ell=0$ and $\ell=1$ modes

We carry out the linear stability analysis for the monopole and dipole to
investigate the symmetry breaking.

$\ell=0$ mode: The monopole

The equations for the monopole $\ell=0$ are given by (see Eqs. C.1):

$\displaystyle\dot{\bm{D}}$ $\displaystyle=$
$\displaystyle\omega\bm{S}\times\bm{B}\ ,$ (C.12) $\displaystyle\dot{\bm{S}}$
$\displaystyle=$ $\displaystyle\omega\bm{D}\times\bm{B}+\mu\bm{D}\times\bm{S}\
,$ (C.13)

where $\int dv\bm{S}_{v}=\bm{S}$ and $\int dv\bm{D}_{v}=\bm{D}$ are the angle-
integrated vectors of the sum and difference total polarization vector,
respectively. We seek to find the eigenfrequency of the system described by
the equations above. The transverse components of $\dot{\bm{D}}$ obey the
following equations:

$\displaystyle\dot{D}^{x}$ $\displaystyle=$ $\displaystyle\omega S^{y}\ ,$
$\displaystyle\dot{D}^{y}$ $\displaystyle=$ $\displaystyle-\omega S^{x}\ ,$
(C.14)

and similarly, the ones for $\dot{\bm{S}}$ are given by

$\displaystyle\dot{S}^{x}$ $\displaystyle=$
$\displaystyle\mu(-D^{z}S^{y}+D^{y}S^{z})+\omega D^{y},$
$\displaystyle\dot{S}^{y}$ $\displaystyle=$
$\displaystyle\mu(D^{z}S^{x}-D^{x}S^{z})-\omega D^{x}\ .$ (C.15)

Using the following linear combinations

$\displaystyle\epsilon_{S}$ $\displaystyle=$ $\displaystyle S^{x}-iS^{y}\ ,$
$\displaystyle\epsilon_{D}$ $\displaystyle=$ $\displaystyle D^{x}-iD^{y}\ ,$
(C.16)

and combining Eqs. C.2 and C.2 we can obtain a system of equations for the
transverse components $\epsilon_{S}$ and $\epsilon_{D}$. The linearization
assumes that the neutrino transverse component $\epsilon_{D,S}$ is much
smaller than $S^{z}$ and $D^{z}$ so that higher order terms such as
$\epsilon_{S,D}^{2}$ can be effectively neglected. Thus, the equations of
motion for the transverse components are

$\displaystyle\dot{\epsilon}_{S}$ $\displaystyle=$ $\displaystyle-i\mu
D^{z}\epsilon_{S}+i(\mu S^{z}+\omega)\epsilon_{D}\ ,$
$\displaystyle\dot{\epsilon}_{D}$ $\displaystyle=$ $\displaystyle
i\omega\epsilon_{S}\ .$ (C.17)

This set of equations can be expressed in matrix form as

$\displaystyle\begin{bmatrix}\dot{\epsilon}_{S}\\\
\dot{\epsilon}_{D}\end{bmatrix}=(-i)M\begin{bmatrix}\epsilon_{S}\\\
\epsilon_{D}\end{bmatrix}\ ,$ (C.18)

where the matrix $M$ is defined as

$\displaystyle M=\begin{bmatrix}\mu D^{z}&-\mu S^{z}-\omega\\\
-\omega&0\end{bmatrix}\ .$ (C.19)

The eigenvalues of $M$ are:

$\displaystyle\Omega^{\pm}_{\ell=0}=\frac{1}{2}\mu
D^{z}\pm\frac{1}{2}\sqrt{(D^{z})^{2}\mu^{2}+4S^{z}\mu\omega+4\omega^{2}}\ .$
(C.20)

According to Eq. C.18, an imaginary $\Omega^{\pm}_{\ell=0}$ would imply that
$\epsilon_{S,D}\propto\exp(\pm\Omega t)$ indicating the existence of an
unstable solution in flavor space. The vacuum oscillation frequency $\omega$
must be negative for $\Omega^{\pm}_{\ell=0}$ to have an imaginary component,
otherwise the argument of the square root is always positive provided that
$\mu>0$. We show the analytical behavior of the root $\Omega^{\pm}_{\ell=0}$
in Fig. 18 (top panels), where one can see that the root displays an unstable
solution only in the inverted mass hierarchy for positive $\mu$.

Figure 18:  Top: Growth rates for the monopole ($\ell=0$) and dipole
($\ell=1$) as analytically computed (Eqs. C.20 and C.21). Bottom:  Comparison
between the prediction of the dipole growth rate (analytical) and the growth
rate computed by solving the EOMs (numerical). All cases assume
$\theta_{V}=0$. This figure shows that solutions that break isotropy (unstable
$\ell=1$) are possible only for $\omega>0$ provided that $\mu>0$. As discussed
in Sec. 3, the MNR provides the conditions for isotropic-breaking solutions to
occur for both mass orderings, i.e. the green dashed lines (left) will extend
to negative $\omega$.

$\ell=1$ mode: The dipole

We can use the growth rate of the monopole (Appendix C.1) to obtain the growth
rate of the dipole by re-defining the neutrino-neutrino interaction strength.
The growth rate of the dipole is given by Eq. C.20 under the transformation
$\mu\rightarrow-\mu/3$, namely

$\displaystyle\Omega^{\pm}_{\ell=1}=-\frac{1}{6}\mu
D^{z}\pm\frac{1}{2}\sqrt{\frac{1}{9}(D^{z})^{2}\mu^{2}-\frac{4}{3}S^{z}\mu\omega+4\omega^{2}}\
.$ (C.21)

Note that the transformation $\mu\rightarrow-\mu/3$ is equivalent to
simultaneously changing $\mu\rightarrow\mu/3$ and $\omega\rightarrow-\omega$.
For this reason, the growth rates of the $\ell=1$ mode are smaller than those
of $\ell=0$ (re-scaled interaction strength) and occur in the $\omega>0$
range, as opposed to the monopole mode.

In Fig. 18 we show show the growth rates $\Omega^{+}_{\ell=0}$ and
$\Omega^{+}_{\ell=1}$. Notably, when solutions are stable in the $\ell=0$
mode, the system could still develop a growing $\ell=1$ mode, therefore
spontaneously breaking the symmetry of the initial angular distributions [52,
53, 38, 35]. Moreover, in the bottom panels of Fig. 18 we compare the
analytical growth rates $\Omega^{\pm}_{\ell=1}$ to the growth rates of our
flavor evolution which have been extracted using the finite difference method;
the agreement between our simulations and the analytical estimates is
excellent.

## References

  * [1] E. Burns, _Neutron Star Mergers and How to Study Them_ , _Living Rev. Rel._ 23 (2020) 4, [1909.06085].
  * [2] H.-T. Janka and A. Bauswein, _Dynamics and Equation of State Dependencies of Relevance for Nucleosynthesis in Supernovae and Neutron Star Mergers_ , pp. 1–98.  Springer Nature Singapore, Singapore, 2023.  2212.07498.  10.1007/978-981-15-8818-1-93-1.
  * [3] B. D. Metzger, _Kilonovae_ , _Living Rev. Rel._ 23 (2020) 1, [1910.01617].
  * [4] L. Wolfenstein, _Neutrino Oscillations in Matter_ , _Phys. Rev._ D17 (1978) 2369–2374.
  * [5] S. P. Mikheyev and A. Yu. Smirnov, _Resonance Amplification of Oscillations in Matter and Spectroscopy of Solar Neutrinos_ , _Sov. J. Nucl. Phys._ 42 (1985) 913–917.
  * [6] H. Duan, G. M. Fuller and Y.-Z. Qian, _Collective Neutrino Oscillations_ , _Ann. Rev. Nucl. Part. Sci._ 60 (2010) 569–594, [1001.2799].
  * [7] A. Mirizzi, I. Tamborra, H.-T. Janka, N. Saviano, K. Scholberg, R. Bollig et al., _Supernova Neutrinos: Production, Oscillations and Detection_ , _Riv. Nuovo Cim._ 39 (2016) 1–112, [1508.00785].
  * [8] I. Tamborra and S. Shalgar, _New Developments in Flavor Evolution of a Dense Neutrino Gas_ , _Ann. Rev. Nucl. Part. Sci._ 71 (2021) 165–188, [2011.01948].
  * [9] S. Richers and M. Sen, _Fast Flavor Transformations_ , pp. 1–17.  Springer Nature Singapore, Singapore.  2207.03561.  10.1007/978-981-15-8818-1-125-1.
  * [10] S. Chakraborty, R. Hansen, I. Izaguirre and G. G. Raffelt, _Collective neutrino flavor conversion: Recent developments_ , _Nucl. Phys._ B908 (2016) 366–381, [1602.02766].
  * [11] A. Malkus, J. P. Kneller, G. C. McLaughlin and R. Surman, _Neutrino oscillations above black hole accretion disks: disks with electron-flavor emission_ , _Phys. Rev. D_ 86 (2012) 085015, [1207.6648].
  * [12] A. Malkus, A. Friedland and G. C. McLaughlin, _Matter-Neutrino Resonance Above Merging Compact Objects_ , 1403.5797.
  * [13] A. Malkus, G. C. McLaughlin and R. Surman, _Symmetric and Standard Matter-Neutrino Resonances Above Merging Compact Objects_ , _Phys. Rev. D_ 93 (2016) 045021, [1507.00946].
  * [14] M.-R. Wu, H. Duan and Y.-Z. Qian, _Physics of neutrino flavor transformation through matter–neutrino resonances_ , _Phys. Lett. B_ 752 (2016) 89–94, [1509.08975].
  * [15] Y.-L. Zhu, A. Perego and G. C. McLaughlin, _Matter Neutrino Resonance Transitions above a Neutron Star Merger Remnant_ , _Phys. Rev. D_ 94 (2016) 105006, [1607.04671].
  * [16] D. Väänänen and G. C. McLaughlin, _Uncovering the Matter-Neutrino Resonance_ , _Phys. Rev. D_ 93 (2016) 105044, [1510.00751].
  * [17] M. Frensel, M.-R. Wu, C. Volpe and A. Perego, _Neutrino Flavor Evolution in Binary Neutron Star Merger Remnants_ , _Phys. Rev. D_ 95 (2017) 023011, [1607.05938].
  * [18] J. Y. Tian, A. V. Patwardhan and G. M. Fuller, _Neutrino Flavor Evolution in Neutron Star Mergers_ , _Phys. Rev. D_ 96 (2017) 043001, [1703.03039].
  * [19] S. Shalgar, _Multi-angle calculation of the matter-neutrino resonance near an accretion disk_ , _JCAP_ 02 (2018) 010, [1707.07692].
  * [20] A. Vlasenko and G. C. McLaughlin, _Matter-neutrino resonance in a multiangle neutrino bulb model_ , _Phys. Rev. D_ 97 (2018) 083011, [1801.07813].
  * [21] G. Sigurdharson, I. Tamborra and M.-R. Wu, _Resonant production of light sterile neutrinos in compact binary merger remnants_ , _Phys. Rev. D_ 106 (2022) 123030, [2209.07544].
  * [22] V. Cirigliano, M. W. Paris and S. Shalgar, _Effect of collisions on neutrino flavor inhomogeneity in a dense neutrino gas_ , _Phys. Lett. B_ 774 (2017) 258–267, [1706.07052].
  * [23] O. Just, A. Bauswein, R. A. Pulpillo, S. Goriely and H.-T. Janka, _Comprehensive nucleosynthesis analysis for ejecta of compact binary mergers_ , _Mon. Not. Roy. Astron. Soc._ 448 (2015) 541–567, [1406.2687].
  * [24] O. Just, M. Obergaulinger, H.-T. Janka, A. Bauswein and N. Schwarz, _Neutron-star merger ejecta as obstacles to neutrino-powered jets of gamma-ray bursts_ , _Astrophys. J. Lett._ 816 (2016) L30, [1510.04288].
  * [25] H.-T. Janka, T. Eberl, M. Ruffert and C. L. Fryer, _Black hole: Neutron star mergers as central engines of gamma-ray bursts_ , _Astrophys. J. Lett._ 527 (1999) L39, [astro-ph/9908290].
  * [26] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, M. Shibata and K. Taniguchi, _Dynamical mass ejection from the merger of asymmetric binary neutron stars: Radiation-hydrodynamics study in general relativity_ , _Phys. Rev. D_ 93 (2016) 124046, [1603.01918].
  * [27] S. Shalgar, I. Padilla-Gay and I. Tamborra, _Neutrino propagation hinders fast pairwise flavor conversions_ , _JCAP_ 06 (2020) 048, [1911.09110].
  * [28] I. Padilla-Gay, S. Shalgar and I. Tamborra, _Multi-Dimensional Solution of Fast Neutrino Conversions in Binary Neutron Star Merger Remnants_ , _JCAP_ 01 (2021) 017, [2009.01843].
  * [29] S. Shalgar and I. Tamborra, _Neutrino decoupling is altered by flavor conversion_ , _Phys. Rev. D_ 108 (2023) 043006, [2206.00676].
  * [30] S. Shalgar and I. Tamborra, _Neutrino flavor conversion, advection, and collisions: Toward the full solution_ , _Phys. Rev. D_ 107 (2023) 063025, [2207.04058].
  * [31] I. Padilla-Gay, I. Tamborra and G. G. Raffelt, _Neutrino fast flavor pendulum. II. Collisional damping_ , _Phys. Rev. D_ 106 (2022) 103031, [2209.11235].
  * [32] D. F. G. Fiorillo, I. Padilla-Gay and G. G. Raffelt, _Collisions and collective flavor conversion: Integrating out the fast dynamics_ , 2312.07612.
  * [33] G. L. Fogli, E. Lisi, A. Marrone and A. Mirizzi, _Collective neutrino flavor transitions in supernovae and the role of trajectory averaging_ , _JCAP_ 12 (2007) 010, [0707.1998].
  * [34] H. Duan, G. M. Fuller, J. Carlson and Y.-Z. Qian, _Simulation of Coherent Non-Linear Neutrino Flavor Transformation in the Supernova Environment. 1. Correlated Neutrino Trajectories_ , _Phys. Rev. D_ 74 (2006) 105014, [astro-ph/0606616].
  * [35] S. Shalgar and I. Tamborra, _Symmetry breaking induced by pairwise conversion of neutrinos in compact sources_ , _Phys. Rev. D_ 105 (2022) 043018, [2106.15622].
  * [36] C. Rackauckas and Q. Nie, _Differentialequations.jl – a performant and feature-rich ecosystem for solving differential equations in julia_ , _The Journal of Open Research Software_ 5 (2017) .
  * [37] J. Bezanson, A. Edelman, S. Karpinski and V. B. Shah, _Julia: A fresh approach to numerical computing_ , _SIAM Review_ 59 (2017) 65–98.
  * [38] H. Duan, _Flavor Oscillation Modes In Dense Neutrino Media_ , _Phys. Rev. D_ 88 (2013) 125008, [1309.7377].
  * [39] O. Just, A. Bauswein, R. A. Pulpillo, S. Goriely and H.-T. Janka, _Comprehensive nucleosynthesis analysis for ejecta of compact binary mergers_ , _Monthly Notices of the Royal Astronomical Society_ 448 (feb, 2015) 541–567.
  * [40] M.-R. Wu, I. Tamborra, O. Just and H.-T. Janka, _Imprints of neutrino-pair flavor conversions on nucleosynthesis in ejecta from neutron-star merger remnants_ , _Phys. Rev. D_ 96 (2017) 123015, [1711.00477].
  * [41] L. M. Murchikova, E. Abdikamalov and T. Urbatsch, _Analytic Closures for M1 Neutrino Transport_ , _Mon. Not. Roy. Astron. Soc._ 469 (2017) 1725–1737, [1701.07027].
  * [42] J. Cernohorsky and S. A. Bludman, _Maximum Entropy Distribution and Closure for Bose-Einstein and Fermi-Dirac Radiation Transport_ , _Astrophys. J._ 433 (Sept., 1994) 250.
  * [43] S. Richers, _Evaluating approximate flavor instability metrics in neutron star mergers_ , _Phys. Rev. D_ 106 (2022) 083005, [2206.08444].
  * [44] I. Izaguirre, G. G. Raffelt and I. Tamborra, _Fast Pairwise Conversion of Supernova Neutrinos: A Dispersion-Relation Approach_ , _Phys. Rev. Lett._ 118 (2017) 021101, [1610.01612].
  * [45] T. Morinaga, _Fast neutrino flavor instability and neutrino flavor lepton number crossings_ , _Phys. Rev. D_ 105 (2022) L101301, [2103.15267].
  * [46] I. Padilla-Gay, I. Tamborra and G. G. Raffelt, _Neutrino Flavor Pendulum Reloaded: The Case of Fast Pairwise Conversion_ , _Phys. Rev. Lett._ 128 (2022) 121102, [2109.14627].
  * [47] Boost, “Boost C++ Libraries.” http://www.boost.org/, 2015.
  * [48] S. Sarikas, D. de Sousa Seixas and G. G. Raffelt, _Spurious instabilities in multi-angle simulations of collective flavor conversion_ , _Phys. Rev. D_ 86 (2012) 125020, [1210.4557].
  * [49] T. Morinaga and S. Yamada, _Linear stability analysis of collective neutrino oscillations without spurious modes_ , _Phys. Rev. D_ 97 (2018) 023024, [1803.05913].
  * [50] H. Duan, G. M. Fuller and Y.-Z. Qian, _Collective neutrino flavor transformation in supernovae_ , _Phys. Rev. D_ 74 (2006) 123004, [astro-ph/0511275].
  * [51] S. Hannestad, G. G. Raffelt, G. Sigl and Y. Y. Y. Wong, _Self-induced conversion in dense neutrino gases: Pendulum in flavour space_ , _Phys. Rev. D_ 74 (2006) 105010, [astro-ph/0608695].
  * [52] G. G. Raffelt and D. de Sousa Seixas, _Neutrino flavor pendulum in both mass hierarchies_ , _Phys. Rev. D_ 88 (2013) 045031, [1307.7625].
  * [53] G. G. Raffelt, S. Sarikas and D. de Sousa Seixas, _Axial Symmetry Breaking in Self-Induced Flavor Conversion of Supernova Neutrino Fluxes_ , _Phys. Rev. Lett._ 111 (2013) 091101, [1305.7140].