Collective fast neutrino flavor conversions in a 1D box: Initial conditions and long-term evolution

We consider a simple neutrino-dense system which has translation symmetry in the $x$- and $y$- directions in space, the same as in Refs. Martin et al. (2020); Bhattacharyya and Dasgupta (2020, 2021). For simplicity, initially we assume a pure system consisting of only $\nu_{e}$ and $\bar{\nu}_{e}$ (before setting small flavor perturbations; see below). Focusing on fast flavor conversions, we neglect the momentum changing collisions and only keep the neutrino–neutrino forward scattering contributions Fuller et al. (1987); Pantaleone (1992); Sigl and Raffelt (1993) in the effective Hamiltonian, i.e., we drop the terms originating from vacuum neutrino mixing and neutrino–matter forward scattering contributions.

Assuming the differential neutrino angular distribution $dn_{\nu_{e}(\bar{\nu}_{e})}/d\Omega$ is uniform in $z$ and taking the two-flavor approximation, the equation of motion, which governs the evolution of the normalized neutrino density matrices (see below) $\varrho$ (for neutrinos) and $\bar{\varrho}$ (for antineutrinos) is given as follows¹¹1When only considering the neutrino–neutrino self-interaction term, one can redefine the antineutrino density matrix by $\bar{\rho}\rightarrow-\sigma_{y}^{\dagger}\bar{\rho}\sigma_{y}$ such that neutrinos and antineutrinos are treated in equal footing Duan et al. (2006a). However, here we choose to treat them separately for the purpose of consistency with follow-up works that may include other terms in the Hamiltonian and the collisions.:

with

in the flavor basis. Since we omit the vacuum mixing contribution, the flavor evolution of $\rho$ and $\bar{\rho}$ do not depend on the neutrino energy.

In Eq. (1), the Hamiltonian $H$ and $\bar{H}$ can be explicitly written down as:

where $\mu=\sqrt{2}G_{F}n_{\nu_{e}}$ with $G_{F}$ being the Fermi coupling constant. The asymmetry parameter $\alpha=n_{\bar{\nu}_{e}}/n_{\nu_{e}}$ quantifies the ratio between the number density of $\bar{\nu}_{e}$ and $\nu_{e}$. The angular distribution function $g_{\nu(\bar{\nu})}(v_{z})$ relates to the physical neutrino phase-space distribution function $f_{\nu_{e}(\bar{\nu}_{e})}$ without any flavor perturbations as

where $E_{\nu}$ is the energy of a neutrino. Note that the azimuthal symmetry in phase space of $f_{\nu_{e}(\bar{\nu}_{e})}$ with respect to the $z$ direction is implicitly taken. In this notation, $\rho_{ee}=\bar{\rho}_{ee}=1$ for our system without any initial flavor perturbations while all other matrix elements are zero.

Since $\mu$ is the only dimensional quantity in Eq. (1), we define $\mu\equiv 1$ and express $t$ and $z$ in dimensionless form (in the units of $\mu^{-1}$) hereafter.

We simulate the flavor evolution of neutrinos inside a box of size $L_{z}=1200$ in $z\in[-600,600]$ with the periodic boundary condition. Similar to Ref. Martin et al. (2020), we parametrize $g_{\nu(\bar{\nu})}$ as

with normalization condition $\int_{-1}^{1}dv_{z}g_{\nu(\bar{\nu})}(v_{z})=1$. Throughout the paper, we fix $\sigma_{\nu}=0.6$ and $\sigma_{\bar{\nu}}=0.5$, such that $\bar{\nu}_{e}$ are more forward-peaked than $\nu_{e}$. For the asymmetry parameter $\alpha$, we take $\alpha=0.9,1.0$, $1.1$, $1.2$, and $1.3$. Figure 1 shows the corresponding $\nu$ELN distributions $G_{\nu}(v_{z})\equiv g_{\nu}(v_{z})-\alpha g_{\bar{\nu}}(v_{z})$. With increasing value of $\alpha$, the $\nu$ELN crossing where $G_{\nu}(v_{z,c})=0$ occurs at smaller $v_{z,c}$ with $|G_{\nu}(v_{z})|$ being larger (smaller) for $v_{z}>v_{z,c}$ ($v_{z}<v_{z,c}$). The values of $v_{z,c}$ for $\alpha=0.9,1.0,1.1,1.2$, and $1.3$ are 0.65, 0.45, 0.33, 0.23, and 0.15, correspondingly.

