On the Occurrence of Crossings Between the Angular Distributions of Electron Neutrinos and Antineutrinos in the Supernova Core

TL;DR

This study models neutrino angular distributions in supernova cores, finding ELN crossings occur mainly near decoupling regions late in supernova evolution, with perturbations unlikely to induce crossings.

Contribution

The paper introduces a spherically symmetric supernova toy-model that predicts conditions for ELN crossings based on neutrino distributions without hydrodynamical instabilities.

Findings

ELN crossings develop near the decoupling region

Crossings occur when electron neutrino and antineutrino densities are comparable

Perturbations are unlikely to induce ELN crossings in this setup

Abstract

Neutrino fast pairwise conversions have been postulated to occur in the dense core of a core-collapse supernova (SN), possibly having dramatic consequences on the SN mechanism and the observable neutrino signal. One crucial condition favoring pairwise conversions is the presence of crossings between the electron neutrino and antineutrino angular distributions (i.e., electron neutrino lepton number crossings, ELN crossings). A stationary and spherically symmetric SN toy-model is constructed to reproduce the development of the neutrino angular distributions in the dense SN core in the absence of perturbations induced by hydrodynamical instabilities. By iteratively solving the neutrino Boltzmann equations including the collisional term, our model predicts that ELN crossings can develop only in the proximity of the decoupling region and for a sharp radial evolution of the baryon density, when the electron neutrino and antineutrino number densities are comparable. Such conditions are likely to occur only in the late SN stages. Interestingly, flavor instabilities induced by spatial or temporal perturbations are unlikely to generate ELN crossings dynamically within our simplified setup.

Figures

Tables

Table 1. Table 1 Ratio of the electron neutrino and antineutrino number densities ( n ν e / n ν ¯ e subscript 𝑛 subscript 𝜈 𝑒 subscript 𝑛 subscript ¯ 𝜈 𝑒 n_{\nu_{e}}/n_{\bar{\nu}_{e}} ) at the radius of the onset of free streaming of ν e subscript 𝜈 𝑒 \nu_{e} for “case A” and “case B,” see text for details.

	$n_{{\nu }_e}/n_{{\overline{\nu }}_e}$
case A	$2.02$
case B	$1.29$

Table 2. Table 2 Ratio of the electron neutrino and antineutrino number densities ( n ν e / n ν ¯ e subscript 𝑛 subscript 𝜈 𝑒 subscript 𝑛 subscript ¯ 𝜈 𝑒 n_{\nu_{e}}/n_{\bar{\nu}_{e}} ) at the radius of the onset of free streaming of ν e subscript 𝜈 𝑒 \nu_{e} for three different post-bounce times t p . b . = 0.25 , 0.5 subscript 𝑡 formulae-sequence p b 0.25 0.5 t_{\mathrm{p.b.}}=0.25,0.5 and 1 1 1  s of our benchmark SN model.

	$n_{{\nu }_e}/n_{{\overline{\nu }}_e}$
$t_{p.b.}=0.25$ s	$1.30$
$t_{p.b.}=0.5$ s	$1.14$
$t_{p.b.}=1$ s	$1.07$

Equations

$$i\left(\frac{\partial \rho (\vec{x},\vec{p},t)}{\partial t}+\vec{v}\cdot \nabla \rho (\vec{x},\vec{p},t)+\vec{F}\cdot {\nabla }_p\rho (\vec{x},\vec{p},t)\right)$$

$$\frac{d{n}_{\nu }}{dE}=\frac{1}{(2\pi {)}^2}\frac{E^2}{e^{(E-{\mu }_{\nu })/T}+1},$$

$$\cos \theta =\sqrt{1-{\left(\frac{r_{\text{max}}}{r}\right)}^2\left(1-{\cos }^2{\theta }_0\right)}.$$

$$\cos {\theta }_i\frac{d{\rho }_i^{\uparrow }(r)}{dr}$$

$$⟨\sigma ⟩=\sum _m\frac{{\int }_{E_{\min }}^{E_{\max }}{\sigma }_m(E)\frac{d{n}_{\nu }}{dE}dE}{{\int }_{E_{\min }}^{E_{\max }}\frac{d{n}_{\nu }}{dE}dE},$$

$$n+\nu (\bar{\nu })$$

$$p+\nu (\bar{\nu })$$

$$\nu (\bar{\nu })+e^{\pm }$$

$$n+e^{+}$$

$$p+e^{-}$$

$$\cos {\theta }_i\frac{d{\rho }_i^{\downarrow }(r)}{dr}$$

$$\cos {\theta }_i\frac{d{\rho }_i^{\uparrow }(r)}{dr}$$

$$n_{{\nu }_e,{\bar{\nu }}_e}={\int }_{-1}^1{n}_{{\nu }_e,{\bar{\nu }}_e}(\cos \theta )d\cos \theta .$$

$${\lambda }_{{\nu }_e,{\bar{\nu }}_e}^{-1}\simeq \sum _l{\left(\frac{1}{⟨\sigma n_{\mathrm{targets}}⟩}\right)}_l^{-1},$$

$${\rho }_{B,\mathrm{caseA}}(r)=10^{14}\exp \left(0.25(5-r)\right)\text{gm/cc},$$

$${\rho }_{B,\mathrm{caseB}}(r)=10^{14}\exp \left(0.5(5-r)\right)\text{gm/cc}.$$

$$N_{{\nu }_e,{\bar{\nu }}_e}=4\pi R_{{\nu }_e,{\bar{\nu }}_e}^2{\int }_0^{\infty }\frac{E^2}{e^{(E-{\mu }_{{\nu }_e,{\bar{\nu }}_e})/T}+1}\sim R_{{\nu }_e,{\bar{\nu }}_e}^2{T}^3;$$

