Spatiotemporal linear instability analysis of collective neutrino flavor conversion in 4-dimensional spacetime

In an environment with high-density neutrinos formed in a core-collapse supernova (CCSN), the neutrinos exhibit nonlinear and complex oscillation behaviors due to their self-interactions. The onset of this nonlinear oscillation can be investigated by linearizing the evolution equation for small perturbations around the flavor eigenstates. While the condition under which the flavor eigenstates are unstable has been investigated in many studies, how the perturbations evolve in spacetime has yet to be elucidated. In this paper, we analytically and correctly derive the asymptotic behaviors of the linear perturbations in 4-dimensional spacetime in the linear regime for a 2-beam neutrino model using the recently proposed Lefschetz thimble formulation. The result suggests that the perturbations grow in the directions between the two neutrino beams. We also briefly discuss the possible effects of neutrino flavor conversion on the explosion mechanism of a CCSN. In particular, the result implies that the flavor instability in the preshock region may propagate into the postshock region, contrary to the previous study focusing on the group velocity in 1-dimensional space. How to treat the case of a more realistic continuous spectrum is also discussed.


I. INTRODUCTION
Neutrino oscillation is a well-known phenomenon in which the survival probability of each neutrino flavor oscillates due to the deviations between the flavor eigenstates and mass eigenstates of neutrinos. When background matter is considered, forward scattering on charged leptons changes the effective mass of the neutrinos and, hence, their oscillation behaviors. This is called the Mikheyev-Smirnov-Wolfenstein (MSW) effect [1,2], and its behavior is well understood. Additionally, neutrino oscillation is similarly affected by neutrino selfinteractions, which, unlike the MSW effect, transform neutrino oscillation into a nonlinear phenomenon. This phenomenon is called collective neutrino oscillation, and its complex behaviors have been intensively studied by many researchers.
Collective neutrino oscillations are important only in environments with large amounts of neutrinos, such as supernovae and the early universe. In particular, flavor conversion may affect the explosion mechanism of a core-collapse supernova (CCSN) [3][4][5][6][7][8][9][10]. It is widely believed that the stagnant shock in a supernova core needs to be revived through the introduction of additional energy for the explosion to be successful. Nucleons in the postshock region can obtain energy by absorbing ν e and ν e emitted from neutrinospheres, while the absorption of the other flavors of neutrinos is almost kinematically prohibited. Therefore, flavor conversion can change the effective heating rate of the shock and may play a key role in determining the success of the explosion.
Although numerical approaches to solving the nonlinear equations for collective neutrino oscillation suffer from enormous computational costs, attempts have been made to address these equations under certain assumptions in many studies. For example, the time evolution of homogeneous monochromatic single-angle neutrinos can be understood as analogous to the motion of a pendulum and a spinning top and exhibits several interesting phenomena, such as synchronized and bipolar oscillations [11,12]. Numerical calculations treating many modes have also been performed and have revealed some new phenomena, such as spectral swaps or splitting . It should be stressed that in all of these studies, the nonlinear differential equations have been solved by imposing certain symmetries to reduce their dimensionality in spacetime.
Linear stability analysis has also been applied to investigate the conditions for the onset of flavor conversion [35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53]. Such analysis involves considering small perturbations around the flavor eigenstates and searching for growing modes from the derived dispersion relations (DRs). Flavor conversion can occur if there is a mode with a real wave vector k and a complex angular frequency ω. Many studies have predicted that flavor conversion occurs when the angular distribution of the electron lepton number (ELN) crosses 0, although this has not been proven mathematically. Moreover, such linear analysis suggests the importance of multidimensionality. Generally, the flavor eigenstates are more unstable in higher dimensions, and the results of flavor oscillation could drastically change in reduced dimensions. Therefore, it is important to consider 4-dimensional spacetime to understand realistic flavor conversion behaviors.
Normal mode analysis can reveal the conditions for the onset of flavor conversion and which plane wave modes grow or decay. In reality, however, perturbations may take the form of wave packets, and how they evolve in spacetime cannot be trivially determined from the DRs. Recently, it has been suggested that the spatiotemporal behaviors of perturbations can be investigated via the method originally proposed by Briggs [54][55][56]. One can derive how perturbations grow from the coalescence features of the analytic continuation of the DR k(ω). Insta-bilities can be classified as either absolute or convective instabilities on the basis of the behaviors of perturbations in the form of wave packets. This analysis, however, assumes that a perturbation is a wave packet in only one direction and homogeneous in the other directions, which is not realistic in supernovae. Recently, we developed a general and powerful method to investigate the spatiotemporal evolution of linear perturbations by using the Lefschetz thimble formulation [57]. This formulation can be used to treat multidimensional perturbations for arbitrary DRs. In this paper, we apply this method to the collective neutrino flavor conversion in a 2-beam neutrino model and reveal the flavor conversion behaviors in 4-dimensional spacetime. This is the first study to treat the spatiotemporal evolution of collective neutrino oscillations in 4-dimensional spacetime, although this is done in the linear regime.
The results for the 2-beam model show an absence of absolute instabilities in 4-dimensional spacetime, although they are present when only 2-dimensional perturbations are considered. Moreover, the results suggest that flavor instabilities grow toward the directions between two neutrino beams and may imply the impact of flavor conversion on shock heating mechanisms in CC-SNe.
This paper is organized as follows. In Sect. II, we linearize the kinetic equations that describe neutrino oscillations to treat perturbations around the flavor eigenstates. In addition, we introduce the 2-beam neutrino model and reduce the linearized equations. In Sect. III, we perform a spatiotemporal analysis of the obtained linearized equations for the 2-beam model and discuss general cases other than the 2-beam model. Sect. IV concludes the paper.

