Collective fast neutrino flavor conversions in a 1D box: Initial condition and long-term evolution

We perform numerical simulations of fast collective neutrino flavor conversions in an one-dimensional box mimicking a system with the periodic boundary condition in one spatial direction and translation symmetry in the other two dimensions. We evolve the system over several thousands of the characteristic timescale (inverse of the interaction strength) with different initial $\bar\nu_e$ to $\nu_e$ number density ratios and different initial seed perturbations. We find that small scale structures are formed due to the interaction of the flavor waves. This results in nearly flavor depolarization in a certain neutrino phase space, when averaged over the entire box. Specifically, systems with initially equal number of $\nu_e$ and $\bar\nu_e$ can reach full flavor depolarization for the entire neutrino electron lepton number ($\nu$ELN) angular spectra. For systems with initially unequal $\nu_e$ and $\bar\nu_e$, flavor depolarization can only be reached in one side of the $\nu$ELN spectra, dictated by the net neutrino $e-x$ lepton number conservation. Quantitatively small differences depending on the initial perturbations are also found when different perturbation seeds are applied. Our numerical study here provides new insights for efforts aiming to include impact of fast flavor conversions in astrophysical simulations while calls for better analytical understanding accounting for the evolution of fast flavor conversions.

We perform numerical simulations of fast collective neutrino flavor conversions in an onedimensional box mimicking a system with the periodic boundary condition in one spatial direction and translation symmetry in the other two dimensions. We evolve the system over several thousands of the characteristic timescale (inverse of the interaction strength) with different initialνe to νe number density ratios and different initial seed perturbations. We find that small scale structures are formed due to the interaction of the flavor waves. This results in nearly flavor depolarization in a certain neutrino phase space, when averaged over the entire box. Specifically, systems with initially equal number of νe andνe can reach full flavor depolarization for the entire neutrino electron lepton number (νELN) angular spectra. For systems with initially unequal νe andνe, flavor depolarization can only be reached in one side of the νELN spectra, dictated by the net neutrino e − x lepton number conservation. Quantitatively small differences depending on the initial perturbations are also found when different perturbation seeds are applied. Our numerical study here provides new insights for efforts aiming to include impact of fast flavor conversions in astrophysical simulations while calls for better analytical understanding accounting for the evolution of fast flavor conversions.

I. INTRODUCTION
The discovery of the flavor oscillations of neutrinos by terrestrial experiments and astrophysical observations is one of the most exciting milestones in neutrino physics. Ongoing and planned experimental projects are expected to further pin down the yet-unknown parameters in neutrino mixing -the mass ordering and the CP -violating phase -while searching for potential signatures of physics beyond the Standard Model [1]. However, despite the success of the theory of neutrino mixing in explaining the majority of experimental data, one aspect that is poorly understood yet is how neutrinos oscillate in astrophysical environments dense in neutrinos. Analytical and numerical works over the past decades have shown that in an environment where the self-interactions between neutrinos cannot be ignored, various collective phenomena can arise due to the nonlinear and strong coupling nature (in the flavor space) of the system; see e.g., [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19] and review articles [20][21][22][23]. Exploratory works also demonstrated that such collective neutrino flavor oscillations may largely affect our understanding of important astrophysical events such as the core-collapse supernovae and the merger of two neutron stars [24][25][26][27][28][29][30][31][32].
In particular, Refs. [37] and [41] examined how fast flavor conversions develop in the nonlinear regime using numerical simulations in a one-dimensional (1D) box which possesses translation symmetry in the x and y directions while periodic in the z direction. Reference [37] showed that coherent and wavelike patterns in flavor space can develop with a point-source like perturbation. On the other hand, Ref. [41] found that the system can quickly settle into a state where flavor depolarization happens when averaged over the z domain, for a major part of neutrino spectrum with either point-source like or random perturbation seeds. These findings seem to contradict each other and require further examinations.
In this work, we systematically investigate the dependence of the outcome of fast neutrino flavor conversions in an 1D box on the initial condition of the perturba-tion seeds. By adopting advanced numerical schemes, we evolve the systems to late times when quasi steady states are achieved. We also explore the dependence of the steady-state outcome on the initialν e to ν e number density ratio.
In Sec. II, we describe our models, the initial and boundary conditions, and the adopted numerical schemes. We also briefly review the notation of the neutrino polarization vectors. In Sec. IV, we discuss our numerical simulation results and the implications. We summarize our main findings in Sec. V. Throughout the paper, we adopt natural units and takeh = c = 1.

