Fast flavor conversions of supernova neutrinos: Classifying instabilities via dispersion relations

TL;DR

This paper introduces a systematic method, based on dispersion relations and plasma physics techniques, to classify and analyze fast flavor instabilities of supernova neutrinos, providing both qualitative insights and quantitative predictions.

Contribution

It applies a plasma physics-inspired classification of flavor instabilities to supernova neutrinos, enabling detailed analysis of their growth and stability properties.

Findings

Instabilities are characterized by complex frequency and wave number behaviors.

Stable modes occur when either both frequencies and wave numbers are real or when only wave numbers are complex.

Supernova conditions with lepton number crossings can lead to absolute flavor instabilities.

Abstract

Supernova neutrinos can exhibit a rich variety of flavor conversion mechanisms. In particular, they can experience "fast" self-induced flavor conversions almost immediately above the core. Very recently, a novel method has been proposed to investigate these phenomena, in terms of the dispersion relation for the complex frequency and wave number ($\omega$,$k$) of disturbances in the mean field of the $\nu_e\nu_x$ flavor coherence. We discuss a systematic approach to such instabilities, originally developed in the context of plasma physics, and based of the time-asymptotic behavior of the Green's function of the system. Instabilities are typically seen to emerge for complex $\omega$, and can be further characterized as convective (moving away faster than they spread) and absolute (growing locally), depending on $k$-dependent features. Stable cases emerge when $k$ (but not $\omega$) is complex, leading to disturbances damped in space, or when both $\omega$ and $k$ are real, corresponding to complete stability. The analytical classification of both unstable and stable modes leads not only to qualitative insights about their features but also to quantitative predictions about the growth rates of instabilities. Representative numerical solutions are discussed in a simple two-beam model of interacting neutrinos. As an application, we argue that supernova and binary neutron star mergers exhibiting a "crossing" in the electron lepton number would lead to an absolute instability in the flavor content of the neutrino gas.

Equations

$$S\sim e^{i(\mathbf{k}\cdot \mathbf{x}-\omega t)}.$$

$${\partial }_t{\varrho }_{\mathbf{p},\mathbf{x},t}+{\mathbf{v}}_{\mathbf{p}}\cdot {\nabla }_{\mathbf{x}}{\varrho }_{\mathbf{p},\mathbf{x},t}=-i[{\Omega }_{\mathbf{p},\mathbf{x},t},{\varrho }_{\mathbf{p},\mathbf{x},t}],$$

$${\Omega }_{\mathbf{p}}={\Omega }_{\mathrm{vac}}+{\Omega }_{\mathrm{MSW}}+{\Omega }_{\nu \nu },$$

$${\Omega }_{\mathrm{MSW}}=\lambda \text{ diag}(1,0),$$

$${\Omega }_{\nu \nu }=\sqrt{2}{G}_F\int \frac{d^3\mathbf{q}}{(2\pi {)}^3}({\varrho }_{\mathbf{q}}-{\bar{\varrho }}_{\mathbf{q}})(1-{\mathbf{v}}_{\mathbf{p}}\cdot {\mathbf{v}}_{\mathbf{q}}),$$

$$\varrho=\frac{f_{\nu_{e}}+f_{\nu_{x}}}{2}\left(\begin{array}[]{cc}1&0\\ 0&1\end{array}\right)+\frac{f_{\nu_{e}}-f_{\nu_{x}}}{2}\left(\begin{array}[]{cc}s&S\\ S^{*}&-s\end{array}\right)\,\ ,$$

$$i({\partial }_t+\mathbf{v}\cdot {\nabla }_{\mathbf{x}})S_{\mathbf{v}}$$

$$G_{\mathbf{v}}=\sqrt{2}{G}_F{\int }_0^{\infty }\frac{dE{E}^2}{2{\pi }^2}\left[f_{{\nu }_e}(E,\mathbf{v})-f_{{\bar{\nu }}_e}(E,\mathbf{v})\right],$$

$$\mu =\sqrt{2}{G}_F(n_{{\nu }_e}-n_{{\bar{\nu }}_e})\simeq 6{\text{m}}^{-1}\frac{(n_{{\nu }_e}-n_{{\bar{\nu }}_e})}{10^{31}{\text{cm}}^{-3}},$$

$${\bf v}=\Big{(}\sqrt{1-v_{z}^{2}}\cos\varphi,\sqrt{1-v_{z}^{2}}\sin\varphi,v_{z}\Big{)}\,\ ,$$

$$S_{\mathbf{v}}(t,\mathbf{x})=Q_{\mathbf{v}}{e}^{i(\mathbf{k}\cdot \mathbf{x}-\omega t)}.$$

$$S_{\mathbf{v}}(t,\mathbf{x})\to S_{\mathbf{v}}(t,\mathbf{x})\exp \left[-it\left(\lambda +\int d{\Gamma }^{\prime }{G}_{{\mathbf{v}}^{\prime }}\right)+i\int d{\Gamma }^{\prime }\mathbf{x}\cdot {\mathbf{v}}^{\prime }{G}_{{\mathbf{v}}^{\prime }}\right],$$

$$(\omega -\mathbf{k}\cdot \mathbf{v})Q_{\mathbf{v}}=-\int d{\Gamma }^{\prime }(1-\mathbf{v}\cdot {\mathbf{v}}^{\prime })G_{{\mathbf{v}}^{\prime }}{Q}_{{\mathbf{v}}^{\prime }}.$$

