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

The collisionless Boltzmann equation, or the Vlasov equation, is an equation that describes the evolution of a distribution function $f$ in an $N$-dimensional phase space, where $N$ is the number of spatial and momentum coordinates being considered. Using vector quantities for velocity, v, position, r, and acceleration, a, it reads:

BGK modes are exact, nonlinear solutions to the time-independent form of this equation where the acceleration is given by the Lorentz force equation so that it becomes, for species $s$,

where $f_{s}$ is the distribution function describing the electron species ($s=e$) or the ion species ($s=i$), and $q_{s}$ and $m_{s}$ are the associated charge and mass for each species, respectively. In this paper, we utilize only a single ion species. Solutions to this equation coupled with Poisson’s equation,

for a charge density $\rho_{q}$, are what is typically needed for BGK mode solutions. As mentioned previously, we will also make use of the time-steady Ampère’s law for initialization, given here as:

Using the zeroth-order moment of the distribution function, $\rho_{q}=\sum_{s}q_{s}\iiint f_{s}(\textbf{r},\textbf{v})\,d^{3}v$, Poisson’s equation becomes:

where we have assumed a uniform background ion density, $n_{0}$. Similarly, we utilize the first-order moment of the distribution function to rewrite Ampère’s law as

where ions are assumed to be stationary and do not contribute to the current density.

Following the recent work of McClung et al., we assume cylindrical symmetry of the electron distribution function, $f_{e}(\textbf{r},\textbf{v})=f_{e}(\rho,v_{\rho},v_{\phi},v_{x})$, the electric scalar potential, $\Psi=\Psi(\rho)$, and the magnetic vector potential, $\textbf{A}=A_{x}(\rho)\hat{x}+A_{\phi}(\rho)\hat{\phi}$, such that $\textbf{B}=\nabla\times\textbf{A}=-\frac{dA_{x}}{d\rho}\hat{\phi}+\frac{1}{\rho}\frac{d}{d\rho}({\rho}A_{\phi})\hat{x}$. Note that we have chosen $\hat{x}$ to be the axial direction instead of the typical $\hat{z}$ to match the convention of PSC. For uniform B, we require that $A_{\phi}=B_{0}\rho/2$, where $B_{0}$ is the background magnetic field parameter, and $A_{x}$ be a constant. This uniform B was what we used in previous simulations, but we can now initialize with a nonuniform B that satisfies Ampère’s law.

For normalization, we use the same scheme given in McClung et al.McClung et al. (2024) We normalize such that the electron thermal velocity, $v_{e}$, is the unit for v, the Debye length, $\lambda_{D}=v_{e}/\omega_{pe}$, is the unit for r, $n_{0}e\lambda_{D}^{2}/\epsilon_{0}$ is the unit for $\Psi$, $n_{0}e\lambda_{D}/\epsilon_{0}v_{e}$ is the unit for B, and $n_{0}/v_{e}^{3}$ is the unit for $f_{e}$. In the simulations, several physical constants are set to the following:

With these assumptions and normalizations, the forms of the Vlasov, Poisson, and Ampère equations in cylindrical coordinates respectively become, for an electron distribution function,

and

where $\beta_{e}=v_{e}/c$ and $c$ is the speed of light.

The BGK solution is then an electron distribution function that satisfies these equations self-consistently. Since it has been shown that a distribution function that depends only on the total normalized energy, $w=v^{2}/2-\Psi$, is not a valid solution because it does not remain localized,Ng and Bhattacharjee (2005) we utilize a distribution function that depends on additional constants of motion, such as the normalized canonical angular momentum, $l=\rho(v_{\phi}-A_{\phi})$, and linear momentum, $p=v_{x}-A_{x}$. One possible solution is

where $h_{0}$, $k$, and $\xi$ are parameters.Ng (2020) In the simulations presented here, $\xi$ is set to 0 for simplicity, so that the solution only depends on electron angular momentum and total energy. The $h_{0}$ parameter can be set to distinguish between two cases regarding the density of electrons in the center of the structure, which is in the center of the simulation domain. For the EDH, $h_{0}$ is positive and less than one, and the central normalized electron density is less than unity. For the EDB, $h_{0}$ is negative and the central normalized electron density is greater than unity. Once this form of the distribution function is assumed and the parameters are assigned, numerical integration is performed to solve for the corresponding electric scalar and magnetic vector potentials to be initialized in the simulation. We also initialize the particle velocities according to the exact method of initialization mentioned in McClung et al. McClung et al. (2024)

Because of the form of the distribution function, we can expect instabilities to form. As we showed in our past work with the EDH, the reduced distribution function for the azimuthal velocity, $f_{e}(v_{\phi})$, has a “bump-on-tail” shape at some radial positions. Velocity distributions of this form are well known to possibly lead to an instability.Chen (2019) The distribution function gives a clockwise flow to electrons around the center. Since the bump on the tail is at a positive value of $v_{\phi}$, we expect an instability that propagates in the counter-clockwise direction.

