In this paper, we provide insight into the role and impact that a positively biased electrode (anode) has on bulk plasma potential. Using two-dimensional Particle-in-Cell simulations, we investigate the plasma potential as an anode transitions from very small (“probe” mode) to large (“locking” mode). Prior theory provides some guidance on when and how this transition takes place. Initial experimental results are also compared. The simulations demonstrate that as the surface area of the anode is increased transitions in plasma potential and sheath polarity occur, consistent with experimental observations and theoretical predictions. It is expected that understanding this basic plasma behavior will be of interest to basic plasma physics communities, diagnostic developers, and plasma processing devices where control of bulk plasma potential is important.

Measuring or controlling a plasma potential is a fundamental component of most plasma-based phenomena and device operation. The use of Langmuir probes is predicated on the assumption that the probe has negligible effect on bulk plasma properties being measured. From Ref. 1, we have “…probes are quite small and under suitable conditions, produce only minor local perturbations of the plasma.” From Ref. 2, we have “The simple theory of electrical probes is based on the tacit assumption that, although the observed probe currents are dependent upon the acceptance of ions and electrons by the probe, this does not affect the plasma in any way.” And from Ref. 3, we have “… it is postulated that the plasma outside the sheath is unperturbed by the presence of the probe.” Further considerations can also be found in Ref. 4.

As a complement to understanding when a probe is small enough to be considered to have negligible impact on the bulk plasma, one would like to know how an immersed electrode is affecting the bulk plasma. That is, using a relatively small electron-collecting electrode, or anode, to modify bulk plasma potential is clearly of value to modifying, say, the bulk plasma electron energy distribution.5,6

Prior work in Ref. 7 provides a theoretical description for the transition in sheath behavior in the vicinity of an anode in a conceptualized chamber where all other surfaces are grounded. In the limit of a vanishingly small anode, the bulk plasma properties are unchanged due to the presence of the anode. Alternatively, for a sufficiently large anode, the bulk plasma potential will be determined by, or locked to, the anode potential instead of being anchored to the grounded walls. This theory is based on balancing current to the grounded walls and current to the anode through particle flux (loss) terms and does not incorporate a number of other plasma phenomena that may be present in real devices (e.g., global drifts, instabilities, discharge behavior, etc.). For example, the theoretical model assumes uniform sheaths along the full extent of anode and wall surfaces, but in reality there are dielectric interfaces at the boundaries modifying the sheath structure locally, and other complications. Nonetheless, it provides an excellent framework for both computational and experimental testing.

Summarizing from Ref. 7, the global balance of current due to electron and ion losses from the plasma determines the plasma potential with respect to the anode and grounded walls. In the case of very small anodes, the extracted electron current becomes negligible and has essentially no impact on the bulk plasma characteristics. In this case, the anode potential exceeds the plasma potential and the anode-plasma interface is an electron sheath. A larger anode will extract more electrons from the plasma. This must be balanced by a reduction in the electron current lost to the chamber wall, which implies that the plasma potential must rise as the anode size rises. However, as long as the anode potential is sufficiently greater than the plasma potential essentially no ions can reach the anode. Increasing the anode size further, we reach a sufficiently large anode such that the bulk plasma potential will approach the anode potential and ions will reach the anode. This anode size corresponds to the beginning of a transition regime in which the anode-plasma interface is neither a classical electron sheath nor a classical ion sheath. For an even larger anode (and therefore larger bulk plasma potential), significant ions will reach the anode and an ion sheath will form. In this large anode regime, there are now ion sheaths at all surfaces, although a smaller potential difference is present at the anode. This behavior is termed “potential locking” or “attachment” when the bulk plasma potential has locked onto and is determined by the anode potential.

Estimates of the qualitative changes in the nature of the anode-plasma interface as a function of the ratio of anode-to-wall collection surfaces are reproduced here from Ref. 7. Letting AA be the anode area, AW be the wall (grounded cathode) area, and μ=2.3me/mi represent a mass ratio factor (me is the mass of an electron and mi is the mass of the ion),

Ratio rangeExpected interface
AA/AW < μ Electron sheath 
μ ≤ AA/AW < 1.7μ Transition 
1.7μ ≤ AA/AW Ion sheath 
Ratio rangeExpected interface
AA/AW < μ Electron sheath 
μ ≤ AA/AW < 1.7μ Transition 
1.7μ ≤ AA/AW Ion sheath 