A. Equation of motion
We consider a simple neutrino-dense system which has translation symmetry in the x-and y-directions in space, the same as in Refs. [37,40,41]. For simplicity, initially we assume a pure system consisting of only ν e and ν e (before setting small flavor perturbations; see below). Focusing on fast flavor conversions, we neglect the momentum changing collisions and only keep the neutrinoneutrino forward scattering contributions [55][56][57] in the effective Hamiltonian, i.e., we drop the terms originating from vacuum neutrino mixing and neutrino-matter forward scattering contributions.
Assuming the differential neutrino angular distribution dn νe(νe) /dΩ is uniform in z and taking the two-flavor approximation, the equation of motion, which governs the evolution of the normalized neutrino density matrices (see below) (for neutrinos) and¯ (for antineutrinos) is given as follows 1 : in the flavor basis. Since we omit the vacuum mixing contribution, the flavor evolution of ρ andρ do not depend on the neutrino energy.
1 When only considering the neutrino-neutrino self-interaction term, one can redefine the antineutrino density matrix byρ → −σ † yρ σy such that neutrinos and antineutrinos are treated in equal footing [2]. However, here we choose to treat them separately for the purpose of consistency with follow-up works that may include other terms in the Hamiltonian and the collisions.
The initial νELN angular distribution G(vz) for systems with differentνe to νe number density ratios α. The νELN crossing where G(vz) = 0 shifts to smaller vz with larger α.
In Eq. (1), the Hamiltonian H andH can be explicitly written down as: where µ = √ 2G F n νe with G F being the Fermi coupling constant. The asymmetry parameter α = nν e /n νe quantifies the ratio between the number density ofν e and ν e . The angular distribution function g ν(ν) (v z ) relates to the physical neutrino phase-space distribution function f νe(νe) without any flavor perturbations as where E ν is the energy of a neutrino. Note that the azimuthal symmetry in phase space of f νe(νe) with respect to the z direction is implicitly taken. In this notation, ρ ee =ρ ee = 1 for our system without any initial flavor perturbations while all other matrix elements are zero.
Since µ is the only dimensional quantity in Eq. (1), we define µ ≡ 1 and express t and z in dimensionless form (in the units of µ −1 ) hereafter.
For a system without the vacuum Hamiltonian, it requires a small perturbation in ρ ex (ρ ex ) to trigger the flavor instabilities. For simplicity, we assume that the perturbations are independent of v z . Specifically, we use the following description for our initial condition at t = 0 as: For (z), we explore the following two different choices: 1. Point-source like perturbations centered at z = 0: 2. Random perturbations: (z) is randomly assigned by a real value between 0 and 0 .

C. Numerical methods
We discretize z and v z with N z and N vz grids, and evolve and¯ defined at the grid centers. In evaluating the advection terms ∂ /∂z and ∂¯ /∂z, we use two different methods to check for consistency. The first method uses the forth-order (central) finite difference method with artificial dissipation using the thirdorder Kreiss-Oliger formulation [58]. For the second method, we use finite volume method plus the seventh order weighted essentially nonoscillatory (WENO) scheme [59,60]. For time evolution, we use fourth-order Runge-Kutta method. The numerical implementation details will be reported in a separate publication together with the public release of our simulation code COSEν [61].
Using the above initial and boundary conditions, we perform simulations with fiducial numerical parameters of N Z = 12000 and N vz = 200, with a fixed time step size ∆t = C CFL ∆z, where ∆z = L z /N Z and the Courant-Friedrichs-Lewy number C CFL = 0.4. In Appendixes A and B, we show the comparison of results obtained with different resolutions and with two different numerical methods for advection. For the rest of the paper, all results are obtained using the finite volume method with seventh order WENO scheme, which provides better numerical accuracy given the same resolution.

