High-Resolution Particle-In-Cell Simulations of Two-Dimensional Bernstein-Greene-Kruskal Modes

We now look at a specific case with $B_{0}=4$ to compare the moment runs with normal and reversed flow in more details. Fig. 3 shows curves of $n_{e}$ vs $y$ along the $y$-axis at different instants within about three $\tau_{c}=2\pi/\omega_{\text{c}}=\pi/2$, for the normal case in panel (a) and the reversed case in panel (b). These instants are chosen when the $n_{e}$ value at the origin are roughly the largest or smallest during one $\tau_{c}$. Therefore, $n_{e}$ curves between two instants are somewhere in between the corresponding two groups of curves that are shown. The black continuous curves on both panels are the input $n_{e}$ profiles for the analytic BGK mode.

PIC simulations are known to suffer from noise due to finite particle number. 10.1063/1.862452 The curves in Fig. 3 demonstrate this, despite being from runs with our largest grid size of $2048^{2}$ cells. The noise is low enough that the pulsation amplitude for both normal and reversed case is resolved. While the noise level is about the same for the normal and reversed cases, the reversed case has a much larger pulsation amplitude. Based on this result and similar results for every other relevant case, we conclude that the difference in behavior is based on true dynamics, not numerical noise.

The fact that the normal flow case has a small pulsation amplitude, as compared with the reversed flow case, as well as the initial depletion of $n_{e}$ at the center, shows that the BGK-like initial condition does represent a more time-steady and stable case. By BGK-like, we mean this configuration only matches the three lowest-order moments (density, flow velocity, and effective temperature) of the true BGK mode solution. Such a configuration cannot represent the non-Maxwellian distribution as shown in Fig. 1, especially for low $B_{0}$. Nevertheless, the above results do show that even just matching the moments can make a configuration approximately time-steady in the strong $B_{0}$ limit.

We have put many simulation movies into the Supplemental Materials, including the two density movies corresponding to Fig. 3. The movies for these two runs provide a better contrast between the different pulsation levels for the two cases.

Generally, the pulsation level is larger for smaller $B_{0}$, and much larger for the reversed case than the normal case. The $n_{e}$ value at the center can be greater than unity at times for smaller $B_{0}$ cases. This value as a function of time turns out to be a simple visualization to show the pulsation of a run. Figs. 4, 5, and 6 show the spatial averaged $n_{e}$ value at the center, $n_{e0}$, as functions of time, for $B_{0}$ in the weak, medium, and strong ranges (with overlaps).

In Figures 4, 5 and 6, we show two sets of curves using $256^{2}$ and $512^{2}$ grids to show the possible dependence on resolution. The figures indicate that the plasma behavior is essentially the same for the two grid resolutions. Therefore, the dynamics of the moment runs are adequately resolved using a resolution of $256^{2}$ or higher.

Angular frequencies of the pulsating behavior of moment runs are given in Table 2. These values were obtained by taking a Fourier transform of the mean electron density in the middle 4 cells of the $256^{2}$ runs, or an equivalent-sized square area for higher-resolution runs. The angular frequency $\Omega$ is then $2\pi$ times the frequency with the largest amplitude. Although minor, pulsation for high-$B_{0}$ moment runs were detectable using this method.

The angular frequencies are very close to the electron gyrofrequency $\omega_{ce}$, especially for reversed runs and for runs with higher $B_{0}$. The discrepancy can in part be explained by the effects of the electron $E\times B$ drift in a cylindrically symmetric radial electric field, which increase the effective gyrofrequency for non-reversed cases more than the reversed cases, since the electric field magnitude oscillates (even with sign reversal) for the latter. This hints that the pulsation is due to gyromotion.

Another look at the velocity-space distributions in Fig. 2 suggests that electrons gyrate between the “spike” region and the central hole. The spike region is present in Eq. (5), as Fig. 7 shows and of which Fig. 1 is a vertical, renormalized cross section. Indeed, the spike region seems to lie on the line where an electron’s gyrodiameter $2\lvert v_{\phi}\rvert/B_{0}$ equals its radial coordinate $\rho$ and gyromotion would take an electron at $\rho$ inward to near the origin, and an electron near the origin outward to this $\rho$ location. This analysis led in part to support in PSC for initialization of the true electron distribution function, as described in section III.4.