This phenomenon has also been investigated experimentally,8 where a number of scenarios are described that produce the transition from “probe” mode to “locking” mode. A number of plasma characteristics were measured that can be used to infer what kind of anode-plasma interface exists (e.g., ion sheath, electron sheath, or something different).

In the experiments, a Gaseous Electronics Conference (GEC) reference cell is filled with helium gas and plasma is generated in a hollow cathode region above the main plasma region through use of an emitting thoriated tungsten filament. Details of the experimental setup are described in Ref. 8. The segmented electrode assembly in Fig. 1 was used to vary the size of the anode at the bottom of the GEC chamber in the experiment by setting the potential for an inner set of rings, and leaving the rest grounded.

FIG. 1.

Anodic surface used in plasma chamber to allow ease of varying anode size.

FIG. 1.

Anodic surface used in plasma chamber to allow ease of varying anode size.

Close modal

A number of techniques often used for interrogating plasma parameters carry the risk of perturbing the plasma (e.g., physical probes). Through the use of simulation, we can avoid this. For example, we can interrogate features of the plasma very close to the anode in simulation results that are extremely difficult (or currently unachievable) to measure in a physical experiment.

The computational study also uses a He plasma for which μ = 0.018 and the area ratio AA/AW thresholds are 1.8% area (transition from electron sheath) and 3.5% area ratio (transition to ion sheath). In the context of device scales, the transition region is quite narrow, occurring over a range of less than 2% of the total surface area. For Ar, the thresholds are even smaller, 0.56% and 1.0% of the total wall area, and the transition regime narrower at ∼0.5%.

Although there were some results from prior one-dimensional simulations,9 qualitatively different mechanisms can take place in two or more dimensions. At a minimum, the basic interpretation of anode area to wall area is very different in one versus higher dimensions. Higher dimensions also allow for more complicated behavior, such as ions preferentially flowing to the sides of the anode surface.

The rest of the paper proceeds as follows. We first introduce the computational model setup, including similarities and differences with the theoretical and experimental systems. We include a discussion of simulation parameters and how they are selected. Next we provide and discuss simulation results. Finally, we conclude with an analysis of the comparisons between computational, theoretical, and experimental results. Overall, this work aims to connect a richer simulation framework to the theory while maintaining consistency with the theoretical and physical descriptions to the extent possible.

Simulations were performed with Aleph, an electrostatic plasma simulation code implementing Particle-In-Cell (PIC),10,11 and Direct Simulation Monte Carlo (DSMC)12 kinetic techniques. The basic algorithm models plasma by tracking computational particles that are coupled through an electrostatic field. Each computational particle represents some number of real ones (e, He, or He+). We also include one interparticle interaction (He-He+ elastic collision). Because this methodology tracks individual particles undergoing specific kinetic events, it is an excellent tool for simulating low temperature non-Maxwellian plasmas. Although the domain is two spatial dimensions, we included three velocity components (a “2d3v” method).

The configuration for the simulations is intended to mimic important aspects of the experiment8 and be as compatible with the theory as possible. Fig. 2 provides a comparison between the two domains. In the experiment, plasma was generated using a thermionic emitter in a source chamber (top). Plasma flowed through a 75 mm diameter hole into a larger chamber in which the experiments were conducted. The experimental features we capture in the simulations include the anode position and dielectric boundary between the anode and the (grounded) wall. The computational model domain is 15 cm wide and 5 cm tall (versus 22.8 cm diameter and 18 cm height for the experiment). We take advantage of centerline symmetry to reduce the computational domain (mesh) to 7.5 cm × 5 cm. The computational domain is initialized with a neutral background of 1 mTorr He at 273 K (nHe is ∼3.54 × 1013 cm−3). To model the experimental setup of plasma flowing from a source region, we generate plasma by identifying a region in the upper portion of the chamber to be continuously repopulated with neutral plasma (Te = 4 eV). This is compatible with the theoretical description (the process maintains charge neutrality), and while not identical to the method used in the experiment, is a simplified model of a physical method to generate plasma. We did not attempt to match the physical plasma generation rate in the experiment's cathode region but instead targeted the steady state plasma densities in the near-anode area, as measured in the experiment (ne is ∼1 × 109 cm−3).

FIG. 2.

Comparison of abstracted experimental geometry to computational model domain. Shaded areas represent plasma generation regions. Solid black surfaces at the bottom are the anodes.

FIG. 2.

Comparison of abstracted experimental geometry to computational model domain. Shaded areas represent plasma generation regions. Solid black surfaces at the bottom are the anodes.