For a system without the vacuum Hamiltonian, it requires a small perturbation in $\rho_{ex}$ ($\bar{\rho}_{ex}$) to trigger the flavor instabilities. For simplicity, we assume that the perturbations are independent of $v_{z}$. Specifically, we use the following description for our initial condition at $t=0$ as:

For $\epsilon(z)$, we explore the following two different choices:

1. 1.

   Point-source like perturbations centered at $z=0$: $\epsilon(z)=\epsilon_{0}\exp[-z^{2}/50]$.

2. 2.

   Random perturbations: $\epsilon(z)$ is randomly assigned by a real value between $0$ and $\epsilon_{0}$.

We take $\epsilon_{0}=10^{-6}$ throughout the paper.

We discretize $z$ and $v_{z}$ with $N_{z}$ and $N_{v_{z}}$ grids, and evolve $\varrho$ and $\bar{\varrho}$ defined at the grid centers. In evaluating the advection terms $\partial\varrho/{\partial z}$ and $\partial\bar{\varrho}/{\partial z}$, we use two different methods to check for consistency. The first method uses the forth-order (central) finite difference method with artificial dissipation using the third-order Kreiss-Oliger formulation Kreiss and Oliger (1973). For the second method, we use finite volume method plus the seventh order weighted essentially nonoscillatory (WENO) scheme Shu (1998, 2003). For time evolution, we use fourth-order Runge-Kutta method. The numerical implementation details will be reported in a separate publication together with the public release of our simulation code COSE$\nu$ George et al. (2021).

Using the above initial and boundary conditions, we perform simulations with fiducial numerical parameters of $N_{Z}=12000$ and $N_{v_{z}}=200$, with a fixed time step size $\Delta t=C_{\rm CFL}\Delta z$, where $\Delta z=L_{z}/N_{Z}$ and the Courant–Friedrichs–Lewy number $C_{\rm CFL}=0.4$. In Appendixes A and B, we show the comparison of results obtained with different resolutions and with two different numerical methods for advection. For the rest of the paper, all results are obtained using the finite volume method with seventh order WENO scheme, which provides better numerical accuracy given the same resolution.

Before we discuss our results, let us define the neutrino and antineutrino polarization vectors $\mathbf{P}$ and $\mathbf{\bar{P}}$ that are often used in literature. The three components of the polarization vector $\mathbf{P}$ are defined as

which satisfy $\rho\equiv(\mathbf{P}\cdot\bm{\sigma}+P_{0}\mathcal{I})/2$, where $\sigma_{i}$ are the Pauli matrices and $\mathcal{I}$ is the identity matrix. For antineutrinos, we have

With this definition, the equation of motion for $\mathbf{P}$ has the same form as for $\mathbf{\bar{P}}$, and we can treat neutrinos and antineutrinos in equal footing using the $\nu$ELN spectrum $G_{\nu}$ defined in Sec. II.2 to account for the contributions from both neutrinos and antineutrinos in the Hamiltonian. Dropping the bar for antineutrinos, we have:

where $\mathbf{H}(t,z,v_{z})=\int_{-1}^{1}dv^{\prime}_{z}G_{\nu}(v^{\prime}_{z})\mathbf{P}(t,z,v^{\prime}_{z})(1-v_{z}v^{\prime}_{z})$. Clearly, in the absence of collisions, the length of the polarization vectors are unity for all $z$ and $v_{z}$, given that our initial condition has $|\mathbf{P}|=1$ uniformly in $z$. Another physical quantity that is conserved globally is the net neutrino $e-x$ lepton number in the entire simulation domain Bhattacharyya and Dasgupta (2020, 2021):

