# Kinetic modeling of three-dimensional electrostatic-solitary and surface
waves in beam neutralization

Nakul Nuwal#<EMAIL_ADDRESS>Deborah A. Levin# Igor D. Kaganovich†
#University of Illinois Urbana-Champaign, IL, USA, †Princeton Plasma Physics
Lab, NJ, USA

###### Abstract

This work studies the fundamental plasma processes involved in the
neutralization of an ion beam’s space-charge by electrons emitted by a
filament using Particle-in-Cell simulations. While filament neutralization is
economical, previous experiments have shown that a variety of waves become
excited in this process that limit the space-charge neutralization. In this
work, the formation and movement of electrostatic solitary waves (ESWs), which
have low dissipation rates, are characterized for 2D planar and 3D cylindrical
beams and are observed to generate waves that survive for a long time and slow
the process of beam neutralization. Further, through a 1D Bernstein-Greene-
Kruskal (BGK) analysis, we find that the non-Maxwellian nature of the beam
electrons gives rise to large-sized ESWs that are not predicted by theory
which assumes that the electrons may be described by a Maxwellian
distribution. Our PIC simulations are sufficiently sensitive to be able to
resolve important three-dimensional effects in a 3D cylindrical geometry that
lead to the excitation of Trivelpiece-Gould surface waves due to high energy
electrons present at the beginning of neutralization.

##  I Introduction

Ion beams are used in various engineering applications such as particle
accelerators, ion-thrusterskorkut2015 ; jambunathan2018chaos ;
nuwal2020kinetic , nanopantographychen2019nearly , and ion-
implantationconrad1989plasma . The electrons are introduced to neutralize the
beam and reduce its divergence along the axis. The source of electrons can be
from a background plasma or through an external source such as filaments or
cathodes. Although better neutralization is achieved using the background cold
electron plasmavay2001intense ; vay1998charge ; kaganovich2001nonlinear ;
berdanier2015intense , this method is expensive as it requires the production
of plasma in the vicinity of the ion beam. Filaments, on the other hand, are
easier to install and control, however, experiments have found multiple
challenges in their use in the ion beam neutralization. For example, the rate
of neutralization is affected by the energy and location of the electron
source relative to the ion beamlan2020neutralization ; kaganovich2009physics .
Additionally, in recent numerical works by Lan and
Kaganovichlan2020neutralization ; lan2020neutralization2 , electrostatic-
solitary-waves (ESWs), also known as electron solitons, were observed to form
when electrons were injected into a 2D ion beam via a filament source. This is
problematic since a solitary wave is a non-linear wave structure which
sustains itself for a long time with low dissipationbernstein1957exact , thus
delaying neutralization.

An ESW forms when enough electrons are trapped in a local electric potential-
well. ESWs can be multidimensional structures, and their movement along the
beam depends on the electric potential profilelan2020neutralization2  and
electron velocity distribution of the beam. ESWs were first described in one-
dimension by Bernstein, Green, and Kruskalbernstein1957exact  where a
spatially symmetric electric potential structure was shown to be stable and
stationary when it trapped a finite fraction of electrons. A number of
theoretical solutions of ESW or electron holes have been found by other
authors, such as Sagdeev, and Schamelschamel1979theory , where the Vlasov-
Poisson equation system was solved for the electric potential and electron
velocity phase distribution, assuming a stationary solitary wave.
Schamelschamel1979theory  showed an analytical solution of a finite (small)
amplitude slow moving electron hole by solving the Vlasov-Poisson system of
equations using the Classical-potential methodschamel1979theory ;
bernstein1957exact ; hutchinson2017electron . Since their theoretical
prediction, these waves have been observed in experimentssaeki1979formation
and in the geomagnetic tailmatsumoto1994electrostatic , which is comprised of
charged particles. While such theoretical analysis provided insights into the
properties of individual solitary waves in a Maxwellian far-field electron
distribution, we will show how they appear with a different shape and size in
a neutralizing beam where the electron distribution is non-Maxwellian.
Further, we will show how the ESWs in a 3D cylindrical beam are different from
a 2D planar beam configuration and demonstrate why the 1D BGK analysis is
inadequate due to more complex electron trajectories in an unmagnetized plasma
beam.

Figure 1: Schematic of physical features that appear along the beam axis
during beam neutralization, where $A_{1},A_{2}$, and $A_{3}$ are the
amplitudes of the electron holes that form near the electron source, solitary
waves, and surface wave, respectively, and $A_{1}>A_{2}>A_{3}$.

ESWs affect the beam neutralization because they create local pockets of
positive potential, as shown in the Fig. 1 schematic, thereby preventing the
beam from achieving a perfect neutralizationlan2020neutralization2 . Figure 1
shows a schematic of the electric potential profile along the axis of a
typical beam in which an electron source is placed along its axis. As will be
shown later, the continuous emission of electrons from the electron source
results in the formation of a virtual cathode and two electrons holes on its
either sides, which can be seen in the electron potential-wells labeled by
amplitude $A_{1}$ in Fig. 1. Further downstream of the electron source, ESWs
appear in the beam, shown by the waves with the amplitude $A_{2}$ in Fig. 1,
and they move along the beam with speeds ranging from zero to
0.92$v_{te}$hutchinson2017electron , where $v_{te}$ is the electron thermal
speed. The difference between an electron hole and an ESW is that an electron
hole does not move while an ESW does. As will be shown later, many of these
wave features depend on the shape and geometry of the beam. Since for a number
of applicationsconrad1989plasma  a beam is cylindrical in shape, performing a
numerical study for a 3D beam helps us to understand the effects of geometry
on solitary and other plasma waves. In this work, we aim to understand the
beam neutralization in 2D planar and 3D cylindrical beams and the formation of
ESWs, and their movement along the beam axis using the Particle-in-Cell (PIC)
method.

Depending on the geometry of the plasma, small disturbances at the plasma-
vacuum interface can result in the excitation of surface waves, as was
theorized by Trivelpiece and Gouldtrivelpiece1959space . These waves were
studied because of their ability to interact with the electromagnetic waves of
certain frequencies that cause interference with the passive measurement
methods used in plasma beam diagnosticskrall1973principles . The amplitude of
these waves, designated as $A_{3}$ in Fig. 1, are lower than ESWs. These waves
are especially likely to appear in systems such as ours, where the electrons
neutralize a high-electric-potential beam and the average electron temperature
of the beam declines with time. Such surface waves, if excited, will not be
damped by high energy electrons. As will be shown later, we found Trivelpiece-
Gouldtrivelpiece1959space  surface waves in our 3D plasma beam simulations.
However, such waves did not appear in the 2D planar beam case because a
higher-phase speed is required to excite them in a planar beam geometry, as
will be discussed further in Sec. IV.

In this paper, we first discuss the numerical approach of the PIC calculations
in Sec. II. Then, in Sec. III, we describe the excitation, movement, and
structure of the ESWs that form along the beam. In Sec. IV, we discuss the
excitation and propagation of surface waves in the beam, where we show a
comparison of our PIC results with the theoretical disperson relations for 2D
and 3D beams. Finally, in Sec. V, we list the key findings and conclusions of
our numerical study.

(a) 2D planar beam

(b) 3D cylindrical beam

Figure 2: Schematic (not drawn to scale) of the (a) 2D and (b) 3D numerical
domains. Here, the channel width is $2b$ and the ion-source width or diameter
for 2D planar and 3D cylindrical beams is $2a$.  Table 1: Numerical case
parameters for 2D and 3D beams1-7

Numerical parameter | 2D planar | 3D cylindrical
---|---|---
|  $a=2.5$ mm |  $a=2.5$ mm |  $5.0$ mm |  $7.5$ mm
$\Delta x$ (mm) | 0.46 | 0.46 | 0.46 | 0.46
Number of GPUs (Nvidia Tesla-K20) | 128 | 128 | 128 | 128
Particles-per-cell | 400-650 | 450-800 | 250-500 | 200-420
$F_{\rm num}$ | 40 | 30 | 60 | 80
Particles per species (Millions) | 33 | 46 | 92 | 154
No. of PIC grid cells (Millions) | 0.9 | 3.7 | 3.7 | 3.7

  * •

1 $\omega_{pe}=7.46\times 10^{8}$ rad/s, $n_{\rm 0}$ = $1.75\times 10^{14}\
{\rm m^{-3}}$, and $\Delta t=80\times 10^{-12}$ s.

  * •

2 $n_{\rm Ar^{+}}$ = $1.75\times 10^{14}\ {\rm m^{-3}}$ at the ion source;
electron current, $I_{e-}=I_{\rm Ar^{+}}/3$, $I_{\rm Ar^{+}}/3$ being the ion
current.

  * •

3 Half-width ion source in 2D planar beam case, $a=2.5$ mm.

  * •

4 Radius of ion source, $r_{\rm Ar^{+}}\equiv a=2.5$, $5.0$, and $7.5$ mm for
different 3D cylindrical cases.

  * •

5 $b=30$ mm for the 2D planar and $b=15$ mm for all 3D cylindrical cases.

  * •

5 $T_{e}=2$ eV at the electron source for all cases.

  * •

6 All ions are monoenergetic and introduced with velocity, $v_{z}=427$ km/s
and no thermal velocity.

  * •

7 Debye length, $\lambda_{\rm D_{0}}=0.78$ mm, at plasma density $1.75\times
10^{14}\ {\rm m^{-3}}$ and $T_{e}=2$ eV.

##  II Numerical approach and general plasma conditions

The numerical study in this paper uses our in-house time-explicit CPU-GPU
based electrostatic Particle-in-Cell (PIC) solver CHAOS (Cuda-based Hybrid
Approach for Octree Simulations)Jambunathan_CnF . We solve Poisson’s equation
to compute the electric potential scalar and electric fields to move charged
particles in the domain. To satisfy the numerical validity of the time-
explicit nature of the coupling between Poisson’s equation and particle
movement, we use a time-step, $\Delta t<0.1\omega_{pe}^{-1}$ and grid cell-
size, $\Delta x<\lambda_{\rm D}$godunov1959difference , where $\omega_{pe}$
and $\lambda_{\rm D}$ are the local electron plasma frequency and Debye
length, respectively. The code CHAOS has previously been used in plume-
plasmajambunathan2018chaos ; nuwal2020kinetic , and plasma-surface
interactionnuwal2021influence  studies, and provides the multi-GPU scalability
necessary for this study.

In this work, we performed calculations of 2D planar and 3D cylindrical beams
in a metal cavity. For the 2D planar beam calculations, we used a numerical
domain of $0.375\times 6\times 40$ cm, where $2b=6$ cm is the channel width,
as shown in Fig. 2(a). The size of our ion source and the domain length is of
the same order as other plasma wave experimentslynov1979observations ;
ikezi1971electron ; stepanov2016dynamics  and nanopantographychen2019nearly
applications. Since CHAOS is a 3D code, the third dimension of the domain was
suppressed for the 2D case by not moving particles in the $x$ direction,
applying a zero potential gradient boundary along the $x$ direction
($\partial\phi/\partial x=0$), and setting the $x$ direction thickness to be
$0.375$ cm, which is much thinner than the other two dimensions of the
simulation domain. All other parameters are given in Tab. 1. We also performed
three cases for 3D cylindrical beams, each with a different ion beam radius.
For all 3D cylindrical beam cases, a domain of size $1.5\times 1.5\times 40$
cm, shown in Fig. 2(b), was used with a uniform grid of cell size, $\Delta
x=0.46$ mm. The domain width was reduced from $2b=$ 6 cm in the 2D planar beam
to $2b=$ 3 cm in the 3D cylindrical beam case to reduce the computational
expense while keeping the same domain length. All walls of the 3D domain were
considered charge-absorbing with $\phi=0$ V, similar to other
numericallan2020neutralization ; lan2020neutralization2  and experimental
workstepanov2016dynamics ; lynov1979observations , where Poisson’s equation
was solved with Dirichlet boundary conditions at all domain walls. The
numerical parameters for all 3D cases are also listed in Tab. 1.

