The role of noise in PIC and Vlasov simulations of the Buneman instability

12/23/2021
by   Arash Tavassoli, et al.
University of Saskatchewan
0

The effects of noise in particle-in-cell (PIC) and Vlasov simulations of the Buneman instability in unmagnetized plasmas are studied. It is found that, in the regime of low drift velocity, the linear stage of the instability in PIC simulations differs significantly from the theoretical predictions, whereas in the Vlasov simulations it does not. A series of highly resolved PIC simulations with increasingly large numbers of macroparticles per cell is performed using a number of different PIC codes. All the simulations predict highly similar growth rates that are several times larger than those calculated from the linear theory. As a result, we find that the true convergence of the PIC simulations in the linear regime is elusive to achieve in practice and can easily be misidentified. The discrepancy between the theoretical and observed growth rates is attributed to the initial noise inherently present in PIC simulations, but not in Vlasov simulations, that causes particle trapping even though the fraction of trapped particles is low. We show analytically that even weak distortions of the electron velocity distribution function (such as flattening due to particle trapping) result in significant modifications of the growth rates. It is also found that the common quiet-start method for PIC simulations leads to more accurate growth rates but only if the maximum growth rate mode is perturbed initially. We demonstrate that the quiet-start method does not completely remedy the noise problem because the simulations generally exhibit inconsistencies with the linear theory.

READ FULL TEXT VIEW PDF

Authors

page 18

10/15/2021

New applications for the Boris Spectral Deferred Correction algorithm for plasma simulations

The paper investigates two new use cases for the Boris Spectral Deferred...
06/09/2011

Predicting growth fluctuation in network economy

This study presents a method to predict the growth fluctuation of firms ...
10/05/2021

Inference and De-Noising of Non-Gaussian Particle Distribution Functions: A Generative Modeling Approach

The particle-in-cell numerical method of plasma physics balances a trade...
04/17/2018

An Efficient SIMD Implementation of Pseudo-Verlet Lists for Neighbour Interactions in Particle-Based Codes

