Inspecting neutrino flavor instabilities during proto-neutron star cooling phase in supernova:
I. Spherically symmetric model

Core-collapse supernovae (CCSNe) are dramatic events at the final stage of the stellar evolution for stars with a zero-age main sequence mass of $\gtrsim 10\,M_{\odot}$. When the central core experiences a core collapse and forms a proto-neutron star (PNS), a huge amount of neutrinos are released, which deleptonize and cools down the CCSN core. Meanwhile, the shock wave is launched by the core bounce and blows accreting matter off. Although the shock wave undergoes stagnation or even retraction during its post-bounce evolution, neutrino transport from PNS to the post-shock region and their interactions with matter behind the shock wave via various weak processes (see, e.g., Langanke & Martínez-Pinedo 2003) can reinvigorate the shock expansion with the aid of multi-dimensional fluid instabilities. This seems to be a primary driving force of CCSNe for a wide range of progenitor masses (see, e.g., Burrows et al. 2020), which is well known as neutrino-heating mechanism (Janka, 2012, 2017; Mezzacappa et al., 2020; Burrows & Vartanyan, 2021; Fischer et al., 2024; Yamada et al., 2024).

In dense neutrino environments such as CCSN cores, refractive effects by neutrino self-interactions, which have been neglected in the standard model of CCSNe, can induce large flavor conversions (Sigl & Raffelt, 1993; Duan et al., 2006; Tamborra & Shalgar, 2021; Capozzi & Saviano, 2022; Richers & Sen, 2022; Volpe, 2024) Among them is the fast flavor instability (FFI) (Sawyer, 2005, 2016), which is triggered by a zero crossing in ELN (electron neutrino-lepton number)-XLN (heavy-lepton neutrino lepton number) angular distribution. Yet another is the collisional flavor instability (CFI) (Johns, 2023), which is induced by the disparity in reaction rates between neutrinos and antineutrinos. Both have attracted significant attention recently due to the potential for instability under the supernova shock. The appearance of ELN-XLN angular crossing has been explored in state-of-the-art CCSN models (Dasgupta et al., 2017; Abbar et al., 2020; Nagakura et al., 2021a; Akaho et al., 2024a), and the occurrence of CFI has been also confirmed similarly using linear stability analysis (Liu et al., 2023a; Akaho et al., 2024a). These surveys for flavor instabilities imply that flavor conversion ubiquitously occurs in the post-shock region. Recent studies also suggest that the neutrino radiation field is significantly changed by flavor instabilities, which have an impact on many processes and outcomes in CCSNe such as fluid dynamics (Nagakura, 2023; Ehring et al., 2023a, b; Xiong et al., 2023b, 2024; Shalgar & Tamborra, 2024), neutron star kick (Nagakura & Sumiyoshi, 2024), and neutrino signal (Wu et al., 2015; Nagakura & Zaizen, 2023; Ehring et al., 2023a; Xiong et al., 2024).

These recent results of flavor conversions stimulate our interest in the post-explosion phase, as neutrinos continue to play important roles in cooling the remnant PNS and determining nucleosynthesis in the ejecta. As suggested by earlier studies (see, e.g., Takahashi et al. 1994; Qian & Woosley 1996), heavy nuclei including (weak) r-process and $\nu$p-process elements could be synthesized in neutrino-driven winds if such winds form (e.g., Witt et al. 2021; Wang & Burrows 2024; Nevins & Roberts 2024). Nucleosynthesis has, however, a delicate sensitivity to the relative difference between electron-type neutrinos ($\nu_{e}$) and their anti-partners ($\bar{\nu}_{e}$) (see, e.g., a recent review of Fischer et al. 2020a and references therein), indicating that it requires accurate treatments of neutrino kinetics that can resolve flavor-dependent features. This suggests that the occurrence of flavor conversions, if any, could be one of the largest uncertainties in these models. PNS convection is expected to be present and generate flavor instabilities in the PNS (e.g., Glas et al. 2020; Abbar et al. 2020; Delfan Azari et al. 2020). Although parameterized models of flavor transformation during PNS cooling cause the amount of ejected mass to increase and more proton-rich conditions, the net effects of flavor instability in self-consistent models are unknown. This study aims to build on these works by determining whether the conditions for the FFI arise in realistic neutrino distributions during the PNS cooling phase.