In both 2D and 3D cases, ions were introduced in the domain with a constant
axial velocity of 427 km/slan2020neutralization ; stepanov2016dynamics  and
zero initial velocity in the transverse direction from an ion source of width
(diameter for 3D cases) $2a=5$ mm at $z=0$ mm, as shown in the Fig. 2
schematic. As discussed before and also shown in Lan et
al.lan2019neutralization , the beam diverges along the z-axis because of self-
repulsive space-charge among ions. Electron injection starts when the ions
reach the end of domain, $z=40$ cm, at $t=1\ \mu s$. Electrons in both 2D
planar and 3D cylindrical beams were injected from the axial location of
$z=10$ cm with a Maxwellian velocity distribution function of electron
temperature, $T_{e}=2$ eV, in all three directions. Similar to Lan et
al.lan2020neutralization , the injected electron current, $I_{e^{-}}$, was one
third of the ion current $I_{\rm Ar^{+}}$ and the electrons were introduced
from a rectangular source of width 2 mm for the 2D planar beam. For all the
cases of a 3D cylindrical beam, the electrons were introduced from a circular
source of radius $r_{e^{-}}$, such that the ratio of ion and electron source
radii, $a/r_{e^{-}}=2.5$, which is consistent with the ion to electron source
width ratio for the 2D planar beam case. Unlike the 2D case, the charged
particles were allowed to move in all three directions for the 3D beam cases.
The number of particles-per-cell ranges from 200 to 800 for all cases, and a
convergence study of the 3D case with $a=5.0$ mm is presented in Appendix A.

In general, regardless of the specific geometry, there are two phases of beam
neutralizationlan2020neutralization : (1) when the rate of neutralization of
the beam, $\partial\phi/\partial t$, is high, and (2) when the beam has
reached a quasi-steady state value and $\partial\phi/\partial t$ is very
small. In all our cases, the first phase lasts from $t=1$ to $3.5\ \mu s$ and
the second phase starts at $t=4.0\ \mu s$, as seen in Fig. 3(a). In a
neutralizing beam, ESWs form in the first phase and then dissipate very slowly
in the second phaselan2020neutralization2  of neutralization. In this work, we
focus only on the first phase of neutralization. As more electrons are
introduced, the beam electric potential decreases, and the electron density
profile in the beam transverse direction widens and assumes a shape that is
similar to the ion density profile, as shown by its evolution in Fig. 3(b).
The electron density profiles at different times, shown in Fig. 3(b), also
compare well with the results of Lan et al.lan2020neutralization  with less
than 1% difference. Because of the conducting walls biased at $\phi=0$ at
$z_{\rm min}$ and $z_{\rm max}$, ion plasma sheaths form on both these
boundaries, as shown in Fig. 3(c). These plasma sheaths reflect low energy
electrons and result in the formation of two electron streams in the beam,
which create plasma waves, such as ESWsomura1996electron .

(a) Electric potential at $z=20$ cm comparison with Lan et
al.lan2019neutralization .

(b) Electron and ion density profiles at $z=20$ cm compared at $t=2$ and $4.3\
\mu s$

(c) Evolution of electric potential profile along the beam axis with time.

Figure 3: Time evolution of (a) probe electric potential, (b) electron density
profiles, and (c) axial electric potential profile for a 2D planar beam.

Since we will be discussing differences in ESWs and surface wave formations
for 2D vs 3D geometries, we provide a comparison of the evolution of their
electric potential and ion density profiles with time in the first phase of
neutralization. Figure 4(a) shows that the normalized electric potential and
ion density profiles have a similar shape for both 2D planar and 3D
cylindrical beams at $t=1.2\ \mu s$, right after the electron injection starts
at $t=1.0\ \mu s$. The electric potential in Fig. 4(a) is normalized by
$\phi_{\rm max}=220$ and 14 V for 2D planar and 3D cylindrical beams,
respectively. Because of higher electric potential in the 2D planar case, the
ion beam is initially more expanded along the beam axis compared to the 3D
cylindrical case at $t=1.2\ \mu s$, as shown in Fig. 4(a). However, with time
the ion beam becomes more coloumnated in the 2D planar beam case, compared to
the 3D cylindrical beam case, as shown in Fig. 4(b), where the ion density
value drops by factors of about 1.1 and 2.82 from $z=0$ cm to $z=30$ cm in the
2D planar and 3D cylindrical beam cases, respectively at $t=4.0\ \mu s$. This
comparison shows that the beam-neutralization process is slower in the 3D
cylindrical beam case than the 2D planar beam case because of the lower beam
potential of the initially un-neutralized beam. The electric potential ‘bumps’
in Fig. 4(b) provide clues about the presence of electrostatic solitary waves
(ESWs) which are excited in both 2D and 3D cases with $a=2.5$ mm. In Fig.
4(b), the amplitude of the ESWs near the end of the first phase of
neutralization is higher for the 2D than 3D case, likely because of the higher
electric potential of the initially un-neutralized 2D planar beam, which
enables it to trap more electrons than the 3D beam.

(a) $t=1.2\ \mu s$.

(b) $t=4.0\ \mu s$.

Figure 4:  Comparison of time evolution of plasma parameters along the beam
axis for 2D and 3D cases. Here, $n_{0}=1.75\times 10^{14}\ {\rm m^{-3}}$ and
$\phi_{\rm max}=220$ and 14 V for 2D and 3D cases, respectively. In (b),
$\Delta\phi$ is the amplitude of the ESW, i.e. the difference between the peak
and base potential.

##  III Excitation of ESWs in beam neutralization

During the neutralization of a beam, a fraction of electrons are trapped in
regions of the beam which results in the formation of non-linear waves called
ESWs, electron holes, or solitonshutchinson2017electron ; bernstein1957exact ,
as mentioned in the previous section. Figure 5(a) shows a time-snapshot of the
axial velocity phase space for electrons at $t=3.12\ \mu s$ where multiple
regions of trapped electrons, i.e. ESWs, appear along the beam. These ESWs
always appear at the local maxima in the electric potential, as shown in Fig.
5(a), and their spatial location can be estimated by numerically
differentiating the electric potential profile (see Appendix A). The $z$
locations of these local maxima are shown by the red bubbles in Fig. 5(a) for
$t=3.12\ \mu s$. The location of these ESWs extracted over multiple discrete
time-steps is shown in Fig. 5(b) where their positions, shown by the red
bubbles, is overlayed on the electric potential variation in time along the
beam axis.

(a) Electron phase distribution, $f_{e-}$, is computed by normalizing number
density with $n_{0}=1.75\times 10^{14}\ {\rm m^{-3}}$, and the electric
potential profile at $t=3.12\ {\rm\mu s}$. The black line shows electric
potential on the right vertical axis. The red dots show the position of the
ESW along the dashed black line in (b).

(b) Movement of ESWs along the beam axis with time overlayed on $\phi(z)$ vs
$t$.

Figure 5:  Time evolution of ESW positions along the beam axis,
$(x,y)=(3.0,3.0)$ cm, for the 2D planar beam case, $a=2.5$ mm and $b=30$ mm.
In (b), the ESW shown by blue ‘X’s merge with the large ESW near the electron
source at $z=10$ cm, and ESWs shown by white and green ‘X’s merge with each
other. The ESW shown by yellow ‘X’ moves with a constant speed from about
$t=2.75$ to $3.75\ \mu s$.

###  III.1 Movement of ESWs along beam axis

Figure 5(b) shows the movement, collisions, and merging of various ESWs with
time and, as will be shown later, can be used to estimate the speed of an ESW
to determine its properties, such as the trapped electron particle population
and their energies. Figure 5(b) also shows a gradual decrease in electron
potential in the beam from about 120 V at $t=2.5\ \mu s$ to about 5 V at
$t=4.5\ \mu s$. In Fig. 5(b), we can see finger-like structures in the $\phi$
contours on the $z$ $vs$ $t$ plot for $t>4\ \mu s$ which coincide with the red
dots. These structures represent the slow moving ESWs that survive in the beam
at the end of the first phase of beam neutralization. These structures then
take a long time to dissipate in the second phase of the beam neutralization,
as was shown in Lan et al.lan2020neutralization2 .  The continuous nature of
these points shows that ESWs travel along the beam with velocities ranging
from nearly zero to $20\ \rm cm/\mu s$ with one of the fastest ESWs shown by
blue ‘X’s in Fig. 5(b). These speeds of ESWs are in the same range of zero to
$0.92\sqrt{(2k_{b}T_{e})/m_{e}}$ shown in Hutchinsonhutchinson2017electron ,
where $k_{b}$ and $m_{e}$ are Boltzmann constant and the mass of electron,
respectively. Since the slope of the curves passing through these ESW
locations gives their velocities, the lines that are nearly parallel to the
$t$ axis show ESWs that are stationary, whereas the lines that are ‘flatter’
represent ESWs moving with very high speeds in Fig. 5(b).

Similar to Lan and Kaganovichlan2020neutralization2 , these ESWs merge with
the other ESWs, as shown by the merging trajectories of a number of ESWs in
Fig. 5(b). An example of this merging can be seen by following two ESWs marked
in blue ‘X’ in Figs. 5(b) and 6(a) from $t=3.04$ to $3.96\ \mu s$, in which
both the ESWs move towards the electron source at $z=10$ cm and finally merge
with the stationary electron hole at $z\approx 11$ cm at $t=3.36\ \mu$s. Both
the ESWs retain their approximate shape, i.e. electric potential and length,
before merging, as shown in Fig. 6(a). Similarly, the ESWs marked in white and
green ‘X’ in Figs. 5(b) and 6(b) collide and merge between $t=4.0\ \mu s$ and
$4.4\ \mu s$. As seen in Fig. 6(b), the smaller ESW, shown by green ‘X’ tries
to pass through the larger ESW, shown by white ‘X’. However, because of a low
relative velocity between these ESWs, they merge together to now move at a
speed lower than either of those individual ESWs, as seen in Fig. 5(b), where
the merged ESW, marked by white ‘X’ at $t=4.40\ \mu s$, has a much smaller
velocity. In both cases, the merging of ESWs results in an ESW of larger
length to accommodate a larger number of trapped particles (as discussed
further in Sec. III.2). Such inelastic collisions between ESWs have also been
previously observed in experiments where they merge when their relative speed
is smallsaeki1979formation . The electron hole with the largest length that
forms next to the electron source does not significantly move and remains
attached to the electron source and therefore cannot be considered an ESW.

(a) Snapshots of ESWs for $3.04<t<3.36\ {\rm\mu s}$.

(b) Snapshots of ESWs for $3.84<t<4.40\ {\rm\mu s}$.

Figure 6:  Movement of ESWs in 2D planar beam case, represented by the X marks
in Fig. 5(b). In (a), two blue ESWs move in the $-z$ direction and merge with
the attached electron hole located at $z\approx 12$ cm. In (b), the two ESWs,
shown by white and green X marks, merge into a single ESW. Here, the electron
distribution $f_{e-}$ is computed by normalizing number density with
$n_{0}=1.75\times 10^{14}\ {\rm m^{-3}}$.