Note that locally this net lepton number can be transferred from one location to another due to the advection. However, globally $L_{e-x}$ must be a conserved quantity in time for the entire system.

In the regime where $|\rho_{ex}|\ll 1$ and $|\bar{\rho}_{ex}|\ll 1$, the evolution of the system can be qualitatively understood by performing the standard technique of linearized instability analysis Banerjee et al. (2011); Izaguirre et al. (2017). By inserting $\rho_{ex}\propto\exp[-i(\omega t-k_{z}z)]$ (same for $\bar{\rho}_{ex}$) and keeping terms in the lowest order of perturbations, Eq. (1) leads to a dispersion relation of the collective mode of the system $\omega(k_{z})$ Izaguirre et al. (2017); Capozzi et al. (2017); Yi et al. (2019). The complex branch of $\omega$ for real $k_{z}$ signifies the existence of flavor instabilities such that small off diagonal perturbations can grow exponentially in time as the system evolves. Note that here we only include the solution that preserves the the axial symmetry. For the symmetry-breaking solutions, we omit them here because our simulation setup does not allow symmetry-breaking solutions.

In Fig. 2, we show in panels (a) and (b) the dispersion relation $\omega(k_{z})$ for systems with $\nu$ELN asymmetry parameter $\alpha=0.9$ and $1.1$. For $\alpha=0.9$, there are two regions in $k_{z}$ where complex $\omega$ exist. The maximal value of Im$(\omega_{\rm max})\simeq 0.044$ locates at $k_{z,{\rm max}}\simeq-0.165$. At the same $k_{z,{\rm max}}$, $(d\omega/dk_{z})_{\rm max}\simeq 0.515$. For $\alpha=1.1$ where the $\nu$ELN crossing is deeper in larger $v_{z}$ (see Fig. 1), Im$(\omega_{\rm max})\simeq 0.107$ at $k_{z,{\rm max}}\simeq-0.225$. The corresponding $(d\omega/dk_{z})_{\rm max}\simeq 0.278$.

In the bottom panels [(c) and (d)] of Fig. 2, we show the evolution of $|P_{\perp}(z)|=\sqrt{P_{1}^{2}+P_{2}^{2}}=2|\rho_{ex}|$ for $v_{z}=0.995$ in the liner regime for the same $\nu$ELNs with point-source like perturbations. Clearly, an initial Gaussian packet of $P_{\perp}$ grows exponentially with time and the peak position of $P_{\perp}$ moves toward positive $z$ direction for both the cases. The growth rate of the maximal value of $P_{\perp,{\rm max}}$ and the velocity $dz(P_{\perp,{\rm max}})/dt$ both agree well with the values of Im$(\omega_{\rm max})$ and $(d\omega/dk_{z})_{\rm max}$ obtained in the dispersion relation above. The growth of perturbations using $\alpha=1.1$ is indeed much faster than that with $\alpha=0.9$. Moreover, the $P_{\perp}(z)$ with $\alpha=1.1$ spreads over both positive and negative $z$ directions. For the case with $\alpha=0.9$, although the width of $P_{\perp}(z)$ also increases with time, it mainly drifts to the positive $z$ direction without spreading over the negative $z$ direction.

When the perturbations grow to the nonlinear regime, wavelike oscillatory features develop Martin et al. (2020). In Fig. 3, we show the snapshots of $P_{3}(z,v_{z})$ at different times for $\alpha=0.9$ and $1.1$ using the same point-source like perturbations as in Sec. III.1. For the case of $\alpha=0.9$, the flavor evolution behavior of the system is in general similar to those reported in Ref. Martin et al. (2020). Flavor waves develop and mainly propagate toward the positive $z$ direction. This is consistent with the growth of perturbations in the linear regime discussed earlier. An interesting feature shown here is that although the flavor oscillations initially can affect all $v_{z}$ modes, illustrated by the vertical stripes in the upper two subpanels in Fig. 3(a), this effect diminishes when time proceeds and flavor conversions are roughly confined in $v_{z}\gtrsim v_{z,c}\simeq 0.65$.