In this paper, we present the first detailed study of the occurrence of FFI and CFI in the PNS cooling phase. We determine the PNS structure and its long-term evolution by our PNS cooling scheme, in which we employ a quasi-static approach and diffusion approximation for neutrino transport under spherical symmetry. For stability analysis of flavor instabilities, multi-angle neutrino transport simulations are carried out for some selected time snapshots. We then survey FFI and CFI in the neutrino data using well-established approximate methods to identify these instabilities.

This paper is organized as follows. We start with introducing essential methodologies to develop models of PNS structure and its long-term evolution in Sec. 2. In Sec. 3, we describe how to model neutrino radiation field under frozen fluid- and metric backgrounds by using a Monte Carlo calculation of neutrino transport. Briefly summarizing the linear stability analysis of flavor conversion in Sec. 4, we show the results in Sec. 5. Finally, in Sec. 6, we conclude the present study and describe the need for further study towards drawing more robust conclusions.

During cooling stages of the PNS, it is nearly in hydrostatic equilibrium and evolves quasistatically while emitting neutrinos. In this study, we model the PNS structure at each instant under a condition that the general relativistic spacetime is stationary and spherically symmetric. The static matter distribution can be, hence, determined by solving Tolman-Oppenheimer-Volkoff (TOV) equation (e.g., Burrows & Lattimer 1986; Pons et al. 1999).

In most of the PNS interior except for the narrow outer envelope, neutrinos interact with the matter on short length scales compared to the system size, and so are effectively trapped. For this reason, their transport of energy and lepton number in space can be approximated as a diffusion processes. Since neutrino interaction rates depend sensitively on neutrino energy and flavor, we employ a multi-group (i.e., multi-energy) flux limited diffusion (MGFLD) scheme to determine the neutrino radiation field and its interaction with matter in our PNS cooling model. We note that a flux limiter is adopted in the region where the diffusion approximation is invalid so as to guarantee that neutrino transport is causal and become accurate in the free-streaming limit. We refer to Appendix A for more details of our scheme.

The neutrino-matter interactions in our MGFLD scheme are computed based on the same approach as that in Bruenn 1985, but with some extensions. We employ the method in Suzuki & Nakamura 1993 for the nucleon bremsstrahlung processes which play important roles for production of non-electron type neutrinos during the PNS cooling. Below we provide a list of weak processes implemented in our code:

where $p,n,A,N,e$ are protons, neutrons, nuclei, nucleons, and electrons, respectively. The plasmon decay process into neutrino pairs is also added to the electron-positron pair annihilation processes by scaling the total neutrino emissivity as calculated by Kohyama et al. 1986.

As the initial condition of the PNS, we take a fluid snapshot of a CCSN model simulated by Nakazato et al. 2018. The CCSN model was developed based on a spherically symmetric general relativistic neutrino-radiation hydrodynamic simulation of a 15 $M_{\odot}$ progenitor model in Woosley & Weaver 1995. The CCSN model employed the Togashi equation-of-state (EOS) (Togashi et al., 2017), which is a nuclear EOS constructed by the variational many-body theory with realistic nuclear force potentials. Following the same procedure as in Nakazato et al. 2018, we remap a fluid profile (i.e.,, entropy and electron fraction profiles as functions of baryon mass) from the center to a radius just behind shock wave at 300ms after core bounce. The corresponding baryon mass covered by the spatial grid is 1.4675$M_{\odot}$, which also represents the baryon mass of the PNS. We construct the hydrostatic PNS structure with steady neutrino flow from this initial condition, and then simulate the following long-term evolution of the PNS. Just as in the CCSN model, we assume that $\mu$-, $\tau$ flavors, and their antipartners are identical. This would be a reasonable condition unless on-shell muons appear (Bollig et al., 2017; Fischer et al., 2020b). We follow the post-explosion evolution of PNS up to 50 seconds from the remapped initial condition; the result is summarized in Fig. 1. In the figure, we display the radial profiles of baryon mass density ($\rho$), electron fraction ($Y_{e}$), temperature ($T$), and lapse function ($\alpha$) at selected time snapshots for our PNS cooling model.

MGFLD transport approximately calculates the flow of energy and lepton number, but does not provide information about the angular structure of the radiation fields needed to determine whether neutrino flavor instabilities present. Based on the PNS model, we carry out another neutrino transport simulation for stability analysis of flavor conversion. The detail is described in the next section.

In general, linear stability analysis of neutrino flavor conversion requires full angular information of neutrino distributions in momentum space (see Sec. 4 for more details). We, hence, carry out another multi-angle neutrino transport simulation on top of the PNS fluid background.

