Effect of collisions on neutrino flavor inhomogeneity in a dense neutrino gas

TL;DR

This study explores how collisions influence the development of flavor inhomogeneities in a dense neutrino gas, revealing that collisions can suppress instabilities at small scales but not at large scales, especially for inverted mass hierarchy.

Contribution

It provides new insights into the role of incoherent collisions in neutrino flavor stability, particularly highlighting scale-dependent effects and hierarchy differences.

Findings

Collisions suppress flavor instability growth at small scales.

Inverted hierarchy instabilities are resistant to suppression by collisions at large scales.

Normal hierarchy shows broader instability range without collisions.

Abstract

We investigate the stability, with respect to spatial inhomogeneity, of a two-dimensional dense neutrino gas. The system exhibits growth of seed inhomogeneity due to nonlinear coherent neutrino self-interactions. In the absence of incoherent collisional effects, we observe a dependence of this instability growth rate on the neutrino mass spectrum: the normal neutrino mass hierarchy exhibits spatial instability over a larger range of neutrino number density compared to that of the inverted case. We further consider the effect of elastic incoherent collisions of the neutrinos with a static background of heavy, nucleon-like scatterers. At small scales, the growth of flavor instability can be suppressed by collisions. At large length scales we find, perhaps surprisingly, that for inverted neutrino mass hierarchy incoherent collisions fail to suppress flavor instabilities, independent of the coupling strength.

Equations

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

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

$$\mathsf{H}$$

$$\bar{\mathsf{H}}$$

$${\mathsf{H}}_0$$

$${\mathsf{H}}_{\nu \nu }$$

$${\rho }^{\theta }(t,\vec{x})$$

$$i\frac{\partial }{\partial t}{\rho }_{\vec{k}}^{\theta }(t)$$

$$+\mu \sum _{k^{\prime }}{\int }_0^{2\pi }d{\theta }^{\prime }\left[({\rho }_{\vec{k}-\overrightarrow{k^{\prime }}}^{{\theta }^{\prime }}(t)-\alpha {\bar{\rho }}_{\vec{k}-\overrightarrow{k^{\prime }}}^{{\theta }^{\prime }}(t)),{\rho }_{\overrightarrow{k^{\prime }}}^{\theta }(t)\right]$$

$$\times (1-\cos (\theta -{\theta }^{\prime })).$$

$$\displaystyle\frac{d}{dt}\rho^{\theta}(\vec{x},t)\Big{|}_{\rm coll}=$$

$$+\frac{1}{2}\left\{{\Pi }^{\text{gain}},\left(1-{\rho }^{\theta }(\vec{x},t)\right)\right\}.$$

$${\Pi }^{\text{loss/gain}}$$

$${\rho }^{\theta }(t,\vec{x})=\left(\begin{array}{cc}1-{\Delta }^{\theta }(t,\vec{x})& {ϵ}^{\theta }(t,\vec{x})\\ {ϵ}^{\theta \ast }(t,\vec{x})& {\Delta }^{\theta }(t,\vec{x})\end{array}\right)$$

$${\bar{\rho }}^{\theta }(t,\vec{x})=\left(\begin{array}{cc}1-{\bar{\Delta }}^{\theta }(t,\vec{x})& {\overline{ϵ}}^{\theta }(t,\vec{x})\\ {\overline{ϵ}}^{\theta \ast }(t,\vec{x})& {\bar{\Delta }}^{\theta }(t,\vec{x})\end{array}\right)$$

$${\rho }_{\vec{k}}^{\theta }(t)=\left(\begin{array}{cc}{\delta }_{\vec{k},0}& {ϵ}_{\vec{k}}^{\theta }(t)\\ {ϵ}_{-\vec{k}}^{\theta \ast }(t)& 0\end{array}\right)$$

$$i\frac{\partial }{\partial t}{ϵ}_{\vec{k}}^{\theta }(t)=\vec{v}\cdot \vec{k}{ϵ}_{\vec{k}}^{\theta }(t)-\omega {ϵ}_{\vec{k}}^{\theta }(t)$$

$$+2\pi \mu (1-\alpha ){ϵ}_k^{\theta }-\mu {\int }_0^{2\pi }d{\theta }^{\prime }({ϵ}_k^{{\theta }^{\prime }}-\alpha {\overline{ϵ}}_k^{{\theta }^{\prime }})(1-\cos (\theta -{\theta }^{\prime }))$$

$$-i2\pi c_1{ϵ}_{\vec{k}}^{\theta }(t)+i{\int }_0^{2\pi }d{\theta }^{\prime }(c_1-c_2\cos (\theta -{\theta }^{\prime })){ϵ}_{\vec{k}}^{{\theta }^{\prime }}(t),$$

$$i\frac{\partial }{\partial t}{\overline{ϵ}}_{\vec{k}}^{\theta }(t)=\vec{v}\cdot \vec{k}{\overline{ϵ}}_{\vec{k}}^{\theta }(t)+\omega {\overline{ϵ}}_{\vec{k}}^{\theta }(t)$$

$$+2\pi \mu (1-\alpha ){\overline{ϵ}}_k^{\theta }-\mu {\int }_0^{2\pi }d{\theta }^{\prime }({ϵ}_k^{{\theta }^{\prime }}-\alpha {\overline{ϵ}}_k^{{\theta }^{\prime }})(1-\cos (\theta -{\theta }^{\prime }))$$