In the 3D cylindrical beam case with beam radius $a=2.5$ mm and channel width
$b=15$ mm, the movement of ESW is shown in Fig. 7. Similar to the 2D planar
beam case discussed above, the electron hole with the largest length, whose
location is shown by the magenta dashed line in the time-snapshot of $t=3.44\
\mu s$ in Fig. 7(b), remains stationary and attached to the electron source at
$z=10$ cm, as shown by the nearly vertical line of magenta color created by
connecting the electron hole locations near $z\approx 12$ cm in Fig. 7(a).
Also, similar to the 2D planar beam, we observe collisions and merging of ESWs
in the 3D cylindrical beam case. An example of a collision can be seen in Fig.
7(b), where the ESW shown by a yellow ‘X’ collides with the attached ESW at
$z\approx 12$ cm and initiates the movement of the previously stationary ESW
at $z=8$ cm. Similarly, the ESW shown by the white ‘X’ in Fig. 7(b) merges
with the attached electron hole at $z\approx 12$ cm. Finally, the two small
ESWs, marked by green ‘X’ in Figs. 7(a) and 7(b), collide, pass-through, and
finally merge with each other. For both 2D planar and 3D cylindrical cases
with $a=2.5$ mm, we find that most ESWs seem to move in the $-z$ direction
with about the same range of ESW speeds between zero and $21$ cm/$\mu s$.

(a) Time evolution of electric potential along beam axis and movement of ESWs.

(b) Time snapshots of ESWs for $3.44<t<4.64\ {\rm\mu s}$, shown by dashed
white box in (a).

Figure 7:  Time evolution of electric potential and electron ESWs along the
beam axis ($(x,y)=(1.5,1.5)$ cm) for 3D cylindrical beam with beam radius,
$a=2.5$ mm. In (a), the magenta colored line near $z=12$ cm shows the electron
hole attached to the electron source at $z=10$ cm. In (b), the yellow and
white ESWs move and merge with the attached electron hole, shown by the dashed
magenta line in phase snapshot of $t=3.44\ \mu s$. The ESWs shown in green
collide with each other and merge. In (b), the electron distribution,
$f_{e-}$, is normalized by $n_{0}=1.75\times 10^{14}\ {\rm m^{-3}}$.

The nature of ESWs, however, dramatically changes for a 3D cylindrical beam
case with a larger beam radius, $a=5.0$ mm and the same channel width, $b=15$
mm, as the previous cylindrical beam case. For this case, the widest ESW forms
in the upstream direction of the electron source because of a higher electric
potential near a larger ion source at $z=0$ cm. This ESW oscillates in the
region $0<z<10$ cm with an amplitude of about 8 cm, as shown in Fig. 8(a), and
this amplitude of oscillation decays with time as the ESW shrinks in shape and
size. Unlike the cases shown before, the widest electron hole for the $a=5.0$
mm case is detached from the electron source. While an attached electron hole
forms downstream of the electron source for the $a=5.0$ mm case, it is much
smaller in size compared to the upstream ESW. A similar ESW movement is seen
for the case with $a=7.5$ mm where we observe a detached ESW oscillating
upstream of the electron source, as shown in Fig. 8(b). Also, the upstream ESW
for the $a=7.5$ mm case is of larger length than the $a=5.0$ mm case because
of the larger number of ions near $z=0$ mm, which attract more beam electrons
in the upstream region of the electron source. The ESWs formed downstream of
the electron source, for both the $a=5.0,$ and $7.5$ cm cases, are initially
very small in size, and they soon disappear with time. This disparity in the
formation of ESWs for different cases shows that the size and position of the
ESWs formed in a 3D beam have a strong dependence on the beam radius and ion
and electron currents. The summary of locations of ESWs for different cases is
listed in Tab. 2.  The effect of trapped electron population on the shape and
size of an ESW in 2D planar and 3D cylindrical beams will be discussed next.

Table 2: Position of ESWs and properties of surfaces waves in all PIC cases

Case no. | Type of beam |  $a^{1}$ (mm) | ESWs positions2 | Surface waves $\left(\frac{\omega}{\omega_{pe}},\kappa a\right)$
---|---|---|---|---
1. | 2D planar | 2.5 mm | All downstream | -
2. | 3D cylindrical | 2.5 mm | All downstream | -
3. | 3D cylindrical | 5.0 mm | One upstream, none downstream | (0.33,0.51)
4. | 3D cylindrical | 7.5 mm | One upstream, none downstream | (0.36,0.79)

  * •

1 $a=$beam half-width/radius for 2D planar/3D cylindrical beam.

  * •

2 $z$ positions w.r.t. $e^{-}$ source, which lies at $z=10$ cm, along the beam
axis.

(a) $a=5.0$ mm.

(b) $a=7.5$ mm.

Figure 8:  Time evolution of electric potential along the beam axis
($(x,y)=(1.5,1.5)$ cm) for 3D cylindrical beam cases with $a=5.0$ and $7.5$
mm.

###  III.2 Shape of ESWs in a beam: trapped vs untrapped electrons

Characterizing ESWs based on their amplitude and size is important in
understanding their role in the neutralization of beams. While ESWs were found
in the numerical study of a 2D beam by Lan et al.lan2020neutralization2 , a
quantitative comparison of their amplitude and size with the BGK theory was
not discussed. In this section, we isolate one ESW each from 2D planar and 3D
cylindrical beams of $a=2.5$ mm, and compare them with the 1D BGK
theorybernstein1957exact . As discussed in the previous section, ESWs move
along the beam as independent entities with a low dissipation, and they have a
population of trapped and untrapped electrons associated with them. In the
work of Bernstein-Greene-Kruskalbernstein1957exact  (BGK), an analytical
relation between electric potential profile and electron velocity distribution
was derived for a stationary solitary wave. In
Hutchinsonhutchinson2017electron , it was shown that the analysis by BGK
applies to a moving ESW as well if the total energy and electron velocity
distribution is evaluated in the frame of the moving ESW. Similar to the
symmetric electron holes of BGK and Hutchinson, the electric potential
profiles of our ESWs are almost symmetric. The electric potential profile of
an ESW traps the low kinetic energy electrons, which are reflected from its
edges and keeps these electrons bouncing in the potential well. These trapped
electrons appear as a circular vortex-like structure in phase space, as was
shown in Figs. 6(a), 6(b) and 7(b).

Figure 9(a) shows a contour plot of the average total electron energy, $E_{\rm
tot}$, per electron for the 2D planar beam case at $t=3.12\ \mu s$,

$\displaystyle E_{\rm tot}(v_{z},z)=\frac{\sum_{i=1}^{N_{\rm
total}}\left[\frac{1}{2}m_{e}|(\vec{v_{e_{i}}}-\vec{v_{s}})|^{2}-e\phi\right]\delta_{v_{e_{i}}v_{z}}\delta_{z_{e_{i}}z}}{\sum_{i=1}^{N_{\rm
total}}\delta_{v_{e_{i}}v_{z}}\delta_{z_{e_{i}}z}},$ (1)

where $\vec{v_{e_{i}}}=v_{e_{iz}}$ and $\vec{v_{s}}$ are the z-directional
electron and ESW velocities, $\delta$ is the Kronecker delta function,
$z_{e_{i}}$ is the $z$ position of the $i^{th}$ electron, and $N_{\rm total}$
is the total number of electron computational particles in the domain. For the
phase-space points $(v_{z},z)$ where there are no electrons, $E_{\rm tot}=0$.
For trapped electrons, the term in square brackets of Eq. 1 is less than zero.
The average total energy per particle, $E_{\rm tot}=0$ eV, shown by the dashed
blue line in Fig. 9(a), marks the approximate boundary in phase space between
the trapped and un-trapped electrons. The discontinuities in the contour plot
of Fig. 9(a), seen at the edges of ESWs at $z=13,$ 19.4, 20.9, 23.9, 25.8, and
28.7 cm, are because of the difference in values of $v_{s}$ for each ESW, as
estimated from Fig. 5(b). Based on the average total energy criteria, the
fractions of trapped ($N_{\rm t}$) and untrapped ($N_{\rm ut}$) computational
electron particles for each of the ESWs numbered in Fig. 9(a) are listed in
Tab. 3. It can be seen that the length ($\lambda_{\rm ESW}$) and amplitude
($\Delta\phi$) of an ESW increases with a higher fraction of trapped
electrons, which is consistent with a 1D BGK modebernstein1957exact . Among
the four ESWs shown in Fig. 9(a), ESW ‘2’ is of interest because it is
relatively isolated from the other ESWs and it moves with a constant speed for
a long duration (see yellow ‘X’ in Fig. 5(b)). The RHS of Fig. 9(a) shows the
electron velocity distribution function (EVDF) extracted at the center of ESW
‘2’ ($z=$ 22.5 cm) in the phase plot on the LHS for the 2D planar beam case.
For ESW ‘2’, the un-trapped particle region of the EVDF has a peak at
$v_{z}/v_{0}\approx 0.5$, which corresponds to the electrons moving with a
velocity close to the ion beam velocity of $v_{z}/v_{0}=0.427$, where
$v_{0}=1\times 10^{6}$ m/s. The trapped region of the EVDF in Fig. 9(a) is
symmetric about the ESW speed, $v_{s}=$ $14.5$ $\rm{cm/\mu s}$. This symmetry
is also consistent with a typical 1D BGK mode where continuous bouncing in the
potential-well produces a symmetric distribution of trapped
electronsbernstein1957exact .

(a) Average total energy of electrons in a 2D planar beam. The RHS figure is
the EVDF at the center of the ESW ‘2’, $z=22.50$ cm, at $t=3.12\ \mu s$.

(b) ESW ‘2’ in a 2D planar beam at $t=3.12\ \mu s$.

(c) Trapped electron distribution from CHAOS (dashed) vs Eq. 2 (solid) for 2D
beam at $t=3.12\ \mu s$.

Figure 9:  Average electron energy distribution for different ESWs in a 2D
planar beam at $t=3.12\ \mu s$. In (b), the electric potential profile and
normalized trapped electron density in phase space are shown. The dashed black
line in the bottom part of (b) are iso-contour lines of $E_{\rm tot}$, with
their values labeled in white boxes. In (c), a comparison of predicted EVDF of
trapped population based on Eq. 2, shown by solid lines and CHAOS PIC result,
shown by dashed lines. Here, $f_{t}$ and $f_{e-}$ are normalized by
$n_{0}=1.75\times 10^{14}\ {\rm m^{-3}}$.

Furthermore, in a 1D BGK mode, all the electrons follow constant energy
trajectoriesbernstein1957exact ; hutchinson2017electron . Therefore, if the
electric potential profile of an ESW’s potential-well, and the trapped
particle EVDF at the center of the ESW are known, the trapped EVDF at
different locations inside the ESW can be analytically
evaluatedbernstein1957exact ; hutchinson2017electron  by,

$\displaystyle f((v_{z}-v_{s}),\Delta
z)=f_{0}\left(\sqrt{0.5(v_{z}-v_{s})^{2}-\frac{e}{m_{e}}(\Delta\phi(\Delta
z)-[\Delta\phi]_{\rm max}})\right),$ (2)

where $\Delta\phi=\phi-\phi_{\rm base}$, $\phi_{\rm base}$ being the lowest
electric potential value in the vicinity of the ESW, $\Delta z=z-z_{\rm
center}$, $z_{\rm center}$ being the center of the ESW, and $f_{0}$ is the
EVDF at the center of the ESW, shown on the RHS of Fig. 9(a), and
$\Delta\phi(z=z_{\rm center})=[\Delta\phi]_{\rm max}$hutchinson2017electron .
For our 2D planar beam, we obtain a 1D electric potential profile, from PIC,
along the beam axis where the value of $\Delta\phi$ ranges from
$[\Delta\phi]_{\rm max}$ at the center of the ESW to 0 V at its edge, as shown
in the top part of Fig. 9(b). Figure 9(b) also shows the distribution of
trapped particles in the phase space (bottom) with iso-energy lines of $E_{\rm
tot}$ for the ESW, labeled ‘2’ in Fig. 9(a). Using the EVDF at the center of
ESW ‘2’, EVDFs at three different locations, shown in the top part of Fig.
9(b), are evaluated using Eq. 2. The EVDFs from Eq. 2, shown by solid lines in
Fig. 9(c), compare well with the PIC result from CHAOS, shown by dashed lines
of the same colors in Fig. 9(c). The minor discrepancy is because of the error
in estimating the ESW velocity, $v_{s}$, from Fig. 5(b). The predicted
untrapped EVDFs using Eq. 2 also compared well with the PIC result. This
comparison shows that the ESWs we see in the 2D beam neutralization case have
electron trajectories that follow the 1D BGK mode type electron distribution
bernstein1957exact .