A. Kinetic equations for neutrinos
We begin with the kinetic equation for the neutrino density matrix f [45,[58][59][60][61]: which describes the evolution of streaming neutrinos in a potential generated by matter and neutrinos. (x µ ) = (t, x) represents the coordinates in spacetime, and Γ = (E, v) denotes the energy and flight direction of the neutrinos. We express antineutrinos by means of a density matrix with negative energy, as follows: The Hamiltonian H is given by where H vac corresponds to vacuum oscillation and H int to the potential. The vacuum oscillation term is where M 2 is the neutrino mass-squared matrix. The potential term is written as where j µ α is a lepton number 4-current of charged lepton α and The collision term C [f] thermalizes the density matrix f, or causes the decoherence of neutrino oscillations in general. The effects of this term on collective neutrino oscillations have recently been discussed [34,62,63]. In this study, we neglect these effects because they are usually much smaller than the flavor instability induced by neutrino potentials in supernovae.

B. Linearization
Here, we consider 2-flavor neutrinos ν e and ν x , because each pair of flavors can be decoupled from the others even if we consider the number of flavors to be more than 2 [53]. Then, the mass-squared matrix can be expressed as with the neutrino mixing matrix where θ is the mixing angle. Additionally, the density matrix f is a 2 × 2 Hermitian matrix and can be decomposed by means of (τ a ) ≡ (I 2 /2, σ/2) with the Pauli matrices σ as follows: The commutator of the Hermitian matrices A = A a τ a and B = B a τ a can be expressed as where c ab is the Levi-Civita symbol if a, b and c are permutations of (1 2 3) and 0 otherwise. Therefore, Eq. (1) can be recast as where the decomposed Hamiltonian components are and with ∆m 2 ≡ m 2 1 − m 2 2 . To investigate the onset of flavor conversion, we focus on the deviation from a flavor eigenstate and express the density matrix as where f 0 = f ee + f xx is the incoherent part, which does not play a role in flavor mixing, and f c = f ee − f xx is the coherent part. Substitution of this expression into Eq. where and s is the term of zeroth order in ε, expressed as This term appears because the flavor eigenstate is not a fixed point of Eq. (1). In the linear analysis we consider, we address the time evolution of ε from the flavor eigenstate ε(0, x, Γ) = 0. Additionally, the spatial domain can be taken to be the open space R 3 , and the boundary condition for ε is given by Eq. (16) with f c = 0 and Λ c = 0 at |x| → ∞. Under these initial and boundary conditions, the absence of s yields ε(x, Γ) = 0 at all times; hence, s is regarded as the seed perturbation.
Here, we define to recast Eq. (16) as up to linear order in ε. The right-hand sides (r.h.s.) of Eqs. (23)- (25) are the source of the linear evolution, and the DR of this system of linear equations depends only on the left-hand sides (l.h.s.). We can confirm that ε 0 is included only in Eq. (23) and is decoupled from the other components of ε. Additionally, Eq. (23) gives the trivial DR for free-streaming massless particles: which has no instability. Therefore, our main target is the DR for Eqs. (24) and (25). If we neglect flavor mixing, then from Eq. (1), we have c = v · ∂ ln f c ∼ C /f c , which is of the same order of magnitude as the inverse of the mean free path of the neutrinos. This quantity is usually much smaller than the growth rate of neutrino flavor instabilities and can be safely neglected when computing the DR. The values of ω c/s are also minimal compared to the growth rate for the typical energy of supernova neutrinos when we focus on a sufficiently small radius in a supernova. As the radius increases and the neutrino density decreases, however, ω c/s becomes comparable to the growth rate, and a gradual transition to vacuum oscillation occurs.
The flavor instability when the ω c/s values are negligible is called a fast instability; when an instability with finite ω c/s vanishes in the limit of ω c/s → 0, it is called a slow instability [48]. We note that Eqs. (24) and (25) are coupled with each other via the terms proportional to ω s . In general, it is not appropriate to neglect ω s when focusing on slow instabilities, although Airen et al. [48] did so.
Here, we focus on fast instabilities; we neglect ω c/s as well as c on the l.h.s. of Eqs. (24) and (25): These equations are no longer coupled with each other, and the DR for ε 3 is also given by Eq. (26). By taking the average over the energy E, Eq. (27) is rewritten as where is the ELN angular distribution and Here, we omit a concrete expression for the source term s and considers to depend on the spacetime position x and the flight direction v, although the r.h.s. of Eq. (27) does not, for the reason discussed below. Equation (21) can be separated into two contributions: (1) the homogeneous part (0, 0, −ω s (E), 0) T and (2) the inhomogeneous part, which is the remainder. The homogeneous part arises from the vacuum oscillation term and induces only modes with the wave vector k = 0. On the other hand, the inhomogeneous part depends on the position and can induce all modes. This study focuses on the latter, and we consider the spatiotemporal behaviors of the linear response to this inhomogeneous source term.
We need to be careful when treating a fast instability. When we omit ω s on the l.h.s. of Eqs. (24) and (25), these equations are decoupled from each other, and the perturbation S seems to be induced only by the homogeneous source −ω s , as in Eq. (27). Even a small ω s , however, can convert the inhomogeneous source −c in Eq. (25) into S, and it then grows exponentially when a flavor eigenstate is unstable. In other words, the effective source term proportional to ω s c should also be contained in Eq. (27) in the fast regime to mimic the contribution from the inhomogeneous source term, and therefore,s should be considered to have an (x, v) dependency.

C. Two-beam model
To investigate the instabilities in detail, we adopt a 2-beam neutrino model in which the ELN angular distribution is expressed as where v 1 and v 2 are the directions of the neutrino beams and G 1 and G 2 correspond to their ELN intensities. The delta function δ(v, v ) on the unit sphere is defined as for an arbitrary continuous function f on the unit sphere.
In addition, we focus on a sufficiently short time and a sufficiently small region such that the temporal and spatial variations in Λ c (x) can be treated as constant. In this model, Eq. (29) can be written as where D andS are defined as andS We note that Λ c is a real 4-vector and does not affect the instability. For instability analysis, the DR given by the zeros of

III. SPATIOTEMPORAL INSTABILITY ANALYSIS
A. General formulation The spatiotemporal behaviors of the perturbationS are obtained by solving Eq. (34) for the appropriate initial and boundary conditions. In our setup, the inhomogeneous sources, which takes the form of a wave packet in realistic situations, drives the perturbation. To treat such perturbations, we simply focus on the Green's function G for the linearized equation defined by because its asymptotic behavior at t → ∞ is essentially the same as that of the perturbations [54,55].
Here, we discuss the necessity of this prescription. Izaguirre et al. [45] stated that a complex wave vector k for a real ω in a DR implies "spatial instability", and several studies based on linear analysis for stationary solutions have actually investigated this kind of "instability". We should, however, be careful about what this really means [49]. The existence of a complex k for a real ω ensures only the existence of spatially exponentially growing/decaying modes that oscillate or are constant in time. Therefore, such a "spatial instability" seems to appear if we consider stationary solutions of Eq. (1). However, it is not guaranteed that stationary solutions, and hence this "spatial instability", will be realized after time evolution. To ensure this, we should at least investigate the time evolution of the perturbations from the stationary solution and confirm the damping of the perturbations. Moreover, even spatial growth is not guaranteed by this "spatial instability". For example, if we consider the diffusion equation (∂ t − ∂ 2 x )f (t, x) = 0, which has a "spatial instability" from the DR k = ± √ iω, a constant source imposed at a certain point decays in space. After all, such a "spatial instability" is not always relevant in realistic situations. On the other hand, a complex ω for a real k does imply instability in time evolution. The explicit form of the DRs ω(k) for k ∈ R 3 directly shows the time evolution of each plane wave mode. However, how perturbations in the form of wave packets, which are superpositions of uncountably infinite numbers of plane waves, actually behave in spacetime cannot be trivially determined from the DRs; hence, we should focus on the behaviors of the Green's function, which is also a superposition of all modes.
The asymptotic behavior of G can be evaluated by using the Lefschetz thimble method [57]. By applying the Laplace transform for time and the Fourier transform for space, G is expressed as where M ≡ L × R 3 is the Laplace-Fourier contour and (u µ ) ≡ (1, u) with the parameter vector u ∈ R 3 . The residue theorem yields where C ≡ D ∩ C × R d and D ≡ k ∈ C d+1 |∆(k) = 0 . The contour of integral C can be deformed to a sum of Lefschetz thimbles {J σ } as follows: where C, K σ is the intersection number of C and K σ and the {K σ } are called the dual thimbles. The Lefschetz thimble J σ and the dual thimble K σ are associated with the critical points of the "height function" defined on D as follows: which corresponds to the real part of the exponent of the dominant factor e −ik·u . The critical points are given by the stationary conditions for h constrained on D, which can be simplified as The Lefschetz thimble associated with the critical point k σ is obtained by solving the upward flow equation with the boundary condition lim s→∞ K(s) = k σ , and the dual thimble K σ is found by solving the same differential equation with the boundary condition lim s→−∞ K(s) = k σ . When |k σ ·x| |k σ ·ut|, the integral over the Lefschetz thimble J σ is dominated by the contribution around the critical point k σ , so the integral is evaluated as follows: where k m is the critical point with the maximum height h(k m ) among all critical points satisfying C, K σ = 0. Now, we can classify instabilities in terms of the behavior of G. If G grows for u = 0, this means that perturbations will grow at every point throughout the whole space; such an instability is called an absolute instability. On the other hand, if G decays for u = 0 but grows for another u, then the perturbations will grow but flow away from the generated point; such an instability is called a convective instability. If G decays for all u, then the system is stable.

B. Two-dimensional perturbation
Now, we can conduct a spatiotemporal instability analysis of the 2-beam model considering 4-dimensional perturbations. In this subsection, however, we first perform a 2-dimensional analysis for two reasons: we would like to reveal the problems with the 2-dimensional analysis, and in the case of the 2-beam model, the results of the 2-dimensional analysis are also helpful for the 4dimensional analysis.
For later use, we introduce some vector notations as follows: • Vector component parallel to b: • Projection of a vector onto the direction parallel to b: • Projection of a vector onto the plane perpendicular to b: In addition, we define which are perpendicular to each other because |v 1 | = |v 2 | = 1. The 2-dimensional Green's function G n can be defined as where the unit vector n gives the direction along which we consider perturbations and x n ≡ x · n as defined above. When we take n to be the z-direction, the situation is the same as in Ref. [56]. The Laplace-Fourier transform for time and the n-direction yields and G (2) n can be expressed as Although we should compute the asymptotic behaviors of G (2) n (t, x + ut) for all u using the method given above, we can achieve an equivalent task by considering only u = 0 in this case. When the variables of integration are transformed such that Eq. (54) becomes This integral can be evaluated by means of the Lefschetz thimble method with height function h(w, κ) ≡ Im w and DR∆ where The critical points for this system are given by This system of equations can be solved analytically, and we obtain two critical points (w + , κ + ) and (w − , κ − ): The dual thimbles K ± can be obtained by solving the upward flow equation with the boundary condition lim s→−∞ K(s) = k σ . Numerically, this is achieved by modifying the boundary condition as follows: K(0) = k ± + (0, H ± δ), where δ ∈ R is a small quantity and The critical points and dual thimbles are shown in Figs. 1, 2, 3 and 4 for several values of and α. All of these figures show the projections onto the Im κ − Im ω plane, and the original contour C lies within the profile of Im κ = 0 (cyan line). What we need to pay attention to is the intersection numbers of C and K ± . For example, in Fig. 1, each of the dual thimbles K ± has a single intersection with C, i.e., | C, K ± | = 1, which means that the Lefschetz thimbles J ± both contribute to the integral. We note that the sign of the intersection number C, K σ can be determined by giving K σ , and consequently J σ , an orientation. For our purposes, however, factors other than the dominant exponential factor are not important, so we do not consider them. In contrast, the dual thimbles in Fig. 4 seem to intersect with C twice. In fact, however, they cross C from opposite sides at each intersection, and the intersection numbers for the two intersections are +1 and −1; consequently, the sum is C, K ± = 0.
The situation is more complicated in Figs. 2 and 3. The dual thimble K − is absorbed into the higher critical point (w + , κ + ) and cannot be well defined. This situation is understood to be related to the fact that these parameter sets lie on Stokes rays [64], and we can avoid it by adding small perturbations to the imaginary parts of the parameters (see the bottom panels in the figures). In Fig. 3, it can be confirmed in the bottom figure that both of the dual thimbles have intersections with C, i.e., | C, K ± | = 1. In Fig. 2, the dual thimble K + apparently does not intersect with C, and K − also does not, in the same manner as in Fig. 4; thus, neither of the Lefschetz thimbles influences the integral. It should be noted that the signs of C, K − can be either positive or negative, depending on the perturbations added. However, this is not important for the asymptotic behavior of the Green's functions because they are not the highest critical points to which the associated Lefschetz thimble contributes to the integral.
Although we show figures only for these 4 parameter sets, we can confirm that the intersection numbers depend only on the signs of and |α| − 1. Therefore, for |α| < 1, the Green's function G (2) n (t, ut) is exactly 0 for all t. For |α| > 1 and > 0, both of the Lefschetz thimbles J ± contribute to the integral, but the growth rate Im w ± = 0 means that the Green's function is damped as ∼ 1/ √ t. For |α| > 1 and < 0, the growth rate is Im w + = σ (2) , where As a result, the asymptotic behavior of G (2) n can be summarized as where M (2) , whose concrete expression is omitted, is a constant matrix independent of t. Based on these results, we can classify the instability for each parameter set as shown in Table I. For ≥ 0, the Green's function decays in time, meaning that the flavor eigenstate is stable. If < 0 and |V n | < |v n | are satisfied, the Green's function grows in the rest frame with u = 0, corresponding to an absolute instability. Finally, if < 0 is satisfied and |V n | < |v n | is not, then the Green's function decays for u = 0 but grows for u ∈ (V n − |v n |, V n + |v n |), corresponding to a convective instability.
The results are the same as those in Ref. [56] if we choose the z-direction as n. We stress, however, that in a 2-dimensional treatment, not only is the space of the ϵ = 1 α= 0.5

Classification Condition Stable
≥ 0 Convectively unstable < 0, |Vn| ≥ |vn| Absolutely unstable < 0, |Vn| < |vn| perturbations restricted, but the interpretation of the results in 4-dimensional spacetime is also ambiguous. The results depend on the direction n, and in fact, it can happen that the system is convectively unstable for one n but absolutely unstable for another n. Needless to say, the nature of the true instability should not depend on which direction we choose as n. What we have done here is to consider a perturbation that takes the form of a wave packet in the direction parallel to n but is homogeneous in the directions perpendicular to n and investigate the corresponding spatiotemporal evolution. Therefore, even if we were to investigate all possible directions n, perturbations with the form of wave packets in 4-dimensional spacetime would not be taken into account. How perturbations truly behave in 4-dimensional spacetime cannot be captured by a 2-dimensional analysis. Nevertheless, a 2-dimensional analysis may be justified for systems that are open in the n direction but bounded in the directions perpendicular to n. Even in such a case, however, we need to conduct the analysis for all plane wave modes perpendicular to n with wave numbers K ⊥n , and the DR will be modified to ∆(ω, k) = det D(ω, kn + K ⊥n ) = 0 (65) in general, although we treat only K ⊥n = 0 in the discussion above.

C. Four-dimensional perturbation
The general methodology for considering a 4dimensional perturbation has already been presented in Sect. III. A. In the 2-beam model, the integral in the calculation of the Green's function can be performed par-tially analytically, and the analysis reduces to the one in 2-dimensional spacetime.
Initially, D satisfies and the inverse matrix is expressed as The integral expression for G is given by Eq. (40). When we decompose the variables of the integral k into the components parallel and perpendicular to v and use Eq. (66), G is written as The variable transformation ω ≡ ω − V · k ⊥v yields Now, we note that the integral over k ⊥v can be performed, and those over ω and k v can be replaced by the 2-dimensional Green's function G (2) v : Therefore, the DR to be considered is again the one expressed in Eq. (57), but α is now given by We can obtain the asymptotic behavior of G v by substitutingv into n in Eq. (64), and the final result is where and M is a constant matrix, whose concrete expression is omitted. A schematic picture of the growth rate for each frame with velocity u is shown in Fig. 5. In the 4-dimensional case, a negative value of , which corresponds to the case in which crossings exist in the ELN angular distribution, is necessary for instability, as in the 2-dimensional case. The most important feature is that for < 0, the delta function in Eq. (72) constrains the velocity of the frame in which the system is absolutely unstable on u − V v. Thus, an absolute instability appears only for V = 0, and the instability is convective otherwise (see Table II). These results suggest that the instabilities grow toward the direction between the two neutrino beams. It may be interpreted intuitively as follows; a perturbation generated at (t, x) = (0, 0) can travel with the velocity v 1 and v 2 at each point; therefore, at time t, all the points the perturbation can reach are x = (sv 1 + (1 − s)v 2 )t with s ∈ [0, 1], which is consistent with the support of Eq. (72).  For CCSNe, it has been pointed out that ELN crossings may occur in decoupling regions [62]. Usually, ν e overwhelmsν e in almost all directions, andν e dominates only for almost radial directions in this region. This situation may be approximated by the 2-beam model with a radial antineutrino beam with G 1 < 0 and v 1 = e r and a neutrino beam in another direction with G 2 > 0 and v 2 = e r , where e r is the unit vector parallel to the radial direction. Therefore, flavor conversion occurs at least toward the radial direction, which may affect the heating of the stalled shock.
In addition, the possibility of ELN crossings in the preshock region was suggested in our previous study [52]. Because the average energy ofν e is larger than that of ν e ,ν e scattered on nuclei in the preshock region can overwhelms ν e while ν e usually doesν e in the radial direction. We concluded that the flavor instability does not propagate into the postshock region from the group velocity of the unstable modes in 1 spatial dimension. The results shown in this study, however, might overturn this conclusion. Radially-peaked ν e distribution seems to be imitated by a beam with G 1 > 0 and v 1 = e r in the 2-beam model while the scatteredν e distribution does by a radially-ingoing beam with G 2 < 0 and v 2 = −e r . As a result, from Eq. (72), perturbations can grow both outward and inward, which leads flavor conversions in the postshock region. In other words, scattered neutrinos seems to convey the instability into the postshock region as the intuitive interpretation mentioned above. Realistically, however, neutrinos do not have a discrete spectrum, as in the 2-beam model, but rather a continuous one. Further investigations should be conducted to determine how instabilities propagate in realistic supernovae.

D. Continuous distribution
Although we should consider a continuous spectrum in general, as mentioned above, it seems difficult to obtain the DR ∆(k) ≡ det D(k) = 0 because D has uncountably infinite dimensions, and it seems necessary to calculate the functional determinant of D. It should be noted that spectral discretization generates many spurious modes, which may affect the features of the instability, and thus, we cannot disregard the continuous nature of the spectrum [65,66]. Fortunately, we can prove that the DR is given by the determinant of a 4 × 4 matrix called the polarization tensor and that the instability depends only on this determinant.
The Green's function for the linearized equation in Eq. (29) is defined as By defining a µ v (k) ≡ , the Green's function G can be expressed as in a locally constant and uniform approximation of the potential Λ c (x). By substituting this expression into the definition of a µ v , we obtain the following linear equation for a µ v : where is called the polarization tensor [45]. By substituting this into Eq. (75), the following expression for the Green's function is derived: Therefore, all of the singularities of G originate from the roots of det Π(k − Λ c ) = 0 and v · (k − Λ c ) = 0 for arbitrary v. The latter does not influence the instability because ω is real for all k ∈ R 3 . Thus, it is sufficient to investigate the instability for the DR det Π(k) = 0.

IV. CONCLUSION
In this paper, we have conducted a spatiotemporal linear instability analysis of collective neutrino flavor conversion in 4-dimensional spacetime. We have analytically derived the asymptotic form of the linear perturbations for a 2-beam model and highlighted the importance of 4-dimensional perturbations for collective neutrino oscillations. In addition, we have shown that collective neutrino oscillations might influence shock heating not only in the decoupling region but also in the preshock region of a CCSN. We have also discussed the application of the method of instability analysis to continuous spectra, which are more realistic in nature.
As many previous studies have pointed out, ELN crossings allow neutrino flavor conversion, and the spatiotemporal behavior of such conversion in 4-dimensional spacetime is revealed in this paper for the first time. When we consider 2-dimensional perturbations, the instabilities can appear as both absolute and convective depending on the directions chosen for the perturbations, and such an analysis cannot elucidate how perturbations truly behave in 4-dimensional spacetime. When 4-dimensional perturbations are taken into account, it is clarified that the instabilities are actually convective and the perturbations grow toward the direction between the neutrino and antineutrino beams.
Needless to say, one of our goals is to gain an understanding of nonlinear spatiotemporal behaviors. However, it seems that more qualitative studies are also needed for guidance, especially in the current situation, in which numerical calculations seem to be extremely difficult. Additionally, more systematic studies are certainly needed to determine how flavor instabilities are actually triggered and evolve in supernovae. It is true that the 2beam model is simply a toy model, but we hope that this first-ever study on the spatiotemporal behaviors of flavor conversion in 4-dimensional spacetime will help develop an understanding of the more realistic flavor conversion behaviors in supernovae.