$$\det \left[\begin{array}{cccc}1+\int d\Gamma \frac{G_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& \int d\Gamma \frac{v_x{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& \int d\Gamma \frac{v_y{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& \int d\Gamma \frac{v_z{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}\\ \int d\Gamma \frac{v_x{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& -1+\int d\Gamma \frac{v_x^2{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& \int d\Gamma \frac{v_x{v}_y{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& \int d\Gamma \frac{v_x{v}_z{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}\\ \int d\Gamma \frac{v_y{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& \int d\Gamma \frac{v_x{v}_y{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& -1+\int d\Gamma \frac{v_y^2{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& \int d\Gamma \frac{v_y{v}_z{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}\\ \int d\Gamma \frac{v_z{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& \int d\Gamma \frac{v_x{v}_z{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& \int d\Gamma \frac{v_y{v}_z{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}& -1+\int d\Gamma \frac{v_z^2{G}_{\mathbf{v}}}{\omega -\mathbf{v}\cdot \mathbf{k}}\end{array}\right]=0.$$

$$\text{det}\left[{\Pi }^{\mu \nu }(\omega ,\mathbf{k})\right]=0,$$

$${\Pi }^{\mu \nu }={\eta }^{\mu \nu }+\int \frac{d\mathbf{v}}{4\pi }{G}_{\mathbf{v}}\frac{v^{\mu }{v}^{\nu }}{\omega -\mathbf{k}\cdot \mathbf{v}},$$

$$\mathbf{x}=z,\mathbf{k}=k,$$

$$D\left(i\frac{\partial }{\partial t},-i\frac{\partial }{\partial z}\right)S(z,t)=0.$$

$$D(\omega ,k)=0,$$

$$\omega =\Omega (k)\in \mathbb{C}\mathrm{with}k\in \mathbb{R},$$

$$k=K(\omega )\in \mathbb{C}\mathrm{with}\omega \in \mathbb{R}.$$

$$S(z,t)=\int dk{g}_k{e}^{i(kz-\Omega (k)t)},$$

$$S(z,t)=\int d\omega f_{\omega }{e}^{i(K(\omega )z-\omega t)},$$

$$D\left(i\frac{\partial }{\partial t},-i\frac{\partial }{\partial z}\right)G(z,t)=\delta (z)\delta (t).$$

$$D(\omega ,k)G(\omega ,k)=1,$$

$$G(z,t)=\frac{1}{(2\pi {)}^2}{\int }_{L_{\omega }}{\int }_{F_k}d\omega dk\frac{e^{i(kz-\omega t)}}{D(\omega ,k)},$$

$$G(z,t)=\frac{-i}{2\pi }{\int }_{F_k}dk\frac{e^{i(kz-\Omega (k)t)}}{{\frac{\partial D}{\partial \omega }|}_k},$$

$$V=\frac{z}{t}=\mathrm{const}$$

$$V={\frac{\partial \Omega }{\partial k}|}_{\hat{k}}.$$

$$G(x,t)\sim -\frac{e^{i\frac{\pi }{4}}}{(2\pi {)}^{\frac{1}{2}}}\frac{e^{i(\hat{k}z-\Omega (\hat{k}))t}}{{\frac{\partial D}{\partial \omega }|}_{\hat{k}}{\left({\frac{{\partial }^2\Omega }{\partial k^2}|}_{\hat{k}}\cdot t\right)}^{\frac{1}{2}}},$$

Full text

I Introduction

Flavor conversions of supernova (SN) neutrinos are a field of intense investigation. Conversions in the deepest SN regions would have a potential impact on the explosion dynamics (e.g., the revival of the shock-wave) and on the nucleosynthesis of heavy elements in the stellar matter. Flavor transition effects, occurring during neutrino propagation from inner to outer SN regions, would also imprint peculiar features on observable spectra in the next Galactic SN neutrino burst (see [1] for a recent review). It is expected that different regimes of flavor conversions $\nu_{e}\to\nu_{x}$ ($x=\mu,\,\tau$) would be encountered by SN neutrino with energy $E$ and squared mass difference $\Delta m^{2}$ while they propagate through background fermions with density $n$. Early studies focussed on the so-called Mikheyev-Smirnov-Wolfenstein (MSW) matter effects [2, 3], occurring when the neutrino oscillation frequency in vacuum $\omega_{\rm vac}=\Delta m^{2}/(2E)$ is comparable to the matter potential $\lambda=\sqrt{2}G_{F}n_{e}$ induced by the net electron density $n_{e}$. This condition, tipically satisfied at “large” distances $r\sim{\mathcal{O}}(10^{4})$ km from the SN core, can produce resonant flavor conversions [4] that carry information on the SN matter profile, including features related to the shock wave [5, 6, 7, 8] and turbulence [9, 10, 11]. Later it was realized that also the SN neutrino density $n_{\nu}$ can be so high to provide a background medium inducing neutrino flavor change [12, 13]. When the neutrino-neutrino interaction potential $\mu\sim\sqrt{2}G_{F}n_{\nu}$ is dominant, the neutrino flavor evolution becomes nonlinear [14]. First large-scale numerical studies of these effects [15, 16], supported by analytical insights [17], stimulated a torrent of still ongoing activities. Through these studies it was realized that neutrino-neutrino interactions can lead to “self-induced” or “collective” flavor conversions at typical radii $r\sim{\mathcal{O}}(10^{2}-10^{3})$ km and frequencies $\sim(\omega_{\rm vac}\mu)^{1/2}$, leading to spectral swaps and splits [16, 18, 19, 20, 21]. The robustness of these effects has been recently questioned, since relaxing some of the symmetry assumptions imposed on simplified emission models [22, 23, 24, 25, 26, 27, 28, 29, 30] can lead to dramatic changes and even to flavor decoherence (see [31] for a recent review). In any case, it was implicitly assumed in the literature that all these effects would occur relatively far from the neutrino emission region, being suppressed by dominant matter terms a closer radii [32, 33, 34, 35].

In contrast, little attention was devoted to the possibility of fast flavor conversions (with frequencies as large $\sim\mu$) at much shorter radii (with $r\sim{\mathcal{O}}(1)$ m from the neutrino sphere) highlighted in [36] and developed in simple toy models [36, 37, 38]. Surprisingly, these conversions would be so fast ($\mu\gg\omega_{\rm vac}$) that they would wipe out any dependence on the vacuum mass-mixing parameters whose role is just to provide a “seed” for the self-induced flavor dynamics. The key insight [36, 37, 38] is that near the neutrino decoupling region the angular distributions of the different neutrino species are rather different and this would cause a speed-up of the flavor conversions. Notably, non-electron species $\nu_{x}$ would decouple earlier (and thus would be more forward peaked) than the electron species $\nu_{e}$ and $\bar{\nu}_{e}$ and. In turn, due to the neutron richness of stellar matter, the $\bar{\nu}_{e}$ would decouple earlier (and thus would be more forward peaked) than ${\nu}_{e}$. On the other hand, all other studies on self-induced flavor conversions had assumed the same (half-isotropic) angular distribution for any species. This approximation is reasonable for slow flavor conversions far from the neutrinosphere111See however [39, 40] for a discussion of the impact of non-trivial angular distributions on slow self-induced conversions. but misses a crucial ingredient of fast conversions, namely, backward-going modes. Recently, fast flavor conversions have been independently studied by some authors [41, 42] finding that these can be possible near the SN core if the neutrino flux ratios and angular asymmetries produce a crossing between the zenith-angle spectra of $\nu_{e}$ and ${\bar{\nu}}_{e}$. In particular, the presence of neutrinos traveling towards the core can make these fast conversions more generic [42]. Fast self-induced flavor conversions are still in an exploratory phase. However, their characterization is crucial in order to have a reliable description of the neutrino flavor evolution during a stellar collapse.

In the emerging literature about fast flavor conversions, a novel and very interesting approach to study these effects was recently proposed in [43]. This approach is based on the dispersion relation for the frequency and wavenumber $(\omega,{\bf k})$ in the mean field of $\nu_{e}\nu_{x}$ coherence, which is essentially the off-diagonal element of the neutrino density matrix $\varrho({\bf p},{\bf x},t)$, that we will call $S$ in the following. One looks for solutions of the linearized equations for the flavor evolution in the form

[TABLE]

Typically such a solution may exist only if $\omega$ and ${\bf k}$ are related by an appropriate equation, called the dispersion relation. If either ${\bf k}$ or $\omega$ develop imaginary parts leading to positive real arguments in the exponential, the solution is expected to blow in space or time, thus signaling an “instability”. Of course, the growth is not exponential forever, as it would break the linear approximation. It is, however, appropriate to diagnose the onset of instabilities by studying linear stability. It should be emphasized that not all modes with complex ${\bf k}$ or $\omega$ correspond to an instability: the appearance of a positive real term in the exponential argument of Eq. (1) would be a naive and not sufficient criterion. The purpose of this paper is to discuss more refined criteria to assess the onset of instabilities and to classify them, building upon the dispersion relation approach proposed in [43] for SN neutrinos, and on formally analogous methods developed for unstable phenomena in other contexts, such as plasma physics and fluid dynamics [44]. The conditions for fast flavor conversions will be elucidated by means of both analytical calculations and numerical simulations.

The outline of our paper is as follows. In Sec. 2 we present the equations of motion of the dense neutrino gas in their complete but linearized form. In Sec. 3 we discuss the dispersion relation for the neutrino ensemble and present a kinematical classification of the possible instabilities of the system. In Sec. 4 we present the theory of the instabilities, based on the asymptotic properties of the Green’s function of the problem. In Sec. 5 we apply this general theory to a simplified two-beam model. We also interpret the possible dispersion relations in terms of “particle-like” vs “tachyon-like’ behaviors. We compare our analytical predictions with numerical solutions of the linearized equations of motion. Finally in Sec. 6 we summarize our results and discuss the potential consequences for flavor conversions in SNe and neutron star mergers.

II Flavor evolution and dispersion equation at high neutrino density

In the absence of external forces and collisions, the dynamics of the space-dependent $\nu$ occupation numbers or Wigner function $\varrho_{{\bf p},{\bf x},t}$ with momentum ${\bf p}$ at position ${\bf x}$ and time $t$ is ruled by the following equations of motion (EoMs) [45, 46]

[TABLE]

where, in the Liouville operator on the left-hand side, the first term accounts for explicit time dependence, while the second term, proportional to the neutrino velocity ${\bf v_{p}}$, encodes the spatial dependence due to particle free streaming. In the absence of oscillations, the right-hand-side would be zero and the EoM would reduce to a Vlasov continuity equation. In the presence of oscillations, the matrix $\Omega_{{\bf p}}$ is the Hamiltonian

[TABLE]

containing the vacuum, matter and self-interaction terms, that leads to the evolution of $\varrho_{{\bf p}}$ over space and time. Here and in the following, to lighten our notation, we drop the subscripts ${\bf x}$ and $t$.

For our purposes, Eq. (3) can be simplified in an effective two-generation scenario, since only two flavors $(\nu_{e},\,\nu_{x})$ are involved and, moreover, the relevant dynamics near the supernova core will be shown to be independent on the vacuum oscillation parameters which just trigger (but do not govern) fast conversions. The vacuum term is thus the matrix $\Omega_{{\rm vac}}=\textrm{ diag}(-\Delta m^{2}/4E,+\Delta m^{2}/4E)$ in the mass basis, where $E=|{\bf p}|$ for ultra-relativistic neutrinos. For antineutrinos, the EoMs are the same but with the replacement $\Omega_{{\rm vac}}\to-\Omega_{{\rm vac}}$, thus it is convenient to think of antineutrinos of energy $E$ as neutrinos of energy $-E$, making their EoMs identical. The matter term in Eq. (3) reads

[TABLE]

in the weak interaction basis, where $\lambda=\sqrt{2}G_{F}n_{e}$. Finally, the term due to $\nu-\nu$ interactions is given by

[TABLE]

where the angular factor $(1-{\bf v}_{\bf p}\cdot{\bf v}_{\bf q})$ leads to multi-angle effects [15], as neutrinos moving on different trajectories experience different self-interaction potentials. Concerning the density matrix $\varrho_{\bf p}$ in two-flavor case we can write it in the weak-interaction basis as

[TABLE]

where $f_{\nu_{e}}$ and $f_{\nu_{x}}$ are the initial occupation numbers. The complex scalar field $S_{\bf p}(t,{\bf x})$ represents the $\nu_{e}\nu_{x}$ flavor coherence for the mode ${\bf p}$, while the real field $s_{\bf p}(t,{\bf x})$ satisfies $s_{\bf p}^{2}+|S_{\bf p}|^{2}=1$. Note that since we are assuming the ${\bar{\nu}}$ have negative energy and negative $\varrho$, the ${\bar{\nu}}$ coefficients are $-(f_{\bar{\nu}_{e}}+f_{\bar{\nu}_{x}})/{2}$ and $-(f_{\bar{\nu}_{e}}-f_{\bar{\nu}_{x}})/{2}$. Neutrinos are produced as flavor eigenstates and no flavor mixing occurs as long as $S_{\bf p}=0$.

While the self-induced flavor evolution described by Eq. (2) is a nonlinear problem, the onset of these conversions can be examined by linearizing the equations, observing that $|S|\ll 1$ initially and that $s=\sqrt{1-|S|^{2}}\simeq 1$ to linear order in $S$. In the context of fast conversions (where $\mu\gg\omega_{\rm vac}$ dominates), the $\Omega_{{\rm vac}}$ term just produces the initial seed $|S|\neq 0$ and then becomes irrelevant, so that we can take $\Omega_{{\rm vac}}=0$ hereafter. In this limit, the energy $E$ disappears from the EoMs, and the modes of $S$ can be labeled by the same unit vector ${\bf v}$ (with $|{\bf v}|\approx c=1$) for both $\nu$ and $\bar{\nu}$.

Assuming that the occupation numbers as well as the matter density are homogeneous and stationary within the test volume, one gets the following linearized equations

[TABLE]

where $d\Gamma=d{\bf v}/4\pi$ and $G_{\bf v}$ is the angle distribution of the electron lepton number (ELN) carried by neutrinos [43], namely

[TABLE]

where the $f$’s represent the neutrino distributions in energy and in emission angle. The corresponding ELN potential in Eq. (7) is given by $\mu=\int d\Gamma G_{\bf v}$, with an associated current ${\mathbf{j}}=\int d\Gamma G_{\bf v}{\bf v}$. For typical SN neutrino densities [34] one numerically finds

[TABLE]

corresponding to a time scale of ${\mathcal{O}}(1)$ ns.

Note that in the previous equation we are assuming $f_{\nu_{x}}(E,{\bf v})=f_{\bar{\nu}_{x}}(E,{\bf v})$, otherwise one should add a term $-[f_{\nu_{x}}(E,{\bf v})-f_{\bar{\nu}_{x}}(E,{\bf v})]$.

A further simplification occurs in studying flavor conversions at distances from the SN neutrinosphere much smaller than its radius of ${\cal{O}}(10)\,{\rm km}$, since the curvature of the emitting surface becomes irrelevant. Then one can simply use Cartesian components for the velocity,

[TABLE]

where $v_{z}=\cos\vartheta$ is the component along the $z$-axis, and $\vartheta$ and $\varphi$ the zenith and azimuthal angles, respectively. Note that $v_{z}$ can take negative values, i.e., the zenith angle $\vartheta$ can take values from [math] up to $\pi$, representing neutrinos with trajectories that range from radially outwards to radially inwards, and not merely up to $\pi/2$ as usually taken in the “bulb” model. With the same logic, the background matter and neutrino density can be locally assumed as constant.

III Dispersion relation and classification of the instabilities

III.1 Dispersion relation

Starting from the linearized equations of motion [Eq. (7)], one can proceed by seeking solutions via the ansatz

[TABLE]

After a global phase transformation

[TABLE]

to remove the matter and neutrino density terms in the square-bracketed term on the r.h.s. of Eq. (7) (see [43]), one obtains

[TABLE]

While this transformation makes the equation look simpler, one must be careful that the modes labelled by ${\bf k}=0$ or $\omega=0$ after the transformation do not correspond to what were homogeneous or stationary modes.

The above equation has solutions only if

[TABLE]

Note that in [43] the previous equation has been written as

[TABLE]

in terms of a “polarization tensor”

[TABLE]

with $v^{\mu}=(1,{\bf v})$ and $\eta^{\mu\nu}=\textrm{diag}(+1,-1,-1,-1)$ being the metric tensor.

Thus, wave-like solutions as in Eq. (11) exist only if $\omega$ and ${\bf k}$ (“normal modes”) are appropriately related as per Eq. (14) [or, equivalently, Eq. (15)], which is called the “dispersion relation” of the system. For the sake of simplicity, hereafter we shall assume that the relevant dynamics involves only one space coordinate $z$ and its conjugate $k$,

[TABLE]

and that instabilities do not grow in $(x,\,y)$ plane orthogonal to the flow along $z$. The study of multidimensional instabilities is beyond the scope of this work.

III.2 Types of Instabilities

For a system described by a generic linear differential equations in space and time,

[TABLE]

The ansatz in Eq. (11) leads to a dispersion relation of the kind,

[TABLE]

which is solved by normal modes, providing an expansion basis to decompose any space-time perturbation $S$.

One can learn a great deal about the stability properties of the system, by studying the conditions which (do not) lead to imaginary parts for $\omega$ or $k$ in the dispersion relation. In general, the solutions can be cast either in the form

[TABLE]

or in the form

[TABLE]

Solutions in terms of Eq. (20) involve a temporal stability analysis where one searches for growing modes in time for a wave propagating with real wavenumber ${k}$. If not only $k$ but also $\Omega({k})$ is real, the wave is stable, while if $\Omega({k})$ develops a positive imaginary part, it can become unstable at large time $t$. On the other hand, solutions in terms of Eq. (21) involve a spatial stability analysis. In this case, if ${K}(\omega)$ develops an imaginary part for some real $\omega$, exponentially growing or decaying solutions would emerge as a function of the distance $z$.

Following the seminal work in [48], the stability of a linearized system can be studied in terms of the long-time behavior of a generic perturbation $S$, expressed either as “space-like” wavepacket (when the perturbation can grow with time, while propagating in space)

[TABLE]

or as a “time-like” wavepacket (when the perturbation can grow along a spatial direction, as time elapses),

[TABLE]

depending on the problem under study, e.g., the response to a space-localized perturbation, or to a forcing harmonic frequency, etc. By means of general arguments about the convergence properties of integrals with oscillating integrands, it was shown [48] that the perturbation asymptotics generally fall into one of the following four categories, which have been almost universally adopted in the literature on linear instabilities:

• $\bullet$

absolute instability, if the perturbation blows up “on site” and around;

• $\bullet$

convective instability, if the perturbation decays locally, but blows up elsewhere as it moves away;

• $\bullet$

complete stability, if the perturbation is neither enhanced nor damped;

• $\bullet$

stability with damping, if the perturbation decays (exponentially) in space.222 For electromagnetic disturbances, this case would lead to the so-called “evanescent waves” and “non-transparency” phenomena [44].

The time-asymptotic behavior of the function $S(z,t)$ for absolute and convective unstable cases is sketched in Fig. 1. In [48] the above (un)stable cases were also set in correspondence with the emergence of complex values of $(\omega,\,k)$ in the dispersion relation:

• $\bullet$

absolute instability, when $\omega$ can be complex for real $k$, but $k$ remains real for all real $\omega$;

• $\bullet$

convective instability, when $\omega$ can be complex for real $k$, and $k$ can be complex for real $\omega$ ;

• $\bullet$

complete stability, when $\omega$ is real for all real $k$, and $k$ is real for all real $\omega$;

• $\bullet$

stability with damping, when $\omega$ is real for all real $k$, but $k$ can be complex for real $\omega$.

We shall see an example of such correspondence in the simple two-mode system of Sec. V. However, the above diagnostics in terms of complex or real $(\omega,\,k)$, despite being useful in simple systems, does not lead to quantitative insights and, most importantly, may not be general enough to capture more complex situations. Already from the above classifications of instabilities we realize that the $z$ and $t$ variables (as well as their conjugates $k$ and $\omega$), do not play symmetrical roles in the dispersion relation approach. Indeed, perturbations causally propagate forward in time (and are absent for $t<0$), while they can propagate both downstream ($z>0$) and upstream ($z<0$) in space, involving Laplace and Fourier transforms for the conjugate variables $(t,\,\omega)$ and $(z,\,k)$, respectively. Thus, while a positive imaginary part of $\omega$ is a clear indicator of a time-growing instability, the sign of the imaginary part of $k$ is not an obvious diagnostic of space-growing (or decaying) perturbations. In addition, for impulsive perturbations one is usually interested in the local $(z=0)$ behavior of $S$ at large time ($t\to\infty$), while for steady harmonic forcing one focuses on space growth. Finally, multi-valued solutions may be more easily tractable in terms of $K(\omega)$ or $\Omega(k)$, depending on the system. Thus, both temporal and spatial analyses may be needed, in order to fully understand the instability conditions of the system, see e.g., [44, 47]. A rather general approach, originally developed in the context of plasma physics [49, 50] and further extended also to fluid dynamics [44, 47], is described in the following section.

IV Theory of Instability

We discuss below the general theory that leads to criteria distinguishing the different types of instabilities presented in the previous section, the so-called Briggs-Bers criterion [49, 50]. We follow the Green’s function formulation given in [44] and expanded in later reviews, see e.g. [51, 52, 53].

IV.1 Green’s Function

The Green’s function ${G}(z,\,t)$ encodes the response of the linear system governed by Eq. (18) to an impulsive forcing (or “source”) term,

[TABLE]

If ${G}(z,\,t)$ is known, the response $S(z,\,t)$ to a generic forcing $f(z,\,t)$ is obtained by a double convolution of $G$ and $S$. In terms of conjugate variables, the Green’s function is simply the inverse of the dispersion function,

[TABLE]

and the response $S$ is given by $S(\omega,\,k)=G(\omega,\,k)f(\omega,\,k)$.

According to Eq. (25), the Green’s function $G(z,\,t)$ admits an integral representation as a double (Laplace-Fourier) inverse transform of $1/D$,

[TABLE]

which is governed by the poles of $G(\omega,\,k)$, or equivalently by the zeros of $D(\omega,k)$, in the complex planes of $k$ and $\omega$. The Fourier integration domain $F_{k}$ is the usual real axis $k\in(-\infty,\,+\infty)$, closed by an upper half-circle for $z>0$ (enclosing a family of poles $k^{+}$) and by a lower half-circle for $z<0$ (enclosing a family of poles $k^{-}$), in order to ensure convergence. The Laplace contour $L_{\omega}$ must be shifted to $\omega\in(-\infty+i\sigma_{0},\,+\infty+i\sigma_{0})$, with $\sigma_{0}$ high enough to lie above all poles $\omega(k)$ with $k\in\mathbb{R}$ (this is often referred to as a Bromwich contour). In this way, $L_{\omega}$ can be closed by an upper half-circle to satisfy causality ($G(z,\,t)$=0) for $t<0$, while for $t>0$ it is closed by a lower half-circle. Figure 2 shows a sketch of such integration contours in the $\omega$ and $k$ planes. Note that poles may form discrete or (possibly multi-valued) continuous sets, depending on the system under study. In particular, we assume that the integral on $F_{k}$ in Eq. (26) is performed first, followed by the integral on the $L_{\omega}$. Then in the figure, the thick curves in the $\omega$-plane correspond to all solutions of the relation $D(\omega_{j},k)=0$ for any real $k$ and $j=1,\ldots N$ correspond to the $j$th branch of the solution. For example, if $D(\omega,k)=0$ is quadratic in $\omega$, as it happens in the two-beam example we discuss later, there are two branches for $\omega$ at every real $k$ as shown here. Meanwhile, the dots in the $k$-plane correspond to solution of the dispersion relation for a given $\omega=\omega_{0}$ on the isocontour of $\textrm{Im}(\omega)$ on $L_{\omega}$, i.e., $D(\omega_{0},k_{l}(\omega_{0}))=0$, $l=1,\ldots M$ and $\omega_{0}$ on $L_{\omega}$. For example, if $D(\omega,k)=0$ is cubic in $k$ there are three roots for a given choice of $\omega_{0}$, as shown.

General features of $G(z,\,t)$ emerge already in the simple case of a single pole $\omega=\Omega(k)$ for real $k$. For $t>0$, by closing the $L_{\omega}$ contour with a lower half-circle (clockwise oriented), the residue at $\Omega(k)$ yields [51, 52]:

[TABLE]

with the overall minus sign due to the contour going clockwise. The behavior of this integral at $t\to\infty$ along a propagation ray with finite $z/t$,

[TABLE]

can be obtained by the steepest descent method around the point where the phase $(kV-\Omega)t$ is stationary. Namely the point with coordinates $\hat{\omega}=\Omega(\hat{k})$ and $\hat{k}$ where

[TABLE]

A second-order expansion of the integrand around this stationary point yields the time-asymptotic result (derivation omitted) [51, 52]:

[TABLE]

which corresponds to a wavepacket with frequencies and wavenumbers related by the “group velocity” in Eq. (29).

The shape of the wave packet changes with the ray $V$ and, in particular, the phase $(\hat{k}V-\hat{\omega})t$ may acquire an imaginary part $\sigma$ for some $V$,

[TABLE]

If $\sigma(V)\leq 0$ for any $V$, then the flow is stable. Conversely, if $\sigma(V)>0$ for some $V$, then the wavepacket amplitude grows as $e^{\sigma t}$. In particular, if $\sigma>0$ for $V=0$, then the amplitude grows exponentially at $z=0$, i.e., locally (case of absolute instability), otherwise it moves away faster than it spreads out (case of convective instability).

One can thus translate the classification of the instabilities discussed in Sec. III.2 in terms of properties of the Green’s function along propagation rays $V=z/t$. In particular, a stable flow implies that

[TABLE]

while an unstable flow arises when the wave packet grows along at least for some $z/t$,

[TABLE]

Convective or absolute instability can be further distinguished by the local response at $z=0$, i.e., along the ray $z/t=0$. An unstable flow displays a convective instability if

[TABLE]

otherwise the instability is absolute,

[TABLE]

Note that, in a temporal stability analysis with $\omega=\Omega(k)$ for real $k$, it is simply $\sigma=\textrm{Im}(\Omega(k))$, and the above criteria would seem to be related only to the position of poles in the complex $\omega$ plane. However, this is often but not always the case, and a more powerful (geometric) criterion can be envisaged to characterize the instabilities.

IV.2 Geometric Criterion for Instability

We consider first the time asymptotics of an initially localized disturbance $\delta(z)\delta(t)$, in order to refine the characterization of absolute and convective instabilities. Then we consider the evolution of a localized harmonic forcing $\delta(z)e^{i\omega_{0}t}$, which helps to refine the discrimination between convective (unstable) and damped (stable) cases.

IV.2.1 Absolute vs Convective Instability

The criterion is based on two key observations. The first is that the Fourier and Laplace contours in Eq. (26) and in Fig. 2 can be continuously deformed in the corresponding complex planes, as far as they do not cross any pole. Studying the large-time behavior of $G$ implies lowering the $L_{\omega}$ contour towards the half-plane with negative $\textrm{Im}(\omega)$, the integral being then dominated by the highest pole in the complex $\omega$ plane. Referring to Fig. 3 we take the contour $L_{\omega}$ and lower it downward (left panel). If we focus our attention on a specific point on $L_{\omega}$, denoted by $\omega_{0}$, we see that as $L_{\omega}$ is lowered that point gets closer to a pole. Since the corresponding poles in the $k$-plane, i.e., $k_{l}=k_{l}(\omega_{0})$, depend on $\omega_{0}$, they will also move around in the $k$-plane as $L_{\omega}$ is moved downward (right panel). In particular, when the point on $L_{\omega}$ that was $\omega_{0}$ eventually hits a pole on the highest branch $\Omega_{j}$, some pole(s) in $k$-plane must also hit the real $k$-axis. This is because $\Omega_{j}$ is defined by $D(\Omega_{j},k_{l})=0$ with $k_{l}$ being real. Once $k_{l}$ becomes real, the Green’s function $G(\omega,z)$ is no longer analytical because $k_{l}$ lies on the $F_{k}$ contour of Eq. (26) which is the real $k$-axis. In order to analytically continue $G(\omega,z)$ below the largest value of ${\rm Im}(\Omega_{j})$ it is then necessary to deform the $k$-integration contour around those poles which approach or cross the real $k$-axis, i.e., choose a new contour ${\tilde{F}}_{k}$ as shown in Fig. 3 (right panel). Figure  4 shows the different possibilities regarding deformation in the complex $k$-plane. redIn the left panel we show the roots of the dispersion relation in the $k$-plane for a point $\omega_{0}$ on the original $L_{\omega}$ contour, i.e., at $\textrm{Im}(\omega)\to\infty$. As $L_{\omega}$ is lowered towards ${\rm Im}(\omega)\to 0$, the poles move in the $k$-plane. Some of them can cross the real $k$-axis but they can be avoided by deforming the ${F}_{k}$ contour (middle panel). However, a singularity occurs if two moving poles come close together and pinch the contour of integration $\tilde{F}_{k}$. This latter becomes “stuck” and cannot be deformed around the two merging poles, signalling an instability.

The second observation is that the modes at zero group velocity $V=0$, obtained for some ($\omega_{c},\,k_{c}$), are double roots of the dispersion relation, since they obey not only $D=0$ but also $\partial D/\partial k=0$ via Eq. (29). The dispersion relation can be locally expanded as

[TABLE]

which implies that $k=k(\omega)$ has two local solutions $k^{\pm}$ and that the pinching point is just $k_{c}$, where the two poles $k^{\pm}$ coalesce. The corresponding $\omega$-pole is located at $\omega_{c}$, with $\textrm{Im}(\omega_{c})$ determining the nature of the instability. If $\textrm{Im}(\omega_{c})<0$, the contour $L_{\omega}$ can be safely lowered below the real axis, so that $G\to 0$ at large time and the instability is convective. If instead $\textrm{Im}(\omega_{c})>0$, the instability necessarily grows in time and is absolute. In a nutshell, absolute instabilities arise when modes with zero group velocity develop a positive imaginary part $\textrm{Im}(\omega_{c})$, in which case both downstream and upstream flows ($k^{\pm}$) mix and blow up.

The “pinching” criterion is consistent with the classification of the previous subsection, but is also more general. In fact, there may be cases where a double root of $D=0$ is found, but the $k$-poles approach the deformed $F_{k}$ contour from the same side rather than from opposite sides, in which case instabilities do not arise. The criterion is also amenable to numerical and graphical implementations which help to identify instabilities: A cartesian grid spanning the complex $\omega$-plane near $\omega=\omega_{c}$ is mapped into a deformed grid in the $k$-plane with a saddle point at $k=k_{c}$; vice versa, a cartesian grid in the complex $k$-plane around $k_{c}$ is mapped into a deformed grid in the $\omega$-plane with a cusp point at $\omega_{c}$ [47].

IV.2.2 Convectively unstable and Damped stable solutions

The previous criterion for instabilities was mainly based on a stability analysis in the time domain. Further insights can be gained by looking at the evolution of disturbances in space. In particular, the damped stable solution introduced in Sec. III B can be better characterized and separated from the convective unstable solution (both decaying away at $z=0$), by considering the upstream and downstream response of the system to a harmonic forcing with real frequency $\omega_{f}$ and localized at $z=0$, namely, $f(z,t)\sim\delta(z)e^{-i\omega_{f}t}$ for $t>0$. In terms of conjugate variables:

[TABLE]

The response $S$ is then governed by $S=G(\omega,k)f(\omega,k)$, which has a pole at $\omega=\omega_{f}$ in addition to the Green’s function poles. If there are no absolute instabilities, then the $\omega$-poles are located at $\textrm{Im}(\omega_{f})=0$ and possibly at $\textrm{Im}(\omega)<0$ if there are convective instabilities.

As already mentioned, in order to get the large-time behavior one lowers the $L_{\omega}$ contour. Initially the $L_{\omega}$ contour is located at high $\textrm{Im}(\omega)>0$, and the real $F_{k}$ contour separates upper poles $K^{+}$ with $\textrm{Im}(K^{+})>0$ from lower poles $K^{-}$ with $\textrm{Im}(K^{-})<0$, where the dispersion relation is solved in terms of $k=K(\omega)$. The $L_{\omega}$ contour can be lowered down to the real axis, where it hits the highest pole at $\omega=\omega_{f}$, which dominates the large-time evolution of both downstream and upstream modes as [44]

[TABLE]

Note that the poles $K^{\pm}(\omega_{f})$ have been displaced from their original position. If the initial sign of $\textrm{Im}(K)$ is preserved for both poles $K^{\pm}$, then the above solutions decay exponentially, both upstream and downstream, leading to a damped stable case. If instead at least one $K^{\pm}$ pole changes sign in $\textrm{Im}(K)$ (i.e., if the pole moves into the “wrong” half-plane for complex $k$), as shown in the middle panel of Fig. 4, then the response is exponentially amplified in at least one direction (either upstream or downstream), corresponding to a convective instability which moves away from the source.

This argument provides a criterion to distinguish the cases of stable damping and of convective instability: if $\textrm{Im}(K(\omega))$ changes sign when $\textrm{Im}(\omega)$ varies from $+\infty$ to [math] at fixed $\textrm{ Re}(\omega)=\omega_{f}$, then the flow is convectively unstable; if instead $\textrm{Im}(K(\omega))$ does not change sign, then the (stable) flow is damped.

V Two-beam case

The general formalism discussed in the previous Section is applied below to a specific case, widely considered in plasma physics (see, e.g., [49, 44]), and also suitable as a toy model for fast flavor conversions, as shown in [43]. Namely we assume a two-beam neutrino model with a ELN spectrum

[TABLE]

Using this angular spectrum, one can rotate away the first term in the right-hand-side in Eq. (7) [as in Eq. (12)]. By assuming azimuthal symmetry with respect to the $z$-direction, and translational invariance with respect to the transverse directions, the problem reduces to the evolution in $z$ and $t$. Under these simplifying assumptions, one gets a system of coupled equations for the two velocity modes

[TABLE]

where $v_{1}$ and $v_{2}$ indicate the projection of $\nu$ velocity ${\bf v}_{1,2}$ along the $z$-direction [see Eq. (10)], and

[TABLE]

Note that these equations are equivalent to those of of a two-level system excited by an electric field (see, e.g., [54, 55]). If one writes explicitly these set of equations for the real and the imaginary parts of $S_{1}$ and $S_{2}$ one gets two sets of coupled differential equations

[TABLE]

where $(f_{1},f_{2})=(\textrm{Im}\,S_{1},\textrm{ Re}\,S_{2})$ or $(\textrm{Im}\,S_{2},\textrm{ Re}\,S_{1})$.

V.1 Stability Analysis

We assume the ansatz

[TABLE]

which implies the following dispersion relation for $\omega$ and $k$:

[TABLE]

where

[TABLE]

and the amplitudes in Eq. (45) are related by

[TABLE]

The dispersion relation in Eq. (46) can be solved either as

[TABLE]

or as

[TABLE]

In the absence of a coupling between the two modes ($\varepsilon=0$), Eq. (49) would give $\Omega_{\pm}(k)=v_{1,2}k$, and the two modes would cross at $k=0$. Switching on the coupling ($\varepsilon\neq 0$), two different situations arise, depending on the sign of $\varepsilon$. If $\varepsilon>0$, then $\Omega(k)$ is real for real $k$ and the system is stable. In contrast, if $\varepsilon<0$, then $\Omega(k)$ is complex for $k$ in the range

[TABLE]

and the system becomes unstable. The transition from stability to instability involves a sign change of $\varepsilon=g_{1}g_{2}$ and thus of the relative sign of $G_{1}$ and $G_{2}$ in Eq. (39). This conclusion extends to more general $G(\theta)$ where one needs a crossing from positive to negative ELN intensities to have an instability.

For $\varepsilon<0$, the nature of the instability (as discussed in the previous section) is determined by the imaginary part of the frequency $\omega_{c}=\Omega(k_{c})$, at the wavenumber $k_{c}$ corresponding to null group velocity $V$. In our two-mode system $V$ reads

[TABLE]

and vanishes for complex $k_{c}$ and $\omega_{c}$ given by

[TABLE]

and

[TABLE]

respectively. It is easy to verify that $k_{c}$ is a double root of $K(\omega)$ for $\omega=\omega_{c}$. The above equations entail an unstable solution with $\textrm{Im}(\omega_{c})>0$ (corresponding to an absolute instability) only for $v_{1}v_{2}<0$. For unstable cases $(\varepsilon<0)$, the convective or absolute nature of the instability is thus determined by the sign of $v_{1}v_{2}$. For stable cases $(\varepsilon>0)$, the function $K(\omega)$ develops an imaginary part only if $v_{1}v_{2}<0$, in which case the stable mode is exponentially damped. The sign of $v_{1}v_{2}$ thus distinguishes also damped and completely stable modes. Summarizing, one can recover the four categories introduced in Sec. III.2 in terms of the sign of $\varepsilon$ and of $v_{1}v_{2}$ as follows (see [44]).

• $\bullet$

Complete stability ($\varepsilon>0$, $v_{1}v_{2}>0$), see Fig. 5. In this case $\Omega(k)$ is real for real $k$ and $K(\omega)$ is real for real $\omega$. No imaginary part is developed, and the system is completely stable.

• $\bullet$

Stability with damping ($\varepsilon>0$, $v_{1}v_{2}<0$), see Fig. 6. The function $\Omega(k)$ is real for real $k$, and the system cannot exhibit convective or absolute instability. However, $K(\omega)$ is complex in the frequency range

[TABLE]

graphically corresponding to a “gap” in $\omega$ , where perturbations decay exponentially in the space coordinate $z$.

• $\bullet$

Convective instability ($\varepsilon<0$, $v_{1}v_{2}>0$), see Fig. 7. In this case the dispersion relation can provide complex $\omega$ for real $k$ and complex $k$ for real $\omega$. Following the criteria outlined in Sec. IV, we consider the path of the complex roots $K(\omega)$ of Eq. (50) from large to small $\rm{Im}(\omega)$, assuming $v_{1,2}>0$ for definiteness. For $\rm{Im}(\omega)\to\infty$, both roots are in the same upper half of the complex plane,

[TABLE]

When $\textrm{Im}(\omega)\to 0$, the roots $K(\omega)$ in the range of of Eq. (55) are instead complex conjugate, and the one for which $\textrm{Im}(K(\omega))<0$ has migrated from the upper to the lower half of the complex plane. Therefore, in the range of Eq. (55) we have a convective instability propagating in the positive $z$ direction [see Eq. (38)].

• $\bullet$

Absolute instability ($\varepsilon<0$, $v_{1}v_{2}<0$), see Fig. 8. In this case the function $K(\omega)$ is real for all real $\omega$, but the function $\Omega(k)$ is complex in the range of Eq. (51), corresponding to a gap in $k$. We note from Eq. (56) that for $\textrm{Im}(\omega)\to\infty$, since $v_{1}v_{2}<0$ the two roots $K(\omega)$ are in opposite half-planes. Such roots coalesce at a single pinching point as $\textrm{Im}(k)$ is gradually lowered, down to imaginary $\omega_{c}$ value derived from Eq. (54).

The behavior of the $K(\omega)$ roots in the complex $k$ plane for convective and absolute instabilities is further illustrated in Figs. 9 and 10, respectively. In Fig. 9, nearly vertical and solid curves correspond to $K(\omega)$ roots of the dispersion relation for the indicated value of $\textrm{Re}(\omega)$, while nearly horizontal and dashed curves correspond to isocontours of $\textrm{Im}(\omega)$. The two colors refer to the two roots of the dispersion relation. Lowering $\textrm{Im}(\omega)$ towards zero, one of the two solutions at given $\textrm{Re}(\omega)$ migrates from upper to lower plane and gives rise to a convective instability. In Fig. 10, nearly vertical lines indicate $K_{+}(\omega)$ and $K_{-}(\omega)$ roots of the dispersion relation for the indicated value of $\textrm{Re}(\omega)$. Horizontal dashed lines indicate isocontours of $\textrm{Im}(\omega)$. By lowering $\mathrm{Im}(\omega)$, the pinching point is eventually found at $\mathrm{Im}(\omega_{c})=0.46$ and $\mathrm{Im}(k_{c})=-0.38$, in agreement with Eqs. (53) and (54) for the chosen two-mode parameters.

We conclude this section with some comments on the growth rate along a ray $V=z/t$ for unstable flows ($\varepsilon<0$), as given in Eq. (31). For the two-mode system we find:

[TABLE]

where

[TABLE]

The maximum growth rate $\sigma=\sqrt{\varepsilon}$ occurs at $\delta=0$, namely, $V=(v_{1}+v_{2})/2$, which applies to both convective and absolute instabilities. Only for $v_{1}v_{2}<0$ (absolutely unstable case) there growth rate is also defined at $V=0$, and one finds then $\sigma=\textrm{Im}(\omega_{c})$ with $\omega^{2}_{c}$ as in Eq. (54). Finally, notice that the growth rate vanishes at $V=v_{1,2}$.

V.2 Particle-like vs Tachyon-like Dispersion

Here, we briefly mention another possible interpretation of the instabilities discussed above. We consider the set of the coupled differential equations [Eq. (44)] for the two-beam problem transforming the time coordinate as [56, 57]

[TABLE]

so that the apparent group velocity with the new time coordinate are the same in magnitude but opposite in sign and given by

[TABLE]

In this reference system one obtains as coupled equations

[TABLE]

Applying to the first equation of the system in Eq. (62) the differential operator of the second equation (or viceversa) one obtains the Klein-Gordon equation 333Note that in the moving frame, one can introduce a spinor wave field $\psi=(f_{1},f_{2})^{T}$, writing Eq. (62) as a one-dimensional Dirac equation [57].

[TABLE]

with

[TABLE]

From Eq. (63) one obtains as dispersion relation

[TABLE]

If $m^{2}>0$, this quantity plays the role of a mass term, then Eq. (65) is the dispersion relation of a particle, having a gap in $\omega$. From the dispersion relation, one realizes that if $\omega^{2}>m^{2}{\bar{c}}^{4}$, Eq. (63) would have oscillatory waves as the solution, while if $\omega<m^{2}{\bar{c}}^{4}$ it would represent damped oscillatory waves. This is consistent with what found for the in the damped case in the previous Section for two counter-propagating modes with $\varepsilon>0$.

If we now move to the case of $m^{2}<0$, Eq. (63) would represent a Klein-Gordon equation with imaginary mass. In this case the dispersion relation of Eq. (65), with a gap in $k$, would be the one expected for “tachyons” [58, 59]. From this dispersion relation if $k^{2}>m^{2}{\bar{c}}^{2}$ one would expect normal oscillatory motion. Conversely, for $k^{2}<m^{2}{\bar{c}}^{2}$ we would have an exponential growing solution. It is intriguing to realize that the absolute instability found for two counter-propagating modes with $\varepsilon<0$ is of a “tachyonic” type.

V.3 Numerical Results

In order to illustrate the predictions of the stability analysis for the four cases discussed in the previous Sections, we work out representative numerical solutions of Eq. (44). We assume a length interval $z\in[0,\,L]$. In the following we will work in the units in which the neutrino potential in Eq. (9) is $\mu=1$. Therefore times and length are expressed in units of $\mu^{-1}$. Moreover, we will assume ${\mathcal{O}}(1)$ initial values for $f_{1}$ and $f_{2}$, since in the linear regime these are just arbitrary normalization factors.

When $v_{1}v_{2}>0$ we assume that the two modes are emitted at $z=0$. Conversely when $v_{1}v_{2}<0$ we assume that the mode with $v_{1}>0$ is emitted at $z=0$, while the mode with $v_{2}<0$ is emitted at $z=L$.

At first, we look for plane-wave solutions, given by Eq. (45). We assume the same numerical parameters as for the four cases in Sec. V.1. In Fig. 11 we plot the orbits of the solution ($f_{1},\,f_{2}$) as a function of time $t$ at some fixed $z$ (the exact value being irrelevant for our discussion). To orbits start at the square point (at $t=0$) and terminate at the star point (at some large $t$). If $\omega$ and $k$ were taken both real and exactly related by the dispersion relation, the orbits would be stable ellipses (not shown). We instead allow small numerical deviations from the dispersion relation, in order to mimic perturbations and to study their evolution. In the completely stable case ($v_{1}v_{2}>0$ and $\varepsilon>0$, upper left panel), perturbations are neither enhanced nor suppressed, and the orbits remain quite stable in time. In the stable case with damping ($v_{1}v_{2}<0$ and $\varepsilon>0$, upper right panel), perturbations eventually lead to a decay of the solution amplitudes, signaled by the shrinking of the orbit towards $f_{1}=f_{2}=0$. In the convectively unstable case ($v_{1}v_{2}>0$ and $\varepsilon<0$, lower left panel), perturbations are amplified and orbits diverge for a transient time, until the unstable flow moves away and decays at any fixed $z$. Finally, in the case of absolute instability ($v_{1}v_{2}<0$ and $\varepsilon<0$, lower right panel), after a transient the orbit parameters tend to grow indefinitely with time.

In order to go beyond monochromatic waves, we perform a second numerical experiment by launching two wavepackets with different velocities $v_{1}$ and $v_{2}$, with an initial shape

[TABLE]

Numerically we fix $k=15$, $\sigma_{f}=0.5$, while all the other two-mode parameters are as given in the previous section. The evolution in $z$ at different $t$ for the four different cases of Sec. V.1 is shown in Figs. 12, 13, 14, and 15. In each case the fastest mode ($v_{1}$) is represented in red, while the slowest ($v_{2}$) is in blue. The two wave packets act as seeds of perturbations, which will generally propagate from some $(z,\,t)=(z_{0},\,t_{0})$. The relevant coordinates, especially for unstable (growing) modes are then $z-z_{0}$ and $t-t_{0}$. In particular, rays are defined by constant values of $V=(z-z_{0})/(t-t_{0})$.

In Fig. 12 (completely stable case), the two packets travels in the same direction ($v_{1}v_{2}>0$) with positive coupling $(\varepsilon>0)$. Besides the wavepackets, small perturbation are created, and all signals propagate without any amplification. In Fig. 13 (damped case), the two packets travel in opposite directions ($v_{1}v_{2}<0$) with positive coupling ($\varepsilon>0$). Perturbations are created in this case too, but are not amplified; actually they are damped, but the size and the numerical observation time are too small to make the damping graphically evident. In Fig. 14 (convective instability) the two packets travel in the same direction ($v_{1}v_{2}>0$) with negative coupling. A disturbance is created at some finite $z_{0}$ and $t_{0}$, and then grows exponentially but also moves away. We have numerically verified that the growth of the disturbance closely follows Eq. (57) along any ray $V=(z-z_{0})/(t-t_{0})$ and, in particular, that the perturbation is contained in the ray interval $V\in[v_{1},\,v_{2}]$, and is maximally enhanced along the intermediate ray $V=(v_{1}+v_{2})/2$ by an exponential growth factor $\propto\varepsilon^{1/2}t$. Finally in Fig. 15 (absolute instability) the two packets travel in opposite directions ($v_{1}v_{2}<0$) with negative coupling ($\varepsilon<0$). The growth of the disturbance embraces the original point $z_{0}$ where it has been generated, with the predicted local amplification factor—see the comments after Eq. (57)–(58).

VI Summary and perspectives

In our work we have discussed in detail a classification of the fast instabilities which may arise in the flavor evolution in space-time $(z,t)$ of a dense and self-interacting neutrino gas, such as close to the neutrino sphere in core-collapse supernovae. This classification is based on the dispersion relation $D(\omega,\,k)$ among the among the $(t,\,z)$-conjugate coordinates ($\omega,\,k$). The dispersion relation has been recently introduced in [43] for SNe and independently elaborated in more general contexts, such as plasma physics and fluid dynamics [44]. If the disturbances in the mean field of the $\nu_{e}\nu_{x}$ flavor coherence grow, propagating away from the point of origin, they are associated to convective instabilities. Conversely, if the disturbances grow in amplitude and extent, embracing the point of origin, they are called absolute instabilities. Cases with no instabilities cases may also be differentiated into completely stable (with neither growth nor decay of disturbances) and damped stable (with decay of disturbances).

Starting from the dispersion relation, at least for simple systems, one finds that if $\omega$ is real for all real $k$ and vice versa the system is completely stable. If $\omega$ is real for all real $k$, but $k$ is complex for some real $\omega$, the flow is both stable and damped. Instead, if $k$ is complex for some real $\omega$ and $\omega$ is also complex for some real $k$, a convective instability arises. Finally, if $k$ is real for all real $\omega$, and $\omega$ is complex for some real $k$, the instability is absolute. Deeper criteria can be envisaged to identify instabilities via poles of $D(\omega,\,k)=0$ in the complex $k$ and $\omega$ planes. In particular, instabilities can emerge from an evaluation of the time-asymptotic behavior of the Green’s function of the system, which is related to the dispersion relation by a double (Laplace-Fourier) integral representation, amenable to complex-calculus techniques. These techniques, extensively developed in the field of plasma physics and fluid dynamics, have been presented and discussed herein in the context of self-interacting neutrinos, with particular attention to a simple two-modes system. For a two-beam model of dense neutrinos that are forward-scattering off each other, we considered four cases leading to the four possibilities described above, and presented a comparison of the predictions from the linear instability theory with the numerical solutions of the linearized equations of motion. The comparison demonstrate that the two-beam system can be fully understood theoretically and is under control numerically.

The results obtained for the simple two-beam neutrino model represent a basis to attempt extensions to more general angular spectra $G_{\bf v}$ [see Eq. (8)], as expected in a realistic SN. As pointed out in [43], one needs a crossing from positive to negative $G_{\bf v}$ (i.e., in the ELN) in order to have an instability. Conversely, an ELN with no crossing would give either a completely stable evolution (if $v_{1}v_{2}>0$) or at most a damped stable one (if $v_{1}v_{2}<0$). The ELN giving fast flavor conversions in SNe would correspond to situations in which there is a significant backward-going neutrino flux, corresponding to $v_{1}v_{2}<0$ in the simple two-beam model. Depending on a possible crossing in the ELN this would lead to damped solutions or to absolute instabilities. In this context, a dedicated investigation of the energy and angle distributions of the neutrino radiation field has been presented in [60] for several spherically symmetric (1D) supernova simulations. In the cases studied, the ELN near the neutrino-sphere has backward going modes but still does not show any crossing. According to our analysis, in first approximation this situation would correspond to the damped stable case, so that no instability (i.e., fast conversion) should show up in such 1D SN models. However, one cannot exclude that these finding may change in 3D models, for example in the presence of LESA (Lepton-Emission Self-sustained Asymmetry) [61]. LESA manifests itself in a pronounced large-scale dipolar pattern in the ELN emission and naturally implies a change of sign in $\nu_{e}-{\bar{\nu}}_{e}$ angular distributions. It is therefore conceivable that, especially in the regions where the ELN changes its sign, crossings in the ELN angular distributions may occur. In this case one could expect the emergence of absolute instabilities. In this situation, the time evolution of the neutrino gas would be dominated by the spectrum of wave-numbers $g_{k}$ in a given region of space close to neutrino emission. At this regard, fast flavor conversions have been shown to arise for neutrino angular distributions inspired by LESA models, assuming space homogeneity [42]. The next logical step would be to remove the homogeneity assumption, thus introducing a spectrum $g_{k}$. A solution of this problem could be obtained by following the strategy developed in [24, 26], i.e., by Fourier transforming the equations of motion with respect to the space coordinate, in terms of coupled equations for different conjugate modes. If the $k=0$ (homogeneous) mode is unstable, the instability can then cascade to smaller scales by the coupling among the different modes.

It has also been recently pointed out that a crossing in ELN could be a generic feature of neutrino emission from binary neutron-star mergers [62]. In this case fast conversions would occur, and the associated instability would be of absolute type. All these cases deserve dedicated studies to gain both a deeper understanding and a broader perspective. Indeed, the analysis of the dispersion relation and the relative classification of the instabilities should be extended to cases presenting many emission angles, both in the zenith and azimuthal ones. As pointed out in [43], this would lead to instabilities breaking the axial symmetries, that definitely need further investigations.

The use of the dispersion relation to classify the neutrino flavor instabilities could also be extended to the study of slow self-induced flavor conversions in SNe. In this case the instabilities would depend also on the vacuum oscillation term. The slow self-induced conversions occuring for free streaming neutrinos far from the neutrino-sphere are dominated by forward modes (i.e., by half-isotropic zenith angle distributions [32]). In the two-beam models, this case would correspond to $v_{1}v_{2}>0$, leading to either a completely stable case or a convective instability. Numerical simulations have been found a transition from a stable situation at low-radii (where the neutrino system would exhibit synchronized oscillations [32]), to an unstable region at larger radii. The boundary between these two regions, called synchronization radius, has been taken as a sort of inner boundary for the subsequent evolution. Using this approach, spatially growing solutions have been found. Recently, the assumption of a stationary emission from the spatial boundary has been questioned [28, 29, 30], on the basis that any frequency would have some amplitude at the boundary that would be amplified in the further evolution. In this regard it may happen that the spatially growing solution found in previous works are not stationary. Indeed, in the presence of convective instability, the perturbations would propagate while they grow. This effect could shift the onset of the flavor conversions towards larger radii, where the neutrino density is smaller and thus the instability is weaker. This important issue also requires dedicated investigations.

From this discussion it appears that the phenomenology of self-induced flavor conversions in SNe would be much richer than previously thought. In the presence of fast conversion, the usual characterization of flavor oscillations in terms of spatial evolution from an inner boundary should be probably revised, in favor of a time evolution within a small region of space near the neutrino decoupling region. This would make the approach to SN somewhat similar to the study of neutrinos in the early Universe. From present studies [36, 37, 38, 42], one might argue that fast conversions could lead to a quick flavor equilibration among different neutrino species, if instabilities are general enough. If flavor equilibration were complete, further oscillation effects would be ineffective. Otherwise, one could characterize different regimes, e.g., fast conversions near SN core followed by spatial slow conversions at larger distances. We believe that our current study can provide valuable and inter-disciplinary tools to predict the possible scenarios occurring of the dense SN neutrino gas. The many open questions call for further numerical and analytical research in this field.

Acknowledgments

We acknowledge useful discussions with Amol Dighe, Georg Raffelt, Günter Sigl and Irene Tamborra during the development of this work. We also thank Georg Raffelt and Irene Tamborra for useful comments on the manuscript. Alessandro Mirizzi acknowledges kind hospitality at TIFR (Mumbai) and at CERN where part of this work was done. The work of Francesco Capozzi is supported by NSF Grant PHY-1404311 to J.F. Beacom. The work of Basudeb Dasgupta is partially supported by the Dept. of Science and Technology of the Govt. of India through a Ramanujam Fellowship and by the Max-Planck-Gesellschaft through a Max-Planck-Partnergroup. The work of Eligio Lisi, Antonio Marrone and Alessandro Mirizzi, is supported by the Italian Istituto Nazionale di Fisica Nucleare (INFN) through the “Theoretical Astroparticle Physics” project and by Ministero dell’Istruzione, Università e Ricerca (MIUR).