In particle-based simulations, neighbour finding (i.e finding pairs of p...
03/18/2020

The influence of initial perturbation power spectra on the growth of a turbulent mixing layer induced by Richtmyer-Meshkov instability

This paper investigates the influence of different broadband perturbatio...
10/06/2020

Inferring Microbial Biomass Yield and Cell Weight using Probabilistic Macrochemical Modeling

Growth rates and biomass yields are key descriptors used in microbiology...
01/16/2019

New particle representations for ergodic McKean-Vlasov SDEs

The aim of this paper is to introduce several new particle representatio...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Kinetic simulations are a powerful tool for studies of the linear and nonlinear behavior of plasmas. The particle-in-cell (PIC) method and the continuum Vlasov method are two widely used simulation approaches. The PIC method, which has been available for several decades, has successfully captured many physical phenomena in various scenarios. The PIC method, however, is also known for relatively large levels of numerical noise introduced by the discretization and the limited number of macroparticles used to resolve the phase space[langdon1979kinetic]. The noise in PIC simulations exists during the first time step (initial noise), but it also evolves during simulations, thus affecting the results. An alternative to the PIC method, the continuum Vlasov method, is well known as a method that is free of statistical noise. The availability of high-performance computational resources has led to a steady increase of interest and applications of continuum simulations for many physical phenomena and situations that are poorly amenable to the PIC approach [juno2018discontinuous, von2014vlasiator].

It is well known that the noise in the PIC simulations may significantly undermine the physical results of the simulations. For example, the PIC simulations of electron temperature gradient modes[lin2005role, lin2005Electron] yielded a level of turbulent heat transport that deviated greatly from results of gyrokinetic Vlasov simulations of Refs. dorland2000electron, jenko2000electron, jenko2002prediction. The origin of this discrepancy is investigated in Ref. nevins2005discrete. It is shown that the discrete particle noise effects in the PIC simulations of Refs. lin2005role, lin2005Electron undermine the dynamics of the instability, strongly modifying the predictions for the heat transport levels [nevins2005discrete, holod2007statistical]. The role of the PIC noise has also been discussed in the study of electrodynamic filamentation instability [palodhi2019counterstreaming]. It is shown that the noise in PIC simulation affects the mechanism of the instability and results in an incorrect instability threshold. In another study, it is shown that the noise of the PIC simulations can lead to significant artificial heating of plasma in the presence of the Monte Carlo collision operator [turner2006kinetic]. In Refs. nevins2006characterizing, juno2020noise, holod2007statistical, it is emphasized that the role of the discrete particle noise in PIC simulations has to be carefully analyzed and evaluated for each physical situation. For this purpose, several approaches have been proposed in the literature[riva2017methodology, holod2007statistical, kesting2015propagation, turner2016verification, turner2013simulation]. One approach is to benchmark the physical results with different simulation methods in order to build confidence in simulation results and determine the roots of discrepancies and possible numerical artifacts. Benchmarking has been successfully used as a verification tool for numerical codes in several publications [turner2013simulation, villafana20212d, charoy20192d]. One feature of benchmarking is that it tests the entire simulation code as opposed to individual units, and it can also be used on the specific problems of interest rather than synthetic test cases [turner2013simulation, turner2016verification]. Benchmarking with different numerical methods, such as PIC and Vlasov methods, provides additional confidence in the reliability of the simulations as well as highlighting the causes of possible discrepancies.

In this study, we use several PIC and Vlasov codes to investigate the impact of noise in PIC simulations of the Buneman instability. We show that the linear growth rates of the instability are significantly affected by the noise inherent to the PIC simulations. We identify the trapping of electrons (a nonlinear effect) in the early noise-driven potential as a source of the inconsistencies with the linear theory. This relationship is confirmed by continuum (Vlasov) simulations for the same parameters and initial states (and respectively the same level of noise) as in the corresponding PIC simulations. It is also supported by analytical calculations that show a high sensitivity of the linear growth rates in this problem to small distortions of the Maxwellian velocity distribution function. Therefore, we propose that early trapping of electrons induces a small plateau in the velocity distribution function, leading to the much higher linear growth rates observed in PIC simulations.

The similarities and differences between the PIC and Vlasov simulations are presented through a number of simulations. We begin with PIC simulations of the cold-plasma limit, when is relatively large. In this case, the simulated growth rates are shown to be consistent with the theoretical ones. A set of simulations is then presented for a relatively low value of the streaming velocity, , where is the thermal velocity of the electrons, is the initial temperature of the electrons, and is the electron mass. In each simulation, we measure the linear growth rates of several modes and compare them with the results of the linear theory. Some Vlasov simulations are started with an extremely small perturbation that is required by this method to excite the instability. We refer to these simulations as “low-noise” Vlasov simulations (VL1 and VL2 in Table 2). The growth rates measured by the low-noise Vlasov simulations are shown to be consistent with the linear theory. On the other hand, some PIC simulations are started with macroparticles randomly distributed in phase space. We refer to these simulations as “random-start” PIC simulations (PIC1, PIC2, PIC3, PIC4, and PIC 5 in Table 2). The growth rates measured using random-start PIC simulations deviate significantly (up to a factor of 3) from the linear theory. This discrepancy in linear growth rates persists in the random-start PIC simulations using up to macroparticles per cell. In addition, we show that starting a Vlasov simulation with the same level of initial noise as the PIC simulations (VL3 in Table 2) leads to a similar discrepancy between the simulated and the theoretical growth rates.

Reflecting its statistical origin, the noise in PIC simulations scales as , where is the number of macroparticles in each grid cell. The initial noise is a result of the random distribution of the macroparticles in phase space before the first time step. To reduce the adverse effects of the initial noise, a “quiet-start” initialization has been proposed[dawson1983particle, byers1970perpendicularly]. In contrast to the random-start method, in the quiet-start method, the macroparticles are distributed regularly or semi-regularly with appropriate weights in phase space. Accordingly, the initial noise level is made much smaller. We show that using the quiet-start method does not completely solve the noise problem in PIC simulations. Another outcome of the current study is to show how the quiet-start method should be used to improve the accuracy of the observed linear growth rates in PIC simulations. We first show that although the quiet-start can improve the results by reducing initial noise, the growth of modes is still subject to statistical noise, making an accurate measurement of the linear growth rates difficult. However, initially perturbing the mode with the maximum growth rate helps to achieve better consistency with the linear theory.

The outline of the next sections is as follows. In section II, we review the linear theory of our problem and introduce the general setup for the simulations. In section III, we show some results for a large () value that show good agreement between the theoretical linear growth rates and the growth rates measured from the theory. In section IV, we show various simulations with the PIC and Vlasov methods. As a result, we show how the initial noise of random-start PIC simulations adversely influences with the linear growth and undermines the accuracy of growth rate measurements, a problem that does not appear in the low-noise Vlasov simulations. In section V, we show that a small flattening in the distribution function can increase the observed linear growth rates by several factors. This provides a hypothesis as to source of the problem in random-start PIC simulations. In section VI, we show that, although it does lead to some improvements, the quiet-start method is unlikely to completely solve the problem of the noise in PIC simulation. In section VII, we summarize the conclusions of this study.

Ii The Buneman instability and the problem setup

The Buneman-type instabilities are driven by the relative drift of electrons with respect to ions in an unmagnetized plasma. The instabilities can be categorized according to the magnitude of [galeev1984current]. The Buneman instability regime occurs for . On the other hand, the ion-sound instability occurs for , where , is the initial temperature of the electrons, and is the ion mass. Streaming Buneman-type instabilities are important in many topical problems of plasma physics. For example, they are considered as candidates for explaining the turbulence and anomalous resistivity in solar plasmas[ChePoP2017, HellingerGRL2004, CheHH_MPLA2016] and hollow cathode plasmas in Hall thrusters [HaraPSST2019] as well as sources of nonlinear effects in ion-beam fusion applications[StartsevPoP2006].

The Buneman instability has been broadly investigated through numerical simulation [lampe1974two, niknam2011simulation, rajawat2017particle, tavassoli2021backward]. Most of the numerical simulations focus on the nonlinear regimes of the instability, assuming that the linear regime is well understood via analytical dispersion relations. However, the comparison of the linear regime in numerical simulations with the linear theory provides a valuable test for the simulation methods, revealing the validity range of the linear approximation for a particular approach. In this study, we focus on the verification of the linear regime of the Buneman instability in PIC and Vlasov simulations.

The considered equations in the setup of our problem are

(1)

where is the distribution function for ions and electrons, respectively, is the electric field, are the ion and electron densities, and is the charge, which is for the ions and for the electrons. The ions are taken to be Hydrogen with mass amu. The initial temperature for both ions and electrons is eV; the initial plasma density is . The initial conditions are

(2)
(3)

The quantities , and parameterize a small initial perturbation. In the low-noise Vlasov and quiet-start PIC simulations, these parameters are required to excite the instability. In the low-noise Vlasov simulations, we take , while in the quiet-start PIC simulations we take . In the random-start PIC simulations, there is no need for this perturbation, and we take . In all the simulations reported, we use periodic boundary conditions in a system of length 6 mm discretized with a grid of 2048 points. This length is large enough to allow excitation of several modes with mode numbers . The time step used in the simulations is ns. All the time-dependent data are collected at intervals of 500 . The relative drift between ions and electrons () drives the instability in several modes identified by the linear dispersion relation

(4)

where with the linear growth rate, the frequency,

is the wave vector,

is the ion plasma frequency, is the electron plasma frequency, and is the plasma dispersion function.

Iii Linear growth rates from PIC simulations for large drift velocity, .

By choosing the drift velocity of , we approach the cold-plasma limit of the Buneman instability. We perform PIC simulations in this limit with macroparticles per cell. Fig. 0(a) shows the growth of some select modes. These modes are chosen for the linear growth analysis of the case and include the maximum growth rate mode . We can see a distinct linear growth region in the early evolution of the modes. By fitting a line to this region, we calculate the growth rate of each mode. In Table 1

, the calculated growth rates are shown to be quite consistent with the results from the linear theory. The standard error (SE) associated with the measurement of the growth rates is also reported in this table. The SEs of the fits are quite small, showing that the growth of the chosen modes is quite linear and not oscillatory. Our other investigations (not reported here) show that even for as few as

macroparticles per cell, PIC simulations with produce accurate linear growth rates.

The mean of the derivative of the spectral growth over the same time period is an equivalent measure of the growth rate. Variations about the mean provide a measure of the power of the noise present in the growth region. The square of the growth rate over the square of the power of noise was calculated as the signal-to-noise ratio (SNR) of the growth rate

[Lynn1989]. The average SNR of the chosen modes is 21.25 dB, which is much greater than 1. This indicates that the power of noise carried in this case is quite low in the linear growth region. In all the simulations reported in this study, we see that the value of the SNR does not vary much among the chosen modes. Therefore, we only report the SNR averaged over the four chosen modes of each simulation.

((a))
((b))
Figure 1: a) The evolution of individual modes of the electric field. The dashed black line shows the fitted line on the mode. b) The evolution of the electrostatic energy. Both figures are from VSim PIC simulations for the case (see below).
(Theory)
s
(Simulation)
s
SE (Simulation)
%
14 7.08 7.09 1.03
15 8.67 8.60 0.30
16 9.46 9.20 0.24
17 8.13 8.12 0.10
Average 8.34 8.25 0.42
Table 1: Comparison of theoretical growth rates with growth rates observed in VSim PIC simulations with .