We calculate the equilibrium neutrino radiation field on each snapshot using a general relativistic Monte Carlo neutrino transport code Sedonu (Richers et al., 2015; Richers, 2022). We assume spherical symmetry and discretize the domain in radius using the same grid present in the MGFLD transport model. We also map the same metric, matter density, temperature, and equilibrium neutrino chemical potentials from the PNS cooling model. At each radial grid cell, the neutrino distribution is discretized in polar angle with 100 zones uniformly distributed over $-1\leq\cos\theta_{\nu}\leq 1$. We use the SFHo EOS (Steiner et al., 2013) and neutrino interaction rates calculated using NuLib (O’Connor, 2015). We treat charged-current absorption and emission from free nucleons and nuclei and elastic isotropic scattering from nucleons and nuclei. Emission from neutrino pair production and nucleon-nucleon bremsstrahlung radiation is included, and the inverse reactions are treated with an effective absorption opacity calculated by assuming detailed balance. The opacities in NuLib are largely the same as those in Bruenn 1985, but include corrections for weak magnetism and recoil, form-factor corrections due to decoherence, electron polarization corrections, and ion-ion correlations (Horowitz, 1997, 2002; Burrows et al., 2006). We employ 18 energy bins with upper bounds logarithmically spaced from 2 to 311 MeV. We treat regions with large scattering optical depths with the Monte Carlo random walk method to accelerate diffusive transport (Fleck & Canfield, 1984).

We simulate $2.1-4.0\times 10^{10}$ Monte Carlo packets in each snapshot, each representing a number of neutrinos proportional to the total luminosity of the grid cell they are created in. The particle positions and four-momenta move along geodesics, while the number of neutrinos contained in each packet decreases exponentially according to the absorption rate. Monte Carlo packets change the direction after a random interval according to the scattering opacity, and are assumed to scatter isotropically in the frame comoving with the fluid. As particles propagate, they contribute to the local distribution function. Particles are destroyed when they leave the domain or pass within a radius of $10\,\mathrm{km}$ (although for the $30\,\mathrm{ms}$ and $50\,\mathrm{ms}$ snapshots we necessarily remove this condition to avoid artifacts associated with the finite inner boundary). Particles are rouletted if their weight decreases below a factor of $10^{-3}$ of their original weight. Once all particles are destroyed, the resulting distribution is read out and then used for stability analysis of neutrino flavor conversions.

It is important to note that neither the EOS nor the neutrino interaction rates are consistent with those used in the simulation of the PNS evolution. Although a more exact match would be desirable, we take a more expedient first approach to understand qualitative features of the neutrino distributions. We choose the baryon mass density, temperature, and chemical potential of $\nu_{e}$ ($\mu_{p}+\mu_{e}-\mu_{n}$) as the primitive thermodynamic variables from which to reconstruct reaction rates, and allow $Y_{e}$ to float to accommodate this choice. As a result, the Monte Carlo calculation has the same equilibrium distributions as in the MGFLD calculation, even if the interaction rates and fluxes differ. In practice, this choice also results in the similar reaction rates for weak processes. The chemical potentials of electron neutrinos and the change in $Y_{e}$ required are shown in Fig. 2.

In Fig. 3, we compare the energy density of neutrinos computed using MGFLD and Sedonu. Their distributions are nearly identical in optically thick region by construction, since the equilibrium chemical potentials are designed to match between the two methods. However, they gradually deviate with increasing radius in the outer envelope of PNS. The deviation is larger for $\bar{\nu}_{e}$ than $\nu_{e}$ due to the lower opacity of $\bar{\nu}_{e}$ in the more neutron rich environment of the MGFLD calculation. One thing we do notice here is that the deviation is not only due to EOS and weak interactions but also the difference of two transport schemes. Despite both the different transport schemes and the different opacities, the differences between the distributions of neutrino energy density are minor and the two methods are in overall good agreement with each other. The Monte Carlo calculation, however, produces a realistic angular distribution that results from self-consistent diffusion and advection of radiation away from the PNS.

By using the neutrino data obtained from Sedonu simulations, we perform linear stability analysis with a dispersion relation approach. Below, we provide essential details of the analysis but we refer readers to Izaguirre et al. 2017; Liu et al. 2023b; Akaho et al. 2024a for more details. Throughout this work we do not simulate neutrino quantum kinetics, but rather use the results of our classical kinetics calculations to ascertain the presence of quantum kinetic instabilities.