Rewriting Eq. (5) in terms of $\rho,v_{\rho},$ and $v_{\phi}$ reveals it to be the difference of two Gaussians (note the normalization is by Eq. (2), not Eq. (8)):

where

The spike is caused by the shifted Gaussian having a negative coefficient, i.e., $1-1/\alpha<0$. The negative Gaussian is centered at $v_{\phi}=kB_{0}\rho^{3}/\gamma^{2}$, which for $\rho>1$ (or in PSC units: $\rho>.001$) gives $v_{\phi}\approx\frac{1}{2}B_{0}\rho$. Ignoring the effect of $v_{\rho}$, which is normally distributed with a mean of 0, this is the condition for an electron to have a gyrodiameter equal to its distance from the origin and to gyrate towards the origin.

The spike itself is not centered at the negative Gaussian, but slightly higher, due to the influence of the positive Gaussian. This shift is such that not only does $v_{\phi}\to\frac{1}{2}B_{0}\rho$, but $\rho v_{\phi}\to\frac{1}{2}B_{0}\rho^{2}$, i.e., $l\to 0$. Note that the mean of the negative Gaussian does not satisfy the latter condition.

With the pulsating behavior gone, we were able to determine whether the analytic form was time-steady or not, and if so, whether it was stable or unstable. We use $B_{0}=4$ as an example again to contrast it with the results discussed in section IV.1.1, and because the magnetic field was strong enough that the configuration was stable.

Figure 9 shows plots using data from four runs with $B_{0}=4$. In addition to the moment and moment reversed runs using $2048^{2}$ grids from Fig. 3, we show a $2048^{2}$ run initialized exactly. Its total duration is again not very long—up to $t\sim 6.3$, or about 4 $\tau_{c}$—since it is computationally expensive to run at such a high resolution. To see the long term behavior, we performed another exact $B_{0}=4$ run using a $256^{2}$ grid and ran it up to $t\sim 320\sim 204\tau_{c}$.

To show the time-steadiness of the exact runs, we plot the radial electric field component $E_{\rho}$, averaged over all azimuthal angles, as a function of $\rho$, the radial position. This averaging removes the random noise due to finite grid size (as seen in Fig. 3) but preserves the actual dynamics of time evolution. The moment and reversed moment runs effectively serve as controls of this process.

From Fig. 9(a) we see that the high-resolution exact run is indeed time-steady over one $\tau_{c}$, to the extent that all five $E_{\rho}$ profiles are on top of each other. One the contrary, profiles of $E_{\rho}$ for the high-resolution moment and moment reversed runs, shown in Fig. 9(b) and (c), have significant variations within one $\tau_{c}$, except that the profiles at the end of one $\tau_{c}$ come back to the initial profiles. To demonstrate that the solution is actually stable, Fig. 9(d) shows 500 $E_{\rho}$ profiles for the $256^{2}$ exact run over 204 $\tau_{c}$. Except for a few small deviations, these 500 curves are also almost on top of each other. The time-steadiness of the two exact runs can be seen more clearly from the movies in the Supplementary Material.

Our simulations with $B_{0}=2$ and $B_{0}=10$ were similar to those with $B_{0}=4$. The $B_{0}=10$ case is extremely steady; even the corresponding moment run only slightly pulsated. The $B_{0}=2$ cases also look stable generally, but the electron density hole is observed to have slight shifts or deformations in a long duration $256^{2}$ exact run.

Figure 10 shows the $n_{e0}$ values for exact runs using these three values of $B_{0}$. Compared to the corresponding curves of Fig. 6, $n_{e0}$ fluctuates with smaller amplitudes in exact runs than moment runs for $B_{0}=2$ and 4 cases. The fluctuations are ostensibly comparable for the $B_{0}=10$ cases, but the $B_{0}=10$ moment run has clear pulsation in the $n_{e}$ movie, whereas the movie for the exact run does not.

To the best of our knowledge, this is the first time 2D BGK mode solutions of a form similar to Eq. (5) have been confirmed with a kinetic simulation. This is remarkable in the sense that such solutions are fully nonlinear, self-consistent localized solutions of the Vlasov-Poisson system of equations, i.e., Eqs. (3) and (4). Moreover, this result also shows that there exists solutions that are stable over some ranges of parameters, at least with perturbations confined in the 2D plane, and up to rather long durations of the simulations. The stability study for more general perturbations using 3D PIC simulations is ongoing, and will be the subject of a later publication.