Close modal

To match with theory we assign the 2D anode an “area” (length) to produce the same AA/AW ratio range we expect to be of interest from theory, and is also captured in the experiment (experimental anode-to-wall ratios varied from 9.2 × 10−5 to 3.8 × 10−2). That is, the actual modeled “area” of the anode is different in a quantitatively measured sense from the physical anode area, but we compare all results through the AA/AW ratio (i.e., the “in depth” length in the simulation is the same for all boundaries and so “area” ratios in the theory and experiment are equivalent to “length” ratios in the simulation). The dielectric spacer in the model has a width of 2 mm. We include the dielectric area in the wall area for purposes of calculating AW (as we shall see, the dielectric current is not zero on our time scales).

No other plasma chemistry beyond the He-He+ elastic collision was included in this work. Although other chemistry is certainly present in the experiment, this choice is compatible with the theoretical description, which is entirely independent of plasma chemistry. We used a standard variable hard sphere (VHS) model for the elastic collision.

We discretize the simulation to meet usual PIC/DSMC simulation requirements. A maximum electron density of ne = 2 × 109 cm−3 and electron temperature of Te = 4 eV were assumed. These values were estimated from experimental measurements, and the theoretically predicted trends are generally independent of actual values assuming certain assumptions are met (e.g., TeTi). This combination yields a Debye length λD of ∼332 μm, which we resolve by using a mesh size of Δx = λD/3 or ∼111 μm. Simulations utilized unstructured meshes of ∼161 000 nearly uniform triangular elements. Because much of the domain has a lower plasma density, we are resolving Debye lengths quite well in the bulk. With regards to the DSMC sizing requirements, we note that the only interaction, He-He+ elastic collision, has a mean free path equal to approximately 500 Debye lengths, and so is sufficiently resolved. Proper temporal resolution requires a timestep constrained by the maximum plasma electron frequency (2.5 GHz), Δt < = 793 ps, which we satisfy by setting Δt = 100 ps. A typical simulation reaches steady state in ∼12 h on 128 cores (physical time is 20 μs) and is extended ∼24 h further to collect statistics for a total of ∼36 h simulation time. All of the following simulation results are based on 30 μs averages after the initial 20 μs (where the solution becomes stationary).

The repopulation domain is the shaded rectangular region in the model domain in Fig. 2. It is 13 cm × 0.25 cm and is 1 cm, or ∼30λD, away from the boundary edges. The repopulation condition creates equal numbers of He+ and e at the rate of 2.35 × 1015 cm−3 s−1, leading to a steady density of nHe+ = ne- = ∼2 × 109 cm−3 in the repopulation zone as the particles lost to the walls (anode, grounded wall, and dielectric) are balanced by the creation of new particles. The ions are sampled from a Maxwellian distribution with zero drift velocity and a temperature of 1000 K, and the electrons are sampled from a Maxwellian distribution with zero drift velocity and a temperature of 4 eV. An advantage of using this repopulation model for generating plasma versus alternatives (e.g., injection from a surface) is that particles may return and transport through the region. It mimics some physical plasma generation methods (e.g., ionizing gas in a small region via rf coupling), and certain injection singularities are avoided (e.g., letting the plasma potential float along an injection surface adjacent to a grounded wall will produce numerical instabilities). Another feature this configuration shares with experiment is that plasma is generated in a source region and expands or flows into the plasma domain.

As indicated above, the anode region is held at the anode potential of V = 20 V, the dielectric surface is modeled as a surface charging material, and the remaining walls are all held at ground, V = 0 V.

We simulated anode sizes that extend across the range of area ratios that we expect to yield electron sheaths, a transition region, and ion sheaths. Recalling that AA/AW = 1μ and AA/AW = 1.7μ are the theoretically determined thresholds for transitioning from an electron sheath and to an ion sheath, respectively, the area ratios we simulated are: 0μ, 0.721μ, 1μ, 1.35μ, 1.46μ, 1.58μ, 1.7μ, 1.89μ, 2.08μ, and 2.25 μ (some area ratios are determined by a physical length given to the anode, e.g., a 0.5 cm anode corresponds to AA/AW = 0.721μ).

Fig. 3 shows a typical solution in the electron sheath regime (AA/AW = 0.721μ). Although the simulation was performed in a half-domain (taking advantage of symmetry), the results presented here are mirrored to provide a fuller representation. The line over the anode in the upper plot indicates where the data for most of the following discussion are extracted. For all area ratios, there is a slight bump in potential and plasma density at the population region. The local bump in density occurs because there is a net outflowing flux across each of the population region surfaces, causing an approximately ½ drop in number densities. Assuming the plasma electrons approximately follow a Boltzmann distribution in this region, an increase in plasma potential where there is an increase in plasma density is to be expected.