Iv Linear growth rates from Vlasov and random-start PIC simulations for low drift velocity, .

In this section, we report on linear growth rates from several PIC and Vlasov simulations for the case of . As we show, this regime of relatively low drift velocity can be problematic for the PIC simulations. Therefore, we investigate this regime more extensively by performing several PIC and Vlasov simulations. Due to the large number of these simulations, we assign a specific name to each one in this regime. These simulations are listed and described in Table 2.

Simulation Numerical code Initial condition macroparticles per cell
VL1 Semi-Lagrangian Vlasov perturbed
VL2 BOUT++ perturbed
PIC1 EDIPIC Random start, no perturbation
PIC2 VSim Random start, no perturbation
PIC3 XES1 Random start, no perturbation
PIC4 EDIPIC Random start, no perturbation
PIC5 VSim Random start, no perturbation
VL3 Semi-Lagrangian Vlasov Identical to PIC2
PIC6 EDIPIC Quiet-start, perturbed
PIC7 VSim Quiet-start, perturbed
PIC8 XES1 Quiet-start, perturbed
PIC9 VSim Quiet-start, perturbed
PIC10 EDIPIC Quiet-start, no perturbation
PIC11 EDIPIC Quiet-start, perturbed
PIC12 EDIPIC Quiet-start, perturbed
Table 2: The list of simulations with .

The first simulation (VL1) is performed by a locally developed semi-Lagrangian code. The semi-Lagrangian Vlasov scheme is a well-known and tested scheme for solving the Vlasov–Poisson equations[cheng1976integration, gagne1977splitting]