Similarly, we can look at the reduced distribution function for $v_{\phi}$ for the EDB case as well to predict instability. In Figure 1, we see the same sort of bump-on-tail shape of the distribution for the lower value of $B_{0}$. The bump is not as obvious for smaller $\rho$, but becomes more apparent as $\rho$ increases in the distribution function. At $\rho$=10$\lambda_{D}$, the bump (or spike, rather) is much narrower and sharper than it was for the EDH, leading to unstable behavior. On the other hand, the reduced distribution function of $v_{\phi}$ is more indicative of a single-humped distribution for larger $B_{0}$, as shown in Figure 2. This leads to more stable behavior in the simulations for larger $B_{0}$.

We initially ran our simulations without solving Ampère’s law for the solution. This was a valid approximation because we set the electron thermal velocity to be a reasonably small fraction of the speed of light. We do this by adjusting the parameter of $\beta_{e}$. This allowed us to approximate the normalized $\nabla\times\textbf{B}=-\beta_{e}^{2}\iiint f\textbf{v}\,d\textbf{v}\approx 0$. It was found that the magnetic field variance was minimal for these parameters. See McClung et al. McClung et al. (2024) for a more thorough description of this approximation. As mentioned earlier, we have now included Ampère’s law in the solution, so we are no longer limited to very small $\beta_{e}$, as long as it is not too close to unity to invalidate the non-relativistic theory. Since the PSC simulations are relativistic, the fact that there are good agreements between simulations and theory shows that the non-relativistic solutions are valid for the $\beta_{e}$ values used in this study.

As mentioned previously, the ions are assumed to form a uniform static background. This approximation can be generalized in the future to more realistic distributions of ions, such as a Boltzmann distribution.

The setup and units of the simulations are the same as they were in McClung et al. McClung et al. (2024) Here are the key points that affect the units of the simulations’ spatial and temporal domains.

The normalizations given in Eqs. 7e and 7f imply that the Debye length $\lambda_{\text{D}}\equiv v_{\text{e}}/\omega_{\text{pe}}=1$. The magnetic field strength parameter $B_{0}$ has units $em_{\text{e}}/n_{\text{e}}\omega_{\text{pe}}$. Note that $B_{0}=\omega_{\text{ce}}/\omega_{\text{pe}}$, where $\omega_{\text{ce}}\equiv eB/m_{\text{e}}$ is the electron gyrofrequency. Thus, the electron gyroperiod $\tau_{\text{ce}}\equiv 2\pi/\omega_{\text{ce}}$ is $2\pi/B_{0}$.

PSC simulation results and figures shown are normalized according to Eq. 7 with the change

Thus, time still has units of $\omega_{\text{pe}}^{-1}$, but distance now has units of the electron skin depth, $c\omega_{\text{pe}}^{-1}$. This multiplies the length scale compared to Eq. 7 by a factor of $\beta_{\text{e}}$, or $\lambda_{D}=\beta_{e}.$

There are two main sets of parameters used for the presented simulations. The first set of parameters, labeled as the EDH case, is $h_{0}=0.9$, $k=0.4$, which is the same set of parameters used in the 2D runs of McClung et al. McClung et al. (2024) Furthermore, within the two cases, we identify two more cases for different values of $B_{0}$. $B_{0}=1$ will be the “quasi-stable” case, and $B_{0}=0.25$ will be the “unstable case.” This identification follows from our findings in McClung et al. that $B_{0}=0.25$ produced an instability quicker than $B_{0}=1$. McClung et al. (2024) The EDB uses two other sets of parameters. This was done because it was found that a similar set of parameters to the EDH produced much more stable results. Therefore, for the EDB quasi-stable case, $h_{0}=-0.9$, and $k=0.4$ were used. For the unstable case, $h_{0}=-6.0$, and $k=4.8$ were used; this second set of parameters produced a qualitatively similar initial electric potential to the unstable EDH run. An interesting feature of a EDB solution is that it has a negative electric potential so that there is no electrostatic electron trapping. This is not possible for a localized 1D BGK mode, since the electron trapping is essential for the existence of solutions.

Within these cases, we include runs with the generalizations we made to the solution and simulation. To distinguish the cases where we change the $\beta_{e}$ parameter, we will explicitly list the $\beta_{e}$ parameter being used in the simulation. We will also list whether the simulation used VP initialization or VM initialization. When addressing the individual simulations, we will use shorthand to denote the run being discussed (e.g., EDH-VP represents a run that is an electron density hole that was initialized using the Vlasov-Poisson method of initialization). A table of cases being discussed in this paper (along with resolutions) is included in Table 1.