D. Polarization vectors
Before we discuss our results, let us define the neutrino and antineutrino polarization vectors P andP that are often used in literature. The three components of the polarization vector P are defined as which satisfy ρ ≡ (P · σ + P 0 I)/2, where σ i are the Pauli matrices and I is the identity matrix. For antineutrinos, we haveP With this definition, the equation of motion for P has the same form as forP, and we can treat neutrinos and antineutrinos in equal footing using the νELN spectrum G ν defined in Sec. II B to account for the contributions from both neutrinos and antineutrinos in the Hamiltonian. Dropping the bar for antineutrinos, we have: . Clearly, in the absence of collisions, the length of the polarization vectors are unity for all z and v z , given that our initial condition has |P| = 1 uniformly in z. Another physical quantity that is conserved globally is the net neutrino e − x lepton number in the entire simulation domain [40,41]: Note that locally this net lepton number can be transferred from one location to another due to the advection. However, globally L e−x must be a conserved quantity in time for the entire system.

III. RESULTS: FLAVOR EVOLUTION OF THE SYSTEM
In this section, we focus on results obtained with α = 0.9 and α = 1.1. In Sec. III A, we examine how the point-source like perturbations in the off diagonal elements of the density matrices grow in the linear regime and compare their evolution with detailed analysis using the dispersion relation. In Secs. III B and III C, we further discuss how flavor evolution proceeds in the nonlinear regime for cases with point-source like perturbations and random perturbations, respectively.

A. Linear regime
In the regime where |ρ ex | 1 and |ρ ex | 1, the evolution of the system can be qualitatively understood by performing the standard technique of linearized instability analysis [15,62]. By inserting ρ ex ∝ exp[−i(ωt−k z z)] (same forρ ex ) and keeping terms in the lowest order of perturbations, Eq. (1) leads to a dispersion relation of the collective mode of the system ω(k z ) [15,34,36]. The complex branch of ω for real k z signifies the existence of flavor instabilities such that small off diagonal perturbations can grow exponentially in time as the system evolves. Note that here we only include the solution that preserves the the axial symmetry. For the symmetrybreaking solutions, we omit them here because our simulation setup does not allow symmetry-breaking solutions.
In Fig. 2, we show in panels (a) and (b) the disper-sion relation ω(k z ) for systems with νELN asymmetry parameter α = 0.9 and 1.1. For α = 0.9, there are two regions in k z where complex ω exist. The maximal value of Im(ω max ) 0.044 locates at k z,max −0.165. At the same k z,max , (dω/dk z ) max 0.515. For α = 1.1 where the νELN crossing is deeper in larger v z (see Fig. 1), Im(ω max ) 0.107 at k z,max −0.225. The corresponding (dω/dk z ) max 0.278.
In the bottom panels [(c) and (d)] of Fig. 2, we show the evolution of |P ⊥ (z)| = P 2 1 + P 2 2 = 2|ρ ex | for v z = 0.995 in the liner regime for the same νELNs with point-source like perturbations. Clearly, an initial Gaussian packet of P ⊥ grows exponentially with time and the peak position of P ⊥ moves toward positive z direction for both the cases. The growth rate of the maximal value of P ⊥,max and the velocity dz(P ⊥,max )/dt both agree well with the values of Im(ω max ) and (dω/dk z ) max obtained in the dispersion relation above. The growth of perturbations using α = 1.1 is indeed much faster than that with α = 0.9. Moreover, the P ⊥ (z) with α = 1.1 spreads over both positive and negative z directions. For the case with α = 0.9, although the width of P ⊥ (z) also increases with time, it mainly drifts to the positive z direction without spreading over the negative z direction. When the perturbations grow to the nonlinear regime, wavelike oscillatory features develop [37]. In Fig. 3, we show the snapshots of P 3 (z, v z ) at different times for α = 0.9 and 1.1 using the same point-source like perturbations as in Sec. III A. For the case of α = 0.9, the flavor evolution behavior of the system is in general similar to those reported in Ref. [37]. Flavor waves develop and mainly propagate toward the positive z direction. This is consistent with the growth of perturbations in the linear regime discussed earlier. An interesting feature shown here is that although the flavor oscillations initially can affect all v z modes, illustrated by the vertical stripes in the upper two subpanels in Fig. 3(a), this effect diminishes when time proceeds and flavor conversions are roughly confined in v z > ∼ v z,c 0.65. Another interesting feature is when the forwardtraveling wave front interacts with the slowly backward propagating part after t > ∼ 1200, it pushes the whole pattern to drift toward positive z direction. Meanwhile, this interaction breaks the coherent wavelike pattern. Substructures in smaller scale develop such that the orientation of the neutrino polarization vectors varies rapidly in z for v z > ∼ v z,c . Consequently, although at each location z, neutrinos with different v z > ∼ v z,c still have |P| = 1, the "average polarization vector" over the z domain can shrink due to their misalignment. Such a flavor state was referred as "flavor depolarization" in Refs. [40,41] and we will discuss its behavior in more detail in the next section. For v z < ∼ v z,c , most neutrinos remain unaffected. For α = 1.1 shown in panel (b), flavor conversions quickly develop toward both the positive and negative z directions and produces coherent and wavelike structure, once again consistent with that indicated by the growth of perturbations in the linear regime. Similarly, when the forward and backward propagating modes interact after t > ∼ 665, smaller structures develop and cause a major part of v z reaching closer to flavor depolarization. One important difference from the previous α = 0.9 case is now flavor depolarization happens mostly in v z < ∼ v z,c 0.45. We will also discuss this feature and its consequences in Sec. IV B.