Another interesting feature is when the forward-traveling wave front interacts with the slowly backward propagating part after $t\gtrsim 1200$, it pushes the whole pattern to drift toward positive $z$ direction. Meanwhile, this interaction breaks the coherent wavelike pattern. Substructures in smaller scale develop such that the orientation of the neutrino polarization vectors varies rapidly in $z$ for $v_{z}\gtrsim v_{z,c}$. Consequently, although at each location $z$, neutrinos with different $v_{z}\gtrsim v_{z,c}$ still have $|\mathbf{P}|=1$, the “average polarization vector” over the $z$ domain can shrink due to their misalignment. Such a flavor state was referred as “flavor depolarization” in Refs. Bhattacharyya and Dasgupta (2020, 2021) and we will discuss its behavior in more detail in the next section. For $v_{z}\lesssim v_{z,c}$, most neutrinos remain unaffected.

For $\alpha=1.1$ shown in panel (b), flavor conversions quickly develop toward both the positive and negative $z$ directions and produces coherent and wavelike structure, once again consistent with that indicated by the growth of perturbations in the linear regime. Similarly, when the forward and backward propagating modes interact after $t\gtrsim 665$, smaller structures develop and cause a major part of $v_{z}$ reaching closer to flavor depolarization. One important difference from the previous $\alpha=0.9$ case is now flavor depolarization happens mostly in $v_{z}\lesssim v_{z,c}\simeq 0.45$. We will also discuss this feature and its consequences in Sec. IV.2.

For the other $\nu$ELN spectra with $\alpha=1.0$, $1.2$, and $1.3$, the behaviors are qualitatively similar to that with $\alpha=1.1$. Full simulation animations are available at https://mengruwuu.github.io/COSEv1dbox/.

Now, let us look at how the flavor systems evolve when we adopt different seed of random perturbations. We show again $P_{3}(z,v_{z})$ at different time snapshots in Fig. 4 for $\alpha=0.9$ and $1.1$. With random perturbations, some locations have larger $P_{\perp}$ initially such that flavor conversions develop faster around those locations (see the top subpanels). Similar to what discussed in Sec. III.2 for single point-source like perturbations, flavor waves can initially develop independently now at different locations. The flavor waves transported in space interact and break the coherent pattern to form small-scale structures. Thus, with perturbations everywhere in space, the interaction of flavor waves happen at much earlier times and no large-scale coherent structure can be formed, differently from cases with point-source like perturbations. The systems can then reach closer to flavor depolarization in $v_{z}\gtrsim v_{z,c}$ ($v_{z}\lesssim v_{z,c}$) for $\alpha=0.9$ ($1.1$) within a much shorter amount of time.

Once again, for the other $\nu$ELN spectra with $\alpha=1.0$, $1.2$, and $1.3$, the behaviors are qualitatively similar to that with $\alpha=1.1$. Full simulation animations are also available at https://mengruwuu.github.io/COSEv1dbox/.

We define the overall flavor survival probabilities of neutrinos and antineutrinos, $\overline{\langle P_{ee}\rangle}$ and $\overline{\langle{\bar{P}}_{ee}\rangle}$ in the simulation domain by

where the brackets and overline denote averaging over $z$ and $v_{z}$, respectively. The integration limits over $z$ and $v_{z}$ are from $-L_{z}/2$ to $L_{z}/2$ and from $-1$ to $1$, respectively.

Fig. 5 shows the time evolution of $\overline{\langle P_{ee}\rangle}$ and $\overline{\langle{\bar{P}}_{ee}\rangle}$ for all five different $\nu$ELN spectra that we considered (see Fig. 1). The cases with point-source like perturbations and random perturbations are shown in panels (a) and (b), respectively. First, we see that cases with larger $\alpha$ reach the final asymptotic states earlier in time for both point-source like and random perturbations. Meanwhile, all cases with random perturbations reach the asymptotic states much earlier than those with point-source like perturbations as discussed in Sec. III.3.

