Self-induced temporal instability from a neutrino antenna

It has been recently shown that the flavor composition of a self-interacting neutrino gas can spontaneously acquire a time-dependent pulsating component during its flavor evolution. In this work, we perform a more detailed study of this effect in a model where neutrinos are assumed to be emitted in a two-dimensional plane from an infinite line that acts as a neutrino antenna. We consider several examples with varying matter and neutrino densities and find that temporal instabilities with various frequencies are excited in a cascade. We compare the numerical calculations of the flavor evolution with the predictions of linearized stability analysis of the equations of motion. The results obtained with these two approaches are in good agreement in the linear regime, while a dramatic speed-up of the flavor conversions occurs in the non-linear regime due to the interactions among the different pulsating modes. We show that large flavor conversions can take place if some of the temporal modes are unstable for long enough, and that this can happen even if the matter and neutrino densities are changing, as long as they vary slowly.

dispersion due to a large matter density. The line model behaves in this case as a neutrino antenna, propagating a pulsating neutrino signal.
In the present work we will study the development of the temporal instability, choosing several representative cases and comparing, for each of them, the results of the numerical solution of the equations of motion with analogous results from a linear stability analysis. In Sec. 2 we describe the features of the non-stationary and inhomogeneous flavor evolution. We discuss the equations of motion for time-dependent two-dimensional case. We show how it is possible to solve this problem by Fourier transforming these equations, obtaining a tower of ordinary differential equations for the different Fourier modes in space and time. We also present the linear stability analysis of these equations, with explicit algebraic equations for the eigenvalues for the line model. In Sec. 3 we discuss the results of our study for different representative cases. We find a good agreement between the results obtained with this two approaches in the linear regime, while a dramatic amplification of the flavor conversions occur in the non-linear case due to the interaction among the different pulsating modes, that get excited in a cascade. Finally, in Sec. 4 we comment about future developments and we conclude.
2 Setup of the flavor evolution

Equations of motion
We consider the situation that neutrinos and antineutrinos are emitted from a surface and subsequently free-stream, but with forward scatterings with other neutrinos and antineutrinos as well as the background matter. The flavor evolution can be characterized in terms of the (anti)neutrino density matrix E,v , for each neutrino mode with energy E and velocity v. In this formalism the diagonal elements of the density matrix represents the occupation numbers, while the off-diagonal terms encode phase relations that allow one to follow flavor oscillations. Then, the equation of motion (EOM) for such a stream of neutrinos is [33][34][35] (2.1) Of course varies with time and space, i.e., = (t, x). This space-time dependence is always understood, as also for other medium-dependent quantities, e.g., the background net electron density n e (t, x). For notational clarity, hereafter we shall omit showing these dependencies wherever there is no scope for confusion. The Hamiltonian for flavor evolution in the collisionless limit is where, in the two flavors limit, where U θ is the mixing matrix. We introduce the mass-squared difference ∆m 2 = m 2 2 − m 2 1 , which is positive for normal mass ordering and negative for inverted ordering. The integral over energy spans −∞ to +∞, with "negative" E corresponding to antineutrinos of energy E, i.e.,¯ E = − −E . Note the overall minus sign, so that equations for neutrinos and antineutrinos are the same with the replacement E → −E. The integral over v corresponds to integrals over the two independent velocity components, that are typically chosen such that E,v (t, x) remains constant in the absence of oscillations.
The vacuum oscillation frequency ω = |∆m 2 |/(2E) is a more convenient energy label than E. Thus, one has ω > 0 for neutrinos and ω < 0 for antineutrinos. The density matrix can be written as where s 2 = 1 − |S| 2 , with S = 0 when no flavor evolution has occurred. The function Φ ν is a suitable normalization of the emission spectrum g ω,v , defined by where dΓ = dω dv and Φν e,x are the unoscillated flavor-dependent neutrino number fluxes at x. The neutrino-antineutrino asymmetry parameter is given by The equations of motion [Eq. (2.1)] can then be written expanding all the quantities on the Pauli matrices σ as is the matter potential and is the ν-ν potential. The vector P ω,v = g ω,v (Re S ω,v , −Im S ω,v , s ω,v ) T , constructed from the density matrix, is a so-called polarization vector that encapsulates the flavor composition, while B = (sin 2θ, 0, cos 2θ) T and L = (0, 0, 1) T are the analogous vectors constructed from H vac and H mat , respectively. In this paper we present our results for a normal neutrino mass ordering, i.e., m 1 < m 2 , which leads to the minus sign in front of ωB in Eq. (2.9), and for the inverted ordering one simply flips that minus sign to get the equivalent equations. The problem boils down to calculating P ω,v (t, x), given their values at the source as a function of time. Typically this boundary condition is taken to be stationary in time and homogeneous over the source, which may naively suggest that the solution ought to respect these symmetries as well. However, such a solution is often unstable to small spatial and temporal fluctuations, and it is important to ascertain the role of spontaneous breaking of these spatial and temporal symmetries. We will study this in both the linear and nonlinear regime, and we introduce the set-up and the notation in the next sections.