The same method of velocity initialization outlined in our previous work is used in these simulations as well. Initialized fields and quantities, such as $\Psi$, B, charge density $\rho_{q}$, are initialized in the same way as McClung, et al. McClung et al. (2024) for the VP cases, such that they were evaluated at discreet values of $\rho$ and interpolated onto the simulation grid. In 3D, each level of the simulation in the axial direction is initialized with different random particle velocities following the same distribution function. For the VM cases, the same methods of initialization are used for all fields and quantities except for B. We now utilize Ampère’s law to solve for B and interpolate it onto the grid, instead of the uniform B used in our previous work. Finally, the simulations are all initialized with 100 particles per cell.

For the EDH, we ran 2D and 3D simulations for two different values of $B_{0}$ that were found to be either quasi-stable or unstable from our previous work in 2D. All of the presented EDH runs use $\beta_{e}=0.001$ unless otherwise noted.

The spatial domain volume for the quasi-stable 3D EDH-VP run was around $8,000\lambda_{D}^{3}$ or (20$\lambda_{D})^{3}$, which was determined following the same scheme as McClung et al. for the 2D simulations. McClung et al. (2024) In this case, we found that the 3D results were qualitatively similar to the 2D results. The instability that was seen in our previous work in 2D takes a long time (around $t\approx 400$) to develop. Because of the length of time required to run the simulation until the onset of the instability, we limited the 3D run to a just a few electron gyroperiods to see if the initial stable equilibrium was altered by the inclusion of the third dimension. Notably, there seems to be no new instability that develops because of the added third dimension. Essentially, the 3D quasi-stable EDH-VP simulation seems to confirm the time-steadiness, at least initially, in the same way as the 2D run of the same name. This can be seen more clearly in Figure 3.

For the unstable 3D EDH-VP simulation, we utilized a simulation spatial domain size of around 70$\lambda_{D}$ in each direction, giving a spatial resolution that is 0.27$\lambda_{D}$. This simulation domain size is again decided according to the convention outlined in McClung et al. In this case, we saw different results between the 2D and 3D runs, as can be seen in Figure 4.

Qualitatively, the two runs are identical at the start. However, the 3D run seems to maintain stability for a longer time than the 2D run. We can see from the mean value of $n_{e}$ at the center of the simulation that the 3D run takes $\approx 1\tau_{ce}$ longer to reach a similar level of electron density hole filling to the 2D case, which can be seen in Figure 5. Like in our previous work, this value is calculated as the average $n_{e}$ across the middle four grid cells for runs with a resolution of $256^{2}$ or an equivalent area for higher resolution runs.McClung et al. (2024)

Once both of the runs reach a similar time in their instabilities, the electrostatic wave that propagates in the azimuthal direction that we saw in our 2D runs becomes visible, as can be seen in Figure 6. One remarkable feature of the 3D run is that it seems that the instability is in-phase between all levels of the simulation. This in-phase behavior is apparent when comparing the almost identical bottom four rows of Figures 4 and 6. It seems that the instability is slightly out of phase initially, but becomes in-phase over a short time scale. As a consequence of the structure evolving in-phase on all levels of the simulation, we see that there is only a slight deviation of the average $n_{e}$ at the center of the simulation from level to level, depicted in Figure 7. It is possible that Landau damping is the cause for this in-phase nature of the simulation; however, more work needs to be done to explain this phenomenon.

Next, we present the results of the EDB-VM runs with $\beta_{e}=0.01$ unless otherwise specified. This means that the spatial domain size being used in these simulations will be a factor of 10 larger than the domain sizes used in the EDH case. We will defer the discussion of the validity of the VM initialization with a larger $\beta_{e}$ in Sub-section C below. For now, we will only look at the results of the quasi-stable and unstable EDB simulations. These cases have a few similarities to the EDH cases, but there are differences.

For the quasi-stable case, we see the same behavior in both 2D and 3D, much like the EDH. A comparison of the quasi-stable case is shown in Figure 8. Similar to the EDH, we found that an instability starts to develop in 2D at around $t\approx 400$, and the added third dimension does not induce an additional instability.

The much more notable runs for the EDB-VM case can be seen in the unstable case. Here we see some similarities to the EDH, but also some key differences. First, we see in Figure 9 that the instability again takes longer to develop in 3D than it does in 2D. Also, we see that there is a significant oscillation of the density bump in both 2D and 3D before the instability develops. Figure 10 shows the delay in instability development as well as the oscillations, and Figure 11 shows the electron density oscillations. This oscillation seems to occur twice over the course of 2$\tau_{ce}$ with a smaller third oscillation in 3D at around 3$\tau_{ce}$. Again, the delay in 3D seems to be about 1$\tau_{ce}$.