Second, systems with $\alpha=1$ (equal number of neutrinos and antineutrinos) have final asymptotic values of $\overline{\langle P_{ee}\rangle}_{f}\simeq\overline{\langle{\bar{P}}_{ee}\rangle}_{f}\simeq 0.5$, i.e., full flavor depolarization is almost reached for the entire system. For systems that are more asymmetric, the asymptotic $\overline{\langle P_{ee}\rangle}_{f}$ and $\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ are larger, i.e., on average less $\nu_{e}$ and $\bar{\nu}_{e}$ are converted. Comparing $\overline{\langle P_{ee}\rangle}_{f}$ to $\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ for a given $\alpha$, one finds that cases with $\alpha>1$ have larger values of asymptotic $\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ than those of $\overline{\langle P_{ee}\rangle}_{f}$ (vice versa for $\alpha<1$). This is related to the fact that for $\alpha>1$ ($\alpha<1$), neutrinos with $v_{z}\lesssim v_{z,c}$ ($v_{z}\gtrsim v_{z,c}$) experience more flavor conversions. Since for our adopted neutrino angular spectra, $g_{\bar{\nu}}(v_{z})$ is more forward peaked in positive $v_{z}$ than $g_{\nu}(v_{z})$, more flavor conversions in larger $v_{z}$ naturally lead to smaller $\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ than $\overline{\langle P_{ee}\rangle}_{f}$.

Comparing the final values of $\overline{\langle P_{ee}\rangle}_{f}$ and $\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ for cases with point-source like and random perturbations, the point-source like ones generally lead to a slightly smaller $\overline{\langle P_{ee}\rangle}_{f}$ and $\overline{\langle{\bar{P}}_{ee}\rangle}_{f}$ than the random ones. This is due to the reason discussed in Sec. III.2: although the interaction of the flavor waves lead to states close to flavor depolarization, for cases with point-source like perturbations, some large-scale coherent structures can remain until the end of the simulations. This can also be seen in Fig. 6, which shows the spatially averaged flavor survival probabilities $\langle P_{ee}(v_{z})\rangle_{f}$ as functions of $v_{z}$ at the end of our simulations, where $\langle P_{ee}(v_{z})\rangle$ is defined by

Note that $\langle{\bar{P}}_{ee}(v_{z})\rangle$ that can be similarly defined are equal to $\langle P_{ee}(v_{z})\rangle$, when only neutrino-neutrino self-interaction terms are included here.

Figure 6 clearly shows several interesting features. First, for cases with either point-source like or random perturbations, those with $\alpha>1$ ($\alpha<1$) lead to nearly full flavor depolarization ($\langle P_{ee}(v_{z})\rangle_{f}\simeq 0.5$) for velocity modes $v_{z}<v_{z,c}$ ($v_{z}>v_{z,c}$). For velocity modes in the other side of the $\nu$ELN, nearly flavor equipartition cannot be achieved and $\langle P_{ee}(v_{z})\rangle_{f}$ gradually increase from $\simeq 0.5$ to larger values for larger (smaller) $v_{z}$ for $\alpha>1$ ($\alpha<1$). For the case with $\alpha=1$, i.e., symmetric in neutrinos and antineutrinos, full flavor depolarization for the entire $\nu$ELN can be achieved. Comparing results with different perturbations, one sees that due to the remaining large-scale coherent structures (see e.g., Fig. 3), systems with point-source like perturbations give rise to smaller $\langle P_{ee}(v_{z})\rangle_{f}\simeq 0.4$ in some regions in $v_{z}<v_{z,c}$ ($v_{z}>v_{z,c}$) for $\alpha>1$ ($\alpha<1$) than those with random perturbations.

Based on our simulation results, we make a conjecture as follows. For a system with the periodic boundary condition in one spatial direction (while having perfect translation symmetry in the other two directions) where interaction of flavor waves can happen, the system can evolve to an asymptotic state where nearly flavor equipartition can be reached for velocity modes in one side of the $\nu$ELN, when averaged over space. In other words, when averaged over space, the $\nu$ELN crossing nearly vanishes. For such a scenario to happen, the net neutrino $e-x$ lepton number conservation $L_{e-x}$ [see Eq. (10)] would enforce that nearly flavor depolarization can only happen for velocity modes in the shallower side of the $\nu$ELN.