Exact runs with $B_{0}\lesssim 1$ were unstable. Unstable modes can appear due to the bimodal velocity distribution as in Fig. 1. The instability manifests as radial arms in $n_{e}$ that rotate counterclockwise around the origin, as shown in Figs. 11; note that the electron flow is clockwise.

The fluctuations are especially visible in the azimuthal component of the electric field $E_{\phi}$, as in 12 and in the movies in the Supplementary Material. We thus identify these rotating fluctuations as electrostatic waves. The propagation appears as roughly co-rotating at first, and then either forms a spiral pattern as it fades away.

Spiral waves are clearly seen in the $n_{e}$ and $E_{\phi}$ movies for $B_{0}=0.25$, 0.5, and 1. For the $B_{0}=0.1$ case, the initial wave has many more arms than other cases (12 or more), but the co-rotation phase quickly ends with the almost total dissipation of the wave, obscuring any possible spiral phase and making analysis difficult.

For the $B_{0}=1$ case, the instability takes a long time ($t>100$) to appear, but the spiral wave lasts a long time, persisting until we had to stop the simulation at $t=433$. This case, as well as the $B_{0}=0.5$ case, has only four arms, while the $B_{0}=0.25$ case has between 6 and 8. Some of these arms are not present over the full range of $\rho$ for which there are unstable waves; rather, there appears to be bifurcations such that one arm can split into two when going radially outwards, as can be seen in Fig. 12. Generally, not all arms are of the same strength, suggesting a mixture of unstable modes instead of one single mode.

Given the spiral structure of the density fluctuations, the waves should also have $E_{\rho}$ fluctuations such that the wave has a component propagating outwards along the $\rho$ direction. The $E_{\rho}$ fluctuations are more difficult to see because $E_{\rho}$ from the background BGK mode is large compared with the fluctuations. We include one $E_{\rho}$ movie in the Supplemental Materials for the $B_{0}=1$ case to illustrate this point. In fact, the wave starts to have more apparent outward propagation when the spiral pattern emerges, and this could be the physical reason behind the formation of the spiral pattern.

It is observed that the electron density hole begins to fill as soon as the instability appears. This can be seen from the $n_{e0}$ vs. $t$ plots, e.g. Fig. 13 for cases with weak $B_{0}$ and Fig. 14 for cases with medium $B_{0}$.

These two figures are significantly different from corresponding figures from the moment runs, i.e., Figs. 4 and 5. The latter exhibit pulsation, while the former appear to follow logistic curves and saturate at a value below 1. This is clearly show in Fig. 15. We see that the growth of $n_{e0}$ starts near $t\sim 0$ for the $B_{0}=0.1$ and $0.25$ cases, while there may be a linear phase in the growth of $n_{e0}$ for the $B_{0}=0.5$ and 1 cases, as it takes some time for the instability to grow.

As far as our simulations can show, $n_{e0}$ apparently stabilizes at a value below 1 for these unstable cases. In other words, there is still an electron density hole after the initial hole decays, albeit one with an $n_{e0}$ closer to 1. The spike also seems to survive the instability. The snapshots of the exact $B_{0}=1$ run in Fig. 8 show such a state in the last two frames, while the second frame captures a glimpse of the radial arms. If an electron density hole truly does exist after the instability, it is unclear whether it would be in a stable time-steady state similar to the solutions described by Eq. (5).

It is unknown what critical value of $B_{0}$ separates stability from instability, since we have not simulated enough cases with different $B_{0}$. The apparent stability of high-$B_{0}$ runs indicates that such a value exists, and our limited number of cases suggest that it is around $B_{0}\sim 2$ for our choice of other parameters $h$ and $k$. The chosen form of the analytic BGK mode distribution shown in Eq. (5) implies that the critical value of $B_{0}$ depends on these parameters. For example, using a smaller value of $h$ for the same $B_{0}$ will make the distribution more Maxwellian, and a Maxwellian distribution is well known to be stable.

Due to limitations of space and scope, we omit a more thorough analysis of the observed instabilities, including numerical results and a comparison to a theoretical model. We will present such discussion in a separate publication instead.