. In this scheme, the Vlasov equation is split into a convection equation and a force equation. Each of these equations is then solved by the method of characteristics using cubic spline interpolation. The Poisson equation is solved by a spectral method, the FFT. The second Vlasov simulation (VL2) is done with the BOUT++ framework. BOUT++ is a modular platform for 3D simulations of an arbitrary number of fluid equations in curvilinear coordinates using finite-difference methods

[dudson2009bout++, dudson2015bout++]

. Time integration of partial differential equations (PDEs) in BOUT++ is based on the method of lines. The time stepping is performed with the CVODE ODE solver from the SUNDIALS package 

[hindmarsh2005sundials] using variable-order, variable-step multistep methods and is suitable for stiff and nonstiff problems. Spatial derivatives are treated with the third-order weighted essentially non-oscillatory (WENO) scheme for upwind terms and a fourth-order central-difference scheme for other first-order derivatives. In the Vlasov simulations, the velocity boundary conditions are open, and the velocity grid consists of 2001 points. This leads to a velocity resolution of and for the ions and electrons, respectively, where is the ion sound velocity. We start the low-noise Vlasov simulations (VL1 and VL2) with an extremely small initial perturbation ().

Fig. 2 shows the evolution of electrostatic (ES) energy in the low-noise Vlasov simulations (VL1 and VL2), and from here, the linear growth and the nonlinear saturation can be seen. The ES energy in the VL2 simulation starts growing from a larger value than the VL1 simulation. This difference is likely due to the Poisson solver used in the BOUT++ code that introduces some initial noise that is not present in the semi-Lagrangian code. The ES energy in VL2 simulation, however, damps to a value close to the starting energy in VL1 after about 100 ns. This damping leads to some phase difference between the two simulations, so that after 350 ns, the ES energy is higher in the VL2 simulation. The linear growth regimes, which come after about 100 ns in VL1 and 125 ns in VL2, are highly similar in both simulations. For the calculation of linear growth rates, we have chosen four individual modes of the electric field. These modes are in all simulations (PIC and Vlasov) for the case . According to the linear theory, mode has the maximum linear growth rate in our setup. In Fig. 3, the linear growth region is clearly seen for each mode. Table 3 shows the values of the linear growth rates calculated from the low-noise Vlasov simulations and the linear theory are quite consistent with each other. The low SEs reported in Table 3 reflect the fact that the growth is essentially linear. The average SNR of the chosen modes in the linear growth region are 49.13 dB in the VL1 simulation and 17.75 dB in the VL2 simulation. Because the SNR in both simulations is much greater than unity, we can say the power of noise carried in the growth region is quite small.

((a))
((b))
Figure 2: The evolution of the electrostatic energy in the low-noise Vlasov simulations (VL1 and VL2). a) Semi-Lagrangian (VL1) and b) BOUT++ (VL2).
((a))
((b))
Figure 3: a) The evolution of individual modes of the electric field in the low-noise Vlasov simulations (VL1 and VL2). a) semi-Lagrangian (VL1) and b) BOUT++ (VL2). The dashed black line shows the fitted line for the mode.
(Theory)
s
(VL1)
s
SE (VL1)
%
(VL2)
s
SE (VL2)
%
30 0.90 0.9 0.029 0.90 0.014
37 1.08 1.09 0.004 1.08 0.006
44 1.17 1.17 0.017 1.16 0.009
51 1.07 1.08 0.004 1.08 0.006
Average 1.06 1.06 0.014 1.06 0.009
Table 3: Comparison of the theoretical growth rates with the growth rates observed in the VL1 and VL2 simulations.

The PIC simulations are performed with the codes EDIPIC, VSim, and XES1. EDIPIC is a locally developed code that uses the direct-implicit method to integrate the Vlasov–Poisson system of equations in a 1D3V (one spatial dimension and three velocity dimensions) geometry [sydorenko_2006]. VSim is a commercial PIC package that uses the VORPAl computation engine [nieter2004vorpal] to simulate plasmas. In addition, we perform some simulations with XES1 [birdsall2004plasma]. For the simulations PIC1, PIC2, and PIC3, the calculation of the linear growth rates is done using macroparticles per cell. Fig. 4 shows the evolution of the electrostatic energy in random-start PIC simulations. At , the electrostatic energy of PIC simulations is very small because at this time the negative and positive charges are distributed uniformly in the system, so that the system is in a quasi-neutral state. This characteristic is embedded in all PIC simulation codes used in this study, independently of their initialization method. However, after , the ES energy jumps to a finite value. This jump, which is absent in the Vlasov simulations, depends on the initial noise in the velocity space of PIC simulations. Therefore, the ES energy at the second collected time () can be seen as a measure of the initial noise in the simulations. Fig. 4 shows that, relative to the Vlasov simulations (Fig. 2), the ES energy is much larger. This indicates that the amount of initial noise in the PIC simulations is much larger than that of the Vlasov simulations. Fig. 5 shows the evolution of the chosen modes separately. The initial growth in these modes is essentially oscillatory instead of being linear, and therefore, the SE of the growth rate measurements is much larger than unity (Table 4

). We can also define the 99% confidence interval of the measured growth rates as