Evolution of Fourier modes in the line model
To study the role of these fluctuations concretely, we choose as toy-model of neutrino emission the line model where the growth of these instabilities is easily calculable. The geometry of the model is shown in Fig. 1. Neutrinos and antineutrinos are emitted in the x − z plane, in two directions v L = cos ϑ Lx + sin ϑ Lẑ and v R = cos ϑ Rx + sin ϑ Rẑ , as shown. The dynamics is confined to this plane, though one may also think of this as a 3 dimensional problem with rotational symmetry around the x-axis. We consider neutrinos and antineutrinos emitted with a single energy E 0 , so that This fixes the frequency scale of flavor evolution. Therefore, the spectrum is given by which implies the normalization Φ ν = (Φν e − Φν µ )/2. Φ ν and λ are constant along x and in t but vary along z. Away from the source, the neutrino flavor composition can depend on z, as well as x and t. The boundary conditions on P ω,v will have to be specified along both x and t at z = 0, and break the corresponding space-time translation symmetries.
The flavor evolution can be studied using Fourier modes of the polarization vectors. Eq. (2.9) can be converted to a tower of ordinary differential equations by decomposing the polarization vectors into their constituent Fourier modes labeled by p and k [23][24][25]. One writes where p, k are the temporal and spatial frequency modes of the polarization vector. Using this decomposition, one gets a tower of equations for the Fourier modes For our model, we have the simplifications, , ω only picks up the modes at ±ω 0 , and v = v L or v R . The velocity components, v x and v z , being different for the L and R modes, produce multi-angle effects in additional to the spatial and temporal symmetry breaking in our set-up. For our numerical studies, we assume that only ν e and ν e are emitted, with a factor of two excess of ν e over ν e , i.e., = 1, chosen to guarantee that flavor conversions do not take place in the homogeneous case for the adopted value of the neutrino density µ on the boundary. We take the L and R modes to have two different angles ϑ R = 5π/18 and ϑ L = 7π/9. In this way, we mimic the multi-angle matter suppression in the presence of a large matter term. We choose θ = 10 −3 and a normal mass ordering, i.e., ∆m 2 > 0, but the result would be qualitatively similar for the inverted ordering, i.e., ωB → −ωB. The overall frequency-scale is set by ω 0 = 1.