$$-i2\pi c_1{\overline{ϵ}}_{\vec{k}}^{\theta }(t)+i{\int }_0^{2\pi }d{\theta }^{\prime }(c_1-c_2\cos (\theta -{\theta }^{\prime })){\overline{ϵ}}_{\vec{k}}^{{\theta }^{\prime }}(t).$$

$${ϵ}_{\vec{k}}^{\theta }(t)=Q_{\vec{k}}^{\theta }\exp (-i{\Omega }_{\vec{k}}t)$$

$${\overline{ϵ}}_{\vec{k}}^{\theta }(t)={\overline{Q}}_{\vec{k}}^{\theta }\exp (-i{\Omega }_{\vec{k}}t).$$

$$\kappa \equiv \text{Im }({\Omega }_{\vec{k}}),$$

$$\left[{\Omega }_{\vec{k}}-\vec{v}\cdot \vec{k}+\omega -2\pi \mu (1-\alpha )+i2\pi c_1\right]Q_{\vec{k}}^{\theta }$$

$$=-\mu {\int }_0^{2\pi }d{\theta }^{\prime }(Q_{\vec{k}}^{{\theta }^{\prime }}-\alpha {\overline{Q}}_{\vec{k}}^{{\theta }^{\prime }})(1-\cos (\theta -{\theta }^{\prime }))$$

$$+i{\int }_0^{2\pi }d{\theta }^{\prime }{Q}_{\vec{k}}^{{\theta }^{\prime }}(c_1-c_2\cos (\theta -{\theta }^{\prime })),$$

$$\left[{\Omega }_{\vec{k}}-\vec{v}\cdot \vec{k}-\omega -2\pi \mu (1-\alpha )+i2\pi c_1\right]{\overline{Q}}_{\vec{k}}^{\theta }$$

$$=-\mu {\int }_0^{2\pi }d{\theta }^{\prime }(Q_{\vec{k}}^{{\theta }^{\prime }}-\alpha {\overline{Q}}_{\vec{k}}^{{\theta }^{\prime }})(1-\cos (\theta -{\theta }^{\prime }))$$

$$+i{\int }_0^{2\pi }d{\theta }^{\prime }{\overline{Q}}_{\vec{k}}^{{\theta }^{\prime }}(c_1-c_2\cos (\theta -{\theta }^{\prime })),$$

$$Q_{\vec{k}}^{\theta }$$

$$+\frac{i}{D_{\vec{k}}}\left[c_1{a}_{\vec{k}}-c_2\sin \theta s_{\vec{k}}-c_2\cos \theta c_{\vec{k}}\right]$$

$${\overline{Q}}_{\vec{k}}^{\theta }$$

$$+\frac{i}{{\overline{D}}_{\vec{k}}}\left[c_1{\bar{a}}_{\vec{k}}-c_2\sin \theta {\bar{s}}_{\vec{k}}-c_2\cos \theta {\bar{c}}_{\vec{k}}\right],$$

$$D_{\vec{k}}$$

$${\overline{D}}_{\vec{k}}$$

$$a_{\vec{k}}={\int }_0^{2\pi }d\theta Q_{\vec{k}}^{\theta }$$

$$s_{\vec{k}}={\int }_0^{2\pi }d\theta \sin \theta Q_{\vec{k}}^{\theta }$$

$$c_{\vec{k}}={\int }_0^{2\pi }d\theta \cos \theta Q_{\vec{k}}^{\theta }$$

$$\left(\begin{array}{cccccc}-J_1[1](c_1)-1& \alpha I_1[1]& J_1[\sin \theta ](c_2)& -\alpha I_1[\sin \theta ]& J_1[\cos \theta ](c_2)& -\alpha I_1[\cos \theta ]\\ -I_2[1]& -J_2[1](c_1)-1& I_2[\sin \theta ]& J_2[\sin \theta ](c_2)& I_2[\cos \theta ]& J_2[\cos \theta ](c_2)\\ -J_1[\sin \theta ](c_1)& \alpha I_1[\sin \theta ]& J_1[{\sin }^2\theta ](c_2)-1& -\alpha I_1[{\sin }^2\theta ]& J_1[\sin \theta \cos \theta ](c_2)& -\alpha I_1[\sin \theta \cos \theta ]\\ -I_2[\sin \theta ]& -J_2[\sin \theta ](c_1)& I_2[{\sin }^2\theta ]& J_2[{\sin }^2\theta ](c_2)-1& I_2[\sin \theta \cos \theta ]& J_2[\sin \theta \cos \theta ](c_2)\\ -J_1[\cos \theta ](c_1)& \alpha I_1[\cos \theta ]& J_1[\sin \theta \cos \theta ](c_2)& -\alpha I_1[\sin \theta \cos \theta ]& J_1[{\cos }^2\theta ](c_2)-1& -\alpha I_1[{\cos }^2\theta ]\\ -I_2[\cos \theta ]& -J_2[\cos \theta ](c_1)& I_2[\sin \theta \cos \theta ]& J_2[\sin \theta \cos \theta ](c_2)& I_2[{\cos }^2\theta ]& J_2[{\cos }^2\theta ](c_2)-1\end{array}\right)\left(\begin{array}{c}{a}_{\vec{k}}\\ {\bar{a}}_{\vec{k}}\\ s_{\vec{k}}\\ {\bar{s}}_{\vec{k}}\\ c_{\vec{k}}\\ {\bar{c}}_{\vec{k}}\end{array}\right)$$