FIG. 3.

Solution for electron sheath case, AA/AW = 0.721μ. From top to bottom, the panes show He+ number density, e number density, and potential. The number density plots are on a logarithmic scale, and the repopulation zone is quite noticeable. In all three plots, the anode is at the bottom, center. The vertical line above the anode in the upper plot shows where data are extracted for discussion.

FIG. 3.

Solution for electron sheath case, AA/AW = 0.721μ. From top to bottom, the panes show He+ number density, e number density, and potential. The number density plots are on a logarithmic scale, and the repopulation zone is quite noticeable. In all three plots, the anode is at the bottom, center. The vertical line above the anode in the upper plot shows where data are extracted for discussion.

Close modal

A close-up of the anode region is given in Fig. 4. The black segment along the bottom of the plot represents the anode, and the purple segments at the sides along the bottom of the plots show where the surface-charging dielectric surface is located. It is clear that its presence has an impact on local fields, and its influence on the anode behavior is enhanced as the anode size shrinks (while the dielectric remains of constant size). Although this effect can be mitigated (and even removed in a computational model where an anode can be “floated” in space), its effect is also present in the physical system (and any system where the anode is adjacent to dielectric material).

FIG. 4.

Close-up solution for electron sheath case, AA/AW = 0.721μ, near the anode. The anode is depicted in black, and the purple surfaces are the surface charging dielectrics. The influence of the charging dielectrics becomes more prominent as the anode size is reduced.

FIG. 4.

Close-up solution for electron sheath case, AA/AW = 0.721μ, near the anode. The anode is depicted in black, and the purple surfaces are the surface charging dielectrics. The influence of the charging dielectrics becomes more prominent as the anode size is reduced.

Close modal

On the timescales investigated here, there is non-zero current to the dielectric but the relative magnitude is extremely small. We expect the current to approach zero on very large timescales. Fig. 5 shows the currents to anode, wall, and dielectric over 100 μs in an extended time simulation. Letting IA be the current to the anode, IW be the current to the wall, and ID be the current to the dielectric, the ratio of dielectric to anode current is ID/IA = 1.6 × 10−4 (averaging from 90 to 100 μs), or less than 0.02%.

FIG. 5.

Currents to the anode, wall, and dielectric surfaces for electron sheath case, AA/AW = 0.721μ. Note the dielectric current is indistinguishable from 0 at this scale (see text for more detail). A continuous averaging window of 1 μs has been applied. This total simulation time was 100 μs instead of the usual 50 μs.

FIG. 5.

Currents to the anode, wall, and dielectric surfaces for electron sheath case, AA/AW = 0.721μ. Note the dielectric current is indistinguishable from 0 at this scale (see text for more detail). A continuous averaging window of 1 μs has been applied. This total simulation time was 100 μs instead of the usual 50 μs.

Close modal

Although we are primarily interested in the steady state behavior in this work, there are clearly interesting transient phenomena present in Fig. 5. The temporal oscillations in the currents are greatest for small anodes where an electron sheath is present. In Ref. 13, these oscillations were shown to be associated with an electron-ion drift instability that arises as the electrons flow toward the anode. In contrast, the currents in Fig. 6 for a larger anode in the ion sheath regime (AA/AW = 2.25μ) lack similar oscillations. The fluctuations in this situation are due to the statistical noise from computational particles.

FIG. 6.

Currents to the anode, wall, and dielectric surfaces for ion sheath case, AA/AW = 2.25μ. Note the dielectric current is indistinguishable from 0 at this scale (see text for more detail). A continuous averaging window of 1 μs has been applied.

FIG. 6.

Currents to the anode, wall, and dielectric surfaces for ion sheath case, AA/AW = 2.25μ. Note the dielectric current is indistinguishable from 0 at this scale (see text for more detail). A continuous averaging window of 1 μs has been applied.

Close modal

The potentials along the first 3 cm of the centerline for all simulations are shown in Fig. 7. From Fig. 7 we see that the bulk potential converges to a single profile as AA/AW approaches and then exceeds 1.70μ. This behavior is what we refer to as potential “locking.” As demonstrated in Ref. 8, at this anode size the bulk plasma potential would increase with increased anode potential. This is not the case for smaller anodes (e.g., AA/AW ≪ 1).