Linearized stability analysis using Fourier modes
The complete flavor evolution even in this rather simplified model is quite complicated, and it is useful to linearize the equations. To linear order in S, the EOMs [Eq. (2.9)] simplify to [27,36] (2.16) One can again express S ω,v using its Fourier transform, The eigenvalue equation for the Fourier amplitudes Q p,k ω,v is then given by In the model we consider, there are only four neutrino momentum modes labeled by ω = ±ω 0 and v = v L and v R , and the above equation becomes an eigenvalue equation for the four modes, with 4 eigenvalues Ω i=1,2,3,4 , as a function of p and k. These eigenvalues are given by the equation where, Unlike in the case of a symmetric set-up with ϑ L = ϑ R , there is no k → −k symmetry. This is however an artifact of explicitly breaking the L ↔ R symmetry in order to mimic a multi-angle matter effect. These artifacts are avoided in a truly multi-angle scenario with L ↔ R symmetric velocity distributions of many modes. Note, however, the characteristic equation is real and will only give real or complex-conjugate eigenvalues.  The eigenvalues for this model can be calculated in a closed form, but these general expressions are not particularly illuminating. So we numerically evaluate these eigenvalues, and check if any of the eigenvalues have a large imaginary part that would signal rapid flavor conversions. Fig. 2, shows the largest imaginary part of any of the eigenvalues Ω, for a range of λ and µ and for specific values of p and k. We see that for certain large values of λ and µ, for which there were no instabilities in the stationary and homogeneous scenario (top-left panel), become unstable to a pulsating mode with large p (right panels) as the "nose-like" feature extends into these new regions. For a small number of modes as we have chosen here, a similar effect is produced by a choosing a large k (bottom-left panel) but this cancellation becomes very fine-tuned when more velocity modes are included.
One important observation concerns the interplay of matter and temporal instability, i.e., λ and p. A crucial feature of the temporal instabilities is that the eigenvalues depend on λ only through the combination λ + µ v − p, and there is always a p that can remove the effect of matter for all modes, allowing an instability for arbitrary λ. The temporal instability effectively undoes the phase dispersion introduced by matter and instabilities can grow unhindered by matter effects. This is seen in the right panels of Fig. 2. Such an effect is possible also for stationary but spatially inhomogeneous modes with k = 0 (bottom-left panel in Fig. 2), but genuinely multi-angle matter effects with many modes will not be completely suppressed by them, unlike what we see here. This interplay remains true at a nonlinear level, which is best seen using Eq. (2.15). If one looks at a particular combination of the transverse components of the polarization vector, i.e., where the ellipsis refers to terms not depending on λ or p. The amount of flavour conversion is determined by ∂ z |T p,k ω,v |, which clearly depends on λ only via λ−p. Thus at a fully nonlinear level, a temporal pulsation p counteracts the effect of matter density λ.