. The theoretical growth rates can be seen to not lie in the 99% confidence interval of the fits, and thus the measured growth rates cannot be seen as equal to the theoretical growth rates to within the measurement error. We note that the applicability of the 99% confidence interval, for this purpose, is limited to the simulations with significant noise in their linear growth regime (i.e., using the confidence interval is not meaningful in simulations where because in such cases the confidence interval nearly vanishes). The average SNR of the chosen modes in the growth region is dB in PIC1, dB in PIC2, and dB in PIC3. Therefore, the EDIPIC code introduces the least noise power, and XES1 introduces the most noise power in the growth region of the three simulations. We will see that this trend of SNR also applies for the three codes in all other simulations of this study. The small SNR in all three simulations indicates the high noise power in the random-start PIC simulations. To investigate the convergence with respect to the spatial resolution, we repeated the PIC2 simulation with 1024 and 4096 spatial grid points, and the level of noise and the reported results remained close to the original PIC2 simulation. In order to study the effect of only changing the spatial resolution, we note that it is important to not change the number density of macroparticles. The PIC2 simulation was also repeated with a doubled time step size, and again, no significant change was observed in the results.

((a))
((b))
((c))
Figure 4: The evolution of the electrostatic energy in random-start PIC using macroparticles per cell from a) EDIPIC (PIC1) b) VSim (PIC2) c) XES1 (PIC3) simulation.
((a))
((b))
((c))
Figure 5: The evolution of individual modes of the electric field in random-start PIC simulations using macroparticles per cell from a) EDIPIC (PIC1), b) VSim (PIC2), and c) XES1 (PIC3) simulation. The dashed black line shows the fitted line on the mode.
(Theory)
(PIC1)
SE (PIC1)
%
(PIC2)
SE (PIC2)
%
(PIC3)
SE (PIC3)
%
30 0.90 2.79 5.80 2.96 1.83 3.37 1.72
37 1.08 3.68 5.30 2.37 3.46 2.02 4.70
44 1.17 3.00 6.00 3.54 1.46 3.44 1.98
51 1.07 4.03 3.00 2.76 4.14 2.38 2.61
Average 1.06 3.38 5.02 2.91 2.72 2.80 2.75
Table 4: Comparison of the theoretical growth rates with the growth rates observed in the PIC1, PIC2, and PIC3 simulations.

The inaccurate growth rates of the PIC simulations suggest that the noise level in these simulations is so high that it severely influences with the linear growth. Therefore, to reduce the statistical noise level, we increase the number of macroparticles per cell to and redo the PIC simulations (PIC4 and PIC5). The initial electrostatic energy in this case is reduced by an approximate factor of , whereas the initial amplitude of individual modes is reduced by an approximate factor of (compare Fig. 6 and Fig. 7 with Fig. 4 and Fig. 5, respectively). This indicates that the initial noise is reduced approximately by a factor of , as expected. The measured growth rates for the random-start PIC simulations with particle per cell is reported in Table 5. The growth rates of PIC4 simulation with macroparticles per cell are smaller than their counterparts in PIC1 with macroparticles per cell. Accordingly, they are closer to the theoretical growth rates. On the other hand, we see a reduction in spurious oscillation in the linear regime, so that the SEs of the PIC4 simulation are less than those of the PIC1 simulation. In Table 5, we can also see that the average growth rate in the PIC5 simulation is closer to the theory than its corresponding PIC2 simulation. However, in a few modes, such as , we see that the measured growth rate in PIC5 is farther from the theory than it is in PIC2. In both PIC4 and PIC5, the measured linear growth rates are still much larger than the theoretical growth rates. The average SNR of the chosen modes is dB in PIC4 and dB in PIC5.

((a))
((b))
Figure 6: The evolution of the electrostatic energy in random-start PIC using macroparticles per cell from a) EDIPIC (PIC4) b) VSim (PIC5) simulations.
((a))
((b))
Figure 7: The evolution of individual modes of the electric field in random-start PIC simulations using macroparticles per cell from a) EDIPIC (PIC4), b) VSim (PIC5) simulations. The dashed black line shows the fitted line on the mode.
(Theory)
(PIC4)
SE (PIC4)
%
(PIC5)
SE (PIC5)
%
30 0.90 2.45 1.60 2.08 2.23
37 1.08 1.81 1.71 2.56 2.01
44 1.17 3.00 1.90 2.55 2.74
51 1.07 2.69 1.11 3.29 1.17
Average 1.06 2.49 1.58 2.62 2.04
Table 5: The comparison of the theoretical growth rates with the growth rates observed in the PIC4 and PIC5 simulations.