$$n_{{\nu }_e,{\bar{\nu }}_e}(\cos \theta )\propto \{\begin{array}{l}B\left[1-\exp \left[g_{{\nu }_e,{\bar{\nu }}_e}(\cos \theta )\right]\right]\\ \text{if}\cos \theta <\sqrt{1-{\left(\frac{r}{R_{{\nu }_e,{\bar{\nu }}_e}}\right)}^2}\\ 0\\ \text{if}\cos \theta \ge \sqrt{1-{\left(\frac{r}{R_{{\nu }_e,{\bar{\nu }}_e}}\right)}^2},\end{array}$$

$$g_{{\nu }_e,{\bar{\nu }}_e}(\cos \theta )=-2{\kappa }_{{\nu }_e,{\bar{\nu }}_e}{R}_{{\nu }_e,{\bar{\nu }}_e}\sqrt{1-{\left(\frac{r}{R_{{\nu }_e,{\bar{\nu }}_e}}\right)}^2(1-{\cos }^2\theta )}.$$

Full text

1. Introduction

Core collapse supernovae (SNe) are among the densest sources in neutrinos (Mirizzi et al., 2016; Janka, 2012). According to our current understanding, neutrinos are emitted as matter accretes on the proto-neutron star; they transport energy to finally revive the stalled shock-wave and power the SN explosion (Janka, 2017). However, a detailed picture of the role of neutrinos remains unclear.

Understanding the impact of neutrino flavor conversions on the SN inner working is one of the crucial unsolved problems. In fact, neutrino flavor conversions have been traditionally neglected in SN neutrino transport, since they were expected to occur at radii larger than the shock radius and assumed to have a negligible impact on the shock revival; see e.g. Dasgupta et al. (2012) where effects coming from $\nu$–$\nu$ interactions were included in hydrodynamical simulations, albeit not in a self-consistent setup.

Recent work suggests that the large density of neutrinos in the proximity of the neutrino decoupling region may favor the development of pairwise conversions (Sawyer, 2016, 2009; Izaguirre et al., 2017). Contrary to our intuition, neutrino pairwise conversions could occur even in the absence of a hierarchy among the neutrino mass eigenstates and are independent on the neutrino energy; however, fast pairwise conversions crucially depend on the exact shape of the neutrino angular distributions and on the local neutrino density. The scale ruling such phenomenon is determined by $(G_{F}|n_{\nu_{e}}-n_{\bar{\nu}_{e}}|)^{-1}\simeq\mathcal{O}(10)$ cm, with $G_{F}$ being the Fermi constant and $n_{\nu_{e}(\bar{\nu}_{e})}$ the local (anti)neutrino number density. As a direct consequence of the “fast” interaction rate, neutrino pairwise conversions could possibly lead to flavor decoherence (Abbar & Volpe, 2019; Dasgupta et al., 2017; Capozzi et al., 2019), affecting the neutrino heating within the gain radius as well as the SN nucleosynthesis, and the neutrino signal observable at Earth.

The criteria leading to the development of fast pairwise conversions were investigated in Refs. (Izaguirre et al., 2017; Capozzi et al., 2018, 2017). A non-negligible flux of neutrinos streaming backwards towards the proto-neutron star core as well as the presence of a crossing between the angular distributions of electron neutrinos and antineutrinos have been pinpointed as crucial ingredients to trigger fast instabilities in the neutrino flavor space. A crossing between the angular distributions of $\nu_{e}$ and $\bar{\nu}_{e}$, i.e. a change of sign of $(n_{\nu_{e}}-n_{\bar{\nu}_{e}})$ within a specific angular range, will be throughout dubbed electron neutrino lepton number crossing (ELN crossing).

While ELN crossings are bound to occur in the case of compact binary mergers because of the particular merger geometry (Wu & Tamborra, 2017; Wu et al., 2017); the situation seems to be more complex in the SN context. In fact, self-consistent spherically symmetric SN simulations in 1D did not find any evidence for ELN crossings over a range of progenitor masses (Tamborra et al., 2017). On the other hand, multi-D hydrodynamical simulations have reported mixed results (Abbar et al., 2019; Delfan Azari et al., 2019) depending on the different degree of approximation adopted in the implementation of the neutrino transport. In particular, the occurrence of LESA, the lepton-emission self-sustained asymmetry (Tamborra et al., 2014a, b) has been recently confirmed by a larger set of independent simulations (Janka et al., 2016; Walk et al., 2019; Vartanyan et al., 2019; O’Connor & Couch, 2018; Glas et al., 2018). LESA has been also invoked as a possible ingredient favoring the development of the ELN crossings, see e.g. Izaguirre et al. (2017); Dasgupta et al. (2017).

The recent developments described above hint towards the need of a full implementation of the neutrino quantum kinetic equations, including collisions and flavor conversions (Volpe et al., 2013; Stirner et al., 2018; Sigl & Raffelt, 1993; Vlasenko et al., 2014; Blaschke & Cirigliano, 2016; Cirigliano et al., 2017). However, given the major numerical complications involved in the problem, this has not been tackled numerically in a self-consistent manner yet; see Richers et al. (2019) for a recent simplified attempt in this direction which does not include fast pairwise conversions.

Hydrodynamical simulations of the core collapse do not include the physics of neutrino flavor conversions. Despite that, the simulations tracking the evolution of the neutrino angular distributions, see e.g. Refs. (Tamborra et al., 2017; Sumiyoshi et al., 2015; Nagakura et al., 2018), are extremely expensive in terms of computational resources. As a consequence, it is challenging to explore the conditions leading to the development of ELN crossings and to track to the behavior of the neutrino angular distributions in the late stages of the core collapse.

In this work, by assuming spherical symmetry and that a stationary solution is reached, we develop a simplified model to mimic the SN core microphysics and iteratively solve the neutrino Boltzmann equations with the inclusion of the collisional term. Our goal is to estimate the neutrino angular distributions of $\nu_{e}$ and $\bar{\nu}_{e}$ in the neutrino trapping region and follow their evolution until neutrinos are fully decoupled in the different SN phases. By neglecting complications coming from the SN hydrodynamics, we aim at identifying the microphysical ingredients leading to the development of ELN crossings and discuss their implications for the physics of flavor conversions.

The outline of the paper is as follows. In Sec. 2, we present our stationary and spherically symmetric SN model and the iterative method employed to solve the neutrino Boltzmann equations by taking into account the collision term. In Sec. 3, we investigate the occurrence of ELN crossings by introducing a simple evaluation criterion; we first test the latter for a simplified scenario and then discuss its implications through a more realistic framework based on inputs from a 1D hydrodynamical SN model. Section 4 focuses on exploring whether ELN crossings can be generated dynamically by the occurrence of neutrino flavor conversions because of spatial or temporal perturbations. A discussion on our findings and conclusions are reported in Sec. 5.

2. Neutrino Equations of Motion in the dense proto-neutron star

In this Section, we first introduce the neutrino kinetic equations and then describe the method employed to solve them in spherical symmetry and derive the neutrino angular distributions.

2.1. Neutrino Kinetic Equations

In the interior of a SN, the neutrino number density as well as the baryon density are very large and neither neutrino-neutrino self-interactions (coherent forward scattering) nor neutrino-matter scatterings (incoherent scattering) can be ignored. The coherent and incoherent parts of the neutrino evolution can have a significant feedback on each other. In fact, the neutrino-matter cross-section is flavor dependent, while the non-linear neutrino flavor evolution depends on the spatial distribution of neutrinos, which is determined by the neutrino matter scattering cross-section. However, the interplay between the non-linearity arising from the neutrino self-interaction and scatterings with the nuclear matter remains enigmatic.

Throughout this paper, for the sake of simplicity, we assume that ($\nu_{e},\nu_{x}$) describes the neutrino flavor basis with $\nu_{x}$ being a mixture between $\nu_{\mu}$ and $\nu_{\tau}$ and similarly for antineutrinos. The equation of motion describing the neutrino flavor evolution, including coherent and incoherent components, can be conveniently expressed in terms of a $2\times 2$ density matrix, $\rho$ (and $\bar{\rho}$ for antineutrinos).

The neutrino kinetic equation, for each neutrino momentum mode $\vec{p}$, is given by

[TABLE]

where $\vec{x}$ and $t$ describe the spatial and temporal coordinates, and $\vec{v}$ is the neutrino velocity. The Hamiltonian $H$ on the right-hand-side of the Eq. 1 contains three components, a term embedding the neutrino mixing in vacuum, a term describing the interactions of neutrinos with the SN matter background, and a term describing the neutrino-neutrino interactions, see e.g. Mirizzi et al. (2016). Since we are interested in exploring the development of the ELN crossings, we will neglect neutrino flavor conversions hereafter, unless otherwise stated. The term $\mathcal{C}[\rho,\bar{\rho}]$ describes the scattering of neutrinos with the medium.

The left-hand-side of the Eq. 1 contains the convective term, which can be reduced to [$\cos\theta{\partial\rho(\vec{x},\vec{p},t)}/{\partial r}$] under the assumption of spherical symmetry and stationarity. It should be noted that, due to the non-linearity of the neutrino-neutrino term, the convective term may lead to spontaneous breaking of spherical symmetry (Raffelt et al., 2013; Duan & Shalgar, 2015). For the sake of simplicity, in the following, we will assume that spherical symmetry is preserved and ignore the convective term. The term $\vec{F}\cdot\nabla_{p}\rho$ tracks the effect of external forces on the neutrino field; we will assume that the latter has a negligible contribution. Hereafter we will also use $\hbar=c=1$.

The time-dependent density matrix, $\int\rho(\vec{x},\vec{p},t)d\Omega$ (with $d\Omega=d\cos\theta d\phi$ being the differential solid angle element), is normalized such that the terms on the diagonal give the total number density of neutrinos for fixed $(\vec{r},t)$. Notably, deep in the SN core the electron (anti-)neutrinos are highly degenerate and are in thermal equilibrium with the SN medium. Therefore, the (anti-)neutrino number densities are well described by a Fermi-Dirac distribution:

[TABLE]

where the $\mu_{\nu}$ is the neutrino chemical potential, such that $\mu_{\nu_{e}}=-\mu_{\bar{\nu}_{e}}$ and $\mu_{\nu_{x}}=\mu_{\bar{\nu}_{x}}=0$, and $T$ is the temperature of the SN matter. The off-diagonal terms of the density matrix, instead, being related to the off-diagonal terms of the evolution operator, provide a measure of quantum coherence that determines the flavor transition probabilities.

Solving Eq. 1 may not appear to be a daunting task at first glance. However, it should be noted that the collision term $\mathcal{C}[\rho,\bar{\rho}]$ (and $\bar{\mathcal{C}}[\rho,\bar{\rho}]$) entails a six dimensional integral for each $(\vec{x},\vec{p},t)$. The partial differential equation cannot be solved using standard discretization techniques in a 7D space due to the limitations of available computational power. Imposing symmetries on the system may be misleading due to spontaneous breaking of symmetries in momentum and real space (Mirizzi et al., 2016).

The exact shape of the neutrino angular distributions in the SN core depends on the interactions described by the collision term. In general, the collision term, $\mathcal{C}[\rho,\bar{\rho}]$, contains two components: the loss and gain terms. The loss term accounts for the reduction of the neutrino number flux for a given direction and momentum due to scattering to another direction and momentum; the gain term, instead, accounts for scattering of neutrinos into a given direction and momentum state. The exact form of the collision term adopted in this paper will be described in the next subsection. As we will discuss later, because of the neutrino-matter interactions described by the collision term, the radial profile of the baryon density strongly affects the development of ELN crossings.

2.2. Stationary and Spherically Symmetric Supernova Model

In this Section, we introduce the stationary and spherically symmetric model that we develop to explore under which conditions the growth of ELN crossings is favored. Our model is based on the following assumptions:

The number of neutrinos is locally conserved. Absorption and emission are thus treated as effective isotropic elastic scatters. 2. 2.

Energy-averaged, flavor-dependent neutrino distributions are adopted. 3. 3.

Only two neutrino flavor eigenstates are considered. 4. 4.

Any dependence on $t$ and on the azimuthal angle $\phi$ is discarded.

Our spherically symmetric model is sketched in Fig. 1. We assume that neutrinos are radiated by an inner surface of radius $r_{\textrm{min}}$ and propagate until they reach an outermost surface of radius $r_{\textrm{max}}$. Each point $P$ in the SN sphere is characterized by the radius $r$ and the angle $\theta_{0}$, where $\theta_{0}$ has been defined with respect to the outermost surface. For each $P$, we further introduce a set of coordinates $\theta$ to characterize the local angular distribution of neutrinos. The angle in the global coordinate system $\theta_{0}$ is related to the local angle $\theta$ by the following relation:

[TABLE]

For fixed $P$, if $\mathrm{Im}(\cos\theta)\neq 0$, then the correspondent neutrino trajectory is discarded (i.e., the evolution equation is not evolved along that particular trajectory for radii such that $\mathrm{Im}(\cos\theta)\neq 0$).

From Eq. 1 and Fig. 1, one can immediately realize that we are not dealing with an initial value problem and need to define boundary conditions instead of initial conditions. To obtain a stationary solution in the presence of collisions, we will adopt an iterative method as described in the following.

Given that the energy-dependent quantities entering the neutrino equation of motion are averaged over the neutrino energy distribution, we will effectively solve Boltzmann equations depending on the neutrino scattering angle, but not on their momentum. In the first step of our iteration method, we ignore the backward flux of neutrinos (i.e., assume a null neutrino flux for $\cos\theta<1$) and evolve the following equation of motion

[TABLE]

where the subscripts ($i$, $j$) denote the indexes of the angular bins. $\Delta\cos\theta_{j}$ is the width of the $j^{\textrm{th}}$ angular bin, which depends on the radius and is calculated at each radial step, as we will see in the following. Note that since we are considering only effective scatterers and no flavor evolution, the collision integral simplifies to $\mathcal{C}[\rho,\bar{\rho}]\sim C\times\rho$ for the loss and the gain term. The $\cos\theta_{i}$-term on the left hand side of Eq. 4 takes into account the dependence of the path length on the zenith angle.

In order to locally conserve the neutrino number, we set the loss coefficient, $\mathcal{C}^{\textrm{loss}}$, equal to the gain coefficient, $\mathcal{C}^{\textrm{gain}}$. The loss and gain coefficients are proportional to the product of the number of effective scatterers and the cross-section for each interaction channel averaged over the neutrino energy distribution, i.e. $\mathcal{C}^{\textrm{loss}},\mathcal{C}^{\textrm{gain}}=F_{B}n_{\mathrm{targets}}\langle\sigma\rangle$, with $\langle\sigma\rangle$ the cross section averaged over the neutrino energy distribution:

[TABLE]

where the sum over $m$ is a sum over all the channels discussed below, and $[E_{\mathrm{min}},E_{\mathrm{max}}]=[0,500]$ MeV. The cross-sections entering $\mathcal{C}^{\textrm{loss}}$ and $\mathcal{C}^{\textrm{gain}}$ have been calculated for the following reactions:

[TABLE]

Each cross-section has been implemented as prescribed in Bowers & Wilson (1982) with the low density approximation adopted for the neutrino-electron cross section. We also employed flavor-dependent neutrino distributions defined along the lines of Eq. 2, and the respective mass-dependent Fermi-Dirac distributions have been adopted for the nucleons. In particular, the local density of nucleons has been defined through the baryon density $\rho_{B}$ and the electron fraction $Y_{e}$ [$n_{n}\simeq\rho_{B}(1-Y_{e})$ for neutrons and $n_{p}\simeq\rho_{B}Y_{e}$ for protons]. Note that, although neutrinos and antineutrinos undergo the same kind of neutral current interactions, a small difference appears in their cross sections. As we will discuss afterwords, this will also affect the eventual development of ELN crossings. For each reaction, the Pauli blocking factor $F_{B}$ of scattering targets has been computed following Refs. (Raffelt, 1996; Bruenn, 1985) and included in the collision term. Notably, in the deepest regions of the proto-neutron star, $\mathcal{C}^{\textrm{loss}}$ and $\mathcal{C}^{\textrm{gain}}$ are suppressed by Pauli blocking effects, however as it will be shown in the following the ELN crossings occur in the proximity of the decoupling region and therefore Pauli blocking effects do not affect our results.

Equation 4 is solved from $r_{\textrm{min}}$, arbitrarily set to 5 km, up to $r_{\textrm{max}}$, fixed at 60 km. Note that the choice of $r_{\textrm{min}}$ guarantees that the collisional rate is large enough to re-distribute the neutrinos in $\theta$ homogeneously, while $r_{\textrm{max}}$ has been fixed outside of the neutrino trapping region and is large enough to not affect the neutrino angular distributions at decoupling. During each iteration, the resulting flux of $(e,x)$ neutrinos and antineutrinos is stored at 100 radial points per km for each angular bin.

The the loss term in Eq. 4 depends on the number of neutrinos in the $i^{\textrm{th}}$ bin and on the phase space. However, we perform the integration numerically to ensure that the numerical error arising from the discretization of the gain term is canceled by that of the loss term.

We bin the neutrino trajectories with respect to the global angle defined at $r_{\textrm{max}}$ denoted by $\theta_{0}$. In the numerical runs, the angular grid has been chosen to be uniform in $\cos^{2}\theta_{0}$. For all the results presented in this paper 200 angular bins in the global coordinate system have been used. We tested the convergence of our results by increasing the number of bins and made sure that the final angular distributions are stable with respect to the resolution adopted in the numerical runs within variations of less than a percent. Since we use an adaptive embedded Prince-Dormand(8,9) method (Contributors, 2010) to solve the differential equations, the components of the density matrix between the 100 radial grid points for which the data have been stored are obtained through linear interpolation.

After the completion of this first step of our iteration method, we evolve the equations of motion for the backward flux from $r_{\textrm{max}}$ to $r_{\textrm{min}}$

[TABLE]

Here, the superscript ‘pr’ is used to denote that the values for the components of the density matrix are interpolated using the solutions of Eq. 4. The width of the angular bins, indexed by $j$ in the forward direction and $j^{\prime}$ in the backward direction, is the same for forward and backward going neutrinos ($\Delta\cos\theta_{j}=\Delta\cos\theta_{j^{\prime}}$). For the sake of clarity, we have divided the right hand side of Eq. 6 in two parts, each accounting for the phase space in the forward and backward direction respectively. In this step, the gain term in Eq. 6 receives a contribution from the forward flux stored during the previous iteration.

Upon completion of the first iteration (forward and backward), we obtain the initial conditions of the equations of motion for the next iteration round. Since $r_{\textrm{min}}$ is chosen to be in a region of extremely large matter density, the angular distribution is essentially uniform and we can set the forward flux at $r_{\textrm{min}}$ equal to the backward flux for all angular bins. After each backward iteration, we rescale the normalization of the backward flux by an amount that is proportional to the flux at $r_{\textrm{min}}$ to compensate for the loss of neutrinos by diffusion at $r_{\textrm{max}}$, and achieve a steady state solution. The relative normalization between $\nu_{e}$ and $\bar{\nu}_{e}$ is thus determined by the dynamics of collisions in our model; hence the number density at $r_{\textrm{min}}$ is determined by the mean free path – a smaller mean free path implies a larger number of neutrinos of a certain flavor. This is a crucial aspect, as the occurrence of ELN crossings (or lack thereof) is determined in large part by the relative normalization of the $\nu_{e}$ and $\bar{\nu}_{e}$ angular distributions.

Using the initial conditions obtained by the former backward iteration, the equations of motion are then evolved in the forward direction, while using the interpolated values from the solution of Eq. 6:

[TABLE]

It should be noted that if we add Eqs. 6 and 7, multiply by $\Delta\cos\theta_{i}$ and sum over the index $i$, the coefficient of $\mathcal{C}^{\textrm{loss}}$ and $\mathcal{C}^{\textrm{gain}}$ is the same apart from the sign. This implies that the change in the total number of neutrinos is zero.

We repeatedly solve Eqs. 6 and 7 for neutrinos and antineutrinos using interpolated values for the components of the density matrix that have the superscript ‘pr’. We find that around 15 iterations are sufficient to achieve numerical convergence of the results. Notably, this procedure guarantees that the steady state solution is independent on the initial conditions used in the first iteration and determined by self-consistency alone.

In order to be certain that our simple SN model reproduces flavor-dependent neutrino angular distributions in agreement with the literature, we tested that our simple model gives numerical results in agreement with the analytical ones presented in Murchikova et al. (2017) for the case of a uniform sphere emitting blackbody radiation, see Appendix for details.

To test a more realistic scenario, we adopted the inputs of the SN models employed in Ref. (Tamborra et al., 2017) and computed the expected angular distributions. Our results are qualitatively in very good agreement with the evolution of the neutrino angular distributions presented in Tamborra et al. (2017), e.g. their Figs. 4 and 6. However, given the approximate modeling of the collision term, we can reproduce well the evolution as a function of the radius of the angular distributions for each flavor as well as the relative normalization among the angular distributions of different flavors, but not their absolute normalization. Hence, in the following we will show the neutrino angular distributions normalized to the total number of neutrinos which still serves our purposes.

Figure 2 shows an illustrative example of the resultant $\nu_{e}$ differential angular distribution, $n_{\nu_{e}}(\cos\theta)$, normalized to the total number density of $\nu_{e}$ ($n_{\nu_{e}}$) at $35$ km, i.e, after the neutrino decoupling. The differential angular distributions $n_{\nu_{e}}(\cos\theta)$, extracted in correspondence of the angular bin centered on $\cos\theta$, has been normalized to the total number density $n_{\nu_{e}}+n_{\bar{\nu}_{e}}$, defined as,

[TABLE]

The angular distribution has been plotted at different radii and has been obtained by using the inputs from a 1D hydrodynamical model of a $18.6\,M_{\odot}$ SN with SFHo nuclear equation of state (Garching Supernova Archive, ) for $t_{\mathrm{p.b.}}=0.25$ s, which we will use in Sec. 3.2. One can see that the neutrino angular distribution is almost isotropic in $\cos\theta$ when neutrinos are trapped and becomes forward peaked as neutrinos approach the free streaming regime. In particular, the $\nu_{e}$ angular distribution becomes more peaked than the $\bar{\nu}_{e}$ as expected by the different interaction rates of electron neutrinos and antineutrinos. We refer the reader to Refs. (Tamborra et al., 2017; Ott et al., 2008; Sarikas et al., 2012) for more details on the radial evolution of the flavor-dependent neutrino angular distributions.

3. Criteria for the appearance of crossings in the electron neutrino lepton number distribution

In this Section, we explore the microphysical conditions leading to the occurrence of ELN crossings. We first focus on a simple toy model in order to pinpoint the main ingredients favoring the development of ELN crossings and present a criterion to predict the crossing occurrence in the absence of SN hydrodynamical instabilities. The development of ELN crossings in a more realistic SN setup is then investigated by employing inputs from a 1D hydrodynamical simulation of the core collapse.

3.1. Toy-model example

In order to explore what are the conditions leading to the occurrence of ELN crossings, we perform simulations by employing the iterative method described in Sec. 2.2. For simplicity, we rely on a simplified setup that does not reproduce the conditions obtained in a realistic SN model; however, we choose to rely on a simplified scheme to explore the main ingredients leading to the development of ELN crossings. We consider a model with medium temperature and $\nu_{e}$ chemical potential not varying with the radius ($T=\mu_{\nu_{e}}=10$ MeV).

In order to characterize the neutrino interaction strength, the total energy-averaged neutrino and antineutrino mean free path is defined as

[TABLE]

where $l$ denotes the various neutrino-matter reaction channels listed in Sec. 2.2; the neutrino number densities and cross-sections have been modeled as described in Sec. 2.2, and the Pauli blocking term $F_{B}$ has been ignored for simplicity. One can easily show that the dominant interaction rates setting the difference between the $\nu_{e}$ and $\bar{\nu}_{e}$ angular distributions, eventually leading to ELN crossings, are the charged current interactions with a subleading contribution coming from the neutral current interactions. Hence, $\lambda_{\nu_{e}}/\lambda_{\bar{\nu}_{e}}$ and the resultant local number density of $\nu_{e}$ and $\bar{\nu}_{e}$ are the most important quantities in setting the relative ratio between the angular distributions of $\nu_{e}$ and $\bar{\nu}_{e}$, as also discussed in the Appendix. We arbitrarily assume a constant $\lambda_{\nu_{e}}/\lambda_{\bar{\nu}_{e}}\simeq 0.3$ for an electron fraction not varying with the radius ($Y_{e}=0.1$).

The baryon density profile is considered to fall exponentially with respect to the radius and distinguish between two cases. A first case involves a shallow baryon density profile (“case A”),

[TABLE]

and a second case that includes a steeply falling baryon density profile (“case B”)

[TABLE]

Figure 3 shows the stationary solution obtained for the $\nu_{e}$ (in red) and $\bar{\nu}_{e}$ (in green) angular distributions as a function of $\cos\theta$ for “case A” on the top and for “case B” on the bottom. The differential angular distributions $n(\cos\theta)$ have been normalized to the total number density $n_{\nu_{e}}+n_{\bar{\nu}_{e}}$. Note that the angular distributions of $\nu_{x}$ and their antineutrinos are also computed, but we do not show them here for simplicity. The plotted angular distributions have been extracted at $r\simeq 40$ km and $r\simeq 30$ km, respectively; i.e., after the neutrino decoupling. One can see that while “case A” does not lead to the formation of an ELN crossing, “case B” does. For completeness, the bottom panels of “case A” and “case B” in Fig. 3 show the resultant ELN distribution ($\nu_{e}-\bar{\nu}_{e}$) and as a function of $\cos\theta$.

The top panel of Fig. 4 shows the radial evolution of ELN$(\cos\theta=1)$ for “case A” (in violet) and “case B” (in blue)111The ELN has been plotted as indicated in the caption of Fig. 3 but for the angular bin containing $\cos\theta=1$. Note that, since we relied on the local system of coordinates defined by $\theta$, the width of the bin that includes $\cos\theta=1$ depends on the radius. This is not true in the global system of coordinates defined by $\theta_{0}$ and the number of particles remains constant after decoupling.. For each case, the onset of the free-streaming regime of $\nu_{e}$ and $\bar{\nu}_{e}$ is marked by two vertical lines to guide the eye; since we only need an approximate estimation of the region where the neutrinos start to stream freely and are fully decoupled from matter, we define the radius of onset for the free-streaming regime as the radius where the forward flux ($\cos\theta>0$) is 10 times larger than the backward flux ($\cos\theta<0$). One can see that the free-streaming regime is reached at larger $r$ for “case A,” as expected. Moreover, the onset radius of free-steaming for $\nu_{e}$ is always larger than the one of $\bar{\nu}_{e}$ given the difference in the interaction strengths. The two radii are separated by a larger distance in “case A.” The bottom panel of Fig. 4 shows the radial evolution of the ELN without integrating over the $\cos\theta$ bin. As expected, the total number of $\nu_{e}$ and $\bar{\nu}_{e}$ is constant at radii larger than the decoupling one.

By looking at Figs. 3 and 4, one can see that ELN crossings originate in the proximity of the neutrino free-streaming region in “case B.” Moreover, because of the different baryon profiles employed in “case A” and “case B,” the decoupling regions of $\nu_{e}$ and $\bar{\nu}_{e}$ occur in very different spatial regions for “case A” while they are close to each other for “case B.” This implies that the $\nu_{e}$ and $\bar{\nu}_{e}$ angular distributions are more similar to each other in the proximity of the decoupling region in “case B” than in “case A.” This is also proved by the fact that the total local number densities of $\nu_{e}$ and $\bar{\nu}_{e}$ are more similar to each other in the proximity of the decoupling region in “case B,” as shown in Table 1.

This toy-model example shows that ELN crossings can develop only in the proximity of the neutrino decoupling region. Moreover, the decoupling of $\nu_{e}$ should occur in a radial region close to where the one of $\bar{\nu}_{e}$ happens and their local number densities should be comparable. If those conditions are fulfilled, ELN crossings are likely to develop.

3.2. Supernova model example

We now extend our findings to a more complex case involving the radial dependence of the main SN quantities. We base our estimations on the inputs from a 1D hydrodynamical model of a $18.6\,M_{\odot}$ SN with SFHo nuclear equation of state and gravitational mass of $1.4\,M_{\odot}$ (Garching Supernova Archive, ) that we adopt as benchmark case, and select post-bounce time snapshots representative of the different SN phases. Note that the hydrodynamical simulation does not provide the neutrino angular distributions that are instead estimated iteratively through our stationary and spherically symmetric model.

The top panel of Fig. 5 shows the baryon density as a function of the radius for three different post-bounce times $t_{\mathrm{p.b.}}=0.25,0.5$ and $1$ s, in violet, blue and cyan respectively. The radii of the onset of free-streaming of $\nu_{e}$ and $\bar{\nu}_{e}$ are marked through the vertical lines to guide the eye. One can see that, as $t_{\mathrm{p.b.}}$ increases, neutrinos start to free-stream at smaller radii, closer to the SN core. Moreover, as time progresses, the baryon density profile becomes steeper in the proximity of the region of the onset of free streaming. Another interesting aspect is that, for fixed $t_{\mathrm{p.b.}}$, the free-streaming radii of $\nu_{e}$ and $\bar{\nu}_{e}$ differ from each other at earlier post-bounce times and become similar to each other at later times.

For earlier times, the electron fraction $Y_{e}$ is larger at radii smaller than the free-streaming one as shown in the middle panel of Fig. 5. The bottom panel of Fig. 5 shows the ratio of the mean free paths $\lambda_{\nu_{e}}/\lambda_{\bar{\nu}_{e}}$ as a function of the radius; in the proximity of the free-streaming radius, $\lambda_{\nu_{e}}/\lambda_{\bar{\nu}_{e}}\leq 1$. The ratio between the $\nu_{e}$ and $\bar{\nu}_{e}$ number densities at the radius of onset of free streaming is reported in Table 2 for the three studied $t_{\mathrm{p.b.}}$. One can see that $n_{\nu_{e}}/n_{\bar{\nu}_{e}}\rightarrow 1$ as $t_{\mathrm{p.b.}}$ increases.

Figure 6 shows the resultant angular distributions for $\nu_{e}$ and $\bar{\nu}_{e}$, each normalized to the total number density of $\nu_{e}$ and $\bar{\nu}_{e}$, as a function of $\cos\theta$ for $t_{\mathrm{p.b.}}=0.25,0.5$ and $1$ s from top to bottom, respectively. Those angular distributions have been obtained with our iterative method by employing the SN inputs shown in Fig. 5, together with the radial profiles of the chemical potentials for neutrinos and nucleons and the medium temperature extracted from the hydrodynamical simulation. We plot the angular distributions at an arbitrary radius of $40$ km, i.e. after the neutrino decoupling occurred.

For $t_{\mathrm{p.b.}}=0.25$ s, no ELN crossing is found as the number density of $\nu_{e}$ remains larger than the one of $\bar{\nu}_{e}$ due to the shallow baryon density profile. At late times ($t_{\mathrm{p.b.}}\gtrsim 0.5$ s), the baryon density profile becomes steeper; as a result the free-streaming radii of $\nu_{e}$ and $\bar{\nu}_{e}$ become closer and the number densities of $\nu_{e}$ and $\bar{\nu}_{e}$ become comparable. Hence, ELN crossings are favored. Note that we investigated the appearance of crossings by using inputs from 5 time snapshots between $0.125$ s and $1$ s in our SN model [$t_{\mathrm{p.b.}}=(0.125,0.25,0.5,0.75,1)$ s], but we only show results for three of them for simplicity. The first post-bounce time for which we find ELN crossings among the analyzed ones is $0.5$ s. ELN crossings appear for any $t_{\mathrm{p.b.}}\gtrsim 0.5$ s of the ones studied.

The top panel of Fig. 7 shows the radial evolution of ELN$(\cos\theta=1)$ for $t_{\mathrm{p.b.}}=0.25,0.5$ and $1$ s. As discussed for the toy-model, ELN crossings only occur in the proximity of the neutrino free-streaming region. In correspondence of the appearance of crossings, we expect that ELN$(\cos\theta=1)$ changes sign; this is what we find for $t_{\mathrm{p.b.}}=0.5$ and $1$ s.

The total number of neutrinos and antineutrinos remains unchanged after decoupling as shown in the bottom panel of Fig. 7.

In conclusion, our stationary and spherically symmetric SN model strongly suggests that ELN crossings can only occur within the spatial regions where neutrinos and antineutrinos decouple and start to free stream. Moreover, a steep drop of the baryon density profile, as typical of the late SN stages, together with comparable number densities of $\nu_{e}$ and $\bar{\nu}_{e}$ and $\lambda_{\nu_{e}}/\lambda_{\bar{\nu}_{e}}\leq 1$ in the decoupling region, would favor the development of angular distributions that are similar to each other; as a consequence, ELN crossings develop given the different interaction rates of $\nu_{e}$ and $\bar{\nu}_{e}$. Our results confirm the findings of Refs. (Tamborra et al., 2017; Delfan Azari et al., 2019) where only the early SN stages were analyzed and no crossing was found. Notably, hydrodynamical simulations investigating the behavior of the neutrino angular distributions presented in the literature do not exhibit successful explosions or have been artificially exploded; as a consequence, none of them investigates the behavior of the neutrino angular distributions in the late accretion and cooling phases. Our model is the first to predict the occurrence of crossings at later post-bounce times, when the baryon density profile becomes steeper.

Our stationary and spherically symmetric SN model does not include any macroscopic asymmetry eventually induced by hydrodynamical instabilities. The occurrence of the LESA instability has been pinpointed as a possible favorable condition leading to ELN crossings (Izaguirre et al., 2017; Dasgupta et al., 2017). Our stationary SN model hints that, in the SN angular regions of net $\nu_{e}$ (or $\bar{\nu}_{e}$) excess, our criteria should still hold. Hence, we expect that ELN crossings should occur only across the SN angular regions where the ELN changes sign because of LESA (see, e.g., Fig. 1 of Tamborra et al. (2014a)). On the other hand, if the LESA instability is sustained until the late SN phases, then ELN crossings may develop even in regions of net $\nu_{e}$ (or $\bar{\nu}_{e}$) excess, as predicted for the cases with $t_{\mathrm{p.b.}}\geq 0.5$ s in our SN model. We stress, however, that our model cannot fully test the conditions leading to crossings in the presence of the LESA instability self-consistently, since this would require a break of the spherical symmetry. Therefore, our conjectures remain to be tested thorough self-consistent 3D hydrodynamical simulations and will be subject of further work.

4. Dynamical generation of crossings in the electron neutrino lepton number distribution

We have seen in the previous Sections that the presence of ELN crossings, or lack thereof, is determined by the radial profile of $\lambda_{\nu_{e}}/\lambda_{\bar{\nu}_{e}}$, the baryon density profile, as well as $n_{\nu_{e}}/n_{\bar{\nu}_{e}}$, which is related to the former two, in the trapping region. However, as discussed in Sec. 3, the presence of ELN crossings is not only determined by the local conditions, but it indirectly feels the effect of distant regions in the SN through collisions. In this Section, we intend to explore whether ELN crossings can be generated dynamically from flavor instabilities occurring in the SN because of local fluctuations (Abbar & Duan, 2015; Dasgupta & Mirizzi, 2015; Capozzi et al., 2016). According to this scenario, the neutrino angular distributions may be modified dynamically; as a consequence, ELN crossings may form and, in turn, they could foster the growth of fast flavor instabilities.

It should be noted that the neutrino interaction rate in matter is dominated by the neutrino-nucleon cross-section, and the ratio $\lambda_{\nu_{e}}/\lambda_{\bar{\nu}_{e}}$ cannot be changed significantly by the neutrino conversions on a global scale. In the deep interior of a SN, the effective average energy of $\nu_{e}$ is larger than that of $\bar{\nu}_{e}$ due to the non-null chemical potential. This contributes to a larger $\nu_{e}$ interaction rate. Neutrino conversions could in principle reduce the average energy of $\nu_{e}$ and bring it closer the one of $\bar{\nu}_{e}$.

We take the inputs of the SN snapshot at $t_{\mathrm{p.b.}}=0.25$ s used in Sec. 3.2 as benchmark case. This case did not exhibit ELN crossings in the absence of pre-existing flavor conversions. We then impose that at a certain radius $r_{\star}$, with $r_{\star}\in[r_{\mathrm{min}}=5~{}\mathrm{km},r_{\mathrm{max}}=60~{}\mathrm{km}]$, flavor conversions are triggered possibly leading to flavor decoherence; the latter is one of the most extreme scenarios that one could expect. This effect can be mimicked by assuming that $\lambda_{\nu_{e}}/\lambda_{\bar{\nu}_{e}}\rightarrow 1$ for $r\geq r_{\star}$. We then change $r_{\star}\in[r_{\mathrm{min}},r_{\mathrm{max}}]$ in order to test whether the radius of the onset of flavor conversions would affect the dynamical development of ELN crossings.

In none of the studied cases, we find a significant modification of the ELN evolution and ELN crossings do not develop within our simplified SN setup. It is worth noticing that the specific choice of the numerical values adopted for $r_{\mathrm{min}}$ and $r_{\mathrm{max}}$ does not affect our conclusions as long as both radii are chosen to be far away from the decoupling region.

We then conclude that, in our model, the neutrino angular distributions cannot be dynamically modified by flavor instabilities due to local fluctuation and induce ELN crossings dynamically. Although a more thorough analysis of the problem is required, we find that the presence of fast oscillations at radii much larger than the decoupling radius, cannot significantly change conditions in the region of decoupling to facilitate the dynamical generation of ELN crossings in a geometry that is spherically symmetric.

5. Discussion and Conclusions

The development of crossings between the angular distributions of $\nu_{e}$ and $\bar{\nu}_{e}$ (ELN crossings) is of relevance because it can possibly lead to fast pairwise conversions of neutrinos deep in the SN core with major consequences on the SN physics. It is vital to gain a qualitative understanding of this phenomenon. In this paper, using a simple yet insightful technique, we have qualitatively addressed this question focusing on the microphysics of neutrino-matter collisions in the SN core.

The highly non-linear nature of the neutrino flavor evolution along with the feedback on the flavor dynamics coming from neutrino-matter collisions makes a general self-consistent analysis impossible within current means. However, we have aimed to provide a rule of thumb for the occurrence of ELN crossings. To that purpose, we have constructed a simplified stationary and spherically symmetric SN model that takes into account the physics of collisions through an iterative approach, but neglects any asymmetry and further complications coming from the SN hydrodynamical instabilities. It is should be noted that our assumption of spherical symmetry may be substantially broken either due to SN hydrodynamics (Tamborra et al., 2014b; Janka et al., 2016), or because the nature of the neutrino flavor evolution (Duan & Shalgar, 2015; Mirizzi et al., 2015; Abbar et al., 2015).

We have shown that the conditions affecting the development of ELN crossings are not local in nature. In particular, the appearance of ELN crossings is determined by the slope of the baryon density profile together with the requirement that the $\nu_{e}$ and $\bar{\nu}_{e}$ number densities are comparable in the proximity of the decoupling region. Our simple spherically symmetric SN model hints that ELN crossings can only occur in the late stages of the accretion phase and in the cooling phase, under the assumption that a stationary configuration is reached. In fact, at earlier post-bounce times, a large excess of the $\nu_{e}$ number density over the $\bar{\nu}_{e}$ one prevents ELN crossings from occurring. The latter effect is determined by a baryon density profile that slowly varies with the radius, disfavoring the $\nu_{e}$ and $\bar{\nu}_{e}$ distributions from becoming similar. However, in the late accretion phase and cooling phase the distributions of $\nu_{e}$ and $\bar{\nu}_{e}$ naturally become more similar to each other, $\nu_{e}$ and $\bar{\nu}_{e}$ decouple in closer spatial regions, and favorable conditions for ELN crossings arise.

Due to the numerical challenges involved in solving the equations of motion that include neutrino-neutrino interactions, most of the focus has been on the linear stability analysis of the conditions under which instabilities in the flavor space can occur. However, if flavor instabilities are triggered in a small localized region of space, is it not clear if and under which conditions the flavor instability would spread, see e.g. Capozzi et al. (2017); Yi et al. (2019). One aspect of the question is whether the flavor evolution changes the neutrino interaction rates therefore leading to a dynamical generation of ELN crossings. Our stationary and spherically symmetric SN model suggests that ELN crossings cannot be generated dynamically, unless favorable conditions already exists in the SN core.

Our model neglects perturbations coming from global asymmetries induced by hydrodynamic instabilities occurring in SNe. However, it still provides good insights on the generation of ELN crossings under stationary conditions.

A concrete list of necessary and sufficient conditions under which fast pairwise conversions of neutrinos can occur in the SN core still remains unsettled. Our work provides new insights on the solution of this intriguing jigsaw.

We acknowledge insightful discussions with Thomas Janka, Georg Raffelt and Anna Suliga, and are grateful to Robert Bollig for granting access to the data of the $18.6M_{\odot}$ SN model adopted in this work. SS and IT acknowledge support from the Villum Foundation (Project No. 13164). The work of IT has also been supported by the Knud Højgaard Foundation and the Deutsche Forschungsgemeinschaft through Sonderforschungsbereich SFB 1258 “Neutrinos and Dark Matter in Astro- and Particle Physics (NDM).

In order to further test the validity of our scheme, in this Appendix we adopt the uniform sphere problem as a test case. We consider a sphere emitting uniform black-body radiation and compare the findings of our simplified numerical scheme with the analytical results reported in Sec. 4.2 of Murchikova et al. (2017). We expect to find a constant intensity across the surface with limb darkening at the edge of the sphere due to the effects of finite optical depth. This is the case of neutrinos emitted from a certain flavor-dependent neutrinosphere.

To reproduce the uniform sphere case, we assume that electron neutrinos are emitted at a uniform rate inside a sphere of radius $R_{\nu_{e}}$ and electron anti-neutrinos are emitted within a sphere of radius of $R_{\bar{\nu}_{e}}$. As neutrinos and anti-neutrinos propagate outwards they are absorbed at a constant rate, $\kappa_{\nu_{e},\bar{\nu}_{e}}$, inside $R_{\nu_{e}}$ and $R_{\bar{\nu}_{e}}$ respectively. Under these assumptions the total number of neutrinos emitted by the neutrinosphere is

[TABLE]

while the angular distribution at a distance $r$ from the center of the neutrino-sphere is (Murchikova et al., 2017)

[TABLE]

where $B\simeq R_{\nu_{e},\bar{\nu}_{e}}^{2}/T^{3}$, and

[TABLE]

The continuous lines shown in Fig. 8 represent the angular distributions obtained from Eq. 2 with temperatures of $\nu_{e}$ and $\bar{\nu}_{e}$ set to $10.0$ and $11.0$ MeV, respectively, and $\kappa_{\nu_{e}}=0.6$ km*-1*, $\kappa_{\bar{\nu}_{e}}=0.2$ km*-1*. For the top panel, we set $R_{\bar{\nu}_{e}}=20$ km and $R_{{\nu}_{e}}=22$ km; while for the case in the bottom panel, we adopt $R_{\bar{\nu}_{e}}=20$ km and $R_{{\nu}_{e}}=25$ km.

We now want to reproduce the same result by adopting the numerical iterative scheme introduced in this paper in order to test the robustness of our method. Note that, while in the rest of the paper we assume that the total number of particles is conserved within the sphere of radius $r_{\mathrm{max}}$, to reproduce the uniform sphere problem we need to take a constant $\mathcal{C}^{\mathrm{gain}}$ term in Eq. 1 which does not depend on $\rho$ and $\bar{\rho}$ for $r\leq R_{{\nu}_{e},\bar{\nu}_{e}}$ and it is zero otherwise. The $\mathcal{C}^{\mathrm{loss}}$ term stays unchanged for $r\leq R_{{\nu}_{e},\bar{\nu}_{e}}$ and it is zero at larger radii. The absorption coefficient introduced above corresponds to $\lambda_{\nu_{e},\bar{\nu}_{e}}\equiv 2\kappa_{\nu_{e},\bar{\nu}_{e}}$.

The dashed lines in Fig. 8 show the correspondent angular distributions obtained by solving Eq. 1 through our iterative scheme with the same inputs adopted for the analytical model described above. One can see that our numerical scheme reproduces the analytical result with less than $1\%$ error. Moreover, a comparison between the top and the bottom panel suggests that we obtain an ELN crossing when the two decoupling radii do not differ much from each other; the ELN crossing disappears when the ratio of the two decoupling radii is significantly larger than unity.