Numerical examples
In this section we present the flavor evolution for the model described in the last section. We calculate the fully nonlinear evolution of the polarization vectors given by Eq. (2.15), and compare it to the predictions of linearized stability analysis given by Eq. (2.19). The off-diagonal component of the density matrix ω,v , is proportional to g ω,v S p,k ω,v = T p,k ω,v . In our numerical calculations, we take |T 0,0 ω,v | = 10 −7 for the p = 0, k = 0 mode, and all components of all other Fourier modes of the polarization vectors initially get seeded by numerical noise related to the accuracy of the numerical solution, i.e., O(10 −12 ). Thus, the initial conditions are |S p,k ω,v | √ 2 × 10 −12 for all nonzero p. Their growth is predicted using the eigenvalues found through the linear analysis. The notation eµ n p ,n k is used to refer to both T p,k ω,v and g ω,v S p,k ω,v , and measures the degree of flavor conversion in a given Fourier mode 1 . In general, one should consider the flavor evolution of different Fourier modes associated with both spatial and temporal inhomogeneities. However, this would be numerically a rather challenging task. Therefore, we select a particular mode in space, i.e., the homogeneous mode with k = 0, and follow only the development of the Fourier modes in the n p space, suppressing the n k index henceforth 2 . In order to make p dimensionless, it is expressed in multiples of the effective matter potentialλ at z = 0, i.e., All other quantities are measured in units of ω 0 . We include 300 p modes, but with the modes n p > 200 being kept empty. This trick avoids "spectral blocking" that leads to a spurious 1 The careful reader will notice that we drop the (ω, v) index here. Indeed, we take the (+ω 0 , v R ) mode from nonlinear calculations and compare it with the fastest mode in linear theory. The different modes are all nonlinearly coupled, and we expect the growth rate to be dominated by the fastest mode. 2 We thank Georg Raffelt for pointing out that one can consistently solve for only the k = 0 homogeneous Fourier mode while ignoring other modes. Considering any other n k requires tracking the cascade in n k .
rise of the Fourier coefficients at large n p due to truncation of the tower of equations [37]. We observed this effect if we do not use this prescription. However, we checked that using different number of modes but taking the last 1/3 of them empty, it does not affect the behavior of the modes which are not close to the largest n p . In the following, the corresponding Fourier modes for the nonlinear computation and the linearized stability analysis calculation are shown and compared.
In Fig. 3, the amplitudes of the off-diagonal components of the different n p modes, • The choice τ λ = ∞, corresponds to a constant matter density, λ = λ 0 . Here, in agreement with the stability analysis, we find that the Fourier modes around n p ∼ 100 are excited first (at z 2). Indeed, for them the self-induced pulsation compensates the phase dispersion associated with the matter term (p λ ).
At larger z the temporal instability propagates in the Fourier space, where scales both larger and smaller than n p ∼ 100 start to get excited due to the non-linear interaction among the different modes [c.f. Eq. (2.9)]. At z = 30 all the modes with n p 50 are sufficiently large to be nonlinear.
The finger-like feature just above n p ≈ 100 arises due to the growth rate being zero for those modes. For those values of p one finds only real eigenvalues.
• For τ λ = 30, the temporal instability starts at n p ∼ 100. However, the finger-like feature just above n p ≈ 100 now sweeps into lower n p due to the decline of λ, and these modes go off-resonance and their growth slows down.
The cancellation of the matter suppression shifts to lower n p at larger z. As a result, the modes that are first to have A eµ −1 are n p ∼ 80 at z 5. Thereafter, the instability diffuses in the Fourier space.
With respect to the previous case of constant λ, lower n p are favored, because they are excited by both the non-linear effect and by the p λ condition that is satisfied at lower p.
• Finally, in the case with τ λ = 10, the rapid decline of λ enhances the features of the previous case. In particular, the reduced adiabaticity of the evolution implies that the instability is reduced. Modes that satisfy the p λ condition go rapidly off-resonance due to the rapid change of λ. As a consequence, their further growth gets inhibited.  The instability develops only at z 10 for n p ∼ 50. Modes with n p 50 grow till A eµ −1, when nonlinearity takes over, while larger n p modes are suppressed by the large p that mimics the matter suppression.
In Fig. 4, we show that the above observations are justified by comparing the results of the nonlinear computation (continuous curves) with the predictions of linearized stability analysis (dashed curves) for the amplitudes A eµ of specific n p modes (n p = 0, 50, and 90).
The key points to be noted are • In the case of τ λ = ∞, the mode at n p = 90 is the most unstable and grows by ∼ 10 orders of magnitude till z 6. In the same range the mode at n p = 50 grows by only 2 orders of magnitude, while the n p = 0 mode is stable. Till this distance from the source there is perfect agreement between the non-linear and the linear evolution for these three modes.
At larger z when the mode with n p = 90 has A eµ −1. Its evolution becomes nonlinear and produces significant deviations with respect to the linear case due to the interactions of different Fourier modes. In particular, the growth of the n p = 90 mode slows down with respect to the linear prediction and instead the mode at n p = 50 has a rapid increase of 5 orders of magnitude around z = 7, which then continues to rise and reaches A eµ −2 at z = 20. The mode at n p = 0 begins to oscillate around its initial position at z 7 and starts growing at z 10. Its final value is only 5 orders of magnitude larger than its initial value, as the cascade in the Fourier space has not reached the n p = 0 mode in the considered z range.
• With τ λ = 30, we realize that until z 3 the behaviour of the three modes is similar to the one found in the case of τ λ = ∞. In particular, the n p = 90 mode is still the most unstable, and grows by 4 orders of magnitude. However, this mode then undergoes rapid oscillations from z = 2 to z = 4 as p is just larger thanλ, where the linear theory correctly predicts no growth (the dip feature here is the same as the finger in the Fig. 3). The agreement between linear and non-linear results is good till z = 4 for the mode at n p = 90.
However, at larger z, as the matter potential becomes smaller than p, this mode grows less rapidly. As soon as the A eµ is affected by the rapid (non-adiabatic) oscillations owing to the dip-like feature in the growth rates, the linear solution does not reproduce this behavior and starts to deviate from the non-linear behavior. As a consequence at larger z the linear solution overestimates the real growth.
On the other hand, the mode n p = 50 steadily grows and is enhanced by non-linear effects at z 10. For the mode n p = 50 the agreement is good till z = 10, when the mode reaches A eµ 10 −2 and non-linearities begin to play an important role. Again, the growth of the n p = 0 mode is severely underestimated by linear theory, as it ignores the possibility nonlinear coupling of modes beyond z = 10 that leads to an increase of 4 order of magnitudes by z = 20.
• Finally, for τ λ = 10 we see that due to the fast decline of λ, the growth of the mode n p = 90 is inhibited already at z = 1, while the mode with n p = 50 steadily increases till z = 6. Then it oscillates (when p is just larger thanλ) before resuming its growth at z 7. The agreement between the linear and non-linear result is good till z = 1 for the mode at n p = 90, and for z = 5 for the mode at n p = 50. In both the cases the discrepancy occurs when the mode deviates from a monotonic growth and starts to show these fast oscillations.
We see that unlike in the two previous cases, the mode with n p = 0 is not linearly stable for z > 7. In this case its linear growth is larger than the non-linear one. In particular, at z = 12 they disagree by 2 orders of magnitude.
From these comparisons we conclude that the linear solution is a good approximation of the non-linear evolution only until nonlinearities remain small, so that the interaction among the different modes can be neglected, as evident from the case with τ λ = ∞. Once nonlinearities kick-in, the rate of growth is severely underestimated by linear theory and the growth of flavor conversions appears to be qualitatively larger. Furthermore, in the case of variable λ the adiabaticity of the flavor evolution plays a significant role. Indeed, when the numerical solution shows non-adiabatic features resulting from resonances (rapid oscillations in the dip-like feature), the deviations with respect to the linear approximation become significant.