FIG. 7.

Potential profiles at the centerline above the anode for all simulated area ratios.

FIG. 7.

Potential profiles at the centerline above the anode for all simulated area ratios.

Close modal

Figs. 8 and 9 show the electron and ion number densities above the anode. Although the transition from electron sheath to ion sheath can be seen here, we present three canonical sheath structures in the following Figs. 10–12. These correspond to an electron sheath (0.721μ), a transition interface (1.35μ), and an ion sheath (2.25μ). The electron sheath and ion sheath are as expected, but the transition sheath contains some interesting properties.

FIG. 8.

Electron number densities along the centerline above the anode for all simulated area ratios.

FIG. 8.

Electron number densities along the centerline above the anode for all simulated area ratios.

Close modal
FIG. 9.

Ion number densities along the centerline above the anode for all simulated area ratios.

FIG. 9.

Ion number densities along the centerline above the anode for all simulated area ratios.

Close modal
FIG. 10.

Ion density, electron density, and potential near the anode for a small anode case.

FIG. 10.

Ion density, electron density, and potential near the anode for a small anode case.

Close modal
FIG. 11.

Ion density, electron density, and potential near the anode for the intermediate case.

FIG. 11.

Ion density, electron density, and potential near the anode for the intermediate case.

Close modal
FIG. 12.

Ion density, electron density, and potential near the anode for the large anode case.

FIG. 12.

Ion density, electron density, and potential near the anode for the large anode case.

Close modal

These results show that the potential profile is flat over most of the bulk plasma. The smallest anode exhibits a strong electron sheath and no ions reach the anode. The bulk plasma potential is ∼8 V, significantly below the anode potential. The largest anode clearly shows an ion sheath structure, and the bulk plasma is ∼1 V above the anode. In this case, the bulk plasma has locked onto the anode potential. The intermediate case includes features from both electron and ion sheaths. The bulk potential is ∼18 V, still less than the anode potential, a feature usually associated with electron sheaths, but some ions do reach the anode, a feature usually associated with ion sheaths. This can be due to oscillations as shown in Fig. 5 for the electron sheath,13 laterally flowing ions influenced by the dielectric region, and/or collisional effects. In the smallest anode case, where ions are created at ∼12 V, no ions reach the anode. There is likely some influence of the dielectric-to-anode and dielectric-to-wall jumps in the electric field on the exact collection behavior of the different sized anodes (reducing with larger anodes). However, there is clearly a transition from electron sheath to ion sheath through a less easily described intermediate regime at the transition ratio predicted by theory.

A comparison between experiment and simulation is provided in Fig. 13 where the bulk plasma potentials are given. Similar trends are present, but there is a ∼5 V discrepancy between them. It is unclear what causes this discrepancy, but as discussed in Ref. 8, there may be a better interpretation of anode size from the experiment and/or simulation that is more consistent with the theoretical notion that it is a current-collecting surface, i.e., AA should be compared to the surface area of the “capture region” above the anode and not the anode itself. If this results in a larger AA for experiments, that curve would shift to the right giving a better comparison. Indeed, even dimensionality may play a role. The simulation has only 2 “edges” whereas the experiment has a full three-dimensional structure. There are also plausible explanations for vertical shifts in the data, for example, differing bulk Te, but additional interpretations are not explored here. It should be noted that uniform vertical shifts in the results would not change the area ratio locations of certain features, i.e., transition from intermediate to ion sheath structure. Nevertheless, we see good agreement between theory and the simulation, and the correct qualitative behavior in the experimental results.

FIG. 13.

Comparison of bulk plasma potential as a function of anode size from simulation and experiment for anode potential of 20 V. The simulation values were taken at a position 1.5 cm above the center of the anode. Theory predicts an electron-to-intermediate transition near 1μ, and an intermediate-to-ion transition at 1.7μ.

FIG. 13.

Comparison of bulk plasma potential as a function of anode size from simulation and experiment for anode potential of 20 V. The simulation values were taken at a position 1.5 cm above the center of the anode. Theory predicts an electron-to-intermediate transition near 1μ, and an intermediate-to-ion transition at 1.7μ.

Close modal

We have demonstrated and compared the behavior of electron-to-ion sheath transitions due to anode size from a theoretical model, an experimental model, and a computational model. The qualitative behavior is similar amongst all three models, and the quantitative transitions of the theoretical and computational models are in strong agreement. In terms of area ratios, the threshold area ratios occur within 3.5% of the wall area (1.0% for Ar), an extremely tight range, and this is reproduced in the computational model.