Additionally, we see the same in-phase behavior we saw in the EDH case. The azimuthal electric field heat map comparisons shown in Figure 12 again show in-phase behavior in sharper detail. One other feature of this case is the existence of what seems to be two concentric spirals: one smaller amplitude spiral in the center, and a larger amplitude spiral towards the edge of the simulation (however, these spirals could just as well be parts of one larger spiral with a sharp discontinuity at $\rho\approx 0.1$). As the simulation continues, the center spiral seems to fade away, leaving only the larger spiral on the outside. It is currently unknown why these two spirals form in this case and not in the EDH case, and more work is required to understand them.

In this section, we present two possible reasons for the observed delay of the instability in 3D. We also further discuss the in-phase characteristic of the 3D runs.

The first possible reason for the delay is a reduction in particle noise in the simulation. We saw in our previous work that a decrease in particle noise (an increase in particles-per-cell and, therefore, total particles) leads to a delay in the start of the instability for runs that were otherwise the same. McClung et al. (2024) Since the 2D and 3D runs presented here all have the same number of particles-per-cell, the total number of particles is greater in the 3D runs and the particle noise is therefore reduced, so we see this effect.

The second reason for the delay is attributed to the in-phase nature of the instability, the other notable feature of the 3D runs. It is interesting that it seems all levels of the simulation “wait” to work out the fastest growing mode and become unstable afterwards. This “waiting” could also be characterized as fluctuations of the electric field and electron density along the axial direction becoming Landau damped. If there were significant changes in electron density between levels of the simulation, the axial electric field, $E_{x}$, would become large. Landau damping in the simulation keeps this from occurring. This leads to the same instability spiral pattern on all levels developing for the 3D simulation, instead of many different unstable modes evolving at different levels. This period of delay between all levels is what contributes to the overall delay of the instability. It is the combined effects of reduced particle noise and Landau damping that leads to a delay in the instability in 3D, as well as the instability evolving in-phase.

Here, we look at the two generalizations we have made to the assumed solution and simulations. Namely, we wish to investigate the use of larger $\beta_{e}$ and VM initialization. If we can make these generalizations to the BGK mode solutions and simulations, it will have positive effects moving forward.

As for the larger $\beta_{e}$, it was found that there was no notable difference between the results from two runs that only differed in their values of $\beta_{e}$, aside from a change in spatial scale and velocity scale. For the sake of space and having fewer figures here, we simply acknowledge this is the case and leave these two sets of results in the supplementary material. Since we have identified that this is the case, we are no longer limited to small values of $\beta_{e}$. This means we can use values of $v_{e}$ that are better suited for our approximation of a collisionless plasma. Also, because of the nature of the Courant–Friedrichs–Lewy (CFL) condition that PSC uses, since these higher $\beta_{e}$ cases result in a simulation that has a larger spatial scale, then they also have a larger time step. Germaschewski et al. (2016) Therefore, as we increase $\beta_{e}$, the computation time required to run a simulation to a desired physical time decreases. Also, a larger $\beta_{e}$ means that we can no longer assume a constant background magnetic field, and Ampère’s law must be solved.

Our generalization of solving for Ampère’s law does seem to have an impact on the time-steadiness of a few field quantities. This is to be expected, since, while the value of $\beta_{e}$ was small in the VP runs, it was nonzero. This leads to a variation of the initially uniform magnetic field on the order of $\beta_{e}^{2}$ from Eq. (10). As can be seen in Figure 13, the VP method of initialization has this small variance and oscillates in space and time, while the VM initialization does not. These plots consist of azimuthal averages of $B_{x}$ as a function of $\rho$ at different time-steps of the simulation. The VM initialization results in an initial equilibrium of $B_{x}$, which then begins to vary over time as the instability develops towards the end of the simulation. The deviation from the theoretical value of $B_{x}$ near $\rho\approx 0$ is due to this instability, which is the azimuthal electrostatic wave mentioned earlier. Additionally, we saw a decrease in the time variation of $E_{\phi}$ when utilizing the VM method. In the VP case, the significantly time varying $B_{x}$ would induce a significant change in $E_{\phi}$ from Maxwell’s equations. For VM, this $B_{x}$ variation is nearly eliminated at larger $\rho$, so we no longer see a significant change in $E_{\phi}$. This difference is seen in Figure 14. Moving forward, we plan to use the VM method of initialization, as it is more exact and leads to more time-steady results for $B_{x}$ and $E_{\phi}$. The VM method of initialization did not appear to alter the overall instability seen in the VP runs.

More work is needed to explain the in-phase nature of the 3D runs, as well as the concentric spirals seen for the EDB. Currently, a paper discussing the nature of the EDH instability, such as its magnitude and growth rate, is in progress. Also, more visualization that can highlight the differences between the 2D and 3D cases is desired.

In the future, we plan to also investigate how this BGK mode evolves when the $\xi$ parameter is nonzero. More work is needed to be able to initialize this case, so we anticipate this to be in our future work.