To investigate the problem of inaccurate growth rates in random-start PIC simulations, we introduce a test simulation with the semi-Lagrangian Vlasov code. In this simulation (VL3), we tabulate the initial condition of macroparticles in PIC2 simulation to find the corresponding distribution function and use it as the initial condition for the semi-Lagrangian Vlasov code. By doing this, we introduce the same initial noise as the PIC simulations into the Vlasov simulation. We then repeat the ES energy and mode growth rate analyses using the results of the VL3 simulation (Figs. 7(a) and 7(b)). As with the PIC simulations, we see that the growth of the chosen modes is oscillatory, and the resultant growth rates are much larger than those predicted from theory (Table 6). This strongly suggests that the influence of the initial noise in the PIC simulations is the cause of inaccurate growth rates.

((a))
((b))
Figure 8: a) The evolution of individual modes of the electric field. The dashed black line shows the fitted line on the mode. b) The evolution of the electrostatic energy. Both figures are from the semi-Lagrangian code (VL3) with the initial condition taken from PIC2 simulation.
(Theory)
s
(VL3)
s
SE (VL3)
%
30 0.90 3.00 1.57
37 1.08 2.30 2.71
44 1.17 3.50 1.23
51 1.07 2.17 2.99
Average 1.06 2.74 2.13
Table 6: Comparison of the theoretical growth rates with the growth rates observed in the VL3 simulation.

V The effect of a small flattening of electron distribution function on linear growth rates

In Fig. 9, the coherent structures (holes) in the electron velocity distribution function (VDF) are shown. These structures appear early in the PIC3 simulation (similar structures are observed in other PIC simulations with random-start and VL3) and are a result of trapping of electrons and reflect a small flattening in their Maxwellian velocity distribution function (Fig. 9(a)). This flattening is in fact a depletion of the electrons in the positive velocity region of electron VDF that leads to an increase in electrons in the negative velocity region. To model the flattened velocity distribution function, we add and subtract two shifted Maxwellians (beams) from the initial electron VDF of Eq. 3:

(5)

where is the drift velocity of the added beams, is their thermal velocity, and is their density fraction. To replicate the flattened electron VDF in the simulations, we take , , and . Fig. 9(b) shows this VDF and compares it with the Maxwellian VDF (). Using this VDF, the linear desperation relation reads:

(6)

where are the ion and electron Debye lengths.

Figure 9: The holes in the electron distribution function at ns of simulation PIC3.
((a))
((b))
Figure 10: a) Electron VDF in PIC3 simulation at . b) Electron VDF in Eq. 5 for , , and . For comparison, the Maxwellian VDF is also shown in blue, in each figure.

Solving this dispersion relation, we find the growth rates as shown in Fig. 11. We see that the small flattening in the electron Maxwellian VDF leads to much larger growth rates.

Figure 11: Growth rate from the modified dispersion equation (Eq. 6), with , , and . For comparison, the growth rates of original dispersion relation are also shown.

Vi Using quiet-start initialization to reduce the effect of noise in PIC simulations

In this section, we report on several PIC simulations (PIC6 to PIC12) that use the quiet-start initialization. The quiet-start initialization, proposed by J. A. Byers [byers1970perpendicularly], employs a smooth loading of macroparticles in phase space to reduce the noise in PIC simulation. In this method, the initial placement of macroparticles in - space starts with desired space and velocity densities, and , respectively. The method for generating the positions and velocities of each particle from density functions is based on inversion of the “cumulative density”,

(7)

where is the density function and can be either or

. This cumulative density calculates the cumulative probability in each component

or . can be a uniform set of numbers or a numerical sequence that generates quasi-random numbers with low discrepancy. Several sequences have been proposed in the literature, including the bit-reversed (or Hammersley) sequence[hammersley1964monte, birdsall2004plasma, gonichon1993quiet], Sobol sequence[bratley1988algorithm], and Fibonacci sequence[parker1993fully]. The inversion of the function, by either analytical or numerical means, produces the position or velocity of macroparticles. The set for velocity and position should be uncorrelated to avoid unwanted bunching in phase space. The quiet-start used in our PIC simulations utilizes the bit-reversed set for assigning particle positions and a uniform set of numbers for for particle velocity. This method of quiet-start has been described and implemented in Ref. birdsall2004plasma. In practice, a particular mode is perturbed with a finite amplitude initially.

Fig. 12 shows the growth of individual modes in simulations PIC6, PIC7, and PIC8 using the quiet-start initialization. In these simulations, we have only perturbed the maximum growth rate mode initially. The growth rates measured by these simulations are reported in Table 7. An obvious improvement, in comparison with the corresponding random-start PIC simulations (PIC1, PIC2, and PIC3), is that here the growth rate of the perturbed mode is close to its theoretical value in all three simulations. The average growth rates measured in PIC6 (EDIPIC) and PIC7 (VSim) are also the same as the theoretical growth rates within the 99% confidence interval. However, for the PIC8 simulation (XES1), the theoretical average growth rate is not in the 99% confidence interval of the measured growth rate, and therefore, the two growth rates cannot be considered to be equal by this measure. This discrepancy is due to the inaccuracy of the growth rates for produced in the PIC8 simulation. In particular, the individual mode in all three PIC simulations is far from the theoretical value. The average SE in the PIC6 simulation is improved in comparison with its corresponding random-start simulation (PIC1). In contrast, the average SE of the modes is larger in PIC7 and PIC8 than in the corresponding random-start PIC simulations (PIC2 and PIC3, respectively). This indicates that, in general, the linearity of the growth has deteriorated in PIC7 and PIC8 simulations. The average SNR in the growth rate of chosen modes is dB in PIC6, dB in PIC7, and dB in PIC8. These values of SNR are much smaller than what is reported in section IV for the corresponding random-start PIC simulations. This is likely because of the high-frequency oscillations observed in the growth region of the quiet-start simulations (Fig. 12) carry a large power of noise.