The implication of this result is that a size estimate can be made for supporting the assumption that an electron-capturing biased probe has negligible behavior on the bulk plasma (a negligible current sink). If the probe area is assumed to be “small” (e.g., negligible impact on global plasma characteristics), a comparison can be made to the grounded wall region in terms of surface area ratios (a more general treatment could apply to more complex arrangements of current-collecting boundaries). If the area ratio is approximately 1μ or greater care should be employed to ensure the probe is not globally impacting the measured plasma.

Furthermore, if one hopes to use an appropriately sized anode to influence the global plasma in the intermediate transition regime between electron sheaths and ion sheaths, there is only a small range of surface areas that are effective, approximately 2% of the total wall area (0.5% for Ar). The range of effective surface areas is even smaller if one wants to exercise control at the upper end of the transition region where the plasma potential is becoming locked to the anode. Put another way, using an anode to increase the global plasma potential slightly above the potential determined by a completely grounded domain is easier than decreasing the global plasma potential slightly below the potential determined by the anode (the limit where the global plasma potential is locked to the anode potential).

Future work will involve improvements to the computational model, in addition to extensions of the theoretical and experimental models. A key aspect of this will be closely examining the ion and electron velocity distribution functions in all three investigational methods. We will also enrich the scope of plasma chemistry in the computational and theoretical models to incorporate electron-neutral interactions (e.g., ionization and excitation).

This work was supported by the Office of Fusion Energy Science at the U.S. Department of Energy under Contract No. DE-AC04-94SL85000.

1.
M. A.
Lieberman
and
A. J.
Lichtenberg
,
Principles of Plasma Discharges and Materials Processing
, 2nd ed. (
Wiley
,
2005
).
2.
S.
Glasstone
and
R. H.
Lovberg
,
Controlled Thermonuclear Reactions Huntington
(
Van Nostrand Reinhold Co.
,
NY
,
1960
).
3.
J. E.
Allen
,
R. L. F.
Boyd
, and
P.
Reynolds
, “
The Collection of Positive Ions by a Probe Immersed in a Plasma
,”
Proc. Phys. Soc. B
70
,
297
304
(
1957
).
4.
B. E.
Cherrington
, “
The use of electrostatic probes for plasma diagnostics—A review
,”
Plasma Chem. Plasma Process.
2
(
2
),
113
140
(
1982
).
5.
C.-S.
Yip
,
J. P.
Sheehan
,
N.
Hershkowitz
, and
G.
Severn
, “
Mackenzie's Demon with instabilities
,”
Plasma Sources Sci. Technol.
22
,
065002
(
2013
).
6.
C.-S.
Yip
and
N.
Hershkowitz
, “
Manipulating the electron distribution through a combination of electron injection and MacKenzie's Maxwell Demon
,”
Plasma Sources Sci. Technol.
24
,
034004
(
2015
).
7.
S. D.
Baalrud
,
N.
Hershkowitz
, and
B.
Longmier
, “
Global nonambipolar flow: Plasma confinement where all electrons are lost to one boundary and all positive ions to another boundary
,”
Phys. Plasmas
14
,
042109
(
2007
).
8.
E. V.
Barnat
,
G. R.
Laity
, and
S. D.
Baalrud
, “
Response of the plasma to the size of an anode electrode biased near the plasma
,”
Phys. Plasmas
21
,
103512
(
2014
).
9.
S. D.
Baalrud
,
T.
Lafleur
,
R. W.
Boswell
, and
C.
Charles
, “
Particle-in-cell simulations of a current-free double layer
,”
Phys. Plasmas
18
,
063502
(
2011
).
10.
C. K.
Birdsall
and
A. B.
Langdon
,
Plasma Physics via Computer Simulation
(
Taylor & Francis
,
New York
,
1991
).
11.
R. W.
Hockney
and
J. W.
Eastwood
,
Computer Simulation Using Particles
(
CRC Press
,
1989
).
12.
G. A.
Bird
,
Molecular Gas Dynamics and the Direct Simulation of Gas Flows
(
Clarendon Press
,
1998
).
13.
B.
Scheiner
,
S. D.
Baalrud
,
B. T.
Yee
,
M. M.
Hopkins
, and
E. V.
Barnat
, “
Theory of the electron sheath and presheath
,”
Phys. Plasmas
22
,
123520
(
2015
).