$$I_1[f]$$

Full text

1 Introduction

Several recent studies, concerned primarily with anisotropic astrophysical environments, have demonstrated that the nonlinear neutrino self-coupling results in spatial instabilities [1, 2, 3, 4]. These works indicate the importance of taking into account deviations from spherical symmetry in describing the evolution of, for example, the supernova neutrino flavor field. Previous studies of the neutrino evolution [5, 6, 7, 8, 9] in the early universe generally assume homogeneity and isotropy obtains at all times from the epochs of weak equilibrium through neutrino decoupling and Big Bang nucleosynthesis (BBN). The objective of the present study is to address the validity of this assumption in an exploratory calculation. To this end, in this Letter we study the stability of a two-dimensional, two-flavor dense neutrino gas with respect to the growth of seed spatial inhomogeneity driven by neutrino self-coupling. In this model, we also investigate the impact of incoherent elastic collisions on the stability properties. We find, perhaps surprisingly, that elastic scattering-angle dependent collisions are not always effective at damping flavor instabilities.

In extreme environments with large neutrino number densities, such as the interior of a core-collapse supernovae (CCSN) and in the early universe, the nonlinearity associated with neutrino self-interactions lead to an array of interesting effects beyond the “standard” MSW [10, 11]. These effects, which have their origin in the non-linear character of the neutrino evolution equations, include collective oscillations in an idealized model of supernova [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22], matter-neutrino resonant behavior [23, 24, 25], spontaneous breaking of the axial symmetry [26, 27], and emergence of spatial inhomogeneity [1, 2], to cite a few examples.

In this context, the early universe is particularly interesting to study since non-equilibrium and neutrino flavor oscillation effects may be important for precision cosmological probes, such as BBN. In the early universe neutrino oscillations can be induced by flavor and/or lepton asymmetries.

Studies of neutrino weak decoupling in a Boltzmann equation approach [5, 28, 9] suggest that neutrino spectra deviate from equilibrium Fermi-Dirac distributions at the percent level and display flavor dependence ($\nu_{e}$ versus $\nu_{\mu}$/$\nu_{\tau}$) at the level of a few percent [5, 9]. The flavor asymmetry triggers flavor oscillations and a full analysis of this problem is still lacking. Significant progress, however, has been made in Refs. [6, 7, 29].

Additionally, and complementary to the question of neutrino flavor asymmetry, neutrino oscillations play a significant role in the presence of a lepton number asymmetry, namely a difference in the spectra for neutrinos and antineutrinos (e.g. $\nu_{e}$ and $\bar{\nu}_{e}$). An electron lepton asymmetry large enough for collective neutrino oscillations to be important in the early universe can be generated in several well motivated models [30, 31, 32, 33, 34, 35, 36, 37, 38]. Moreover, a recent study [39], examining a range of lepton number flavor asymmetries over six orders of magnitude, identified a variety of qualitatively distinct regimes and non-trivial behavior at very small lepton asymmetry.

The departure point for this study, in light of the complex and sundry array of behavior driven by the nonlinear neutrino self-interaction, is the question of spatial inhomogeneity. We ask the question of whether the assumption of spatial homogeneity and isotropy, unexamined in previous works on the flavor field evolution in the early universe [6, 7], is sustained in the presence of the nonlinear neutrino self-interaction. Our finding is that, for a two-dimensional, two-flavor exploratory model of a neutrino gas, the system is indeed unstable with respect to spatial inhomogeneity. Further, we ask the question of whether incoherent collisions, with a static background field of elastic scatterers, can serve to drive the neutrino field back to homogeneity. We investigate inhomogeneity for both normal and inverted neutrino mass hierarchies. We find that for the inverted neutrino spectrum such collisions (irrespective of the coupling strength) are incapable of stabilizing the system for a significant range of length scales.

Instability does not imply, of course, the persistence of inhomogeneity for all times. Eventually, incoherent collisions may in fact drive the system toward homogeneity. In this case, the associated transient inhomogeneity is associated with entropy generation. And, if sufficiently large, this added entropy may effect the evolution of the light element nucleosynthesis in early universe. This is a question that can only be addressed with detailed numerical simulations; such a code is currently under development [9, 40, 41].

Our consideration of a simplified model of the early universe is driven mainly by pragmatic concerns. The rich treasure of effects that are closely related to nonlinearity in the collective neutrino oscillations makes prediction difficult for realistic physical systems. Any such prediction involving collective neutrino oscillations in, for example, the interior of the proto-neutron star formed in CCSN requires the solution of the neutrino-flavor density matrix in a large dimensional space,111The general case of inhomogeneous and anisotropic environments, relevant for supernovae, corresponds to an equation of motion for the neutrino density matrix in a space of seven dimensions corresponding to three spatial, three momentum coordinates and time. which is beyond the capability of even modern large-scale, multiprocessor computational platforms. However, it is possible to solve these equations for simplified models that assume deviations from highly symmetrical spatial geometries may be neglected thereby reducing the dimensionality of the solution space. It is hoped that, with this approach, insight may be gained into the development of methods that will allow the solution of the full problem; or, at least, to obtain results of phenomenological relevance.