((a))
((b))
((c))
Figure 12: The evolution of individual modes of the electric field in quiet-start PIC simulations, using macroparticles per cell, from a) EDIPIC (PIC6) b) VSim (PIC7) c) XES1 (PIC8) simulation. The dashed black line shows the fitted line on the mode.
(Theory)
(PIC6)
SE (PIC6)
%
(PIC7)
SE (PIC7)
%
(PIC8)
SE (PIC8)
%
30 0.90 0.61 5.61 0.77 8.36 0.66 7.28
37 1.08 1.05 2.52 1.14 2.18 0.61 5.91
44 1.17 1.15 2.70 1.17 2.21 1.18 1.15
51 1.07 1.12 2.62 1.23 1.43 0.62 6.21
Average 1.06 0.98 3.36 1.08 3.55 0.77 5.14
Table 7: The comparison of the theoretical growth rates with the growth rates observed in PIC6, PIC7, and PIC8 simulations.

In the PIC9 simulation, we have perturbed the group of the modes . Table 8 shows that this method of initialization leads to a much smaller average standard error and indicates an improvement in the linearity of the growth compared to the PIC6, PIC7, and PIC8 simulations (see also Fig. 13). This improved linearity, however, leads to a smaller 99% confidence interval, and accordingly, the theoretical average growth rate lies outside the 99% confidence interval of the measurement. Nevertheless, the measured growth rate of PIC9 simulation are much closer to the theoretical values than those from the corresponding random-start simulation PIC2.

(Theory)
(PIC9)
SE (PIC9)
%
30 0.90 1.32 3.63
37 1.08 1.18 1.95
44 1.17 1.26 1.23
51 1.07 1.22 1.08
Average 1.06 1.25 1.97
Table 8: Comparison of the theoretical growth rates with the growth rates observed in the PIC9 simulation.
Figure 13: The evolution of individual modes of the electric field, from PIC9 simulation. The dashed black line shows the fitted line on the mode.

To further investigate the role of the initial perturbation in the quiet-start simulations, we introduce three more PIC simulations (PIC10, PIC11, and PIC12). In these three simulations, the maximum growth rate is not perturbed initially. In the PIC10 simulation, we initialize the simulation with quiet-start but without exciting any mode. In this case, we expect the inherent noise of the PIC simulation to excite the instability. In PIC11 and PIC12, the simulations are initialized with perturbations in the modes and , respectively. Although the individual modes start growing from much lower amplitudes than the random-start PIC simulations, their growth is quite oscillatory at the beginning (see Fig. 14). Also, the measured growth rates are mostly far from the theoretical values. However, in contrast to the random-start PIC simulations, most of the growth rates are underestimated by these simulations (see Table 9). The average SNR of the growth rate is dB in PIC10, dB in PIC11, and dB in PIC12. Therefore, the SNRs are lower than the quiet-start simulations PIC6, PIC7, and PIC8, where the maximum growth rate is perturbed initially.

((a))
((b))
((c))
Figure 14: The evolution of individual modes of the electric field in quiet-start PIC simulations using macroparticles per cell from a) PIC10 b) PIC11 c) PIC12 simulation. All figures are generated from the EDIPIC simulation code. The dashed black line shows the fitted line on the mode.
(Theory)
s
(PIC10)
s
SE (PIC10)
%
(PIC11)
s
SE (PIC11)
%
(PIC12)
s
SE (PIC12)
%
30 0.90 0.62 4.19 0.62 4.14 1.21 2.59
37 1.08 0.59 6.11 0.61 5.67 1.26 1.52
44 1.17 0.34 7.71 0.40 5.95 0.43 6.73
51 1.07 0.63 4.52 0.65 4.88 0.92 3.13
Average 1.06 0.55 5.63 0.57 5.16 0.96 3.49
/
Table 9: The comparison of the theoretical growth rates with the growth rates observed in PIC10, PIC11, and PIC12 simulations.

Vii Conclusion

In this study, we investigated the linear regime of the Buneman instability with several PIC and Vlasov simulations. The different PIC codes show good consistency of their results, and two different implementations of Vlasov simulations are also consistent with each other; the results between the PIC and Vlasov simulations, however, differ significantly. We show that for a relatively small streaming velocity, , the random-start PIC simulations do not reproduce the theoretical linear growth rates, whereas the low-noise Vlasov simulations reproduce them quite accurately. We show that the reason for the discrepancy is the discrete particle noise inherent to PIC simulations. This is demonstrated by initializing Vlasov simulations with the initial conditions of the random-start PIC simulations, which, in the latter case, show a discrepancy similar to PIC results.

