Effects of nonstandard neutrino self-interactions and magnetic moment on collective Majorana neutrino oscillations

We derive the effective Hamiltonian describing collective oscillations of Majorana neutrinos with a transition magnetic moment, allowing for the presence of scalar and pseudoscalar nonstandard neutrino self-interactions (NSSIs). Using this Hamiltonian, we analyze new flavor instability channels of collective oscillations in a core-collapse supernova environment that are opened up by a small but nonzero neutrino magnetic moment. It turns out that, contrary to certain claims in the literature, within the minimally extended Standard Model (i.e., without NSSIs), no new instabilities arise within the linear order, nor do they produce any observable signatures in the neutrino flavor-energy spectra, at least for magnetic moments up to 10 μB. On the other hand, in the presence of NSSIs, new fast and slow instabilities mixing neutrinos and antineutrinos appear, which show up in the spectra even for tiny magnetic moments of the order of 10 μB, leading to considerable distortions of the spectra and nonstandard spectral splits. We study sensitivity of collective oscillations to these, NSSI-induced instabilities in detail and discuss the observability of the NSSI couplings triggering them.


I. INTRODUCTION
During the most violent phase of a core-collapse supernova explosion, emitted neutrinos carry away the lion's share of the explosion energy and turn out to be the first signal from a newborn supernova coming from its innermost regions [1,2]. An observation of supernova neutrinos could thus yield priceless information both on the physical processes at work during the explosion and on the physics of the neutrino itself in such an extreme environment, possibly unveiling signatures of certain Beyond-Standard-Model (BSM) phenomena [3][4][5][6][7]. Indeed, thanks to the progress in experimental techniques of neutrino detection made since the SN 1987A event, modern neutrino observatories, as well as those planned to start operating in the coming years, should be able to yield a number of events quite sufficient to draw the neutrino energy spectra. For instance, for a 10 kpc-away explosion, the JUNO detector would yield as many as several thousand events in both the neutrino and the antineutrino channels, including the inverse-beta-decay detection channel with an unprecedented energy resolution about 1% [8]. Inevitably though, the energy spectra thus observed would be a result of neutrino oscillations on the way from the supernova core (the protoneutron star) to the detector, so, in order to 'decypher' the spectra, one should be able to 'rewind' these oscillations back to the stellar interior.
For a typical core-collapse supernova, the neutrino densities turn out to be so high next to their last scattering surface (also referred to as the neutrino sphere) that neutrino-neutrino forward scattering processes become important for the flavor evolution, ushering into the physics of collective neutrino oscillations [9]. This self-induced, nonlinear process shapes the flavor-energy spectra in the region deep under the Mikheev-Smirnov-Wolfenstein (MSW) resonance surface, i.e., where noncollective oscillations would be suppressed by a gigantic matter potential [10,11]. One of the main drivers of collective flavor transformations is the instabilities, which are intrinsically present in the nonlinear evolution equations on collective oscillations [9,[12][13][14][15][16]. It has been identified, indeed, that a hierarchy of these instabilities, both in the linearized and fully nonlinear regimes, can lead to specific neutrino spectral features, such as the so-called spectral swaps (splits) and the behavior resembling turbulence [17][18][19].
Instabilities often arise when new degrees of freedom come into play. For example, fast instabilities show up beyond the so-called single-angle scheme [12,14,15,20]; lifting the assumption of translational invariance of the solution leads to turbulent flavor patterns breaking this invariance on a wide range of spatial scales [18,19,21]. In view of this, it seems appealing to analyze the effect of the neutrino anomalous magnetic moment [22], which mixes the two helicity states, and thus doubles the number of nontrivially interacting degrees of freedom, on the spectrum of instabilities. This also sounds natural because of superstrong magnetic fields typical for collapsing stars. This problem has been studied in a number of papers, including the derivation of the effective Hamiltonian [6,23,24] and the impact on the flavor evolution within the single-angle [6,25] and multiangle frameworks [26]. Notably, the effective Hamiltonians derived/used in Refs. [23][24][25][26] and Ref. [6] are different in the self-interaction term and lead to drastically different flavor evolutions: in the latter case, new types of instabilities arise that strongly deform the neutrino flavor-energy spectra even for tiny values of the magnetic moment [6].
Certainly, it seems interesting not only to settle the question of the correct Hamiltonian, but rather to ask a more general question: are there any other factors producing neutrino flavor signatures in an interplay with the effects of a (tiny) neutrino magnetic moment? At least one answer to this question is discussed in the present paper. Namely, we study the effect of the so-called NonStandard neutrino Self-Interactions (NSSIs, see Refs. [27][28][29] for a review) on the collective flavor evolution of Majorana neutrinos triggered by the nonzero neutrino transition magnetic moment. Recently, NSSIs, i.e., four-fermion neutrino-neutrino interactions that are absent in the electroweak sector of the Standard Model, have attracted a lot of attention for being able to affect the observable neutrino spectra [30][31][32][33]. One can distinguish two classes of NSSIs, those with flavor-dependent V-A interactions and those with a non-V-A tensor structure of the four-fermion interaction. The NSSIs of the second class, specifically, those involving scalar (νν) 2 and pseudoscalar (νγ 5 ν) 2 terms, are especially interesting to us in our context, since, as we show below, their presence opens up a new instability channel in oscillations of Majorana neutrinos with a nonzero magnetic moment. It is worth admitting here that another instability channel due to scalar-pseudoscalar NSSIs has already been studied in Ref. [33]; however, the corresponding unstable modes do not mix the neutrino helicities and the nonvanishing magnetic moment is not crucial for their development.
According to our analysis which follows, it turns out that presence of a scalar-pseudoscalar NSSI violently deforms the neutrino spectra even for minuscule values of the magnetic moment, at least down to 10 −24 µ B , while in the absence of NSSIs (i.e., for the neutrino interactions within the Standard Model), the effect of the magnetic moment is virtually unobservable even for much greater magnetic moments. This means that supernova neutrino spectra can be used to probe the presence of (pseudo)scalar NSSIs, possibly stemming from an exchange of (pseudo)scalar BSM particles. The steps we take within our analysis of the NSSI-induced instabilities are presented in the following sections. Namely, in Sec. II, we rederive the Hamiltonian for collective Majorana neutrino oscillations with NSSIs and a nonzero neutrino magnetic moment and compare the no-NSSI case with the results of previous derivations. Further, in Sec. III, we carry out a linear stability analysis of our flavor evolution equations, also showing why in the absence of NSSIs, the effect of the magnetic moment should be small. Focusing then on the most interesting NSSI-induced unstable modes mixing the two neutrino helicities, we determine their growth rates, revealing fast and slow branches, as well as a specific intermediate, matter density dependent one. The linear analysis is accompanied by a full numerical simulation in Sec. IV, to study what happens beyond the linear stability regime and to determine the sensitivities of the spectra to the neutrino magnetic moment and the NSSI couplings. The final section V summarizes the results obtained and discusses their possible generalizations. In the Appendix, several identities are proved or listed that we use in the derivation of the effective Hamiltonian in the presence of NSSIs.

II. EVOLUTION EQUATION FOR COLLECTIVE OSCILLATIONS IN THE PRESENCE OF NSSIS
Let us now settle the question of the effective Hamiltonian for collective oscillations of Majorana neutrinos with nonzero magnetic moment we mentioned above and, more importantly, derive the terms in this Hamiltonian that describe the scalar and pseudoscalar NSSIs. For that, let us derive the evolution equations for the neutrino flavor density matrix accounting for the forward scattering processes, starting from field theory. We use relativistic units = c = 1 throughout the paper. We consider Majorana neutrinos in a small representative volume V , interacting with background electrons, protons, and neutrons, plus an external electromagnetic field F µν (x). After the electroweak symmetry is broken, the terms in the Lagrangian that contribute to neutrino forward scattering read [2,22,33,34] L ν = L (2) vac + L Here ν a (x), a = 1, . . . , N f , are the neutrino fields with Majorana masses m a and e(x), p(x), n(x) are the fields describing electrons, protons, and neutrons, respectively. Further, the conjugate of the Majorana field isν a ≡ ν T a C, where C = −iγ 2 γ 0 , while for the Dirac fields ℓ = e, p, n, the conjugate isl = ℓ † γ 0 . The Dirac matrices γ µ  [34,35], where the summation over the mass indices a, b, . . . and flavor indices f, f ′ , . . . is assumed, the latter take on e, µ, τ values in the three-flavor case and e, x values in the two-flavor case. The electromagnetic interaction term (4) contains the neutrino magnetic moment matrix m ab , which is real and antisymmetric in the Majorana case, so that only transition magnetic moments are allowed [22]. Finally, the quartic terms (5), (6) describe the electroweak (V-A) and the nonstandard (scalar/pseudoscalar, SP) neutrino interactions, and the latter is parameterized by two dimensionless NSSI couplings g S,P . Within the forward-scattering approximation, i.e., while momentum-changing processes of the form ν a (p) → ν b (p ′ ), ν a (p) →ν b (p ′ ) with p ′ = p play a minor role, we can ignore coherent superpositions of different neutrino momentum states (while keeping track of coherent superpositions of mass states due to oscillations). In other words, the neutrino state(s) |Φ we consider can be factorized into a infinite tensor product of state vectors |Φ p , each of them describing neutrinos with a fixed momentum p and belonging to the corresponding Fock space, Here, for the sake of the calculations which will follow, we adopt a notation A = (aα), B = (bβ), . . . for the creation operatorsâ † Ap , with the Latin indices a, b, . . . = 1, . . . , N f denoting the neutrino mass states and the Greek ones α, β, . . . = ±1 standing for the neutrino helicities (multiplied by two). As mentioned above, the state vector |Φ p describes neutrinos with momentum p, whose total number N p ∈ {0, 1, . . . , 2N f } is limited due to the Pauli exclusion principle; moreover, in the forward scattering regime, these numbers are conserved for every individual p and Thus, the time evolution of the neutrino state affects only the set of c-number-valued coefficients C A1...AN p (p). It is much more convenient, however, to work with a 2N f × 2N f neutrino flavor density matrix instead of these coefficients whereĤ is the Hamiltonian corresponding to the Lagrangian (1). Our task now is to transform Ehrenfest equation (12) into an effective von Neumann equation for the density matrix where h(p) is the desired 2N f × 2N f c-number matrix of the effective Hamiltonian, possibly depending on the density matrix ρ. In the derivation of the effective Hamiltonian, which follows, let us consider the contributions of the quadratic and quartic parts ofĤ separately.
In the ultrarelativistic approximation we are using, the neutrino field operators ν a (x) read where the 3-momentum p ∈ (2π/L)Z 3 is quantized in the normalization volume V ≡ L × L × L, the corresponding 4-momentum being p µ a ≈ (|p| + m 2 a /2|p|, p) in the ath mass state. The neutrino annihilation operatorsâ aαp enter together with the plane-wave solutions u α (p)e −ipa·x / √ V describing particles with helicity α/2 = ±1/2. In what follows, for brevity, we will refer to negative-and positive-helicity states as neutrino and antineutrino states, respectively, and also refer to α = ±1 as the helicity values, instead of the rigorous helicities α/2 = ±1/2. Some obvious expressions for the polarization bispinors u α (p), as well as for their bilinear combinations we will need below, are listed in Appendix A.
The quadratic part of the HamiltonianĤ in the Ehrenfest equation (12), up to a insignificant multiple of the identity operator, can obviously be written aŝ Cqâ Dr +ââ terms +â †â † terms.
The second and third terms produce vanishing contributions to the Ehrenfest equation, since |Φ possesses definite particle numbers N p in all momentum states. As for the number-conservingâ †â term, the corresponding contribution to the time derivative of the density matrix is evaluated straightforwardly and has the desired von Neumann structure (13) (the ⊃ sign means that quartic terms have been omitted on the r.h.s.). Forward scattering clearly manifests itself in only the diagonal (q = r) entries of the coefficient function λ(q, r) affecting the evolution of the density matrix. Thus, to find the noncollective part of the desired effective Hamiltonian h(p) in Eq. (13), one should simply list the threeâ †â terms inĤ (2) coming from the three Lagrangians L (2) vac , L mat , and L (2) AMM describing the vacuum neutrino mixing and their interaction with matter and magnetic field. The first term gives nothing but a free Hamiltonian, where (M 2 ) cd = m 2 c δ cd is the neutrino mass-squared matrix (in the mass basis here) and, as usual, the term proportional to the identity matrix does not contribute to the commutator in Eq. (16) and can be omitted. The matrix notation we have adopted in Eq. (18) and to be used hereinafter is as follows: the two lines/columns of a 2N f × 2N f matrix λ vac correspond to the two helicities γ, δ = −1, +1, following a pattern In other words, the first and the second diagonal blocks describe neutrinos and antineutrinos, respectively, while the two off-diagonal blocks describe coherent mixtures of neutrinos and antineutrinos.
The matter (MSW) term inĤ (2) can be retrieved by treating background matter and the corresponding e, p, n operators in the mean-field fashion The background matter current J µ ab above is directly extracted from Eq. (3) and after averaging over neutral nonmoving nonmagnetized matter gives where the projector onto the electron neutrino state (P e ) ab ≡ U * f a · δ f,e δ f ′ ,e · U f ′ b and n e,n are the electron and neutron number densities, respectively. Indeed, expectations of all axial vectors vanish, while for polar vectors, l γ µ ℓ = δ µ 0 n ℓ , ℓ = e, p, n. It remains now to directly substitute the neutrino operators (14) into (20) and extract the terms contributing to forward scattering (for the bilinear expressions of the formū(. . .)u, see Appendix A) Finally, the magnetic moment term is evaluated in the same way givinĝ where complex vectors ζ ± (p) = 1 √ 2 χ † ∓ (p)σχ ± (p) are matrix elements of the Pauli matrices between states with opposite helicities (see Appendix A for the definition details and properties of these vectors).
Thus, we have found the noncollective part of the effective neutrino Hamiltonian to be λ vac (p, p) + λ mat (p, p) + λ AMM (p, p), and we can now resort to the quartic terms in the neutrino HamiltonianĤ. These terms, after commutation withâ † Bpâ Ap in the Ehrenfest equation (12), lead to quartic expectations of the form Φ|â †â †ââ |Φ . In order to close our system of equations on the density matrix, namely, to obtain (13), we will now apply a kind of Wick's theorem and express such quartic expectations in terms of quadratic ones Φ|â † Bpâ Ap |Φ = ρ AB (p). Let us start from the contribution of the electroweak Lagrangian (5) having a V-A structure and insert the corresponding Hamiltonian into Eq. (12) i First of all, commutation withâ † Bpâ Ap transforms creation and annihilation operators into creation and annihilation ones, respectively, i.e., it does not spoil the normal ordering. Thus, the commutator can safely be put inside the normal ordering, leading to an expectation of Next, due to symmetries of bilinear expressions in Majorana fields (see, e.g., Ref.
the V-A currents can be replaced by axial vector ones, . Moreover, by virtue of the same symmetry properties, the latter one is a symmetric quadratic form in ν d , and its commutator withâ † Bpâ Ap leads to two identical terms, yielding −ν d (x)γ µ γ 5 ν d (x),â † Bpâ Ap . As a result, the V-A contribution to the evolution equation for the density matrix takes the form where we have adopted an approximation p µ d ≈ p µ ≡ (|p|, p) and also symbolically denoted A quartic expectation value we have encountered can now be (approximately) reduced to products of quadratic expectations, or contractions, by virtue of the Wick's theorem (see Appendix B for details) where a contraction of two spinor fields ϕ i χ j ≡ Φ| : ϕ i χ j : |Φ , with i, j = 1, 2, 3, 4 numbering the components of a bispinor. To avoid using the index notation, however, we note that the second and the third terms in the above expression are equal due to (28) and apply the Fierz identity to them [36], arriving at Φ| :φγ µ γ 5 χ ·ψγ µ γ 5 ω : |Φ ≈φγ µ γ 5 χ ·ψγ µ γ 5 ω + 2φω ·ψχ + 2φγ 5 ω ·ψγ 5 χ +φγ µ ω ·ψγ µ χ +φγ µ γ 5 ω ·ψγ µ γ 5 χ. (31) Now, every of these contractions is a sum of expectation values of the form Φ|â †â |Φ , which is nothing but a specific entry of the density matrix ρ, namely, The above identities are obtained in a straightforward way from the neutrino operator (14) and the bilinears (A6)-(A11). In the first two identities, we have introduced a matrix distinguishing neutrinos and antineutrinos while another matrix featuring in identities (35) and (36) is a difference between the density matrix and the charge conjugate of its transpose The charge conjugation, as defined here, simply swaps the neutrino lines/columns with the antineutrino lines/columns of the density matrix. Finally, after substituting the contractions listed above into the Fierz transformed Wick's theorem (31) and using the latter expression in the Ehrenfest equation (29), we obtain a contribution to the evolution of neutrino density matrix due to electroweak neutrino-neutrino interactions The ℘ diag (q) matrix above is a block-diagonal part of ℘(q) It is now time to find the term in the evolution equation coming from the Scalar (S) and Pseudoscalar (P) nonstandard neutrino interactions. This is done following exactly the same steps as in the V-A interaction case, but starting from Lagrangian (6). Namely, the S/P contribution to the time derivative of the density matrix reads and after applying the Wick's theorem, we arrive at Note that the first contractions in each of the two parentheses vanish due to Eq. (A6). For the second contractions, we use the Fierz identity [36], omitting the vanishing scalar and pseudoscalar terms in it (34), and arrive at where g ± ≡ (g S ± g P )/2. In fact, the vector and the axial vector terms have already been evaluated above, so that it remains only to evaluate the tensor one Again, working in the same way as we did in the case of V-A interaction, we transform Eq. (47) into where we have expressed the scalar products ζ ± (p) · ζ ± (p) = e ±iΓ(p,q) (1 −p ·q)/2 in terms of a single complex phase Γ(p,q) (see Appendix A). The block-off-diagonal part of the matrix above is defined analogously to Eq. (44) In fact, this block-off-diagonal part anticommutes with G, which underpins hermiticity of the g + part of the Hamiltonian. Further, the complex phase Γ(p,q) comes from the phases included in the helicity eigenstates χ ± (p), χ ± (q); one would therefore like to rephase helicity eigenstates, thus getting rid of the e iΓG term in the NSSI Hamiltonian. However, this is impossible to achieve on the whole Bloch spherep,q ∈ S 2 , since the parallel transport of the helicity eigenstates across this sphere is characterized with a nonzero Berry curvature [37].
As we observe, both in the electroweak and the nonstandard scalar-pseudoscalar case, the time derivative of the neutrino flavor density matrix ρ(p) can be cast into a von Neumann form (13), i.e., a form of a commutator with a hermitian matrix h(p). Writing all the contributions to this time derivative together (see Eqs. (18), (23), (25), (43), (52)), we finally arrive at the desired evolution equation on ρ where one can show that is indeed a complex number with the absolute value equal to the strength of the magnetic field across the neutrino momentum (in accordance with the notation we use). The equations of motion possess a gauge freedom with respect to rephasing the helicity eigenstates where γ(p) is a real function on a sphere describing the gauge transformation. This transformation virtually tunes the relative phases of neutrino and antineutrino states and becomes trivial for a block-diagonal density matrix, i.e., if neutrinos do not coherently mix with antineutrinos. It is also instructive here to emphasize that the off-diagonal blocks of ρ(p) do not behave as scalars under rotations. Indeed, the spinor representation of a rotation R(nϑ) ∈ SO(3) around the n axis on ϑ radians has the form from which one readily obtains the transformation of the neutrino annihilation operatorŝ where the phase Ξ is connected with the complex phases chosen in the helicity eigenstates, Just like in the case of a gauge transformation, the diagonal blocks are left intact, while the off-diagonal ones get rephased. In particular, one can demonstrate that for the neutrino-antineutrino block, i.e., its product with the ζ − vector transforms as a vector field without additional phases. As a result, for an isotropic neutrino gas (ρ ′ (p) = ρ(p)), the above equation tells us that ρ νν (p)ζ − (p) is a spherically-symmetric vector field. However, because p · ζ − (p) = 0 (see Appendix A), this requires ρ νν (p) = 0. That is, nontrivial neutrino-antineutrino mixing violates isotropy, so the NSSI Hamiltonian (57) loses its g + term when applied to an isotropic neutrino gas.
In the next two sections, we will study the effect of the nontrivial neutrino NSSIs on the flavor evolution in the neutrino bulb model using the so-called single-angle scheme [9]. To obtain the corresponding evolution equations, we focus on stationary processes with propagating neutrinos and make a replacement transforming the von Neumann equation (13) into a quantum Liouville equation Next, construction of the single-angle scheme should include introduction of the so-called geometric factor D(r/R ν ), (partially) accounting for the non-spherically symmetric distribution of the neutrino numbers far from the neutrino sphere; even in the Standard-Model case, there exist several versions of it [9,17,38]. In our case, since the same factor 1 −p ·q appears in front of both the electroweak and NSSI g − terms in the collective part of the Hamiltonian (57), it is natural to equip their single-angle counterparts with the same geometric factor. Regarding the g + term having an extra e iΓ(p,q)G factor, the same may not be the case, however, we will keep the same geometric factor, qualitatively relying upon the argument that far from the neutrino sphere the solid angle spanned by the neutrino momenta q is small and one can rephase the helicity eigenstates within it to approximately eliminate the phase factor e iΓG . A more rigorous discussion of the geometric factor for the g + term of Hamiltonian is to be given elsewhere.
With these points in mind, we can now write down the evolution equations of the single-angle scheme. Namely, in the flavor basis and the two-flavor approximation, a neutrino with energy E, escaping the protoneutron star radially from its neutrino sphere r = R ν outwards, is described by equations where ∆m 2 is the mass squared difference, η = ±1 marks the normal/inverted mass hierarchy, θ is the vacuum mixing angle, and µ 12 is the transition magnetic moment (the diagonal entries µ 11 = µ 22 = 0 for Majorana neutrinos [22]).
To get rid of additional factors, we have renormalized the density matrix so that ∞ 0 tr ρ E (r)dE = 1; after such a renormalization, the total neutrino number density n ν (r) appears in the nonlinear collective term. The initial condition for the above system of equations is usually placed at the neutrino sphere, where different neutrino flavors/helicities are assumed to be thermalized with well-defined energy spectra with the normalization convention (69) requiring that Equations similar to the above ones for collective oscillations of Majorana neutrinos have been derived earlier in a number of contexts. First of all, V-A neutrino-neutrino interactions were included into the effective Hamiltonian in Ref. [23], allowing for a nontrivial structure of this interaction in the flavor space. The result we obtained here agrees with Ref. [23] in the absence of nonstandard interactions (g ± = 0). In a recent paper [33], the collective NSSI Hamiltonian was derived in the absence of the neutrino magnetic moment. One also observes our agreement with this paper for µ 12 = 0, however, our derivation goes qualitatively beyond this special case. Namely, interaction with the external magnetic field via the magnetic moment introduces mixing between neutrinos and antineutrinos, so that their states cannot be described by two N f × N f density matrices anymore, but a single 2N f × 2N f matrix is necessary to refer to both neutrino and antineutrino flavors and their coherent mixtures. In fact, in the µ 12 = 0 (or B = 0) case, our density matrix becomes block-diagonal and its ρ νν and ρνν blocks representing neutrinos and antineutrinos, respectively, do obey the evolution equations derived in Ref. [33].
We should also mention two papers [6] by A. de Gouvêa and S. Shalgar, in which the single-angle scheme for collective oscillations of Majorana neutrinos was analyzed in the µ 12 = 0 case. The papers analyzed the electroweak (V-A) four-fermion neutrino interaction, arriving at the neutrino self-interaction Hamiltonian of the form Note that such an interaction, in principle, includes off-diagonal blocks mixing neutrino and antineutrino states. In contrast, our collective Hamiltonian (68) does not have such blocks in the Standard-Model regime g ± = 0, even if the density matrix ρ E (r) contains them. As we will see below, this makes a considerable difference in the instability spectrum of the two Hamiltonians and, as a result, in the evolution of the neutrino spectra in magnetic field. In a nutshell, collective oscillations governed by interaction Hamiltonian (72) favor neutrino-antineutrino mixing due to the presence of nontrivial off-diagonal blocks, thus, a tiny mixing introduced into the system by the magnetic moment term is followed by its exponential growth due to instabilities; however, nothing like occurs in the Standard-Model regime of our evolution equations (65), since the only source of neutrino-antineutrino mixing is the (linear, noncollective) magnetic-moment term.
In any case, despite the disagreement of the results of Ref. [6] with ours, as well as with the earlier Ref. [23], the Hamiltonian presented in the former paper can be treated as yet another type of NSSI containing nontrivial off-diagonal blocks. Moreover, it is spectacular that both the equations from Ref. [6] and our equations (65) conserve the block-diagonality property of the density matrix ρ E (r) along the trajectory, provided that µ 12 = 0 and that the initial condition ρ E (R ν ) is block-diagonal (e.g., the one we chose (70)). In this case, ρ = ρ diag = diag(ρ νν , ρνν) and Hamiltonian (72) can be replaced by which is nothing but the well-known expression within the Standard Model, also valid for Dirac neutrinos [9]. Our Hamiltonian (65) also takes the above conventional form, when additionally, g − = 0, i.e., the scalar and the pseudoscalar couplings are equal in the NSSI Lagrangian (6). If, in contrast, g − = 0, then scalar-pseudoscalar NSSI terms are present in the diagonal blocks of the density matrix and survive in the µ 12 = 0 regime. Their effect has recently been analyzed in Ref. [33], revealing considerable deviations of the neutrino spectra from the predictions of the electroweak theory for nonvanishing g − . As a result, constraints can be placed on the g − coupling from the possible measurements of the neutrino spectra from a supernova, regardless of the magnitude of the neutrino magnetic moment. Our Hamiltonian, however, contains another, 'hidden' sector with the coupling g + , which produces no effect in the zero magnetic moment case. It is this coupling that is able to introduce block-off-diagonal terms into the selfinteraction Hamiltonian (68), possibly leading to nontrivial evolution of neutrino-antineutrino mixing. In the special case g − = 0, g + = −1, the corresponding interaction term takes a form resembling the Hamiltonian (72) from Ref. [6] h ), yet, it differs in the transposition of the off-diagonal part of the density matrix. As we will see below, however, the transposition leads to qualitatively different effects, and we have to conclude that the interaction (72) introduced in Ref. [6] cannot originate from scalar/pseudoscalar four-fermion interactions of Majorana neutrinos.
Finally, we would like to mention that we briefly analyzed the case of Dirac neutrinos in our recent paper [25], where coherent mixing between two (anti)neutrino helicity states was introduced by the magnetic moment. It turns out that the corresponding interaction Hamiltonian is also block-diagonal in the absence of NSSIs and, as a result, the effect of the magnetic moment on the neutrino spectra is suppressed, not being subject to instabilities. In the present paper, we are not focusing on the Dirac neutrinos, and the analysis which follows will be restricted to the Majorana case.