C. Nonlinear regime: Random perturbations
Now, let us look at how the flavor systems evolve when we adopt different seed of random perturbations. We show again P 3 (z, v z ) at different time snapshots in Fig. 4 for α = 0.9 and 1.1. With random perturbations, some locations have larger P ⊥ initially such that flavor conversions develop faster around those locations (see the ∼ v z,c (v z < ∼ v z,c ) for α = 0.9 (1.1) within a much shorter amount of time.

IV. RESULTS: EVOLUTION OF AVERAGED QUANTITIES
In the previous Sec. III, we have discussed how fast flavor conversions of neutrinos can develop in space and time with different seed perturbations. In this section, we further examine the time evolution of relevant quantities averaged over z and/or v z in Sec. IV A and Sec. IV B.

A. Overall flavor conversion probability in the box
We define the overall flavor survival probabilities of neutrinos and antineutrinos, P ee and P ee in the simulation domain by where the brackets and overline denote averaging over z and v z , respectively. The integration limits over z and v z are from −L z /2 to L z /2 and from −1 to 1, respectively. Fig. 5 shows the time evolution of P ee and P ee for all five different νELN spectra that we considered (see Fig. 1). The cases with point-source like perturbations and random perturbations are shown in panels (a) and (b), respectively. First, we see that cases with larger α reach the final asymptotic states earlier in time for both point-source like and random perturbations. Meanwhile, all cases with random perturbations reach the asymptotic states much earlier than those with point-source like perturbations as discussed in Sec. III C.
Second, systems with α = 1 (equal number of neutrinos and antineutrinos) have final asymptotic values of P ee f P ee f 0.5, i.e., full flavor depolarization is almost reached for the entire system. For systems that are more asymmetric, the asymptotic P ee f and P ee f are larger, i.e., on average less ν e andν e are converted.
Comparing P ee f to P ee f for a given α, one finds that cases with α > 1 have larger values of asymptotic P ee f than those of P ee f (vice versa for α < 1). This is related to the fact that for α > 1 (α < 1), neutrinos with v z < ∼ v z,c (v z > ∼ v z,c ) experience more flavor conversions. Since for our adopted neutrino angular spectra, gν(v z ) is more forward peaked in positive v z than g ν (v z ), more flavor conversions in larger v z naturally lead to smaller P ee f than P ee f .
Comparing the final values of P ee f and P ee f for cases with point-source like and random perturbations, the point-source like ones generally lead to a slightly smaller P ee f and P ee f than the random ones. This is due to the reason discussed in Sec. III B: although the interaction of the flavor waves lead to states close to flavor depolarization, for cases with point-source like perturbations, some large-scale coherent structures can remain until the end of the simulations. This can also be seen in Fig. 6, which shows the spatially averaged flavor survival probabilities P ee (v z ) f as functions of v z at the end of our simulations, where P ee (v z ) is defined by Note that P ee (v z ) that can be similarly defined are equal to P ee (v z ) , when only neutrino-neutrino selfinteraction terms are included here. Nearly flavor depolarization can be reached in one side of the νELN spectra for cases with α = 1. For α = 1, the entire νELN spectra can achieve nearly flavor depolarization. Figure 6 clearly shows several interesting features. First, for cases with either point-source like or random perturbations, those with α > 1 (α < 1) lead to nearly full flavor depolarization ( P ee (v z ) f 0.5) for velocity modes v z < v z,c (v z > v z,c ). For velocity modes in the other side of the νELN, nearly flavor equipartition cannot be achieved and P ee (v z ) f gradually increase from 0.5 to larger values for larger (smaller) v z for α > 1 (α < 1). For the case with α = 1, i.e., symmetric in neutrinos and antineutrinos, full flavor depolarization for the entire νELN can be achieved. Comparing results with different perturbations, one sees that due to the remaining large-scale coherent structures (see e.g., Fig. 3), systems with point-source like perturbations give rise to smaller Based on our simulation results, we make a conjecture as follows. For a system with the periodic boundary condition in one spatial direction (while having perfect translation symmetry in the other two directions) where interaction of flavor waves can happen, the system can evolve to an asymptotic state where nearly flavor equipartition can be reached for velocity modes in one side of the νELN, when averaged over space. In other words, when averaged over space, the νELN crossing nearly vanishes. For such a scenario to happen, the net neutrino e − x lepton number conservation L e−x [see Eq. (10)] would enforce that nearly flavor depolarization can only happen for velocity modes in the shallower side of the νELN.
An additional remark is: comparing Fig. 5 and Fig. 6, we also point out that although Fig. 6 seems to suggest that a narrower range of velocity modes in νELN distribution reaches nearly flavor depolarization for α = 0.9 than those with α > 1, on average, more neutrinos and antineutrinos are being converted. Once again, this is because the spectra g ν (v z ) and gν(v z ) are more forward peaked. This means that although a system with α < 1 may appears to be affected in less phase space volume in v z , flavor conversions there can actually have a larger impact on physical processes that are flavor dependent in realistic astrophysical environments.

B. Evolution of different velocity modes
Refs. [40,41] proposed that the behavior of the system can be understood by examining the time evolution of P(v z ) . From Eq. (9), one obtains The authors of Refs. [40,41] argued that the spatial averaging of the cross products in Eq. (13) can be approximated as the cross products of two vectors that are spatially averaged separately: By doing so, the time evolution of P(v z , t) may be understood as a vector precessing around an effective Hamiltonian vector H(v z ) = M 0 − v z M 1 (t) . Moreover, the depolarization of P(v z , t) happens when the , where H 3 differs from H 3 by a factor related to M 2 when going to a co-rotating frame (see Ref. [40] for details).
We carefully examine the validity of these claims. In Fig. 7, we show in panel (a) the time evolution of | H 3 | and H ⊥ for α = 0.9 of v z = 1 and α = 1.1 of v z = −1 using random perturbations as examples. Despite these modes undergo nearly flavor depolarization, Fig. 7(a) shows that the values of H ⊥ remain much smaller than | H 3 | nearly during the whole time, in contrast to what discussed in Ref. [41].
In panels (b) and (c) of Fig. 7, we show the evolution of | P(v z ) | and the relative angles θ r between the two vectors on the right-hand sides of Eqs. (13) and (14) for ten v z modes for the case with α = 1.1 with random perturbations. Clearly, for modes with v z < ∼ v z,c that reach flavor depolarization (darker colors), i.e., | P(v z ) | → 0, their cos θ r oscillate rapidly between -1 and 1. This indicates that Eq. (14) is in fact not a good approximation of Eq. (13).
In sum, our findings here suggest that the picture which explains the late-time behavior and the depolarization of the system proposed in Refs. [40,41] using the spatially averaged polarization vectors needs to be revisited.

V. SUMMARY AND DISCUSSIONS
We have performed long-term numerical simulations of fast collective neutrino flavor conversions for systems with translation symmetry in two spatial (x and y) directions using a newly developed code COSEν. Our numerical calculations were conducted based on the finite volume method with seventh order WENO scheme and the finite difference method with Kreiss-Oliger dissipation. Both methods showed good capability of numerical error suppression which allows simulations being carried to longer than a few thousand characteristic timescale of the system. Assuming uniform νELN distributions in z and taking periodic boundary conditions, we have investigated how flavor conversions happen with different νELN spectra, controlled by the initialν e to ν e asymmetric ratio parameter α, as well as how the results depend on the chosen flavor perturbations. Our main findings can be summarized as follows.
First, we found that for systems with point-source like initial perturbations, flavor waves with coherent structures can develop and propagate. With periodic boundary conditions, these flavor waves can then interact and break into smaller scale structures. When adopting perturbation seed randomly placed in z, the flavor waves originated from different locations can interact much faster such that no coherent structure can be formed.
The interactions of the flavor waves lead the system to a final state where part of the velocity space (v z ) is close to flavor depolarization when averaging over the space. For any asymmetric system with α = 1, only one side of the νELN spectrum relative to the νELN crossing point can reach close to an averaged flavor depolarization while neutrinos in the other side of νELN experience less flavor conversions, as constrained by the conservation of the net neutrino e − x number. Specifically, for α > 1 (α < 1), nearly flavor depolarization can be reached for ∼ v z,c ) when the antineutrinos angular distribution are more forward peaked than that of neutrinos. On the other hand, for systems with α = 1, the entire neutrino spectra can reach close to averaged flavor depolarization. This phenomenon is qualitatively similar for systems with either point-source like or random perturbations. Quantitatively, the developed large-scale coherent pattern in cases with point-source like perturbations allow part of the velocity space to have on average more flavor conversions than perfect flavor depolarization.
Comparing our results with those reported in Refs. [37,40,41], our results with point-source like perturbations agree with [37], which, however, only evolves the system for a shorter period of time without allowing flavor waves to interact. On the other hand, Refs. [40,41] obtained results nearly independent of whether the initial perturbations are being point-source like or random. In fact, behavior of random perturbations emerged in their simulations using point-source like perturbations, different from what we obtained here. One potential reason is that Refs. [40,41] used fast Fourier transform to evalu-ate the derivative terms, which might artificially generate errors of random nature in the spatial domain. Moreover, we have tried to verify the mechanisms proposed in Refs. [40,41] in explaining the occurrence of the flavor depolarization. However, our results do not support the proposed mechanisms and suggest that better understanding is needed.
Practically, our numerical findings may provide insights for efforts which attempt to include the impact of fast neutrino flavor conversions in hydrodynamical simulations; e.g., [32]. For instance, if the adopted neutrino transport scheme can provide the antineutrino-toneutrino asymmetry ratio and the crossing points in the νELN spectrum, partial flavor depolarization to one end of the νELN together with the net e − x neutrino lepton number may be applied to the neutrino distribution functions at where νELN crossings are identified.
Needless to say, there are still several improvements to be made in future. For example, our simulations were performed in reduced dimensions. How the symmetrybreaking solutions may develop and affect our conclusions need to be further examined. In our simulations, we have completely omitted the vacuum and the matter terms in the Hamiltonian, as well as the collision terms. Adding these terms and incorporating full treatment with three neutrino flavors may introduce new effects recently investigated in works without advection [63][64][65][66]. Full inclusion of them are to be implemented in future. Last but not least, the potential impact of the many-body nature of the problem remain to be further elucidated [67][68][69][70].   in the main text, our fiducial resolution is with N z = 12000 and N vz = 200. In panel (a) of Fig. 8, we compare the absolute difference of P ee between runs with different resolutions and our fiducial case, taking α = 0.9 with point-source like perturbations. It shows that the difference of P ee amounts to ∼ 10 −3 and ∼ 10 −2 , when we reduce N vz by a factor of 2 and 4, respectively while keeping N z being the same. For runs with fixed N vz = 200 but with reduced N z = 6000 and 3000, the relative differences are smaller than 10 −6 and are thus not shown in the figure. This can also be inferred by looking at the nearly identical black and red curves, which represent the differences with (N z , N vz ) = (12000, 100) and (6000, 100), respectively. However, this does not mean that we can adopt a much reduced resolution in N z . In panels (b) and (c) of Fig. 8, we show the deviation of the length of the polarization vector from unity, averaged spatially and spectrally over g ν (v z ): and the maximal deviation of the length of the polarization vectors in the entire simulation domain (δ|P|) max . Both panels clearly show that the errors are mostly depending on the resolution in N z . For our fiducial case with (N z , N vz ) = (12000, 200), the associated errors in δ|P| and (δ|P|) max are smaller than 10 −8 and 10 −5 , respectively, for α = 0.9 with point-source like perturbations. For all of our simulations with different α and different perturbation seeds, we obtain δ|P| < 10 −4 . For (δ|P|) max , all cases have (δ|P|) max < 10 −2 except for α = 1.2 and 1.3 with point-source like perturbations, for which their (δ|P|) max reach O(1) by the end of the simulation. We note, however, that demanding (δ|P|) max 1 may not be a practical convergence criterion because (δ|P|) max is usually associated with the grids with largest |v z |, which has minor contribution to the Hamiltonian and the averaged quantities. In fact, we have compared the errors in the averaged quantities for all cases with lowered resolutions (N z = 6000) and confirmed that all the averaged quantities agree within 5 × 10 −3 .