This Letter is organized as follows. In Sec. 2, we detail the model of neutrino oscillations employed in this exploratory study of neutrino flavor stability in the early universe. Section 3 discusses the neutrino-nucleus model of the elastic, angle-dependent collision term adopted for the model. In Sec. 4, we perform linear stability analysis of the model and present results in Sec. 5. Finally, in Sec. 6, we explore the implications and limitations of the present study in the context of more realistic models. In particular, we comment on effects that inelastic contributions to the collisions may have on collective neutrino oscillations.

2 Model setup

The simplification of the spacetime and four-momentum dependence of the equations of motion of a dense neutrino gas in conditions similar to that of the early universe is effected in three stages: simplification of the geometry; restriction to two neutrino flavors; and approximation of the collision terms. The equations of motion that govern the evolution of the neutrino density matrix in the early universe depend on spacetime position and the neutrino three-momentum and energy. The equations are simplified by reducing the spatial dimensions of the universe from three to two. Since we are performing the linearized stability analysis at a given instant in the evolution of the system, we may further neglect the spacetime curvature, without loss of generality. Our assumption of elastic collisions with a background static array of nucleon-like scatterers allows the further assumption that neutrino energies are decoupled; energy transport is neglected. Thus we consider the neutrinos to be mono-energetic. This reduces the momentum dependence to that of a single, angular coordinate $\theta$. The resulting equations of motion are four-dimensional – two spatial dimensions, one momentum direction and time. We work with two active neutrino flavors, instead of three.

The purpose of restricting the dimensionality of spacetime to a planar spatial surface is to limit the size of the space of Fourier modes (see Eq. (6) below) that must be considered in the stability analysis. In fact, we expect that this approximation is not too severe in the sense that, all other approximations taken equally, a treatment considering the universe as three spatial dimensions would simply result in a larger number of Fourier modes. The stability behaviors of at least some of the modes in this case would be similar to those of the two-dimensional case since the spatial mode wave number couples explicitly only to the local velocity of the neutrino field. The two-flavor approximation also reduces the size of the space in which the linearized stability analysis is performed. The purpose of restricting the collision term to have only angular dependence in elastic scatterings of the type $\nu_{\alpha}N\to\nu_{\alpha}N$, where the field $N$ is represented by an infinitely massive, immobile object, is to allow consideration of just a single neutrino energy, which is conserved by this process. With these simplifications the universe can be viewed as a square of length $L$ on each side; we may take periodic boundary conditions without loss of generality.

The dense neutrino gas is described in terms of the density matrix in flavor space with elements $f_{\alpha\beta}(x,p)$, where $x$ and $p$ are position and momentum four-vectors, formally related to the Wigner transform of the neutrino correlation function $\langle\nu_{\alpha}(x)\bar{\nu}_{\beta}(y)\rangle$ in medium [42, 43, 44]. Factorization of the total population (flux) $n(\vec{p})$ in momentum bin $\vec{p}$, a scalar quantity, from the correlation function yields the $2\times 2$ density matrices $\rho^{\theta}(\vec{x},t)$ and $\bar{\rho}^{\theta}(\vec{x},t)$ for neutrinos and anti-neutrinos respectively, which describe the flavor content of the system; they satisfy ${\rm Tr}\rho^{\theta}=1={\rm Tr}\bar{\rho}^{\theta}$. Here $\vec{x}$ is a two-dimensional vector specifying the position, $\theta$ is the polar angle characterizing the direction of the neutrino momentum $\vec{p}$, with respect to a given, arbitrary direction (taken as the $y$-axis); $t$ is the time. The diagonal components of the density matrix represent the relative populations for the two neutrino flavors, denoted $\nu_{e}$ and $\nu_{\mu}$. The off-diagonal terms express correlations between the two flavors. The spacetime evolution equations of the density matrices are the quantum kinetic equations (QKEs) [45, 46, 47, 48, 49, 50, 51, 43, 52, 42, 53, 44], which are given, in the limit that we neglect incoherent scattering, by:

[TABLE]

where $\vec{v}=\dot{\vec{x}}$, $\vec{F}=\dot{\vec{p}}$. Here the Hamiltonian $\mathsf{H}$ for neutrino evolution and that for antineutrino evolution $\bar{\mathsf{H}}$ is written as a sum of terms

[TABLE]

where $\mathsf{H}_{0}$ is that part of the time-evolution operator independent of the density matrix, which is comprised of the vacuum and matter parts of the Hamiltonian. We ignore the term $\vec{F}\cdot\nabla_{p}$, which contributes when external forces, such as those induced by the Hubble expansion, are present. It does not contribute significantly in the early universe to the linearized stability analysis. The neutrino-neutrino evolution part of the Hamiltonian, which depends linearly on the density matrix, is denoted by $\mathsf{H}_{\nu\nu}$. These components are given explicitly as:

[TABLE]