Declining µ and λ
As a mock-up for a more realistic SN-like scenario, it is interesting to consider situations in which both µ and λ decline as a function of z. Fig. 2 provides some guidance as to how λ and µ should vary with z, so that there is no instability if p and k are zero but an instability develops when they are nonzero. Using this guidance, we consider that with the following parameters in units of ω 0 , The functions are chosen such that λ and µ for model A varies as per the path shown as the dashed white line in Fig. 2. However, how fast the system evolves along the curve is determined by τ λ and τ µ . It is important to emphasize that while model B also follows almost the same path on the λ − µ plane, it passes through the instability much quicker than model A, owing to its smaller scale heights. Note that no instabilities would be encountered for either model if k = 0 and p = 0.
In Fig. 5 the amplitudes of flavor conversions, A eµ , are shown at distance z for various Fourier modes n p for the models A (left) and B (right). For the model A, we realize that the instability starts at z 2 from modes at n p 100 and develops at smaller scales, reaching the mode at n p = 0 at z = 30. On the other hand for model B, due to the reduced adiabaticity in the flavor evolution the instability never grows significantly. This result tells that situations that present the same instability footprint, may have a rather different flavor evolution depending on net time spent in the unstable region and on the adiabaticity of evolution. If the profiles decline too fast the instability is inhibited. Fig. 6 is the analogue of Fig. 4, obtained by considering models A (left panel) and B (right panel). The results here are similar to what is presented in the previous section in the case of constant µ. In particular, for model A we find that indeed the p = 0 stationary flavor conversion mode can be nonlinearly excited by its coupling to the fastest p = 0 modes by z = 11, which would not be predicted in a purely linear theory. That aside, agreement between linear and non-linear evolution is good till the non-linear evolution presents nonadiabatic features (fast oscillations in the dip-like feature) which are not accounted in the linear theory.

Conclusions
In our work we have studied the development of the temporal instability in self-induced flavor evolution in the line model. We have found that space and time inhomogeneities can  excite flavor conversions at small spatial scales and can spontaneously generate a pulsating component in the flavor composition. In particular, the pulsation can compensate the phase dispersion associated with a large matter term that would otherwise suppress the flavor conversions.
We have considered several scenarios, with varying matter and neutrino densities, finding that the time instability propagates in the Fourier space exciting in cascade different Fourier modes. The development of this cascade crucially depends on the adiabaticity of the matter and neutrino potential. We have shown that the linearized solution of the equations of motion can be a valuable tool to find the growing modes and to predict the flavor evolution till non-linearities and non-adiabatic features become important. However, the linear theory does not account for the interaction among different Fourier modes that can enhance or reduce the development of the pulsating modes.
This raises the possibility that a specific pulsating mode may experience large growth and eventually lead to flavor decoherence, especially if there is a slowly varying or relatively flat region of the matter density profile that extends for ∼ 10 2 km, as sometimes found to be the case behind the shock in SN simulations (see, e.g., [38,39]), that allows long-term adiabatic growth of a single mode. Our model of course is too simple to gain a definite answer on the development of such effects in the flavor conversions on SN neutrinos. In particular, in SN neutrino context one has to take into account true multi-angle effects that would make the simulation of the non-linear flavor evolution extremely challenging. However, one crucial point that is illustrated by our study is that one may need to account for these details of the matter and neutrino profiles in a SN in order to get a definitive answer on the nature of neutrino flavor conversions in SNe.