An additional remark is: comparing Fig. 5 and Fig. 6, we also point out that although Fig. 6 seems to suggest that a narrower range of velocity modes in $\nu$ELN distribution reaches nearly flavor depolarization for $\alpha=0.9$ than those with $\alpha>1$, on average, more neutrinos and antineutrinos are being converted. Once again, this is because the spectra $g_{\nu}(v_{z})$ and $g_{\bar{\nu}}(v_{z})$ are more forward peaked. This means that although a system with $\alpha<1$ may appears to be affected in less phase space volume in $v_{z}$, flavor conversions there can actually have a larger impact on physical processes that are flavor dependent in realistic astrophysical environments.

Refs. Bhattacharyya and Dasgupta (2020, 2021) proposed that the behavior of the system can be understood by examining the time evolution of $\langle\mathbf{P}(v_{z})\rangle$. From Eq. (9), one obtains

where $\mathbf{M}_{i}(z,t)\equiv\int dv_{z}L_{i}(v_{z})G_{\nu}(v_{z})\mathbf{P}(z,v_{z},t)$ with $L_{i}(v_{z})$ the Legendre Polynomials. The authors of Refs. Bhattacharyya and Dasgupta (2020, 2021) argued that the spatial averaging of the cross products in Eq. (13) can be approximated as the cross products of two vectors that are spatially averaged separately:

By doing so, the time evolution of $\langle\mathbf{P}(v_{z},t)\rangle$ may be understood as a vector precessing around an effective Hamiltonian vector $\langle\mathbf{H}(v_{z})\rangle=\langle\mathbf{M}_{0}\rangle-v_{z}\langle\mathbf{M}_{1}(t)\rangle$. Moreover, the depolarization of $\langle\mathbf{P}(v_{z},t)\rangle$ happens when the $\langle H\rangle_{\perp}=\sqrt{\langle H\rangle_{1}^{2}+\langle H\rangle_{2}^{2}}$ roughly exceeds $|\langle{\tilde{H}}\rangle_{3}|$ Bhattacharyya and Dasgupta (2021), where $\langle{\tilde{H}}\rangle_{3}$ differs from $\langle H\rangle_{3}$ by a factor related to $\langle\mathbf{M}_{2}\rangle$ when going to a co-rotating frame (see Ref. Bhattacharyya and Dasgupta (2020) for details).

We carefully examine the validity of these claims. In Fig. 7, we show in panel (a) the time evolution of $|\langle\tilde{H}\rangle_{3}|$ and $\langle H_{\perp}\rangle$ for $\alpha=0.9$ of $v_{z}=1$ and $\alpha=1.1$ of $v_{z}=-1$ using random perturbations as examples. Despite these modes undergo nearly flavor depolarization, Fig. 7(a) shows that the values of $\langle H_{\perp}\rangle$ remain much smaller than $|\langle\tilde{H}\rangle_{3}|$ nearly during the whole time, in contrast to what discussed in Ref. Bhattacharyya and Dasgupta (2021).

In panels (b) and (c) of Fig. 7, we show the evolution of $|\langle\mathbf{P}(v_{z})\rangle|$ and the relative angles $\theta_{r}$ between the two vectors on the right-hand sides of Eqs. (13) and (14) for ten $v_{z}$ modes for the case with $\alpha=1.1$ with random perturbations. Clearly, for modes with $v_{z}\lesssim v_{z,c}$ that reach flavor depolarization (darker colors), i.e., $|\langle\mathbf{P}(v_{z})\rangle|\rightarrow 0$, their $\cos\theta_{r}$ oscillate rapidly between -1 and 1. This indicates that Eq. (14) is in fact not a good approximation of Eq. (13).

In sum, our findings here suggest that the picture which explains the late-time behavior and the depolarization of the system proposed in Refs. Bhattacharyya and Dasgupta (2020, 2021) using the spatially averaged polarization vectors needs to be revisited.