The fundamental quantity in neutrino quantum kinetics is neutrino density matrix ($\rho$)¹¹1The reader should not confuse $\rho$ with the baryon mass density., and the neutrino transport equation can be expressed in terms of $\rho$. In our stability analysis, we work with the quantum kinetic equation (QKE) derived under the mean-field approximation, which can be described as

where the flavor-transforming Hamiltonian is defined in an orthonormal tetrad and is composed of vacuum, matter, and self-interaction terms:

In the first term, $U$ and $M^{2}$ denote the Pontecorvo-Maki-Nakagawa-Sakata matrix and the mass-squared matrix for neutrino mass eigenstates, respectively. The second term induces matter oscillation, where $v^{\mu}=(1,\bm{v})$ is the neutrino’s unit four velocity and $\Lambda_{\mu}=\sqrt{2}G_{\mathrm{F}}\,\mathrm{diag}[\{j_{\mu}^{l}\}]$ is the charged lepton current $j_{\mu}^{l}$ of charged lepton $l$. The third part corresponds to the neutrino self-interaction, and the integral over phase space is defined by

We employ the flavor-isospin convention, i.e., $\bar{\rho}(E)\equiv-\rho(-E)$ for the anti-neutrino density matrix. Regarding the collision term (the second term in the right hand side of Eq. (1)), we employ the same approach as in Liu et al. 2023a; Akaho et al. 2024a, in which the emission- and absorption processes are treated in the following form;

where $j$ and $\kappa$ correspond to emissivity and absorptivity, respectively, but we ignore scattering processes for expediency. In the expression, the subscripts $a$ and $b$ denote a neutrino flavor, and the bracket is defined by

It is worth noting that pair processes are incorporated in our stability analysis under the form of Eq. (4). We compute their effective emissivity $j_{\mathrm{pair}}$ and absorptivity $\kappa_{\mathrm{pair}}$ by carrying out momentum-integration of neutrino distributions with their reaction kernels (see also Liu et al. 2023a; Akaho et al. 2024a) based on the opacities used in the Monte Carlo calculations. As a result, the collision term for heavy-leptonic neutrinos is non-zero in our stability analysis, despite no charged-current reactions. This is in line with the approximations made in both the MGFLD and Monte Carlo calculations.

The resultant QKE for the off-diagonal component of the density matrix can be written as

where emission, absorption, and pair process are summarized as

Here, the vacuum and matter terms are omitted in this equation because they do not affect FFI or CFI. We linearize Eq. (6) under the condition that neutrinos are nearly in flavor eigenstates. We then derive the dispersion relation under the assumption that the wavelength of unstable modes, if any, are sufficiently smaller than the scale height of background neutrinos. By adopting a plane-wave ansatz, $\rho_{ex}\propto\tilde{Q}\exp[-\mathrm{i}k_{\mu}x^{\mu}]$ where $k_{\mu}=(\omega,\bm{k})$, the resultant dispersion relation can be written as,

where

If there is a positive imaginary value of $\omega$ that satisfies this condition, the neutrino distribution is unstable to flavor conversions.

Although it is straightforward to search for positive imaginary values of $\omega$ by solving Eq. (8), we take an alternative way in the present study. The major reason to avoid the direct search is its computational cost. As is well known, spurious modes, corresponding to unphysical eigenmodes, appear in the stability analysis when we handle momentum-integration in Eq. (9) by using discretized angular distributions of neutrinos (Sarikas et al., 2012; Morinaga & Yamada, 2018). We can suppress the artifact by increasing resolutions of neutrino momentum space, but neutrino distributions need to be scanned at each spatial point for multiple snapshots in the present study, exhibiting that direct method is not appropriate. Below, we describe alternative way to identify FFI and CFI by following the procedure as in Akaho et al. 2024a.

FFI is dictated by anisotropy of ELN and XLN. Under the assumption that all heavy-leptonic neutrinos ($\mu$, $\tau$, and their antipartners) are identical, it can be determined only by ELN angular distribution, which is defined as,

where $f_{\nu_{a}}$ is a neutrino occupation distribution for a flavor $a$, which corresponds to the diagonal component in neutrino density matrix. As proven by Morinaga 2022, its zero crossings (or angular crossings) give a sufficient and necessary condition for FFI. For this reason, we inspect ELN angular crossings in neutrino radiation fields, and then judge whether FFI can occur.