III. LINEAR STABILITY ANALYSIS
Before making numerical simulations of the single-angle scheme in various scenarios, it is instructive to take a look at the linear stability of the evolution equations. Our specific interest here regards the effect of a small neutrino magnetic moment; the crucial question is whether it can trigger new types of neutrino flavor instabilities that are hidden in the µ 12 = 0 regime. As we mentioned above, when µ 12 = 0 and the initial condition is of the form (70), the density matrix ρ (0) E (r) retains its block-diagonal form during the whole evolution r ≥ R ν . For such a matrix, the effect of a nonzero g + NSSI term is absent, while the effect of the block-diagonal g − term was studied in a recent paper [33]. Namely, in that paper, it was shown that the g − term is able to trigger instabilities even in the absence of neutrino magnetic moment. Now, on top of possibly nonzero g ± couplings, let us switch on an infinitesimal magnetic moment µ 12 → 0 and write down the equations of motion linearized in this small µ 12 quantity self is a block-diagonal matrix and so commutation with it does not mix the diagonal with off-diagonal contributions; the same applies to commutation with ρ (0) E ; finally, the 'source' term J E is purely blockoff-diagonal because of the form of h AMM . Note also that block-diagonal and block-off-diagonal parts of δh self are determined by the corresponding parts of δρ E , so that evolutions of δρ diag E (r) and δρ offdiag E (r), defined by the two above equations, do indeed decouple in the linear regime. Technically, this means, that one can search for block-diagonal and block-off-diagonal unstable modes separately.
If one sets g + to zero, i.e., turns off the block-off-diagonal part of interaction, then the second equation becomes which features a fixed Hamiltonian h vac,E +h mat +h (0) self . In other words, this equation virtually describes a noncollective flavor evolution and it is straightforward to show that where M F = tr(M † M) is the Frobenius norm of a matrix and J E (r ′ ) is the spectral norm. The latter being equal to |µ 12 B ⊥ (r ′ )| · ρ (0) E (r ′ ) up to a factor of order unity, we conclude that block-off-diagonal perturbations δρ offdiag grow linearly with distance, not exhibiting exponential growth. As for the block-diagonal part δρ diag of the perturbation, it obeys the same equation (78) as in the µ 12 = 0 case, when the density matrix is exactly block-diagonal at all r. Instabilities are well-known to exist here, both in the absence of NSSIs, g + = g − = 0 [9], and in their presence, g + = 0, g − = 0 [33], but they contain no signatures of a nonzero magnetic moment in the linear regime.
We are thus able to conclude that in the absence of a scalar-pseudoscalar NSSI, i.e., within the electroweak theory, a tiny magnetic moment of a Majorana neutrino does not trigger new unstable modes that could result in observable signatures in the neutrino spectra. In order to leave such a signature, the transition magnetic moment should be at least of the order of 1/|B ⊥ |L, where L is the distance covered by the neutrino (see Eq. (81)); for L ∼ R ν ∼ 50 km and |B ⊥ | ∼ 10 12 Gauss, this amounts to about 0.5 × 10 −15 µ B . In the next section, we will show numerically that this statement also holds true beyond the linear stability regime. Note also that the argument that has led to such a conclusion was based only on the block structure of the Hamiltonian matrix, thus, it can be repeated beyond the single-angle scheme, i.e., for the Hamiltonian (57).
Let us now turn on the block-off-diagonal NSSI, g + = 0. Even in this case, Eq. (78) tells us that the fate of blockdiagonal perturbations is determined solely by g − and provides no signatures of a nonzero µ 12 . The block-off-diagonal modes follow Eq. (79), which, in general, has to be treated numerically. However, to probe the possibility of unstable modes, it is usually instructive to carry out a Lyapunov stability analysis. For that, firstly, the source term in the evolution equation (79) is omitted but is replaced by a nontrivial initial perturbation δρ offdiag E (R ν ) on top of a diagonal reference matrix ρ (0) E (R ν ) of the form (70). Secondly, assuming that the mixing angle is virtually zero and (for our toy model) that neutrinos are monochromatic with energy E 0 , we arrive at a constant nonperturbed solution Finally, let us assume that inhomogeneities of the matter and neutrino number densities are negligible at the typical growth scale of unstable perturbations. Under these assumptions, we arrive at a linear homogeneous equation with constant coefficients governing the growth of perturbation δρ offdiag where µ = G F √ 2D(r/R ν )n ν (r) is the neutrino-neutrino coupling strength. Now the question of linear stability reduces to existence of non-real eigenvalues of a linear map L. Parameterizing the (block-off-diagonal) perturbation matrix as one writes the eigenvalue equation for L(δ̺) = λδ̺ in a vectorized form , −n e + n n 2 , n n 2 +2µ diag ∆s e + ∆s x , ∆s e + 2∆s x , −2∆s e + ∆s x , −∆s e − 2∆s x , where ξ = (ξ 1 , ξ 2 , . . . , ξ 8 ) T , ∆s e,x ≡ s νe,x − sν e,x , Ω ± = G F √ 2(n e − n n ) ± ηω, ω = ∆m 2 /2E 0 is the vacuum oscillation frequency, and η = ±1 is the mass hierarchy. The diagonal block L B has real eigenvalues, while the eigenvalues of a non-hermitian 4 × 4 matrix L A may be complex if the initial neutrinos possess a flavor imbalance, s νe,x = sν x,e . Further we will plot the instability rates κ max , i.e., the maximum imaginary parts of the eigenvalues of L A , for s νe = 1, sν e = α ≥ 0, and s νx = sν x = 0, changing the neutrino intensity parameter µ and the NSSI coupling g + . In fact, for the matter potential V e,n ≡ G F √ 2n e,n set to zero, the eigenvalues of L A can be found analytically, From this equation, one observes that for a strong NSSI coupling |g + | > 3, instabilities survive even in the ultradense neutrino gas regime ω = V e,n = 0, their growth rate being κ max = g 2 + − 9 µ|1 − α|. This corresponds to a so-called fast unstable mode (see, e.g., Refs. [14,15]), whose growth rate for ω, V e,n ≪ µ is of the order of µ. It is well known that there are no such modes in the single-angle scheme with purely electroweak (V-A) neutrino-neutrino interactions: all instabilities disappear as ω → 0 [9]. In the case of a scalar-pseudoscalar NSSI, as we see, such a mode appears even for monochromatic neutrinos and its growth rate is plotted in Fig. 1. However, we will not discuss it here, since we are more interested in unstable modes present for a small NSSI coupling. These come from an eigenvalue branch with lim ω→0 Im λ(ω) = 0, thus, these instabilities are slow, i.e., suppressed in the ω ≪ µ regime.
The instability rates for the slow modes are plotted in Fig. 2 for different neutrino/antineutrino ratios α and both hierarchies. One observes that presence of antineutrinos considerably enhances the instability region; for normal hierarchy (η = +1), in fact, the widest instability 'sector' in the (µ/ω, g + ) plane corresponds to α = 1. This sector touches the g + = 0 line at a resonant neutrino density µ = ω/3|1 − α|, i.e., for such a density, the studied type of instabilities arises for arbitrarily small g + . The pattern for inverted mass hierarchy is quite different, also presenting another, high-neutrino-density instability region. In both cases, as one observes, these slow instabilities should grow at typical scales smaller than the vacuum oscillation length, i.e., at several kilometers or even smaller.
Let us now discuss the effect of background matter, i.e., nonzero V e,n potentials, on the stability properties. Interestingly, in the presence of neutrino-antineutrino NSSI coupling, a transformation to a corotating frame in the flavor space does not let one eliminate the MSW term, as it does for g ± = 0 [9]. Namely, a substitution into the equation of motion (83) for the instability does not totally 'hide' the MSW term into the unitary transformation because of the matrix structure of the block-off-diagonal NSSI interaction (84). Note that both the electroweak (g ± = 0) interaction Hamiltonian and the Hamiltonian (72) from Ref. [6] feature a combination of the form ρ − ρ cT instead of a scalar-pseudoscalar ρ T − ρ c in Eq. (84), which lets one eliminate the MSW term by a unitary transformation to a corotating frame. Physically, this results in the fact that in the presence of NSSI, background matter modifies the instability rates; in particular, one observes that even when ω = 0 and V e,n = 0, there are 'intermediate' (neither fast, nor slow) unstable modes whose rates κ are proportional to V e,n . Indeed, in a typical supernova situation, µ V e,n ≫ ω near the neutrino sphere, thus, these should be considerably faster than the slow ones. Correspondingly, the two possibly non-real eigenvalues in the ω = 0 case read In particular, from this expression, it follows that for a small NSSI coupling, instabilities take place in a narrow resonance region around µ = (V n − V e )/3(1 − α). In Fig. 3, we demonstrate their growth rates for an idealized situation ω = V e = 0; the pattern looks quite similar to Fig. 2, but now the rates are measured in V n instead of ω, i.e., they are higher.
As we see, the block-off-diagonal NSSI coming from the scalar and pseudoscalar terms in the Lagrangian may lead to fast, slow, and intermediate instabilities that remain hidden in the absence of neutrino magnetic moment (or external magnetic field). Within the purely electroweak interaction model, these modes are also absent. limiting the status of its conclusions. Moreover, in reality, one is interested in non-negligible NSSI/magnetic moment signatures in the observable neutrino energy spectra, which is beyond the regime of linear perturbations on top of a stationary nonperturbed solution. To explore the issue and estimate the sensitivities of the neutrino spectra to the magnetic moment, in the present section we carry out a numerical analysis of collective oscillations with the block-off-diagonal NSSI within the single-angle scheme.
For our simulation, we work with two neutrino flavors with ∆m 2 = 2.4 × 10 −3 eV 2 and θ = 9 • [34]. The supernova setup is analogous to the one used, e.g., in [6]. Namely, a neutrino is leaving the neutrino sphere with the radius R ν = 50 km radially in the equatorial plane, so that the transversal component of the magnetic field decays as an inverse power law with r, The neutrino density profile and the geometric factor together form an effective neutrino-neutrino coupling where L is the neutrino luminosity parameter, namely, the number of emitted neutrinos per second (we follow the convention of Ref. [9] dividing it by 2πR 2 ν ). By default we take L = 10 55 sec −1 corresponding to the total radiation power L × 10 MeV ∼ 1.5 × 10 50 erg/sec. The electron density is chosen following the single-angle simulations in Ref. [6], while we also add neutrons with the density n n = 1.5n e . The potentials V e,n (r) = G F √ 2n e,n (r) and the above neutrino-neutrino coupling are shown in Fig. 4a, together with the energy scale ω = ∆m 2 /(2 × 10MeV) of vacuum oscillations. At the neutrino sphere, the flavor density matrix is diagonal (see Eq. (70)), with the energy spectra taken from a simulation [39] (see Fig. 4b); these spectra are nothing but Fermi distributions for the three different decoupling temperatures of ν e ,ν e , and ν x /ν x . Upon evolution (65), the flavor/energy spectra are extracted from the density matrix with the help of a flavor/helicity projector P f We study the region R ν ≤ r 250 km far from the MSW resonance (V e,n ≫ ∆m 2 cos 2θ/2E), in which the oscillations are mainly driven by collective effects. Moreover, in reality, the neutrino self-coupling falls off quite rapidly, so that the oscillations virtually freeze around r = 200 − 250 km producing well-defined 'final' neutrino spectra on exit from the region. It is these spectra produced by nonlinear, collective effects that we are interested in within our analysis of NSSI-induced flavor instabilities. Let us first discuss the flavor evolution without (pseudo)scalar NSSIs, i.e., within the (minimally extended) Standard Model, comparing the effect of the self-interaction Hamiltonian (68) derived by us and that of the Hamiltonian (72) claimed in Ref. [6]. As mentioned in Sec. II, these two Hamiltonians lead to identical flavor evolutions in the µ 12 = 0 case. Fig. 5 shows the results of the simulation for µ 12 = 0, in the cases where the effect is visually noticeable. It turns out that with our Hamiltonian the final spectra are virtually insensitive to the neutrino magnetic moment up to at least |µ 12 | ∼ 10 −15 µ B , while the evolution with self-interaction (72) contains considerable magnetic moment signatures already for µ 12 = 10 −19 µ B , especially in the normal hierarchy (the latter was, in fact, stated in [6]). For a luminosity L = 10 54 sec −1 one order of magnitude lower than that in Fig. 5, the magnetic moment signatures virtually disappear, obviously because the instabilities causing them get suppressed. Note that the low sensitivity of the evolution equations (65) to the magnetic moment of a Majorana neutrino agrees with the linear stability analysis in Sec. III, so that we have to conclude that within the Standard Model without NSSIs, the effect of a small magnetic moment should be quite hard to observe. This has also been noted in our previous analysis [25]; the case of a large AMM has been recently studied in [26], including the angular neutrino distributions, and the analysis demonstrates that for a smaller AMM/weaker magnetic field, their signatures are suppressed. Regarding the Hamiltonian (72) that was claimed to hold within the Standard Model in Ref. [6], we have to consider it a nonstandard neutrino self-interaction instead, which leads to a pronounced effect of the neutrino transition magnetic moment. The upper and the lower rows correspond to the normal and the inverted hierarchies, respectively. The first column (a, d) shows the no-AMM case µ12 = 0, the second (b, e) the µ12 = 0 case with self-interaction Hamiltonian (68), and the last one also to µ12 = 0, but based on the Hamiltonian (72) from Ref. [6]. In the µ12 = 0 case, the two Hamiltonians produce the same result. Note that virtual identity of (a,c) and (b,d) plots, respectively, holds up to at least |µ12| ∼ 10 −15 µB Now let us switch to the main issue of the present paper, namely, to the effect of nontrivial scalar/pseudoscalar NSSIs on the collective neutrino flavor evolution. Namely, for the same luminosity L = 10 55 sec −1 , we set the block-diagonal NSSI coupling g − to zero (as mentioned earlier, this coupling has been analyzed in other papers [33]), keeping only the block-off-diagonal term with the g + coupling, and analyze the effect of this small NSSI for a small transition magnetic moment µ 12 = 10 −19 µ B . Again, when the magnetic moment is zero, the density matrix becomes block-diagonal and the g + coupling produces zero effect on the oscillations. For a nonzero magnetic moment, in contrast, nonstandard interactions radically change the neutrino spectra, see Fig. 6.
Indeed, neutrino-antineutrino instability caused by the block-off-diagonal NSSI and triggered by the µB interaction shows up in both hierarchies provided that the NSSI coupling is not too small, |g + | 0.3, and rapidly leads to a complicated, if not a chaotic pattern. This is clearly manifested in the wiggly patterns in Fig. 6 resulting from probability exchange between ν e andν x flavors. For large g + couplings, this exchange even results in a ν e −ν x spectral split, which replaces the ν e − ν x split in the inverted hierarchy case. In fact, animations of the flavor evolution we made also reveal that g + = 0 flavor evolution starts from rapid synchronized oscillations of the spectra with typical lengths of 10 − 100 m (in contrast to the no-NSSI case, where the oscillation lengths are of the order of 10 km), which, in principle, agrees with our conclusions on the fast/intermediate instabilities in Sec. III.
The effects do not drastically depend on the magnetic moment, since it acts just as a seed for an exponentially growing instability; deformation of the spectra can be observed even for µ 12 ∼ 10 −24 µ B (see Fig. 7 hand, the rates of instabilities resulting in nontrivial NSSI signatures strongly depend on the degree of nonlinearity of the equations, i.e., on the luminosity, and these signatures are virtually absent for L = 10 54 sec −1 . Let us now take a closer look at the development of NSSI-induced instabilities for L = 10 55 sec −1 and address the issue of their potential observability. First of all, a natural parameter controlling these neutrino-antineutrino instabilities is the ratio of the total neutrino and antineutrino numbers n ν (r)/nν (r) obtained after integration over the whole energy spectrum. This ratio is conserved in the absence of magnetic moment µ 12 ; however, as we saw above, it is also conserved to a high accuracy even when µ 12 is nonzero, in both noncollective oscillations and collective oscillations without NSSIs, neither of which contain neutrino-antineutrino instabilities. Interestingly, Fig. 8a demonstrates that when NSSIs come into play, the total neutrino and antineutrino numbers rapidly reach an equilibrium plateau. This phenomenon resembles a sort of equilibration discussed in Ref. [26], but for a value of µ 12 B that is many orders of magnitude smaller. The equilibrium neutrino-antineutrino ratio is then conserved up to the neutrino detector, so that antineutrino excess could, in principle, serve as a signature of NSSIs mixing neutrinos with antineutrinos, such as scalar and pseudoscalar interactions. Fig. 8a also reveals a fast character of the instability, as opposed to slow instabilities with growth scales of the order of 10 km. To quantify this instability and the resulting impact of the neutrino spectra, we introduce a spectral residual as the relative integral deviation of the neutrino flavor spectra from the no-AMM case, in which the effect of NSSIs is absent where n (0) f (E; r) correspond to the flavor-energy spectra for µ 12 = 0 and the summation is performed over both −i[h AMM , ρ] to it does not round the latter one down to zero: the 'large' blocks are not added to 'small' blocks of ρ. During the later stages, it is mainly the self-interaction that drives the evolution of the density matrix. neutrino and antineutrino flavors. Now, Fig. 8b clearly demonstrates a drastic difference between the evolutions of neutrino spectra for g + ≤ 0.3 and g + 0.5: in the latter case, the residuals ∆(r) rapidly saturate to quite observable values of the order of 10 −1 and keep this level during further evolution. Note that the level ∆ ∼ 0.01−0.1 corresponds, in fact, to quite a large spectral distortion, such as those shown in Fig. 6. Finally, the spectral residuals at the 'final' point r = 250 km are shown in Fig. 9 versus the luminosity L and the NSSI coupling g + . One observes that at least for L ≤ 10 56 sec −1 , one can probe the NSSI coupling with the sensitivity around 0.25 − 0.4 in both mass hierarchies. It is important, however, that this order sensitivity can be achieved for extremely small values of the neutrino transition magnetic moments, comparable with the figures predicted within the minimally extended Standard Model with the standard, V-A structure of weak interactions.

V. DISCUSSION
The analysis given in the above sections has led us to a number of interesting conclusions regarding the effect of nonstandard, scalar and/or pseudoscalar four-fermion neutrino interactions on the collective oscillations taking place during a supernova explosion. Firstly, it turned out that such NSSIs include the terms that mix neutrinos and antineutrinos, and these are able to open up a new channel of fast neutrino-antineutrino instabilities considerably deforming the neutrino flavor/energy spectra. In fact, these spectral features, including nonstandard splits and chaoticlooking patterns, need a seed to develop; however, a minuscule transition neutrino magnetic moment µ 12 ∼ 10 −24 µ B proves to be enough to trigger the effect via neutrino-antineutrino (helicity flip) transitions in the stellar magnetic field. Note that such a transition magnetic moment could be generated even at the one-loop level of the Standard Model, in which a Glashow-Iliopoulos-Maiani (GIM) cancelation takes place [22,40]. Secondly, the discussed type of ν −ν instability is also quite nonstandard in that the electron/neutron background nontrivially affects the instability rates even for monochromatic neutrinos (see Sec. III devoted to linear stability analysis). From this analysis we also observed that, quite as usual, a neutrino-antineutrino ratio close to unity favors the development of instabilities.
Thirdly, the sensitivity of (anti)neutrino spectra to the NSSI coupling stays at the level of g + ∼ 0.3 for neutrino luminosities L 10 55 sec −1 , notably, virtually regardless of the (nonzero) value of the magnetic moment. Note that the spectral residual ∆(250 km) chosen by us to quantify the NSSI signatures in the neutrino spectra (see Eq. (94)) can change outside the r = 250 km sphere, where collective oscillations are not important, whereas other effects, such as the MSW effect in a turbulent medium, may come into play [41]. Importantly, however, there is another natural measure of scalar/pseudoscalar NSSI signatures in question, that does not change under (conventional) noncollective oscillations: the neutrino-antineutrino ratio n ν /nν. Both this ratio for neutrinos of a given energy E and the one integrated over the whole energy spectrum are almost conserved in noncollective oscillations (the magnetic moment leaves a negligible effect once the collective oscillations are over), thus, the neutrino-antineutrino ratio (see Fig. 8) should bring the information on the NSSI in the innermost supernova layers to the neutrino detector. Moreover, even though in the present paper we have discussed the NSSI effect on collective oscillations within a simplified picture not including angular degrees of freedom, it is highly likely that the neutrino-antineutrino ratio should depend on the 'latitude', i.e., on the angle between the directions to the stellar magnetic pole and the observer. Therefore, helicityflipping scalar/pseudoscalar NSSIs should manifest themselves in a periodic variation of the neutrino-antineutrino flux ratio from an explosion of a (rotating) supernova.
Finally, from our analysis, both numerical and analytical, it inevitably turns out that within the framework with a pure V-A interaction, the effective neutrino oscillation Hamiltonian takes the form derived in Ref. [23], which is unable to exponentially catalyze the growth of coherent neutrino-antineutrino mixing initially introduced via the transition magnetic moment. Note that the linear stability analysis in Sec. III that led to this conclusion was not a Lyapunov stability analysis of a necessarily stationary solution, so it does not rely upon a 'qualitatively right' assumption of vanishing vacuum mixing angle. As a result, within the Standard Model extended with nonzero Majorana neutrino masses, the effect of a small nonzero magnetic moment is very unlikely to be observable. A couple of words should also be said on what was left beyond the present analysis. First of all, Dirac neutrinos should probably have very similar block-off-diagonal terms in the interaction Hamiltonian, if one introduces scalar and/or pseudoscalar terms into the Lagrangian, while in this case these blocks will describe superpositions of active, left-handed neutrinos with 'sterile', right-handed neutrinos rather than neutrino-antineutrino mixing. In the absence of NSSIs, these blocks vanish, so nothing is able to catalyze the helicity-mixing effect introduced by a small magnetic moment, analogously to the Majorana case [25]. Next, we have tested Majorana neutrinos with scalar/pseudoscalar interactions for instabilities initiated by the so-called induced magnetic moment [42], which physically is a helicitydependent Wolfenstein potential in a magnetized medium, and found no new instabilities. Among other possible reasons, this is definitely related to the fact that at least for a nonmoving background medium, the induced magnetic moment term in the Hamiltonian does not flip the neutrino helicity. However, this, as well as other issues regarding the effect of scalar and pseudoscalar nonstandard neutrino self-interactions on the evolution of dense neutrino flows require further study beyond the present paper.