(a) Electron density from PIC result for 2D planar beam at $t=3.12\ \mu s$.

(b) Fit of $\Delta\phi(\Delta z)$ with the PIC result

(c) Far-field EVDF, $f_{\rm FF}$, compared for Maxwellian and non-Maxwellian
cases.

Figure 10:  (a) Normalized electron density with perturbations at the
locations of numbered ESWs from Fig. 9(a). Comparison of assumed and
CHAOS(PIC) $\Delta\phi$ profile (b) and far-field EVDFs for the ESW ‘2’ (c).

(a) Untrapped particle population for $f_{\rm FF}$ as a non-Maxwellian EVDF
obtained from PIC results.

(b) Untrapped particle population for $f_{\rm FF}$ as a Maxwellian EVDF.

Figure 11:  Comparison of normalized $f_{ut}$ (i.e. $f_{\rm ut}/f_{\rm max}$
here) for the cases of a Maxwellian vs non-Maxwellian far-field distribution.
Here, $f_{\rm max}$ is the maximum of $f_{ut}(v_{z},\Delta z)$, and the dashed
red line represents the velocity of ESW ‘2’, $v_{s}=14.5\ {\rm cm/\mu s}$.

We also observe that the ESWs listed in Tab. 3 are about 30-45 Debye lengths
long for the 2D planar beam case. While the typical ESWs analyzed in theory
are about 2$\lambda_{\rm D}$ to 20$\lambda_{\rm D}$ in length and small in
amplitudehutchinson2017electron ; chen2002bgk2 ; lefebvre2010laboratory , in
our case, the ESWs move with low relative speeds, which allows them to merge
with each other to create longer ESWs with larger trapped particle
populations. This merging was also shown in Lan et al.lan2020neutralization2 ,
and it was identified as one of the properties of electron holes in
Schamelschamel1979theory . Such long ESWs usually have a large electron
density perturbation when they form in a Maxwellian far-field background
plasma. However, Fig. 10(a) shows that all the ESWs in the 2D planar beam PIC
case have surprisingly low electron density perturbations between $3$ to
$5\%$. The electron density profile in Fig. 10(a) is computed by sampling all
the electrons at multiple z planes along the beam axis at $t=3.12\ \mu s$.
With further study, we find that the non-Maxwellian EVDF, found in our beam,
requires a smaller electron density perturbation to form an ESW of the same
length than one where the EVDF is Maxwellian in the far-field. We demonstrate
this by using the 1D BGK integral approach, described in Bernstein et
al.bernstein1957exact  and Hutchinsonhutchinson2017electron , to predict the
untrapped and trapped electron density profiles inside ESW ‘2’. In this
approach, the background or far-field EVDF, $f_{\rm FF}$, of an ESW is known
(in our case from the PIC simulations), and the ions are assumed to be of
constant and uniform density in the ESW. This approach has the following four
steps: (1) assume an electric potential profile for the ESW, (2) using
conservation of total electron energy with the assumed potential profile of
the ESW, calculate the untrapped EVDF, $f_{ut}(v_{z},\Delta z)$, from the
known far-field EVDF, $f_{\rm FF}$, (3) compute the untrapped electron density
($n_{ut}(\Delta z)$) of the ESW from $f_{ut}(v_{z},\Delta z)$, and, (4) using
Poisson equation, calculate the trapped electron density ($n_{t}(\Delta z)$).
See Appendix B for details regarding how this approach is applied to the PIC
particle data.

With respect to step (1), we fit the electric potential profile of ESW ‘2’
from the PIC simulation to an analytical electric potential profile suggested
by Schamelschamel1979theory ,

$\displaystyle\Delta\phi={\phi_{0}}{\rm sech^{4}}\left(\frac{b_{s}\Delta
z}{\lambda_{\rm D_{0}}}\right),$ (3)

as shown in Fig. 10(b), with $\phi_{0}=0.876$ V and $b_{s}=0.05$ as the fit
parameters, and $\lambda_{\rm D_{0}}=0.78$ mm is the Debye length near the
electron source where $n_{e}=1.75\times 10^{14}\ {\rm m^{-3}}$ and $T_{e}=2$
eV. This analytically fit electric potential profile makes it easier for us to
apply the 1D BGK approach. Since the total energy (i.e. kinetic + potential
energy) is conserved along the electron trajectories in 1D, we can compute the
untrapped electron phase distribution, $f_{ut}\left(v_{z},\Delta z\right)$,
from a known far-field EVDF, $f_{\rm FF}$, i.e. step (2),

$\displaystyle f_{ut}\left(v_{z},\Delta z\right)=\begin{cases}f_{\rm
FF}\left(\sqrt{v_{z}^{2}-2\frac{e}{m_{e}}\Delta\phi(\Delta z)}\right)&\text{if
$v_{z}>0$}\\\ f_{\rm
FF}\left(-\sqrt{v_{z}^{2}-2\frac{e}{m_{e}}\Delta\phi(\Delta
z)}\right)&\text{if $v_{z}<0$}.\end{cases}$ (4)

(a) The untrapped and total electron density for two cases of $f_{\rm FF}$.

(b) Comparison of trapped and untrapped populations predicted by Eq. 7 and
obtained from PIC results.

Figure 12:  Comparison of untrapped and total electron density for a
Maxwellian vs non-Maxwellian far-field distribution(PIC) using 1D BGK
approach. In (b), the PIC result is compared with the 1D BGK approach and the
densities are normalized by the electron density in the far-field of ESW ‘2’,
i.e. $1.5\times 10^{14}\ {\rm m^{-3}}$.

The two electron streams in the beam results in a non-Maxwellian far-field
EVDF, $f_{\rm FF}$, as shown in Fig. 10(c), where the highest peak represents
the electrons moving in the direction (+z) of the ion beam and the smaller
peak represents the electrons reflected from the plasma sheaths at $z_{\rm
min}$ and $z_{\rm max}$ boundaries and the virtual cathode that forms near the
electron source at $z=10$ mm. The electron phase distribution,
$f_{ut}(v_{z},\Delta z)$, computed using Eq. 4 with the non-Maxwellian far-
field EVDF, $f_{\rm FF}$, of Fig. 10(c), is shown in Fig. 11(a). Figure 11(a)
also shows that, similar to the far-field EVDF, $f_{\rm FF}$, the
$f_{ut}(v_{z},\Delta z)$ predicted by Eq. 4 is asymmetric about the ESW’s
velocity, $v_{s}=-14.5\ {\rm cm/\mu s}$. In contrast, if we assume a
Maxwellian EVDF for the far-field,

$\displaystyle f_{\rm
FF}(v_{z})=C\exp\left[{-\frac{v_{z}^{2}}{v_{te}^{2}}}\right],$ (5)

where $C=350$ is a fit parameter based on the comparison with the PIC result
shown in Fig. 10(c) and $v_{te}$ is the thermal speed, Eq. 4 gives an
untrapped electron phase distribution, $f_{ut}(v_{z},\Delta z)$ that is
symmetric about the ESW ‘2’ velocity, $v_{s}=-14.5\ {\rm cm/\mu s}$, as shown
in Fig. 11(b). In step (3), we use the untrapped electron phase distribution,
$f_{ut}(v_{z},\Delta z)$, to calculate the number density profile of untrapped
electrons, $n_{ut}(\Delta z)$, given by,

$\displaystyle n_{ut}(\Delta z)=\int_{-\infty}^{+\infty}f_{ut}(v_{z},\Delta
z)dv_{z}.$ (6)

A comparison between untrapped electron density profiles, $n_{ut}(\Delta z)$,
for a Maxwellian vs non-Maxwellian (PIC result) far-field EVDFs is shown in
Fig. 12(a), where a Maxwellian $f_{\rm FF}$ predicts a lower untrapped
electron density than the non-Maxwellian $f_{\rm FF}$, obtained from PIC,
inside the ESW. In the usual 1D BGK approach, we invoke Poisson equation to
calculate the trapped electron density profile ($n_{t}(\Delta z)$) inside the
ESW, i.e.,

$\displaystyle n_{t}(\Delta z)$ $\displaystyle=n_{i}(\Delta
z)+\frac{\epsilon_{0}}{e}\frac{d^{2}[\Delta\phi(\Delta z)]}{d[\Delta
z]^{2}}-n_{ut}(\Delta z),$ (7)

where $n_{i}(\Delta z)$ is the ion number density profile, which is assumed to
be constant and $n_{i}=n_{ut}(\Delta z\to\infty)$. For the particular case of
ESW ‘2’, the second term on the RHS of Eq. 7 would be much smaller than the
other two because of its large length, which also makes the trapped and
untrapped electron density profiles nearly complementary and their difference
is equivalent to $4\%$ electron density perturbation in the ESW. Finally,
using Eq. 7, we calculate the trapped electron density profile using the non-
Maxwellian $f_{\rm FF}$ that we obtained from our PIC simulation result, i.e.
step (4). The normalized $n_{t}$ and $n_{ut}$ profiles for the non-Maxwellian
$f_{\rm FF}$ compare well with the PIC result in Fig. 12(b), which reaffirms
our assertion that ESW ‘2’ in the 2D planar beam case is a 1D BGK mode formed
in a non-Maxwellian far-field plasma.

However, using Eq. 7 for a Maxwellian $f_{\rm FF}$ may result in a fictitious
trapped EVDF where it can have negative values for the number of
electronshutchinson2017electron . To avoid that, we instead calculate the
trapped electron phase distribution by assuming that the trapped EVDF at the
center of the ESW ($\Delta z=0$) is known from our PIC result. This treatment
of trapped electrons results in identical trapped electron density profiles,
$n_{t}(\Delta z)$, for both Maxwellian and non-Maxwellian $f_{\rm FF}$’s. So,
the only difference in the total electron density is due to the untrapped
electrons inside the ESW. Figure 12(a) shows a comparison of the total
electron density drop in the Maxwellian and non-Maxwellian far-field cases,
which shows that the Maxwellian $f_{\rm FF}$ would have a larger electron
density perturbation, given that the trapped electron distribution at the
center of the ESW and the potential profile is the same as it is for the non-
Maxwellian $f_{\rm FF}$. This means that a Maxwellian background EVDF, $f_{\rm
FF}$, would likely result in a larger electron density and electric potential
perturbation than a non-Maxwellian one and this is the reason why such long
ESWs with small density and potential perturbations appear in our ion beam,
which has a non-Maxwellian $f_{\rm FF}$.

Table 3: Trapped and un-trapped computational particles of each ESW in the 2D
planar beam at $t=3.12\ \mu s$, and 3D cylindrical beam at $t=3.44\ \mu s$,
shown in Figs. 9(a) and 14, respectively, with $a=2.5$ mm.