For CFI, an efficient approximate method was developed by Liu et al. 2023b, in which the growth rate of $\bm{k}=\bm{0}$ (or homogeneous) mode can be analytically estimated under an assumption that neutrino angular distributions are isotropic. We note that homogeneous mode usually provides the dominant unstable mode (Liu et al., 2023b), and that the neutrino distributions inside the PNS are nearly isotropic. It is also worth noting that the collision term of iso-energetic scattering processes become exactly zero for isotropic neutrino distributions, which also supports our neglect of scattering processes in the stability analysis.

Under these assumptions, two analytical solutions for oscillation frequency $\omega$, so-called isotropy-preserving and isotropy-breaking branches, can be given as,

and

Here each quantity is defined as the following notation:

where $\mathfrak{g}=\sqrt{2}G_{\mathrm{F}}(n_{\nu_{e}}-n_{\nu_{x}})$. The neutrino number density and mean collision rates are given by

We are primarily interested in determining whether Eqs. (11) and (12) have positive imaginary parts. Following Liu et al. 2023b, we can approximate Eqs. (11) and  (12) by taking a limit of $A$ in its relation to $\left\lvert G\alpha\right\rvert$, which can be written as,

for isotropy-preserving mode and

for isotropy-breaking mode.

Eqs. (15) and (16) exhibit that the second term in the right hand side of these equations needs to overwhelm the first one ($-\gamma$) to make the system unstable. This indicates that large neutrino-self interactions are necessary for the occurrences of CFI. Another important remark is that, in realistic CCSN environments, the former condition ($A^{2}\gg\left\lvert G\alpha\right\rvert$) is usually satisfied. On the other hand, the latter case ($A^{2}\ll\left\lvert G\alpha\right\rvert$) can emerge in regions with $A\sim 0$, indicating that so-called resonance-like CFI occurs (Xiong et al., 2023a; Liu et al., 2023b). In fact, a detailed inspection of the neutrino radiation field based on multi-dimensional CCSN simulations suggest that such a resonance-like CFI can appear in the envelope of the PNS (Akaho et al., 2024a). Since the growth of flavor instability could be very fast and its non-linear evolution would lead to flavor swap (Kato et al., 2024), the CFI may play an important role on CCSN and PNS dynamics, albeit very narrow spatial region.

In this study, we find no positive sign of neutrino-flavor instabilities for either FFI or CFI in any time snapshot. Below, we discuss the physical reasons and provide some figures displaying key quantities. It should be noted that we focus on the spatial region at $r>10\,\mathrm{km}$ hereafter. This is because in very optically thick regions ($r<10\,\mathrm{km}$), all neutrinos are in thermal equilibrium, and both FFI and CFI can not occur in these regions²²2This statement would be understandable for FFI, since it is driven by anisotropic angular distributions, indicating that it can not occur in isotropic distributions as in thermal equilibrium states. For CFIs, the argument was given in our previous papers (Liu et al., 2023a; Akaho et al., 2024a), while we shall provide a more formal mathematical proof in another paper (Liu et al., 2024).. We confirmed that $r=10\,\mathrm{km}$ corresponds to a safe radius where all neutrinos are in thermal equilibrium states at all time snapshots (see, e.g., Fig. 3).

We begin with analyses of FFI. Figure 4 displays the angular distributions for each species as one of the representative examples in our results. The plotted distributions are measured at the outermost cell at $1\,\mathrm{s}$ time snapshot in our PNS cooling model, and there is some noise due to the random nature of the Monte Carlo calculations. As shown in the figure, the $\nu_{e}$-distribution substantially dominates over those of the other species for all angles, indicating that there are no ELN angular crossings. Although the particulars of the distribution vary with location and time, the same qualitative behavior is the same at all other locations and times.

That generation of ELN crossings is difficult in this environment can be also understood from the radial profiles for the ratio of the number density of $\bar{\nu}_{e}$ to $\nu_{e}$, shown in Fig. 5. In general, since the PNS is a neutron-rich environment, $\bar{\nu}_{e}$ is decoupled at smaller radii than $\nu_{e}$, resulting in more forward-peak angular distributions. The process can lead to ELN angular crossing, but the large disparity in the number density between $\nu_{e}$ and $\bar{\nu}_{e}$ implies that a large disparity in angular distributions is required to make up for the disparity in number density to create crossings (Nagakura et al., 2019; Abbar et al., 2019, 2020). This never happens in our calculations.