Here $\mu\propto G_{F}n_{\nu}$ denotes the effective self-coupling strength, proportional to the neutrino number density, and $\alpha$ denotes the ratio $n_{\bar{\nu}}/n_{\nu}$. The linear, vacuum/matter term $\mathsf{H}_{0}$ is written for zero vacuum mixing angle ($\theta_{V}=0$) since $\theta_{V}$ does not alter the linearized stability analysis we perform in Sec. 4.222Including $\theta_{V}\neq 0$ would have two effects: (i) it would rescale the vacuum/matter frequency $\omega\to\omega\times\cos 2\theta_{V}$; (ii) it would introduce a source term $\sim\omega\times\sin 2\theta_{V}$ in the evolution of the homogeneous perturbations, which would lead to a linear growth with time, thus not impacting the stability analysis concerned with exponential growth. Physically, the effect (ii) corresponds to the onset of vacuum oscillations. We have carried out the integration over energy assuming a $\delta$-function distribution for energy. Unlike the case of supernova there is no matter suppression of flavor oscillations in a homogeneous and isotropic gas. We note, for later discussion, that the normal hierarchy case corresponds to $\omega>0$ and $\omega<0$ is the inverted case for this equivalent two-flavor system; we may refer to $\omega=\Delta m^{2}/(2E)$ as the ‘scaled’ vacuum frequency.

We are now in a position to investigate the role of collective neutrino oscillations in the early universe and whether the equations of motion of the neutrino density matrix are stable with respect to seed inhomogeneous perturbations. The analysis is carried out through the decomposition of the equations of motion of the density matrices [Eq. (1)] in terms of their Fourier modes,

[TABLE]

Here the sum is over $k_{i}=2\pi n_{i}/L$ for $n_{i}\in\mathbb{Z}$ and $i=x,y$. We note that previous studies have exclusively considered only the homogeneous mode $\vec{k}=(0,0)$ (in two spatial dimensions); all other modes $\vec{k}\neq\vec{0}$ are inhomogeneous. In the remainder of the present work we write the equations of motion [Eq. (1)] in terms of the Fourier modes $\rho^{\theta}_{\vec{k}}$. The objective of the linearized stability analysis is to identify those modes that grow exponentially with time and are therefore unstable to seed perturbations. Substituting Eq. (6) in Eq. (3) and taking the projection to the mode $\vec{k}$ gives

[TABLE]

Pairs of modes that appear in the last term of Eq. (7) sum to $\vec{k}$. Before considering the stability properties of the above collisionless equations in Sec. 4 we first discuss, in the following section, the general form of the collision term and the elastic, angular dependent approximation we make for the purposes of the present study.

3 The collision term

The QKEs presented in the previous section assumed that the neutrinos experience only the effects of forward scattering from other neutrinos. We must, however, take into account the effect of incoherent, direction changing scattering off other neutrinos, electrons and nucleons present in the early universe [45, 5, 9, 40, 41]. In this section, we briefly review the incorporation of momentum changing scattering to the neutrino flavor evolution. We continue to assume the reduced two-dimensional configuration space with monoenergetic neutrinos, which greatly simplifies the scattering matrices. The stipulation that neutrinos remain monoenergetic throughout the scattering requires that the scattering kernel only changes the direction of the momentum and not magnitude. We note that energy transport among neutrinos in scattering events is expected to contribute to additional suppression of the collective effects (see discussion in Sec. 6).

Previous studies have included the effect of incoherent scattering on neutrino oscillations in an approximate manner. Some references have included an off-diagonal damping term with Boltzmann collision integrals in the diagonal terms [54, 55, 56, 57]. (A recent study [29] has considered collision terms off-diagonal in flavor but neglected terms where neutrinos incoherently scatter from neutrinos and antineutrinos.) These treatments, which assume homogeneity and isotropy, have not considered the dependence of their results on the neutrino mass hierarchy. The mass hierarchy, however, plays a role when inhomogeneity and anisotropy are present.

In order to calculate the effect of collisions on collective neutrino oscillations we use the formalism developed in [44], albeit with a modified notation. The general form for the rate of change of the density matrix due to incoherent collisions is written as:

[TABLE]

The first term on the right-hand side encodes the flux loss for neutrinos traveling in the $\theta$ direction and the second term encodes the flux gain for a neutrino being scattered into the direction $\theta$. We note that both $\rho$ and $\Pi^{\rm gain/loss}[\rho]$ are $2\times 2$ matrices and gain and loss terms are expressed in term of anti-commutators denoted by $\left\{\ ,\ \right\}$. The factor $\left(1-\rho^{\theta}\right)$ in the second term is due to Pauli-blocking, a consequence of the fermionic nature of the neutrinos. Qualitatively speaking, this factor arises due to the degeneracy pressure of filled fermionic levels that disfavor scattering of neutrinos to the direction $\theta$ when $\rho^{\theta}$ is sufficiently close to 1.

The general form of the collision rate [Eq. (3)] is too complex to consider in the linearized stability analysis we pursue here. In order to make the problem tractable and to include the effect of angle dependent scattering in a basic way, we approximate the collision integral in two ways. First, we ignore the effect of Pauli blocking in the second term of Eq. (3). This will be problematic in regimes of high neutrino degeneracy but is not too drastic an approximation during the epoch of weak-decoupling/nucleosynthesis that we’re ultimately interested in understanding. As a further simplifying approximation, we consider only the neutrino-nucleon-like scattering processes $\nu_{\alpha}N\to\nu_{\alpha}N$, which is flavor $\alpha$ independent, neutral current scattering. Incorporating the aforementioned assumptions allows the reduction of the general form of scattering matrices given in Ref. [44]. Using Eqs. (3.21)-(3.22) (Eqs. (60)-(61)) in the e-print archive (published) version of Ref. [44], the form of the collision integral can be calculated by explicitly taking the inner product in 2+1 dimensional space. We find that the matrix-valued gain and loss collision contributions are proportional to the flavor-space identity matrix