Appendix B: Results using different numerical schemes
In this appendix, we compare our results obtained with two different numerical methods described in Sec. II C. Panel (a) of Fig. 9 shows the difference of P ee between the run using the finite volume method with seventh order WENO scheme (FV+WENO) and that with the finite difference method supplied by the third-order Kreiss-Oliger dissipation term (FD+KO3), for α = 0.9 with point-source like perturbations and our fiducial resolution of (N z , N vz ) = (12000, 200). One sees that the differences are smaller than 10 −5 throughout the whole time. Panels (b) and (c) of the same figure show once again δ|P| and (δ|P|) max . Here we see that although both quantities remain much smaller than one with either the FV+WENO or the FD+KO3 scheme, the FV+WENO scheme gives rise to errors much smaller than those using FD+KO3 scheme. In addition, we show values of δ|P| and (δ|P|) max using simply the finite difference method without applying any error suppression mechanisms (FD only) in panels (b) and (c) of Fig. 9. It shows clearly that both our FV+WENO and FD+KO3 schemes help suppress numerical errors in δ|P| by more than two orders of magnitudes for this particular parameter set. for simulations with α = 0.9 and point-source like perturbations, using the finite volume method with seventh order WENO scheme (FV+WENO) and the finite difference method supplied by the third-order Kreiss-Oliger dissipation term (FD+KO3). Note that in panels (b) and (c), we show additionally the quantities derived using the finite difference method only without dissipation (FD only).