Let us turn our attention to the CFI. The rationale behind the stable neutrino radiation field with respect to CFI is interpretable with analytical formulae (see Eqs. (15) and (16)). The condition for unstable modes is that the second term $\left\lvert G\alpha\right\rvert/\left\lvert A\right\rvert$ overwhelms the first one, $\gamma$, in Eqs. (15) and  (16) for non-resonant cases ($A^{2}\gg\left\lvert G\alpha\right\rvert$)³³3We also note that there are no regions with $A^{2}\ll\left\lvert G\alpha\right\rvert$ in this study. This exhibits that it does not meet a condition for resonance-like CFI.. The condition can be simplified into $\left\lvert G\right\rvert/\left\lvert A\right\rvert>\gamma/\left\lvert\alpha\right\rvert$. Since $\gamma$ is larger than $\left\lvert\alpha\right\rvert$, we further simplify the inequality as $\left\lvert G\right\rvert/\left\lvert A\right\rvert>1$, which corresponds to a necessary (though not sufficient) condition for the occurrences of CFIs. In Fig. 6, we portray the radial profiles of $G$ and $A$ in the top panel, while those for $\left\lvert G\right\rvert/\left\lvert A\right\rvert$ and $\gamma/\left\lvert\alpha\right\rvert$ are displayed in the bottom. As shown in the bottom panel, $\left\lvert G\right\rvert/\left\lvert A\right\rvert$ is always lower than unity, exhibiting that no CFI occurs in the entire spatial region at all time snapshots.

It is worth noting that the small $G$ is attributed to the large populations of $n_{\nu_{x}}$. In fact, the inequality $\left\lvert G\right\rvert/\left\lvert A\right\rvert>1$ is satisfied when both $\nu_{e}$ and $\bar{\nu}_{e}$ are larger than $\nu_{x}$⁴⁴4We note that the inequality is also satisfied in cases where $\nu_{x}$ is much larger than $\nu_{e}$ and $\bar{\nu}_{e}$. However, such a situation does not arise in CCSN environments.. It means that $\nu_{x}$ suppresses the appearance of CFI, as reported by our previous paper (Liu et al., 2023a). To see the trend clearly, we provide Fig. 7, in which the radial profiles of neutrino number density for each species are portrayed at two time snapshots $t=3$ and $10\,\mathrm{s}$, as examples. As shown in the figure, the $\nu_{x}$ number densities always dominate over those of $\bar{\nu}_{e}$, which results in $\left\lvert G\right\rvert/\left\lvert A\right\rvert<1$. To understand the trend of relatively low densities of $\bar{\nu}_{e}$ compared to other species of neutrinos, we depict the radial profiles of degeneracy of electron neutrinos, which is defined as $\eta_{\nu_{e}}\equiv\mu_{\nu_{e}}/T$, and also $Y_{e}$ in Fig. 8. The positive value of $\eta$ suggests that Fermi degeneracy of $\nu_{e}$ accounts for the suppression of $\bar{\nu}_{e}$.

We remark that the peak profile of $\eta_{\nu_{e}}$ around $r=12\,\mathrm{km}$ at $t=10\,\mathrm{s}$ might be a numerical artifact in the PNS cooling model. In fact, the diffusion approximation for neutrino transport is not reasonable in the region. It should be also mentioned that the matter distributions around the surface of PNS is sensitive to fall-back or long-lasting accretion flows onto the PNS (Nagakura et al., 2021b; Akaho et al., 2024b). We should keep in mind these uncertainties, and the result might be changed in more self-consistent PNS cooling models.

We provide another piece of evidence that $\nu_{x}$ is responsible for suppressing CFI. Figure 9 shows the growth rates and the ratios without the contributions from $\nu_{x}$ (more specifically, we set $n_{\nu_{x}}=n_{\bar{\nu}_{x}}=0$ in computing $G$ and $A$). As shown clearly in the top panel, positive growth rates appear. Also, in the bottom panel of the figure, $G$ becomes substantially larger than that in the cases with $\nu_{x}$ (see also Fig. 6). On the other hand, $A$ remains the same as the case with $\nu_{x}$. This is simply because $A$ does not depend on $\nu_{x}$, as long as $\nu_{x}=\bar{\nu}_{x}$ is satisfied.