[TABLE]

where, $c_{1}$ and $c_{2}$ are constants. In three dimensional space one has $c_{1}\propto G_{F}^{2}n_{\textrm{nuc}}p^{2}$ with proportionality constant of order unity. For equal number of protons and neutrons in the early Universe the ratio $c_{1}/c_{2}$ is fixed and is equal to 4.83, a value which we use for rest of the paper.

We note that the above form for the approximate collision rate in Eq. (9) is energy-local; that is, there is no energy transport. And it is angle dependent, being a function of the angle of scattering $\theta^{\prime}$ relative to the incoming neutrino angle $\theta$.

Combining Eqs. (3) and (9) we obtain the complete equation describing the evolution of flavor density matrices in presence of collective effects and collisions.

4 Linear stability analysis

Nonlinear neutrino self-interaction in the equations of motion for the density matrix [Eqs. (1)] presents a computationally challenging scenario for the general case, particularly when $\vec{k}\neq 0$. In the presence of seed spatial inhomogeneity perturbations, the reduction of the number of independent variables by symmetry is no longer possible. It is, however, possible to gain insight into the stability properties of Eqs. (1) through a linear stability analysis [58]. In this section we describe the formalism and later, in Sec. 6, we discuss the reliability and limitations of linear stability analyses in the context of collective neutrino oscillations.

The linear stability analysis of collective neutrino oscillations is predicated on the basis that, for significant neutrino flavor oscillations to occur, the off-diagonal elements of the Hamiltonian should be comparable in magnitude to the diagonal elements. Since $\mathsf{H}_{0}$ is diagonal (neglecting vacuum mixing), this can occur if the off-diagonal component of density matrices increase rapidly with time.

In order to investigate the conditions (parameter ranges) under which exponential growth of the Fourier modes may occur, we assume that the density matrices have a form describing the presence of seed inhomogeneity perturbations ($\Delta^{\theta}(t,\vec{x})$ and $\epsilon^{\theta}(t,\vec{x})$ for neutrinos and $\bar{\Delta}^{\theta}(t,\vec{x})$ and $\bar{\epsilon}^{\theta}(t,\vec{x})$ for antineutrinos):

[TABLE]

We leave the origin of these seed inhomogeneity perturbations unspecified. They will be present in generic theories of the early universe unless expressly forbidden by some exact symmetry.

The perturbations in the diagonal terms, $\Delta^{\theta}(t,\vec{x})$ and $\bar{\Delta}^{\theta}(t,\vec{x})$ can be ignored upon linearization of the equations of motion, because they contribute only at second order. In terms of the Fourier modes, Eq. (10) reads

[TABLE]

The equations for $\epsilon_{\vec{k}}$ for the various Fourier modes $\vec{k}$ are decoupled at linear order in the perturbations $\epsilon^{\theta}_{\vec{k}}$, $\bar{}\epsilon^{\theta}_{\vec{k}}$:

[TABLE]

For a given set of parameters, exponential growth of $\epsilon$ and $\bar{\epsilon}$ implies flavor instability. The dependence of collective neutrino oscillations on hierarchy can be viewed as being a result of the relative sign between the second term on right-hand side of Eq. (4), which is proportional to the vacuum term, and the third and forth terms, which is proportional to the self-interaction strength $\mu$. However, it should be noted that the relative sign between the vacuum and self-interaction term can also be changed by going from $\alpha<1$ to $\alpha>1$. All the results therefore for $1-\alpha>0$ with normal hierarchy are identical to those for $1-\alpha<0$ with inverted hierarchy. There are two different ways in which we can investigate the flavor instability of the system. One is the brute-force method in which we discretize the angular dependence of $\epsilon$ and $\bar{\epsilon}$ and solve the resulting system of coupled differential equations. A major drawback of this approach is that sometimes we get ‘spurious’ flavor instabilities as a result of discretization which reduce as we increase the number of angle-bins[2].

The second method is a bit more formal but does not have the problem of spurious instabilities. Based on the fact that we are investigating collective modes, we seek solutions to Eqs. (4) of the form

[TABLE]

Existence of a solution with $\kappa>0$, where

[TABLE]

implies the existence of flavor instabilities. Combing Eq. (4) and Eq. (13) we obtain

[TABLE]

for the neutrino modes, and

[TABLE]

for the antineutrino modes.

Recognizing that the solutions to the above equations for the collective mode amplitudes $Q^{\theta}_{\vec{k}}$ and $\overline{Q}^{\theta}_{\vec{k}}$ are spanned by the functions $\{1,\cos\theta,\sin\theta\}$ leads to the ansatz

[TABLE]

where we have defined the collective mode functions $a_{\vec{k}}$, corresponding to modes independent of the angle $\theta$ and the angular-dependent modes $s_{\vec{k}}$ and $c_{\vec{k}}$ (and their barred counterparts). We have abbreviated the denominators as:

[TABLE]

The collective mode functions are projected from Eqs. (4) and (4) as

[TABLE]

Consistency among these relations yields linear equations for the collective modes, given by:

[TABLE]

in terms of the functionals