This discrepancy is further confirmed by the study of the growth-rate sensitivity to a small flattening in the electron VDF. In Section V, we show that a small flattening significantly increases the linear growth rates. In random-start PIC simulations, the flattening of the electron distribution function occurs as a result of the early trapping of electrons in the noise-driven potential, as can be seen in Fig. 9. In Section III, we show that for large streaming velocity, , the random-start PIC simulations can reproduce the linear growth rates with a reasonable accuracy. This limit is close to the cold-plasma limit ( and ), so that the effects of electron VDF are not important, and the linear growth rates are close to their maximum cold-plasma values.

The noise of the PIC simulations can be further reduced by increasing the number of macroparticles. We show that the results up to macroparticles per cell (PIC4 and PIC5 simulations) will reduce the discrepancy, though the growth rates measured in PIC simulations remain far from the linear theory. Computational resource constraints made it impractical to increase the number of macroparticles per cell much beyond ; in principle, such an increase would increase the accuracy. On the other hand, the Vlasov simulations were able to reproduce the linear growth rates accurately within these constraints.

The effect of the noise resulting from the random sampling of the phase space with a limited number of macroparticles can be partially mitigated by the quiet-start[dawson1983particle]. The accuracy of the growth rates from PIC simulations is greatly increased using the quiet-start initialization method but only if the maximum growth rate mode is perturbed initially. However, the noise is still significant in the linear growth of the quiet-start PIC simulations, and some of the growth rates remain inconsistent with their theoretical values. Therefore, the quiet-start method is not likely to completely remedy the problem of excessive noise in practice.

In the PIC simulations of this study, the quiet-start method uses a bit-reversed sequence in the spatial subspace and a uniform sequence in the velocity subspace[birdsall2004plasma]. These sequences are chosen because of their relative popularity and regularity, which greatly decreases the initial noise level. We also tried a bit-reversed sequence in velocity space, but the resulting growth rates were not as accurate as from the uniform sequence, and therefore, we did not report the results. In practice, other sequences may give different accuracy; however, a systematic comparison of the various proposed sequences is beyond the scope of this study.

The issue of increased noise in PIC simulations is also related to a more general discussion as to what degree the PIC method, which works in between the exact Klimontovich equation and the asymptotic Vlasov equation, describes reality. In the PIC approach, the finite-sized charged clouds may still experience some binary interactions absent in the Vlasov equation but to some degree resembling Coulomb particle collisions [scheiner2019particle, okuda1970collisions, palodhi2019counterstreaming, dupree1963kinetic, dawson1983particle].

For the simulations reported in this study, we have calculated the standard error associated with the measurement of the growth rates. For our purpose, the SE also provides a measure of deviation from linear growth, i.e., the extent to which the observed growth is linear as theory suggests. In the random-start PIC simulations, we observe highly oscillatory growth, and therefore, the SE is the largest for these simulations. However, the inaccuracy of the growth rates is so large that they fall outside the 99% confidence interval of the theoretical growth rates. On the other hand, the growth of the low-noise Vlasov simulation is clearly linear, and therefore, the SE is much less than unity for them. Depending on the simulation code, the SE in the quiet-start PIC simulations can be larger or smaller than the SE in the random-start PIC simulations. Another quantity that we calculate for our simulations is the signal-to-noise ratio (SNR) during the growth of unstable modes. The SNR is largest for the Vlasov simulations, indicating the relatively low power of the noise carried in these simulations. In contrast to the SE, we show that the SNR in the quiet-start PIC simulations is less than that using random-start.

Some modern methods of PIC simulation are proposed to reduce the noise level for a given number of macroparticles. Among these methods, the remapping and the delta-f methods have gained special attention recently. In the remapping method, the microparticles of the PIC simulation are frequently interpolated to a grid in phase space [wang2011particle, myers20174th]. In this way, the noise level is decreased, but the computational cost is increased. In the delta-f method, the known part of the distribution function is separated from its variation (i.e., , where is the known distribution function and is its variation)[aydemir1994unified, brunner1999collisional, allfrey2003revised]. The existing macroparticles are only used to resolve the variation part instead of the distribution function, and therefore, the computational resources are allocated efficiently to reduce noise. However, to keep the noise level small, the condition should be satisfied in the simulation. This condition is usually satisfied in the linear regime of instabilities (as in this study), but it may not quite be satisfied in the nonlinear regime if the distribution function significantly deviates from its initial shape. It is expected that the remapping or delta-f method may ameliorate the noise problem of PIC simulations reported in this study. However, confirming this expectation would require other experiments that are beyond the scope of this study.

Acknowledgments

This work was supported in part by the U.S. Air Force Office of Scientific Research FA9550-21-1-0031, NSERC Canada, and by computational facilities of Compute Canada. A.S. acknowledges illuminating discussions with S. Janhunen.

Author Declarations

Conflict of interest

The authors have no conflicts to disclose.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References