Our result shows the importance of accurate modeling of $\nu_{x}$ radiation fields to assess the occurrence of CFIs. It should be mentioned that there are large uncertainties in weak processes of $\nu_{x}$ even for modern CCSN models, which are, for instance, many-body corrections in both neutral- and charged-current reactions (Burrows & Sawyer, 1998, 1999; Horowitz et al., 2017), nucleon-nucleon bremsstrahlung (Friman & Maxwell, 1979; Bartl et al., 2016; Guo & Martínez-Pinedo, 2019), and on-shell muons (Fischer et al., 2020b; Guo et al., 2020; Sugiura et al., 2022), etc. We, thus, encourage careful consideration of how these rates can have an influence on CFIs, whose detailed investigations are deferred to future work.

In this paper, we investigate the possibilities of flavor instabilities during the PNS cooling phase. To this end, we construct a spherically symmetric PNS model by employing a quasi-static approximation and a flux-limited diffusion treatment for neutrino transport. For an accurate assessment of neutrino-flavor instabilities, we run another multi-angle neutrino transport simulation by Sedonu on top of fluid distributions obtained from the PNS cooling model. We note that full angular information on neutrinos is necessary to assess flavor instabilities, particularly for the FFI. According to our analyses, we find no positive sign of either FFI or CFI in our model.

We also discuss the physical reasons for the absence of FFI and CFI. For both FFI and CFI, high degeneracy of $\nu_{e}$ weakens the emission of $\bar{\nu}_{e}$ and hampers the appearance of flavor instabilities. We also find that CFI is sensitive to $\nu_{x}$. As an evidence, we show that CFI can occur if there is no contribution from $\nu_{x}$. We, thus, conclude that the combined effects of $\nu_{x}$ and the weak emission of $\bar{\nu}_{e}$ account for suppressing CFI.

Before we close this paper, we describe some important limitations in the present study. Most importantly, we neglect multi-dimensional effects such as PNS convection. As shown in previous studies, the results of flavor instabilities can be very sensitive to the dimensionality; in fact, FFIs tend to be suppressed in spherically symmetric models, but multi-dimensional fluid instabilities offer environments more conducive to development of the FFI. This is mainly because PNS convection can facilitate deleptonization of CCSN core (Dessart et al., 2006; Buras et al., 2006; Nagakura et al., 2020), that reduces the disparity between $\nu_{e}$ and $\bar{\nu}_{e}$. As a result, ELN angular crossings can occur much more easily than in spherically symmetric models (Delfan Azari et al., 2020; Glas et al., 2020; Nagakura et al., 2021a; Akaho et al., 2023). We also note that reducing the lepton number in the core would facilitate resonance-like CFI (see Akaho et al. 2024a). On the other hand, PNS convection can increase $\nu_{x}$ diffusion (Nagakura et al., 2020), which may work to suppress CFIs. This suggests that PNS convection has two competing effects regarding CFIs, exhibiting requirements of detailed investigations.

Another issue in the present study is neutrino-matter interactions. As described in Sec. 5, neutrino-matter interactions are one of the most important ingredients to determine neutrino radiation fields, but our current treatments contain an incomplete and imperfect set of neutrino interaction processes, and the dynamical diffusion transport uses interaction rates that are different from the post-processed Monte Carlo calculation. In the future, it will be important to do similar analyses with high-fidelity input physics in order to draw more robust conclusions about whether flavor instabilities can occur and have the ability to impact neutrino signals, nucleosynthesis, and PNS cooling. The results in this study serve as an important reference for future work.

M.Z. is supported by the Japan Society for Promotion of Science (JSPS) Grant-in-Aid for JSPS Fellows (Grants No. 22KJ2906) from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) in Japan. H.N. is supported by Grant-in-Aid for Scientific Research (23K03468) and also by the NINS International Research Exchange Support Program. S.R. was supported by a National Science Foundation Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2001760. H.S. is supported by JSPS Grant-in-Aid for Scientific Research on Innovative Areas “Unraveling the History of the Universe and Matter Evolution with Underground Physics” (No. 19H05802 and No. 19H05811) and for Transformative Research Areas “The creation of multi-messenger astrophysics” (No. 24H01817). C.K. is supported by JSPS KAKENHI Grant Numbers JP20K14457, JP22H04577 and JP24H02245. The numerical computations were partly carried out on Cray XC50 at the Center for Computational Astrophysics in National Astronomical Observatory of Japan. This work is also supported by the HPCI System Research Project (Project ID:240079).