[TABLE]

with $J_{1,2}$ defined as

[TABLE]

We determine the values of $\Omega_{\vec{k}}$ for which Eq. (23) is satisfied by requiring the determinant of the $6\times 6$ matrix vanishes. As the determinant of the matrix is a $6^{\textrm{th}}$-order polynomial we obtain up to six different values of $\Omega_{\vec{k}}$ for each value of $\vec{k}$. While multiple solutions are possible, the instability is dominated by the $\Omega_{\vec{k}}$ with largest positive imaginary part. We have numerically solved for $\Omega_{\vec{k}}$ by using standard root finding numerical algorithms. For the results presented in the next section we use the open source GNU Scientific Library (GSL) [59].

Note that for $\vec{k}=0$ (homogeneous) modes the matrix in Eq. (23) has a block diagonal form with three $2\times 2$ blocks corresponding to subspaces spanned by the pairs of mode functions $\{a_{\vec{k}}$, $\bar{a}_{\vec{k}}\}$, $\{s_{\vec{k}}$, $\bar{s}_{\vec{k}}\}$, and $\{c_{\vec{k}}$, $\bar{c}_{\vec{k}}\}$; these subspaces are decoupled. The latter two subspaces are identical (those for the $s_{\vec{k}}$, $\bar{s}_{\vec{k}}$ and $c_{\vec{k}}$, $\bar{c}_{\vec{k}}$ blocks). In this case, the determinant can be calculated analytically and the solutions for $\Omega_{0}=\Omega_{\vec{k}=0}$ are given by

[TABLE]

for the $a_{\vec{k}=0}=a_{0},\bar{a}_{\vec{k}=0}=\bar{a}_{0}$ block, and

[TABLE]

for the $s_{0}$, $\bar{s}_{0}$ and $c_{0}$, $\bar{c}_{0}$ blocks. Despite the fact that the $\vec{k}=0$ modes are homogeneous, it is instructive to study these results for this $\vec{k}=0$ mode since similar behavior is observed for the inhomogeneous $\vec{k}\neq 0$ modes.

For the inverted hierarchy ($\omega<0$), in absence of collisions the angle-dependent collective frequencies $\Omega_{0,\pm}^{(s,c)}$ are real and therefore do not give rise to exponential growth. On the other hand, Eqs.(27) indicate that for the inverted hierarchy ($\omega<0$), the argument of the square-root $\kappa^{2}=\omega^{2}-2\pi\mu(1+\alpha)|\omega|+\pi^{2}\mu^{2}(1-\alpha)^{2}$ becomes negative for values

[TABLE]

where we have defined $\alpha=1-\delta$ and expanded for $\delta\ll 1$ (in our numerical studies we take $\delta=10^{-2}$). Figure 1 shows $\kappa/|\omega|$ plotted against the neutrino effective self-coupling $\mu/|\omega|$; the left panel shows two $\vec{k}$ modes (red and green curves) for the normal (left panel) and inverted (right panel) hierarchies. The effective coupling inequalities for the inverted hierarchy described by Eqs. (29) are clearly exhibited in the right panel of Fig. 1 by the curve corresponding to $\vec{k}=0$ (the red curve). The left-most inequality in Eqs. (29) corresponds to the region $\mu/|\omega|\approx 0$; the right-most inequality corresponds to where the $|\vec{k}|=0$ intersects the abscissa near $\mu/|\omega|\simeq 1.27\times 10^{4}$.

Equation (27) demonstrates that for inverted hierarchy ($\omega<0$) only $\theta$-independent instabilities can develop (proportional to $a_{\vec{k}}$ and $\bar{a}_{\vec{k}}$). Moreover, and perhaps not too surprisingly, the direction-changing form for the elastic collision rate we employ [Eq. (9)] in this exploratory calculation has no effect on such instabilities: the strength of the collisions, as measured by the quantities $c_{1,2}$, do not feature in the expression for $\Omega_{0,\pm}^{(a)}$. This is potentially an important observation since, if this behavior obtains for the full, unapproximated collision rate then nonlinear mode coupling effects, present in the full QKEs, could propagate the instability to other, inhomogeneous modes. This effect could result in a continuous “feeding” of inhomogeneity that may persist for times long compared to other dynamical time scales, despite collisional damping in the inhomogeneous modes.

For the normal hierarchy ($\omega>0$), the angle-independent collective frequencies $\Omega_{0,\pm}^{(a)}$ are real and therefore do not give rise to exponential growth. On the other hand, $\theta$-dependent instabilities (proportional to $s_{0},c_{0}$ and $\bar{s}_{0},\bar{c}_{0}$) described by Eq. (4) may develop. In this case, the growth rate due to the imaginary part of $\Omega^{(s,c)}_{0,\pm}$ is indeed damped by the direction-changing collisions, as shown by the term proportional to $2c_{1}+c_{2}$ in Eq. (4). The numerical results discussed in the next section match the analytic results at $\vec{k}=0$ and show that a similar behavior persists for $\vec{k}\neq 0$: for normal hierarchy $\theta$-dependent instabilities ($s_{\vec{k}}$, $\bar{s}_{\vec{k}}$ and $c_{\vec{k}}$, $\bar{c}_{\vec{k}}$) can develop, while for the inverted hierarchy only the $\theta$-independent modes $a_{\vec{k}}$ and $\bar{a}_{\vec{k}}$ are subject to unstable, exponential growth.