| ESW No. | $N_{\rm t}^{\dagger}/N_{0}^{\\#}$ | $N_{\rm ut}^{\dagger}/N_{0}$ |  $\Delta\phi^{*}$ (V) |  $[\lambda_{\rm D}]_{\rm FF}^{d}$ (mm) | $\lambda_{\rm ESW}^{c}/[\lambda_{\rm D}]_{\rm FF}$ | $v_{s}\ (\rm{cm/\mu s})$
---|---|---|---|---|---|---|---
| 1. | 0.14 | 0.86 | 0.45 | 0.74 | 29.47 | 13.42
2D planar | 2. | 0.25 | 0.75 | 0.87 | 0.69 | 39.88 | 14.50
| 3. | 0.23 | 0.77 | 0.80 | 0.69 | 26.66 | 20.67
| 4. | 0.30 | 0.70 | 1.19 | 0.69 | 43.88 | 18.08
3D cylindrical | 1. | 0.049 | 0.951 | 0.52 | 0.59 | 73.57 | 15.96
| 2. | 0.081 | 0.919 | 1.04 | 0.63 | 114.2 | 21.00

  * •

† $N_{\rm t},$ $N_{\rm ut}$ are number of trapped and untrapped computational
electron particles, respectively.

  * •

# Total electron particles, $N_{0}=N_{\rm t}+N_{\rm ut}$.

  * •

∗ $\Delta\phi=$ amplitude of the ESW.

  * •

c $\lambda_{\rm ESW}=$ length of ESW based on $E_{\rm tot}=0$ eV.

  * •

d $[\lambda_{\rm D}]_{\rm FF}$ is local Debye length in the far-field of the
ESW.

(a) Iso potential surfaces of a 3D beam with $a=2.5$ mm.

(b) 2D vs 3D beam potential profiles at $z=20$ cm.

Figure 13: Iso-potential surfaces of 3D beam in (a) show the oblong shaped
electron holes next to the electron source and solitary waves at its
downstream. In (b) the transverse direction $\phi$ profiles are compared for
2D planar and 3D cylindrical beams at $t=3.12$ and $3.44\ \mu s$,
respectively.  Figure 14:  Trapped and untrapped number of particles in
($r,z,v_{z}$) phase space based on the average total energy per particle,
$E_{\rm tot}$, of electrons. For trapped particles, $E_{\rm tot}<$0 and
$v_{0}=1\times 10^{6}$ m/s. In the inset images, the white space signifies the
absence of computational particles in the phase space. In the RHS inset
figure, ‘A’ and ‘B’ indicate the difference in electron particle densities in
the two opposing streams, which indicates an asymmetric non-Maxwellian EVDF.
Note that the number of computational particles are divided into the 800 bins
in $z$ direction, 30 bins in $r$ direction and 60 bins in $v_{z}$ direction.

For the 3D cylindrical beam case witth $a=2.5$ mm, however, both the electron
holes near the electron source and the ESWs (i.e., solitary waves marked as
‘1’ and ‘2’) have 3D oblong shapes that are elongated along the beam axis, as
shown in Fig. 13(a). Similar to the 2D planar beam, the amplitude
($\Delta\phi$) of ESWs in the 3D beam is directly correlated with the fraction
of electrons that are trapped. However, the electric potential structure is
very different for the ESWs in the 3D cylindrical beam. Figure 13(a) shows
iso-potential surfaces in 3D where the beam has been cut in half along the
$x=1.5$ cm mid-plane. This electric potential is different from the 2D planar
beam case where the electric potential profile was wide and flat near the beam
axis (see the 2D curve in Fig. 13(b)), thus enabling a 1D BGK analysis in that
case. However, for the 3D cylindrical beam case, the electric potential
rapidly declines with radial distance, as shown in Fig. 13(b). Since the
electron source was placed along the axis of the ion beam, the beam structure
in the 3D cylindrical beam case is perfectly axisymmetric, as shown by the
comparison of the electric potential profiles at $z=20$ cm along the $x$ and
$y$ transverse directions in Fig. 13(b). Therefore the ESWs in the 3D
cylindrical beam case can be analyzed in a $(r,z)$ cylindrical coordinates
system.

Similar to the 2D planar beam analysis discussed above, we can evaluate the
trapped and untrapped electron particle populations along the beam by
calculating the total energy of the electrons using all three velocity
components. The average total electron energy, $E_{\rm tot}$, for an
axisymmetric ESW in the 3D beam is given by,

$\displaystyle E_{\rm tot}(r,z,v_{z})=\frac{\sum_{i=1}^{N_{\rm
total}}\left[\frac{1}{2}m_{e}|(\vec{v_{e_{i}}}-\vec{v_{s}})|^{2}-e\phi\right]\delta_{v_{e_{i}}v_{z}}\delta_{z_{e_{i}}z}\delta_{r_{e_{i}}r}}{\sum_{i=1}^{N_{\rm
total}}\delta_{v_{e_{i}}v_{z}}\delta_{z_{e_{i}}z}\delta_{r_{e_{i}}r}},$ (8)
$\displaystyle\text{where
}|(\vec{v_{e_{i}}}-\vec{v_{s}})|=\sqrt{v_{e_{ix}}^{2}+v_{e_{iy}}^{2}+(v_{e_{iz}}-v_{s})^{2}},$

and $\delta$ is the Kroeneker delta function and $r_{e_{i}}$ is the position
of the $i^{th}$ electron in the domain.

Figure 14 shows the average total energy of electrons, $E_{\rm tot}$, per
electron particle in $(r,z,v_{z})$ phase space where the average electron
energy is negative in the regions where an ESW is present, i.e. electrons are
trapped. In the inset of Fig. 14, the trapped and untrapped particle
populations are shown in the same $(r,z,v_{z})$ space based on the total
electron energy, i.e. term in square brackets of Eq. 8, being negative and
positive, respectively for the ESW labeled as ‘2’. In the untrapped population
shown in the RHS inset of Fig. 14, we can see that a larger population of
electrons have velocities $v_{z}>0$ (shown by ‘A’) than $v_{z}<0$ (shown by
‘B’), because of the non-Maxwellian $z$-direction electron distribution,
similar to the 2D planar beam EVDF that is shown in Fig. 9(a).  The ESW ‘2’ in
Fig. 14 for the 3D beam is about $73\lambda_{\rm D}$ in length, which is
slightly longer than the ESWs in 2D beam, as listed in Tab. 3. The ESW ‘1’ in
Fig. 14 has a larger length than ESW ‘2’ because of the ongoing merging
between the two ESWs shown in Fig. 7(b) at $t=3.44\ \mu s$.

(a) EVDF at far-field region and at the center of ESW 2.

(b) Trapped and untrapped electrons in ESW 2.

Figure 15: (a) Electron velocity distribution of ESW 2 from PIC result and (b)
trapped $vs$ untrapped electron population along the ESW. In (a), there is a
non-zero population of untrapped electrons between $v_{z}=-1\times 10^{6}$ and
$5\times 10^{5}$ m/s as these electrons have high $v_{x}$ and $v_{y}$ values.
In (a), the EVDF is normalized by $n_{0}=1.75\times 10^{14}\ {\rm m^{-3}}$. In
(b), the densities are normalized by the electron density in the far-field of
ESW ‘2’ of 3D beam, i.e. $1.4\times 10^{14}\ {\rm m^{-3}}$.

Applying the same 1D analysis to the electrons in the 3D-cylindrical beam, we
use the non-Maxwellian far-field distribution, $f_{\rm FF}$, for the
electrons, shown in Fig. 15(a), obtained from the PIC result of the $a=2.5$ mm
3D beam case. The $f_{\rm FF}$ (far-field) shown in Fig. 15(a) was obtained by
binning all the electrons that lie in the $z=20$ cm plane, i.e. outside the
ESW ‘2’ of Fig. 13(a). Using the electric potential profile along the beam
axis, we obtain the untrapped electron population profile, $n_{ut}(\Delta z)$,
using Eq. 4, as shown in Fig. 15(b), which can be seen to not agree with the
PIC result. The trapped electron population also shows disagreement, and its
value is over-predicted by the 1D BGK analysis when compared to the PIC
result, as shown in Fig. 15(b). This difference is due to the fact that the 1D
BGK analysis makes the assumption that for all the electron trajectories in
the ESW, the electric potential profile along the trajectory path of an
electron is the same as it is along the beam axis of the 3D cylindrical beam.
However, even though the electric potential of ESW ‘2’ is axisymmetric, the
electron trajectories are quite complicated, as can be seen in Fig. 16. For
example, in Fig. 16(a), the trajectories of two electrons show that they cross
several different iso-potential lines when trapped in an axisymmetric
potential-well and are not parallel to the beam axis. The same is true for the
untrapped electrons, shown by three example trajectories in Fig. 16(b), where
all three labeled electrons cross multiple iso-potential lines. In Fig. 16(b),
the electrons labeled ‘1’, and ‘3’ are reflected from the virtual cathode
which forms at the electron source at $z=10$ cm while electron ‘2’ has enough
kinetic energy to eventually leave from the $z=40$ cm boundary. Note that
since these trajectories in Fig. 16 are taken in a static potential field at
$t=3.44\ \mu s$, they represent the constant total energy lines of electrons
in the 3D beam at that time. Other approaches such as that of Chen et
al.chen2002bgk  and Eliasson et al.eliasson2007theory , which take into
account the radial changes in electric potential profiles, also assume that
the electron trajectories are parallel to the beam axis, which, as we show in
Fig. 16, is not true for a non-magnetized beam. For electrostatic plasma
cases, such as ours, a fully 3D-3V BGK approach would have to be used to
account for the 3D effects in a beam. Although spherical electrostatic ESWs
have been studied theoretically in Ng et al.ng2005bernstein  for a canonical
case with a Maxwellian far-field EVDF, a general theory for 3D is yet to be
developed to understand the ESWs that are seen in our 3D beam PIC simulations.

(a) Trajectory of an electron trapped in ESW 2.

(b) Trajectories of electrons not trapped in ESW 2

Figure 16: Examples of electron trajectories which are (a) trapped and (b)
untrapped in ESW 2 (15 $<z<$ 20 cm) superimposed over the electric potential
in cylindrical coordinates at $t=3.44\ \mu s$ for the 3D cylindrical beam case
with $a=2.5$ mm. In (a) the white arrow indicates the starting of the traced
electron trajectory and the red arrows indicate its direction when reflected
from the ESW potential-well boundaries at $z=17$ and 18.5 cm. In (b), traced
electrons ‘1’ and ‘3’ are reflected from the virtual cathode, where ‘1’ leaves
from conducting boundary at $r=1.5$ cm, while ‘3’ keeps bouncing between
virtual cathode and the ion sheath at $z=40$ cm. Traced electron ‘2’ leaves
the domain from the $z=40$ cm conducting boundary condition.

##  IV Excitation of surfaces waves in a plasma beam

###  IV.1 Trivelpiece-Gould surface waves in 3D cylindrical beams

As was shown in the previous section, the geometry of the plasma has important
effects on the types of electron-induced instabilities that may arise in a
plasma beam. For the 3D cylindrical beam cases with beam radius $a=5.0$ and
$7.5$ mm, as the beam potential decreases and the electrons fill in the
cylindrical ion beam, long-wavelength waves appear, which travel along the
beam axis with high phase speeds. Figure 17(a) shows a wave of long
wavelength, $\lambda_{w}=$ 6 cm, traveling at a phase speed of
$v_{\phi}=2.42\times 10^{6}$ m/s in the $+z$ direction with a frequency of
$\omega=\kappa v_{\phi}=2.49\times 10^{8}$ rad/s, where
$\kappa=2\pi/\lambda_{w}$ is the wavenumber, observed in the simulations. The
geometry of the beam, the low frequency ($\omega/\omega_{pe}=0.33$), and the
long wavelength of this wave are typical of a Trivelpiece-Gould (TG) surface
wavetrivelpiece1959space , commonly found in electron beams in cylindrical
cavitiesmoisan1982experimental .

The $\omega-\kappa$ dispersion relation for a surface wave in a cylindrical
cold electron beam with no external magnetic field, first derived by
Trivelpiece and Gould (TG)trivelpiece1959space , is given by,

$\displaystyle\epsilon_{1}\frac{J_{n}^{\prime}({\rm i}\kappa a)}{J_{n}({\rm
i}\kappa a)}=\epsilon_{0}\frac{I_{n}^{\prime}(\kappa a)K_{n}(\kappa
b)-I_{n}(\kappa b)K_{n}^{\prime}(\kappa a)}{I_{n}(\kappa a)K_{n}(\kappa
b)-I_{n}(\kappa b)K_{n}(\kappa a)},$ (9)

where $J_{n}$ is the $n^{\rm th}$ order Bessel function of the first kind and
$I_{n}$ and $K_{n}$ are the modified Bessel functions of the first and second
kind, respectively, $\epsilon_{1}=\epsilon_{0}(1-\omega_{pe}^{2}/\omega^{2})$
is the dielectric constant of the plasma medium, $a$ and $b$ are radii of the
plasma beam and the cylindrical conductor channel, respectively, and $n$ is
the azimuthal ($\theta$ direction) index of the electric potential
perturbationtrivelpiece1959space . Figure 17(b) shows a comparison of the
axisymmetrical ($n=0$) TG wave dispersion relation of Eq. 9 with the long-
wavelength wave for the $a=5.0$ mm case. Similarly for the case with $a=7.5$
mm, these long-wavelength waves show up near $t=4\ \mu s$, as shown by high
frequency fluctuations in Fig. 18(a). For the $a=7.5$ mm case, the long-
wavelength waves appear at $\omega/\omega_{pe}=0.37$ and $\lambda_{w}=5.94$ cm
(see blue ‘X’ in Fig. 17(b)). The surface wave properties of these two cases
are listed in Tab. 2. The TG waves have also been observed in experiments by
Lynov et al.lynov1979observations , and Moisan et al.moisan1982experimental ,
where these waves were excited by pulsing the plasma column in a cylindrical
TG wave-guide cavity. However, no experimental evidence which shows their
spontaneous excitation during beam neutralization has been reported in the
literature, to the best of our knowledge.

The surface wave phase speed, $v_{\phi}$, is higher than the local electron
velocities, as shown in Fig. 18(b) by the EVDF in the region $z>25$ cm for the
cases with $a=5.0$, and 7.5 mm. From the EVDFs shown in Fig. 18(b), it appears
that these waves become excited by the high energy electrons that are present
at the beginning of beam neutralization, starting from $t=1.0\ \mu s$, when
the beam potential is high, causing the electrons to accelerate to high
velocities. This acceleration is enhanced by the large sized ESWs, which
oscillate upstream of the electron source for both the $a=5.0$ and 7.5 mm
cases, as was shown in Fig. 8.  This oscillation of upstream ESWs is also
shown in phase space in Figs. 19(a) and 19(b) between $t=1.52$ and $3.68\ \mu
s$ and $t=1.52$ and $3.28\ \mu s$ for the $a=5.0$ mm and 7.5 mm cases,
respectively. For both cases, we can see that the large ESWs accelerate the
electrons coming from the electron source at $z=10$ cm. These electrons then
gain velocities close to $2.42\times 10^{6}$ and $2.62\times 10^{6}$ m/s in
the downstream region of the electron source at $t=2.40\ \mu s$ and $2.16\ \mu
s$, respectively for the $a=5.0$ and 7.5 mm cases, as shown in Figs. 19(a) and
19(b). These high-speed electrons then excite the surface waves seen in both
these cases at later times with the respective phase speeds of $2.42\times
10^{6}$ and $2.62\times 10^{6}$ m/s for the $a=5.0$ and 7.5 mm cases. After
the high energy electrons create surface waves, these electrons leave the beam
as the beam potential decreases with time. This leaves the beam with
insufficient number of electrons required to resonate and dampen these surface
waves by Landau damping, as shown by the EVDF and $v_{\phi}$ comparisons in
Fig. 18(b). Therefore, these waves survive and can be seen as the beam
potential falls below 3.5 and 4.0 V, respectively, for the $a=5.0$ and 7.5 mm
cases. However, these waves do not seem to be excited for the case with
$a=2.5$ mm because the phase speed of the TG wave corresponding to a
wavelength of $\lambda_{w}=6$ cm is about $1.33\times 10^{6}$ m/s, which is
sufficiently small to be Landau damped by the particles in the tail of the
electron distribution function, as shown by the EVDF comparison with
$v_{\phi}$ in the inset of Fig. 18(b).

(a) TG waves in a 3D beam with $a=5.0$ mm.

(b) Dispersion curve for TGtrivelpiece1959space  waves.

Figure 17: Long wavelength surface waves traveling along the beam axis in a
cylindrical beam. (a) The time evolution of the electric potential profile
($\Delta t=0$ corresponds to $t=4.72\ \mu s$ from the start of the PIC
simulation, at the end of rapid neutralization) shows waves with
$v_{\phi}=2.42\times 10^{6}$ m/s, and wavelength, $\lambda_{w}$ = 6.1 cm. In
(b), the curves show the theoretically calculated phase speed from dispersion
relation of axisymmetric TG surface waves (Eq. 9); ‘X’ marks shows the PIC
simulation ($\omega-\kappa$) results for the phase velocity for different beam
radii.

(a) Time evolution of electric potential at $z=20$ cm along beam axis.

(b) EVDF vs $v_{\phi}$ for different beam radius.

Figure 18:  (a) Beam neutralization vs beam radius, $a$, for the 3D
cylindrical beam cases. In (b), EVDFs are shown at $z=30$ cm and $t=4.72\ \mu
s$ for all cases. In (b), the inset figure has a logarithmic $y$ axis.

###  IV.2 Dispersion relation of surface waves in 2D planar vs 3D cylindrical
electron beams

The question arises why these long wavelength surface waves do not appear in
the 2D planar beam simulation. A fully neutralized 2D planar electron beam in
a channel has a geometry similar to a plasma slab which lies between two
conducting surfaces, as shown in Fig. 20(a), where, $a$ and $b$ are the half-
widths of the plasma beam and the channel, respectively. Following a procedure
similar to Vedenovvedenov1965solid , Krallkrall1973principles , and
Trivelpiece-Gouldtrivelpiece1959space , it can be shown that a cold-electron
beam assumption gives,

$\displaystyle\frac{\omega}{\omega_{pe}}=\frac{1}{\left[1+\frac{\tanh{[(b-a)\kappa]}}{\tanh{a\kappa}}\right]^{1/2}},$
(10)

where $\omega_{pe}=\sqrt{n_{e}e^{2}/m_{e}\epsilon_{0}}$ is the electron plasma
frequency of the beam. Further details are discussed in Appendix C. As seen
from Eq. 10, the 2D surface waves depend on the width of the beam and channel,
where all the geometric parameters converge to $\omega/\omega_{pe}=1/\sqrt{2}$
for high values of wavenumber, $\kappa$, as shown in Fig. 20(b). Most notably,
$\omega/\omega_{pe}\propto\sqrt{\kappa}$ for small values of $\kappa$ and
$a<<b$, i.e. a thin film surface wave, and $\omega/\omega_{pe}\to 1/\sqrt{2}$
for $a\to\infty,\ b\to\infty$, and $a<<b$, i.e. surface wave with no
conducting wallsvedenov1965solid .

A comparison of the 2D planar and 3D cylindrical surface wave dispersion
relation is shown in Fig. 21(a), where it can be seen that the 2D surface
waves have a much higher frequency than the 3D cylindrical waves of the same
wavenumber. Hence, for surface waves to be excited in a 2D planar beam, the
electrons would need to reach much higher velocity values than they would be
required to excite a 3D cylindrical surface wave. For example, to excite a
surface wave of $\lambda_{w}=6$ cm in a 2D planar beam of beam half-width
$a=2.5$ mm, the required phase speed is $v_{\phi}=3.23\times 10^{6}$ m/s,
which is about 2.5 times that required for the corresponding 3D cylindrical
beam of beam radius $a=2.5$ mm. Since electrons are injected with $T_{e}=2$ eV
they do not reach such high velocities in our 2D planar beam case at the time
when they acquire the shape of a planar beam, and therefore we do not observe
surface waves in our 2D numerical simulation results.

(a) $a=5.0$ mm.

(b) $a=7.5$ mm.

Figure 19:  Time evolution of electrons in phase-space extracted in the
$x=1.5$ cm mid-section plane of a 3D cylindrical beam with $a=5.0$ and $7.5$
mm. In (a) and (b), the yellow dashed line represents the phase speed of the
surface waves found in their respective cases. Here, $v_{0}=1\times 10^{6}$
m/s.

(a) 2D plasma beam of width $2a$ in a channel of width $2b$.

(b) 2D Slab Dispersion relation for different values of $a$ and $b$.

Figure 20:  2D surface wave dispersion relation for plasma beam with different
plasma widths, $2a$, and channel widths, $2b$.

###  IV.3 Excitation of surface plasma waves by electrons in the tail of the
EVDF

Since the surface waves for the 3D cylindrical beams were found to be excited
by high energy electrons, we were curious to see whether we could excite them
in the 2D planar beam case by introducing a small fraction of 3% of electrons
with a specific velocity. To do this, we restart our 2D planar beam PIC
simulation from $3.92\ \mu s$, when the beam is near the end of the rapid
neutralization phase (see Fig. 3(a)) and randomly initialize 3% of the
electrons entering from the source at $z=10$ cm with a specific high velocity
of $2.01\times 10^{6}$ m/s. Figure 21(b) shows a comparison of the electric
potential profiles of the 2D planar beam cases with and without the artificial
excitation of surface waves by the high energy electron tail of
$v_{e}=2.01\times 10^{6}$ m/s at $t=4.08\ \mu s$. The excited surface wave
(dashed red curve) has a wavelength of 2.0 cm and a phase speed of
$v_{\phi}\approx 2.01\times 10^{6}$ m/s. From this, $\omega=\kappa
v_{\phi}=6.31\times 10^{8}$ rad/s, and $(\omega,\kappa)$ of the wave compares
well with the 2D surface wave dispersion curve, shown by the black ‘X’ in Fig.
21(a). Similarly, the surface waves, which were not originally self-excited in
the 3D cylindrical beam case with $a=2.5$ mm, were excited by the introduction
of 10% high energy electrons with $v_{e}=1.33\times 10^{6}$ m/s at $t=5.84\
\mu s$, as shown by the comparison of the original and artificially excited
potential profiles in Fig. 21(c) and a green ‘X’ in Fig. 21(a). Note that the
surface waves are only seen in the region of the beam where there are no ESWs
that can dominate the local electric potential and phase space, as they would
make it hard to identify the small amplitude surface waves. This artificial
excitation shows that a small number of high energy electrons can excite a
surface wave in either a planar or cylindrical beam. This also confirms our
assertion that the TG surface waves found in the $a=5.0$ and $7.5$ mm 3D
cylindrical beam cases are excited by a small number of high energy electrons
in the first phase of beam neutralization.

(a) 2D $vs$ 3D dispersion relation.

(b) Artificial excitation in 2D planar beam.

(c) Artificial excitation in 3D cylindrical.

Figure 21:  (a) 2D planar (Eq. 10) vs 3D cylindrical (Eq. 9) surface wave
dispersion relation and excitation of surface waves by tail electrons in (b)
2D and (c) 3D beams with $a=2.5$ mm. In (b) and (c), original result is shown
in black and the artificially perturbed result in red line.

##  V Conclusions

In this work, we demonstrated the formation of ESWs in the rapid
neutralization phase of 2D planar and 3D cylindrical beams by an external
electron emission source. Similar to Lan et al.lan2020neutralization ;
lan2020neutralization2 , ESWs are formed due to the two streams of
electronsomura1996electron  created by the electron reflection from the
sheaths at the axial boundaries (i.e., $z_{\rm min}$ and $z_{\rm max}$) and
the virtual cathode at $z=10$ cm. We demonstrated the tracking of ESWs along
the beam and their collisions, which results in their merging together to form
slow moving ESWs of larger lengths. Since the slow moving (or nearly
stationary) ESWs have a lower probability of colliding with the domain walls
or other ESWs, they will dissipate slowly and hence, will slow down the
process of beam neutralization. By comparing the individual ESW structures, we
identified that the ESWs in the 2D planar beam case corresponded to those of a
1D BGK mode, while the ESWs in the 3D cylindrical beam case did not. Since the
electron trajectories in the latter case were found to be highly three
dimensional, a BGK analysis with higher dimensionality would be required to
analyze 3D ESWs. The 1D BGK analysis of the 2D beam ESWs found in our
simulation led to the conclusion that they were of a larger length than usual
because of the non-Maxwellian nature of the beam electrons. This shows that
such ESWs are possible with a smaller total electron density perturbation as
compared to the case where the beam electrons are Maxwellian.

Furthermore, we demonstrated the spontaneous excitation of Trivepiece-
Gouldtrivelpiece1959space  surface waves in a 3D cylindrical beam with beam
radius, $a\geq 5.0$ mm, in which the initial high energy electrons were found
to be responsible for the excitation of the wave. By deriving the theoretical
dispersion relation for a cold 2D planar beam, we showed that the high phase
velocity requirement was why surface waves did not appear in our 2D planar
beam results. Finally, we demonstrated that by introducing a small fraction of
electrons with a velocity higher than the usual Maxwellian of $T_{e}=2$ eV,
surface waves could be excited even in the 2D planar beam case.

## Acknowledgments

This project at the University of Illinois Urbana-Champaign is supported by
the Department of Energy Fusion Energy Sciences DOE - Award DE-SC00021348. Dr.
I. D. Kaganovich is supported by the US DOE under contract DE-AC02-09CH11466
as a part of the Princeton Collaborative Research Facility (PCRF). We would
also like to thank NCSA for providing us with the computational resources on
Bluewaters petascale computing facilityBode2013 ; Kramer2015 . We would also
like to thank XSEDE for providing us with GPU computing resources on San Diego
Supercomputing Cluster (SDSC).

## Data availability statement

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

## Appendix

##  Appendix A Numerical convergence of PIC simulations and method of ESW
position location

To show numerical convergence, we performed a case with $a=5.0$ mm with twice
as many charged particles per cell as in the original case in which $F_{\rm
num}$ was 60. We performed this case with 128 GPUs until $t=2.2\ \mu s$ to
show convergence in the ion only ($t<1.0\ \mu s$) and ion-electron ($t>1.0\
\mu s$) regime of the beam neutralization. The electric potential at the
$z=20$ cm probe along the beam axis was found to agree within 2% for cases
with $F_{\rm num}=30$ and 60. Also, the time evolution of the electric
potential profile along the beam is shown in Fig. 22, and it can be seen that
the iso-potential lines compare well for both cases.

Figure 22:  Numerical convergence of the cylindrical beam electric potential
for $a=5.0$ mm case. Here, comparison of iso-potential lines with $F_{\rm
num}=60$(black), and $30$(red) is shown.

To locate the ESWs along the beam axis, we use the fact that the position of
the ESW is accompanied with a local maxima in the electric potential ($\phi$)
profile, as shown in Fig. 23. We extract 800 points spanned over the range
$z\in(0,40)$ cm along the beam axis. To smooth the $\phi(z)$, we subdivide 800
points along the beam axis into 160 $z$-location points,

$\displaystyle z_{\rm window}=\frac{1}{N_{w}}\sum_{n=1}^{N_{w}}z_{n},$ (11)

where $z_{\rm window}$ is the average $z$ location of the window of size
$N_{w}$. For all cases in this work, we use $N_{w}=5$. The five points around
those 160 points, $N_{w}$, are then averaged to assign the electric potential
value to those points, i.e.,

$\displaystyle\phi(z_{\rm
window})=\frac{1}{N_{w}}\sum_{n=1}^{N_{w}}\phi(z_{n}),$ (12)

This results in a smoother electric potential profile, $\phi(z_{\rm window})$,
as shown by orange line in Fig. 23. The profile, $\phi(z_{\rm window})$, is
then differentiated twice to satisfy, $d\phi(z_{\rm window})/dz_{\rm
window}\approx 0$ and $d^{2}\phi(z_{\rm window})/dz_{\rm window}^{2}<0$ to
obtain local maxima, shown by orange X in Fig. 23. This process is performed
for various time profiles to generate the locations of ESWs, shown in Figs.
5(b) and 7(a).

Figure 23:  Window averaging of the electric potential along the 2D planar
beam axis at $t=3.12\ \mu s$. The orange X shows the location of the local
maxima of electric potential peak.

##  Appendix B Applying 1D BGK theory to obtain $f_{ut}$ from the PIC EVDF
data

Here, we discuss how we computed the untrapped electron phase distribution,
$f_{ut}$, from the far-field EVDF $f_{\rm FF}$ obtained from the PIC result.
The analytical fit profile from Eq. 3, shown in Fig. 24, makes it easier for
us to perform the step 2 of 1D BGK approach. Note that the far-field in our
PIC result occurs near $\Delta z\approx 15$ mm, however, for this theoretical
analysis, we assume the far-field to be at $\Delta z=45$ mm to ensure that the
Eq. 3 fit decreases to nearly 0 V. From this point on, we will consider Eq. 3
and Fig. 24 as the electric potential profile of the ESW.

Figure 24:  (a) Electric potential profile of the ESW. (b) Schematic of
discretized phase grid of $(v_{z},\Delta z)$ space. In (b), $\delta z=0.15$ mm
and $\delta v_{z}=12,000$ m/s. (c) A schematic of trapped (red) and untrapped
(black) electron trajectories (i.e. constant total energy lines) for an ESW in
($v_{z},\Delta z$) phase space. In (b), since the magenta colored constant
$E_{\rm tot}$ line (also shown in (c)) passes through points 1, 2, and 3 in
the $(v_{z},\Delta z)$ phase space, the number of electrons are identical in
the highlighted magenta colored grid-boxes in (b).

Since all electrons outside of the ESW (i.e. $\Delta z\geq 45$ mm in Fig. 24)
are untrapped electrons, the untrapped distribution, $f_{ut}(v_{z},\Delta
z\to\infty)=f_{\rm FF}$. Numerically, we discretize the z-velocity space
($v_{z}$) from $-3\times 10^{6}$ to $3\times 10^{6}$ m/s into 500 bins
($\delta v_{z}=12000$ m/s), where all the electrons in a bin are considered to
be of the same velocity. We also discretize the $\Delta z$ space into 300
points between 0 to 45 mm. A schematic of the generated ($v_{z},\Delta z$)
phase-space grid is shown in Fig. 24. Now, we want to fill each cell of this
discretized grid with the value of number of electrons. We first fill the
$\Delta z=45$ mm column (i.e. column designated to far-field) of the phase-
space grid with the electron particle velocity data from PIC divided into 500
bins of discrete velocity space to generate the far-field EVDF, $f_{\rm FF}$,
that was shown in Fig. 10(c). Although the distribution is initially shown in
units of number of computational particles, it could easily be changed to the
units of number density, $\rm m^{-3}$ by multiplication by the factor $F_{\rm
num}/{\rm Vol}$, where $\rm Vol$ is the volume of cell of size $\delta z$.

The electrons in the rest of the phase-space of ESW are either trapped or
untrapped. At steady state of the ESW (i.e. in the reference frame of the
moving ESW), the electron trajectories in the phase space look like those
shown in Fig. 24, where along each electron trajectory, the total energy, i.e.
kinetic + potential energy, is conserved. This means that the number of
electrons are the same along each contour line shown in Fig. 24 and,
therefore, if we know the EVDF along any of the $\Delta z$ locations along
with the electric potential profile of ESW, we can compute the rest of the
phase space distribution. Moreover, since we know the EVDF in the far-field,
$f_{\rm FF}$, we can directly compute the untrapped electron phase
distribution from that, by following the dashed black lines in Fig. 24 from
the far-field EVDF, $f_{\rm FF}$, to any $\Delta z$ location inside the ESW.
Because the trapped electrons never leave the ESW and they have no electron
trajectories, shown by dashed red lines in Fig. 24, entering the far-field, we
cannot estimate the trapped EVDF from $f_{\rm FF}$ alone.

We will now apply the conservation of total energy to compute the phase-space
distribution of the untrapped electron particles based on the known $f_{\rm
FF}$. The electrons at the location $\Delta z$ and in the velocity bin
$v_{z}$, have a total energy of $\frac{1}{2}m_{e}v_{z}^{2}-e\Delta\phi$, and
will reach the far-field with a kinetic energy of $\frac{1}{2}m_{e}v_{z0}^{2}$
such that,

$\displaystyle\frac{1}{2}m_{e}v_{z0}^{2}=\frac{1}{2}m_{e}v_{z}^{2}-e\Delta\phi,$
$\displaystyle
v_{z0}=\text{sgn($v_{z}$)}\sqrt{v_{z}^{2}-2\frac{e}{m_{e}}\Delta\phi}.$ (13)

An example of such a trajectory is shown in Figs. 24 and 24, where the magenta
colored iso-contour line travels from the center ($\Delta z=0$) to the far-
field of the ESW. Since all the electrons which had velocity $v_{z}$ at
$\Delta z$ location reach the far-field with the velocity $v_{z0}$, the number
of electrons with velocity $v_{z}$ at $\Delta z$, i.e. $f_{ut}(v_{z},\Delta
z)$, should be equal to the number of electrons in the far-field distribution
velocity bin $v_{z0}$. In other words,

$\displaystyle f_{ut}\left(v_{z},\Delta z\right)=f_{\rm
FF}\left(v_{z0}\right).$ (14)

If we apply this to multiple points along the $\Delta z$ direction, we can
fill the phase grid of $f_{ut}$. After computing the EVDFs for all 300 $\Delta
z$ locations, the untrapped electron distribution in the phase-space was shown
in Fig. 11(a).

##  Appendix C Dispersion relation for 2D surface waves

(a) Schematic for 2D slab and 2D surface.

(b) Comparison of 2D Slab dispersion relation with cylindrical and planar
surface waves.

Figure 25:  2D surface wave dispersion relation for plasma beam with different
plasma widths, and channel widths.

Equation (4.8.6) of Krall et all.krall1973principles  is the general wave
equation,

$\displaystyle\nabla(\nabla.\vec{E_{1}})-\nabla^{2}\vec{E_{1}}=\frac{\omega^{2}}{c^{2}}\epsilon.\vec{E_{1}},$
(15)

used to derive a differential equation for determining the perturbed electric
field, $\vec{E_{1}}$, in the plasma medium of dielectric tensor, $\epsilon$.
However, Trivelpiece et al.trivelpiece1959space  used a simpler approach in
which, $\vec{E_{1}}=-\nabla\phi_{1}$, where $\phi_{1}$ is the perturbed scalar
electric potential. Then, Poisson’s equation in a dielectric medium is given
by,

$\displaystyle\nabla.\epsilon.\vec{E_{1}}=0,$ (16)
$\displaystyle\nabla.(\epsilon.(\nabla\phi_{1}))$ $\displaystyle=0,$

where the dielectric tensor, $\epsilon$, in 2D Cartesian coordinates is given
by,

$\displaystyle\epsilon$
$\displaystyle=\begin{bmatrix}\left(1-\frac{\omega_{pe}^{2}}{\omega^{2}}\right)&0\\\
0&\left(1-\frac{\omega_{pe}^{2}}{\omega^{2}}\right)\end{bmatrix},$ (17)

and $\omega_{pe}=\sqrt{\frac{e^{2}n_{0}}{m\epsilon_{0}}}$, is the electron
plasma frequency. After inserting Eq. 17 into Eq. 16, we obtain,

$\displaystyle\left[1-\frac{\omega_{pe}^{2}}{\omega^{2}}\right]\frac{\partial^{2}\phi_{1}}{\partial
y^{2}}+\left[1-\frac{\omega_{pe}^{2}}{\omega^{2}}\right]\frac{\partial^{2}\phi_{1}}{\partial
z^{2}}$ $\displaystyle=0,\text{for $-a<y<a$}$ (18)

inside the plasma beam. Outside the plasma beam,

$\displaystyle\frac{\partial^{2}\phi_{1}}{\partial
y^{2}}+\frac{\partial^{2}\phi_{1}}{\partial z^{2}}$ $\displaystyle=0,\text{
for $|y|>a$}$ (19)

because the dielectric tensor is the identity matrix in a vacuum. Both Eqs. 18
and 19 have to be solved to obtain the perturbation field $\phi_{1}$ and the
corresponding dispersion relation. These solutions should satisfy the
following boundary conditions,

$\displaystyle\phi_{1}|_{y\to a^{-}}$ $\displaystyle=\phi_{1}|_{y\to a^{+}},$
(20) $\displaystyle\left[\frac{\partial\phi_{1}}{\partial y}\right]_{y=0}$
$\displaystyle=0,$ (21) $\displaystyle\phi_{1}|_{y=b}$ $\displaystyle=0,$ (22)
$\displaystyle\left[1-\frac{\omega_{pe}^{2}}{\omega^{2}}\right]\left[\frac{\partial\phi_{1}}{\partial
y}\right]_{y\to a^{-}}$ $\displaystyle=\left[\frac{\partial\phi_{1}}{\partial
y}\right]_{y\to a^{+}}.$ (23)

The analytical solution of Eq. 18, obtained by separation of variables, is of
the form,

$\displaystyle\phi_{1}$ $\displaystyle=\left[A_{y}\exp{(\kappa
y)}+B_{y}\exp{(-\kappa y)}\right]\exp{\left[i(\kappa z-\omega t)\right]},$
(24)

where $A_{y}$ and $B_{y}$ are the constants of integration, $\kappa$ is the
wave-number in the $z$ direction, and $\omega$ is the frequency for the wave
along the axial (i.e. $+z$) direction. Similarly, the analytical solution of
Eq. 19 is,

$\displaystyle\phi_{1}$
$\displaystyle=\left[C_{y}\exp{(\kappa^{\prime}y)}+D_{y}\exp{(-\kappa^{\prime}y)}\right]\exp{\left[i(\kappa^{\prime}z-\omega
t)\right]},$ (25)

where $C_{y}$ and $D_{y}$ are also the constants of integration, and
$\kappa^{\prime}$ is the wave-number along the axial beam direction. Since the
surface wave produced has to be continuous at the plasma-vacuum interface,
$\kappa^{\prime}=\kappa$. Finally, when we substitute the boundary conditions
in Eqs. 20-23 into Eqs. 24 and 25, we obtain,

$\displaystyle\frac{\omega}{\omega_{pe}}=\frac{1}{\left[1+\frac{\tanh{[(b-a)\kappa]}}{\tanh{a\kappa}}\right]^{1/2}},$
(26)

which is the dispersion relation for the surface waves in a 2D planar beam of
half-width $a$ inside a channel of half-width $b$. Equation 26 converges to
$\omega/\omega_{pe}=1/\sqrt{2}$ for large values of $\kappa$, as shown in Fig.
20(b), which is the dispersion relation for an infinite cold-plasma and vacuum
interface surface wave. The value of $\omega/\omega_{pe}$ varies from 0.3 to
0.707 for long-wavelength waves (i.e. low values of $\kappa$) for $a=2.5$,
$5.0$, and $7.5$ mm, and $b=30$ mm, as shown in Fig. 20(b). Also, the 2D
dispersion relation changes its shape when the beam width, $2a$ is nearly
equal to the channel width $2b$, in which the value of $\omega/\omega_{pe}$
decreases to $1/\sqrt{2}$, shown in Fig. 20(b), with increase in wavenumber,
$\kappa$.

## References

  * (1) Korkut, B., Li, Z., Levin, D., et al., “3-D Simulation of Ion Thruster Plumes Using Octree Adaptive Mesh Refinement,” IEEE,Plasma Sciences, 2015\.
  * (2) Jambunathan, R. and Levin, D. A., “CHAOS: An octree-based PIC-DSMC code for modeling of electron kinetic properties in a plasma plume using MPI-CUDA parallelization,” Journal of Computational Physics, Vol. 373, 2018, pp. 571–604.
  * (3) Nuwal, N., Jambunathan, R., and Levin, D. A., “Kinetic Modeling of Spacecraft Surfaces in a Plume Backflow Region,” IEEE Transactions on Plasma Science, Vol. 48, No. 12, 2020, pp. 4305–4325.
  * (4) Chen, Y.-M., Sawadichai, R., Tian, S., Donnelly, V. M., Ruchhoeft, P., and Economou, D. J., “Nearly monoenergetic positive ion beam with self-neutralized space charge extracted from a pulsed plasma,” Journal of Physics D: Applied Physics, Vol. 52, No. 35, 2019, pp. 355205.
  * (5) Conrad, J., “Plasma source ion implantation: A new approach to ion beam modification of materials,” Materials Science and Engineering: A, Vol. 116, 1989, pp. 197–203.
  * (6) Vay, J.-L. and Deutsch, C., “Intense ion beam propagation in a reactor sized chamber,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, Vol. 464, No. 1-3, 2001, pp. 293–298.
  * (7) Vay, J. and Deutsch, C., “Charge compensated ion beam propagation in a reactor sized chamber,” Physics of Plasmas, Vol. 5, No. 4, 1998, pp. 1190–1197.
  * (8) Kaganovich, I. D., Shvets, G., Startsev, E., and Davidson, R. C., “Nonlinear charge and current neutralization of an ion beam pulse in a pre-formed plasma,” Physics of Plasmas, Vol. 8, No. 9, 2001, pp. 4180–4192.
  * (9) Berdanier, W., Roy, P. K., and Kaganovich, I., “Intense ion beam neutralization using underdense background plasma,” Physics of Plasmas, Vol. 22, No. 1, 2015, pp. 013104.
  * (10) Lan, C. and Kaganovich, I. D., “Neutralization of ion beam by electron injection: Accumulation of cold electrons,” Physics of Plasmas, Vol. 27, No. 4, 2020, pp. 043108.
  * (11) Kaganovich, I., Davidson, R., Dorf, M., Startsev, E., Sefkow, A., Friedman, A., and Lee, E., “Physics of neutralization of intense charged particle beam pulses by a background plasma,” Tech. rep., Lawrence Berkeley National Lab.(LBNL), Berkeley, CA (United States), 2009.
  * (12) Lan, C. and Kaganovich, I. D., “Neutralization of ion beam by electron injection: Excitation and propagation of electrostatic solitary waves,” Physics of Plasmas, Vol. 27, No. 4, 2020, pp. 043104.
  * (13) Bernstein, I. B., Greene, J. M., and Kruskal, M. D., “Exact nonlinear plasma oscillations,” Physical Review, Vol. 108, No. 3, 1957, pp. 546.
  * (14) Schamel, H., “Theory of electron holes,” Physica Scripta, Vol. 20, No. 3-4, 1979, pp. 336.
  * (15) Hutchinson, I. H., “Electron holes in phase space: What they are and why they matter,” Physics of Plasmas, Vol. 24, No. 5, 2017, pp. 055601.
  * (16) Saeki, K., Michelsen, P., Pécseli, H., and Rasmussen, J. J., “Formation and coalescence of electron solitary holes,” Physical Review Letters, Vol. 42, No. 8, 1979, pp. 501.
  * (17) Matsumoto, H., Kojima, H., Miyatake, T., Omura, Y., Okada, M., Nagano, I., and Tsutsui, M., “Electrostatic solitary waves (ESW) in the magnetotail: BEN wave forms observed by GEOTAIL,” Geophysical Research Letters, Vol. 21, No. 25, 1994, pp. 2915–2918.
  * (18) Trivelpiece, A. and Gould, R., “Space charge waves in cylindrical plasma columns,” Journal of Applied Physics, Vol. 30, No. 11, 1959, pp. 1784–1793.
  * (19) Krall, N. A. and Trivelpiece, A. W., Principles of plasma physics, McGraw-Hill, New York, 1973.
  * (20) Jambunathan, R. and Levin, D. A., “Advanced parallelization strategies using hybrid MPI-CUDA octree DSMC method for modeling flow through porous media,” Computers & Fluids, Vol. 149, 2017, pp. 70–87.
  * (21) Godunov, S. K., “A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics,” Matematicheskii Sbornik, Vol. 89, No. 3, 1959, pp. 271–306.
  * (22) Nuwal, N. and Levin, D. A., “Influence of secondary electron emission on plasma-surface interactions in the low earth orbit environment,” Plasma Sources Science and Technology, Vol. 30, No. 3, 2021, pp. 035015.
  * (23) Lynov, J.-P., Michelsen, P., Pécseli, H., Rasmussen, J. J., Saeki, K., and Turikov, V., “Observations of solitary structures in a magnetized, plasma loaded waveguide,” Physica Scripta, Vol. 20, No. 3-4, 1979, pp. 328.
  * (24) Ikezi, H., Barrett, P., White, R., and Wong, A., “Electron plasma waves and free-streaming electron bursts,” The Physics of Fluids, Vol. 14, No. 9, 1971, pp. 1997–2005.
  * (25) Stepanov, A. D., Gilson, E. P., Grisham, L. R., Kaganovich, I. D., and Davidson, R. C., “Dynamics of ion beam charge neutralization by ferroelectric plasma sources,” Physics of Plasmas, Vol. 23, No. 4, 2016, pp. 043113.
  * (26) Lan, C. and Kaganovich, I., “Neutralization of ion beam by electron injection: Accumulation of cold electrons,” Physics of Plasmas, Vol. 27, No. 4, 2020, pp. 043108.
  * (27) Omura, Y., Matsumoto, H., Miyake, T., and Kojima, H., “Electron beam instabilities as generation mechanism of electrostatic solitary waves in the magnetotail,” Journal of Geophysical Research: Space Physics, Vol. 101, No. A2, 1996, pp. 2685–2697.
  * (28) Chen, L.-J. and Parks, G., “BGK electron solitary waves: 1D and 3D,” Nonlinear Processes in Geophysics, Vol. 9, No. 2, 2002, pp. 111–119.
  * (29) Lefebvre, B., Chen, L.-J., Gekelman, W., Kintner, P., Pickett, J., Pribyl, P., Vincena, S., Chiang, F., and Judy, J., “Laboratory measurements of electrostatic solitary structures generated by beam injection,” Physical review letters, Vol. 105, No. 11, 2010, pp. 115001.
  * (30) Chen, L.-J. and Parks, G. K., “BGK electron solitary waves in 3D magnetized plasma,” Geophysical research letters, Vol. 29, No. 9, 2002, pp. 45–1.
  * (31) Eliasson, B. and Shukla, P. K., “Theory for two-dimensional electron and ion Bernstein–Greene–Kruskal modes in a magnetized plasma,” Journal of plasma physics, Vol. 73, No. 5, 2007, pp. 715–721.
  * (32) Ng, C. and Bhattacharjee, A., “Bernstein-Greene-Kruskal modes in a three-dimensional plasma,” Physical review letters, Vol. 95, No. 24, 2005, pp. 245004.
  * (33) Moisan, M., Shivarova, A., and Trivelpiece, A., “Experimental investigations of the propagation of surface waves along a plasma column,” Plasma physics, Vol. 24, No. 11, 1982, pp. 1331.
  * (34) Vedenov, A., “Solid state plasma,” Soviet Physics Uspekhi, Vol. 7, No. 6, 1965, pp. 809.
  * (35) Bode, B., Butler, M., Dunning, T., Hoefler, T., Kramer, W., Gropp, W., and Hwu, W.-m., “The Blue Waters Super-System for Super-Science,” Contemporary High Performance Computing, Chapman & Hall/CRC Computational Science, Chapman and Hall/CRC, April 2013, pp. 339–366.
  * (36) Kramer, W., Butler, M., Bauer, G., Chadalavada, K., and Mendes, C., “Blue Waters Parallel I/O Storage Sub-system,” High Performance Parallel I/O, edited by Prabhat and Q. Koziol, CRC Publications, Taylor and Francis Group, 2015, pp. 17–32.