5 Results

In this section we present numerical results for the stability analysis of collective modes with $\vec{k}\neq 0$ corresponding to inhomogeneous flavor distributions in configuration space. We continue to take the ratio of electron anti-neutrinos to electron neutrinos, $\alpha=0.99$, as in the previous section. As shown in Fig. 2, for the collisionless case ($c_{1}=0=c_{2}$), the inhomogeneous modes with $|\vec{k}|=1000\omega$ (we use $n=10^{4}$ and $2\pi/L=0.1\omega$) are unstable in both the normal (left panel of the figure) and inverted (right panel) mass hierarchies. The effective neutrino self-coupling, $\mu$ over which instabilities exist, we note, is different for normal and inverted hierarchies. However, for both mass hierarchies the small scale inhomogeneous modes are unstable over a larger range of $\mu$ than for the homogeneous modes.

In the presence of collisions we would generally expect a suppression of instabilities for sufficiently strong interaction rates. This expectation is driven principally by the observation that incoherent scattering of the form of Eq. (9), due to its angular dependence, will tend to redistribute regions of neutrino density toward the mean. The expectation of the suppression is supported by the form of Eq. (4), specifically the term $-i\pi(2c_{1}+c_{2})$. Since $c_{1,2}>0$, this term contributes an exponentially damped factor to the collective mode frequency $\Omega^{(s,c)}_{0,\pm}$.

Indeed, for large values of $k\sim 1000\omega$ we find suppression of instabilities for both mass hierarchies, as shown in Fig. 2. As in Fig. 1, the left and right panels give results for the values of $\kappa/\omega$ for the normal and inverted hierarchies, respectively. In each hierarchy, the effect of varying the strength of the collision term, measured by the magnitude of $c_{1}>0$, is shown by the various curves, as detailed in the figure caption. In each panel, we observe a reduction in the width of the region of instability with increasing $c_{1}$.

However, as previously indicated by Eq. (27), for smaller values of $k$ ($k\sim$ few $\times 100\omega$) including zero, we find a suppression due to collisions only in the case of normal mass hierarchy, not in the inverted mass hierarchy. Results for the $\vec{k}=0$ homogeneous modes have been obtained both analytically, through Eq. (27) and numerically via diagonalization of the $6\times 6$ matrix appearing in Eq. (23). Figure 3 shows the numerical results for the homogeneous mode for various values of the incoherent scattering strength parameter $c_{1}$; it agrees with Eq. (27) at high precision for all values of $c_{1}$ considered in this work.

6 Conclusion

We have studied the stability of a two-dimensional dense neutrino gas with respect to spatial inhomogeneity. In the absence of incoherent collisional effects we have found that the system exhibits growth of seed inhomogeneity due to nonlinear coherent neutrino self-interactions and that the effect depends on the neutrino mass spectrum. For normal hierarchy spatial instability exists over a larger range of neutrino number density compared to the inverted case.

We have further considered the effect of elastic (i.e. energy conserving) incoherent (direction-changing) collisions of the neutrinos with a static background of heavy, nucleon-like scatterers, in order to validate the intuitive expectation that incoherent collisions tend to drive the neutrino flavor field of the early universe toward homogeneous isotropic distributions. Our results suggest that this heuristic picture may be oversimplified.

At small scales, corresponding to Fourier modes $|\vec{k}|\gtrsim 100\omega$, the growth of flavor instability can be suppressed by collisions, irrespective of the neutrino mass spectrum. At large length scales ($|\vec{k}|\lesssim 100\omega$) we find, perhaps surprisingly, that for inverted neutrino mass hierarchy direction-changing collisions fail to suppress flavor instabilities, irrespective of the coupling strength. In the extreme long-wavelength case ($\vec{k}=0$) our numerical results are substantiated by an analytic understanding of the problem (see Section 4). The key to understanding the puzzling behavior is that for the simplified two-dimensional model, at large scale ($|\vec{k}|\sim 0$) only $\theta$-independent modes grow unstable for inverted hierarchy, while only $\theta$-dependent modes grow unstable for normal hierarchy.

The appearance of inhomogeneity due to collective neutrino oscillations in the early universe could have implications for the entropy evolution in the neutrino sector, which in turn could affect Big Bang nucleosynthesis. Moreover, the non-trivial dependence of the instability (and its collisional suppression) on the neutrino mass hierarchy is very intriguing and deserves attention in future studies involving more realistic models.

One obvious step towards a more realistic setup would involve generalizing the collision kernels in order to account for energy-changing collisions in our formalism. This would allow one to model neutrino interactions with electrons / positrons in the early universe, which is beyond the scope of this Letter. However, even in absence of a complete calculation, since the collisional suppression of neutrino flavor oscillations gets contributions from direction change as well as energy change, we speculate that the total effect of collisions on the suppression of neutrino inhomogeneities should still be dependent on the neutrino mass hierarchy.

Acknowledgments

The authors would like to thank George Fuller for offering valuable insights. The authors would also like to thank Sajad Abbar, Daniel Blaschke, Joe Carlson, Huaiyu Duan, Evan Grohs, Luke Johns, Lei Ma and Joshua Martin for numerous discussions. The authors acknowledge support by the LDRD program at Los Alamos National Laboratory.