We develop an analytical model for estimating the equilibrium quantities, such as electron temperature and number density, in an electron beam–plasma interaction system. This model provides a convenient way to calculate the effective electron temperature and density by considering the energy balance of the bulk cold electrons. Six energy sources/losses terms relevant to the cold electrons are accounted for, where quasi-linear theory is applied for estimating wave heating at equilibrium. We compare this calculation with the particle-in-cell (PIC) simulation results and find good agreement. Based on these results, we then consider two situations where we can simplify our model. The first is dominated by the balance between electron–electron Coulomb collisions and loss to the anode, which is mostly relevant to the conduction phase of plasma switches. The second is dominated by wave heating balanced by the anode loss, relevant to the electron beam–plasma discharge systems. We then couple our simplified energy balance model with the ion diffusion model and solve both the number density and the electron temperature as functions of the current density, electrode distance, pressure, and applied voltage, where a nice agreement is also obtained when comparing to PIC simulations.

## I. INTRODUCTION

Electron beam–plasma interactions have a wide range of applications in both space plasmas^{1–11} and low-temperature plasma devices^{12–16} and have been continuously studied both theoretically^{17,18} and experimentally.^{19–22} In low-temperature plasma sources, for example, certain hollow cathode discharges contain a process of beam electron emission from the cathode, ionizing the background plasma and generating plasma.^{11–13} High-energy electron-beam generated plasma is also frequently used in the plasma processing systems, with the advantage of a relatively simplified production of species while providing comparably high ion-to-radical production rates, increasing the processing rate.^{14,23–26} The plasma switch is another important example, where the estimate of electron density and temperature in its conduction phase is a prerequisite for the simulations of current interruption.^{27–29} Recently, attempts were made to combine the capacitively coupled (CCP) discharge with a fast electron beam to obtain a better control of an ion energy distribution at the rf-biased electrode.^{30} Therefore, an efficient estimate of electron temperature and plasma number density in the equilibrium state is of great importance in these practical electron beam–plasma devices. The current trend in the plasma processing systems is to operate at low gas pressure (<100 mTorr) to achieve a higher precision level in, for example, the etching width and depth of the wafer in semiconductor systems.^{12,13,24,25,31} At this pressure, the electron kinetic motions become crucial to understand such that the traditional glow discharge model no longer works.^{32} Despite the numerous simulation results for electron beam discharge systems,^{10,13,33–39} a simple physics-based analytical model for estimating equilibrium state plasma quantities at low pressure in the electron beam–plasma interaction systems, which could be directly implemented by experimentalists, is absent at present to the best of our knowledge.

In this paper, we propose an analytical model for estimating the equilibrium quantities, such as electron density and temperature for an ubiquitous beam–plasma system bounded by a negatively biased cathode and a grounded anode and confirm its reliability using Particle-In-Cell (PIC) simulations. We consider the energy balance of cold bulk electrons, taking into account six energy sources/losses. We also show that the model could be simplified under certain conditions to give simpler estimates of plasma density and electron temperature as functions of the current density, electrode distance, gas pressure, and applied voltage to the cathode. We consider two different situations. In the first situation, the system is balanced by Coulomb collision heating and loss to the anode, most relevant to high-density plasmas, such as the conduction phase of a plasma switch,^{27} whereas the second situation has a lower density, which is achieved when the system is balanced by the wave heating and loss to the anode, most relevant to normal electron beam–plasma discharge systems.^{27–29}

This paper is organized as follows: The proposed model for estimating the plasma equilibrium quantities is described in Sec. II, followed by the simulation setup and results in Sec. III. The way to simplify our model and to estimate equilibrium density and temperature are described in Sec. IV, with a comparison to several more PIC simulations. Finally, summary and discussion are provided in Sec. V.

## II. THEORETICAL BACKGROUND

In this section, we show the analytical model that we use to obtain the electron temperature and plasma density estimates. We focus on the energy balance of cold bulk electrons. The beam electrons being scattered by elastic and inelastic collisions are treated as bulk electrons. The bulk electron density has already taken this fraction of electrons into consideration. Thus, our purpose is to estimate the temperature of cold electrons and plasma density from a theoretical standpoint. When such a current-carrying plasma achieves a steady state, the energy that transfers to the cold electrons (say, heating sources) should be balanced by the energy dissipation (say, cooling sources). In our case, three heating sources of cold electrons should be considered:

Heating from the beam electrons through the electron–electron (e–e) Coulomb collisions.

Heating due to the ionization by the electron beam.

Wave heating through the wave–particle interaction.

Three cooling sources can be important:

Cooling due to electron–neutral (e–n) elastic collisions.

Cooling due to electron–neutral (e–n) inelastic collisions.

Cooling due to the loss of energetic population to the anode.

Note that the ionization term may also be a cooling source, depending on the average energy of secondary electrons and cold electron temperature. The energy balance of cold electrons is illustrated by Fig. 1.

^{40}$\u27e8F\u27e9= n e e 4 l n \Lambda 2 \pi \epsilon 0 2 T e \Phi ( v e b v t e )$, where $ n e$ is the number density of cold electrons, $ T e$ is the cold electron temperature, $ln \Lambda $ is the Coulomb logarithm, $ \epsilon 0$ is the vacuum permittivity, $ v t e= 2 T e / m e $ is the thermal speed of cold electrons, $ v eb= 2 E e b / m e$ is the velocity of the beam electrons, $ E e b$ is the beam electron energy, and $ \Phi (x)= 1 \pi x 2 ( \u222b 0 x e \u2212 \xi 2 d \xi \u2212 x e \u2212 x 2 )$. When $x\u2192+\u221e$, $ \Phi (x)\u22481/2 x 2$. The energy transferred to the cold electrons can be estimated by $ W e e= n e b v e b\u27e8F\u27e9$, where $ n e b$ is the beam electron density. Therefore, this term can be estimated by

The unit of $ W e e$ is $ W / m 3$.

*R*is a random number in [0, 1], describing the randomness of electron energy produced by the ionization, and

*I*is the ionization energy of the particles being considered. By integrating over

*R*, the average energy of the secondary electrons can be calculated and is denoted by $ E mean secondary$. Comparing this value with the effective electron temperature and multiplying the ionization frequency $ f ioniz= \sigma ioniz( E e b)N n e b 2 E e b m e$, the total energy transfer rate due to the ionization can be calculated by

^{12,13}

*N*is the density of the background gas, $ \sigma ioniz$ is the ionization cross section, and the $1.5 T e$ comes from the fact that our velocity space has three dimensions.

^{12,13}since the system should be dominated by weak turbulence. Therefore, the wave heating can be roughly estimated by the quasi-linear theory, which describes the interaction between waves and particles by assuming the stochastic motion of electrons in the electric wave field:

^{12,25,27,33,43}$ \u2202 f \u2202 t= \u2202 \u2202 v ( D \u2202 f \u2202 v )$, where

*f*is the electron velocity distribution function and $D= e 2 \pi 2 m e 2 v x| E k ( k = \omega / v x ) | 2$ is the quasi-linear diffusion coefficient, where $ E k$ is the Fourier component of the electric field associated with a Langmuir wave with wavenumber

*k*. This diffusion coefficient can be easily calculated from the Vlasov equation by first linearizing the equation and then taking a spatial average. After applying a spatial Fourier transform to the electric field, $ E k$ in the above expression can be obtained and

*D*can be readily calculated. To get an analytical expression of the stochastic heating, we need to assume a certain distribution of the field spectrum $ E k$. We assume that the field spectrum takes the form of $I(k)=1/2 \epsilon 0 E k 2\u223c I 0 ( k / k 0 ) \u2212 2$, with $ k 0= \omega p e/ v e b$ ( $ \omega p e$ is the electron plasma frequency) and $ I 0=1/2 \epsilon 0 E k 0 2$. This assumption will be further justified by the simulation results in Sec. III. Therefore, by following the derivations in Ref. 32 and integrating the electron energy over the velocity space, we can estimate the stochastic heating due to the wave as follows:

*L*is the system length. This expression will be implemented as a wave heating term of our analytical model of beam–plasma systems. Note that the

*k*in Eq. (3) refers to the wave number of Langmuir waves instead of the Boltzmann constant.

*M*is the mass of the neutral molecules (H

_{2}) or atoms (Ar).

^{34}For every single inelastic e–n collision, we have $ \Delta \epsilon elas H 2( E e \u2265 E excit H 2 = 0.55 eV)=0.55 eV$, and the total energy exchange due to the inelastic e–n collisions can be calculated by the formula

^{44,45}

*erf*is the error function, and $Z(x)=\u2212 exp( \u2212 x)$ $( x 3 + 3 x 2 + 6 x + 6)$ is obtained by fitting the cross section of an inelastic collision as

Note that in the calculation of Eq. (5), we have assumed a Maxwellian electron velocity distribution function.

*x*is the direction perpendicular to the anode. Therefore, we can follow the similar logic as in Ref. 46 and integrate the distribution close to the sheath edge,

## III. SIMULATION RESULTS

### A. Simulation setup

We perform one-dimensional PIC simulations of an electron beam–plasma interaction system using the PIC code, EDIPIC,^{47} to better understand the derived energy balance model in Sec. II. EDIPIC is a parallelized implicit electrostatic PIC code resolving one spatial coordinate and three velocity components (1D3 V), which is a thoroughly benchmarked PIC code developed by Princeton University.^{48} EDIPIC incorporates the Boris scheme for particle pushing, and the field is solved using the Thomas algorithm.^{49} It can perform full kinetic modeling as well as Monte-Carlo Collision (MCC) algorithms for electron–neutral and ion–neutral collisions. A Langevin model for the Coulomb collisions and a surface emission model for the plasma–wall interactions are also included in the code. All the simulation parameters are listed in Table I. The initial plasma density and temperature are $ n e 0= 10 19 m \u2212 3$, $ T e 0=1 eV$ for Coulomb collision heating-dominated cases and $ n e 0= 10 17 m \u2212 3$, $ T e 0=0.2 eV$ for wave heating-dominated cases. The time steps and cell sizes of the Coulomb collision heating-dominated cases are $ \Delta t=1\xd7 10 \u2212 3 ns\u22481.12\xd7 10 \u2212 3 ns=0.2 \omega p e 0 \u2212 1$, $ \Delta x=1.13\xd7 10 \u2212 3 mm$ $=1/2 \lambda D e 0$, while the time steps and cell sizes of wave heating-dominated cases are $ \Delta t=1\xd7 10 \u2212 2 ns\u22481.12\xd7 10 \u2212 2 ns=0.2 \omega p e 0 \u2212 1$, $ \Delta x=7.84\xd7 10 \u2212 3 mm=1/2 \lambda D e 0$. Initially, the number of macro-particles per cell is set to be 400. In the steady state, however, the number of macro-particles per cell depends on the equilibrium plasma density. For wave heating-dominated cases, it is at least 3000, whereas for Coulomb heating-dominated cases, the value varies from 50 to 800. A schematic diagram of the simulation setup is plotted in Fig. 2. Two electrodes are placed at the two boundaries. The cathode is negatively biased, and anode is grounded. In order to prove the validity of our theory, we used two types of gases, hydrogen and argon, for the two types of cases, e–e Coulomb heating-dominated cases and wave heating-dominated cases that we study here. The neutral component has a uniform density in the simulations. As mentioned in Sec. I, the wave heating-dominated cases are relevant to a beam–plasma discharge system and the Coulomb heating-dominated cases are related to the conduction phase of a plasma switch. Note that for the parameters that we are considering here, the Secondary Electron Emissions (SEEs) are not expected to be important since the steady state ion flux to the cathode is around $3\xd7 10 22 s \u2212 1 m \u2212 2$ (see Fig. 5). Assuming the SEE coefficient is $ \gamma S E E=0.1$ (reasonable assumption, see, for instance, Ref. 50) and suppose the SEE flux $ F S E E=1.5\xd7 10 18 s \u2212 1 m \u2212 2$. Comparing to the electron beam flux $ F e b= n e b v e b=1.559\xd7 10 23 s \u2212 1 m \u2212 2$, the SEE flux is much smaller.

Case number . | Case type . | $ n e b/ n e 0$ . | V (eV)
. | $ I e b( A / m 2)$ . | n_{e0} (m^{−3})
. | T_{e0} (eV)
. | P (mTorr)
. | Gas type . | L_{x} (cm)
. |
---|---|---|---|---|---|---|---|---|---|

1 | e–e | 0.0048 | −18 | 15 000 | 10^{19} | 1 | 76.2 | H | 2 |

2 | e–e | 0.0048 | −30 | 25 000 | 10^{19} | 1 | 76.2 | H | 2 |

3 | e–e | 0.0048 | −42 | 35 000 | 10^{19} | 1 | 76.2 | H | 2 |

4 | e–e | 0.0029 | −30 | 15 000 | 10^{19} | 1 | 76.2 | H | 2 |

5 | e-e | 0.0067 | −30 | 35 000 | 10^{19} | 1 | 76.2 | H | 2 |

6 | Wave | 0.048 | −18 | 1 500 | 10^{17} | 0.2 | 7.62 | Ar | 6 |

7 | Wave | 0.048 | −30 | 2 500 | 10^{17} | 0.2 | 7.62 | Ar | 6 |

8 | Wave | 0.048 | −42 | 3 500 | 10^{17} | 0.2 | 7.62 | Ar | 6 |

9 | Wave | 0.029 | −30 | 1 500 | 10^{17} | 0.2 | 7.62 | Ar | 6 |

10 | Wave | 0.067 | −30 | 3 500 | 10^{17} | 0.2 | 7.62 | Ar | 6 |

Case number . | Case type . | $ n e b/ n e 0$ . | V (eV)
. | $ I e b( A / m 2)$ . | n_{e0} (m^{−3})
. | T_{e0} (eV)
. | P (mTorr)
. | Gas type . | L_{x} (cm)
. |
---|---|---|---|---|---|---|---|---|---|

1 | e–e | 0.0048 | −18 | 15 000 | 10^{19} | 1 | 76.2 | H | 2 |

2 | e–e | 0.0048 | −30 | 25 000 | 10^{19} | 1 | 76.2 | H | 2 |

3 | e–e | 0.0048 | −42 | 35 000 | 10^{19} | 1 | 76.2 | H | 2 |

4 | e–e | 0.0029 | −30 | 15 000 | 10^{19} | 1 | 76.2 | H | 2 |

5 | e-e | 0.0067 | −30 | 35 000 | 10^{19} | 1 | 76.2 | H | 2 |

6 | Wave | 0.048 | −18 | 1 500 | 10^{17} | 0.2 | 7.62 | Ar | 6 |

7 | Wave | 0.048 | −30 | 2 500 | 10^{17} | 0.2 | 7.62 | Ar | 6 |

8 | Wave | 0.048 | −42 | 3 500 | 10^{17} | 0.2 | 7.62 | Ar | 6 |

9 | Wave | 0.029 | −30 | 1 500 | 10^{17} | 0.2 | 7.62 | Ar | 6 |

10 | Wave | 0.067 | −30 | 3 500 | 10^{17} | 0.2 | 7.62 | Ar | 6 |

In the simulations, e–n collisions and e–e Coulomb collisions are turned on. Thermionic electrons can be accelerated by the cathode sheath and form an electron beam, interacting with the bulk cold plasma until a steady state is achieved.

### B. Time evolution of physical parameters before reaching a steady state

Figure 3 shows the time evolution of the electron temperature, ion density, and electron loss to the right wall for the e–e Coulomb heating-dominated cases. At 30 000 ns, all the simulations reach a steady state. We see that the steady state electron temperature and density depends primarily on the beam energy (controlled by the cathode sheath potential). Comparing Case 2 $( E e b=30 eV )$ and Case 3 $( E e b=42 eV)$, we see that an increase of beam energy results in a decrease of steady state electron temperature, whereas the plasma density increases with the increase of beam current. This trend can be explained when coupling the heating model described in Sec. II with the ion diffusion model [see Eq. (13) in Sec. IV]. Figure 4 gives the same plot for the wave heating-dominated cases. We see that the different cases have nearly the same steady state electron temperature (note that the initial electron temperature for the cases in Fig. 4 is the same $ T e 0=0.2 eV$, not shown in the figure), whereas the steady state electron density depends strongly on the beam current. These trends can also be explained by our analytical model [see Eq. (15)].

At a steady state, our analytical theory can be directly applied. Before we show any definite results, we need to verify that the system has reached a steady state by inspecting the particle balance of the beam–plasma system. For a steady state, the ion loss should be balanced by the production of ions (i.e., $ \u2202 \Gamma i \u2202 x=S$, where Γ_{i} is the ion flux and *S* is the ionization rate). In Fig. 5, two example cases (Case 2 and Case 7) are considered. We can see that the integral of the ionization rate over *x* is roughly the same as the ion flux for both cases, which indicates that the production and the loss of the ions are balanced. Therefore, we can say that our system achieves a steady state.

### C. Comparison of the analytical model prediction against PIC simulation at a steady state

An important assumption of our model for wave heating is the $I\u223c I 0 ( k / k 0 ) \u2212 2$ spectrum for the quasi-linear model. This power law is shallower compared to the power law of strong Langmuir turbulence,^{12,13,18,20,51–53} indicating that the energy transfer is more efficient. We verify this relation by plotting the steady state energy spectrum for all the ten cases, which is shown in Fig. 6. A power law of approximately −2 is observed in all the cases, thereby justifying our assumption.

In the following part of this section, we analyze the simulation results of the two representative cases, Case 2 and Case 7 in depth, including the profiles, the electron velocity distribution function, and the energy balance of different terms given by Eq. (7) in a steady state to verify our analytical model. In Sec. IV, we further simplify our theory under two different limits, the e–e Coulomb collision dominated heating and wave dominated heating and then coupling with ion diffusion model, to obtain estimates of steady state electron temperature and plasma density. The simplified analytical theory is then compared with all simulation cases (Cases 1–10), and a good agreement is found.

Profiles of the critical parameters for Case 2 and Case 7 are plotted in Fig. 7. Subfigures (a) and (d) show the profile of the electric potential. Some fluctuations near the cathode can be seen, which indicates the generation of electrostatic waves due to the beam injection. Subfigures (b) and (e) show the profile of the ion density, a parabolic-like distribution over *x* is observed. For our case, due to the relatively low background pressure of 76.2 mTorr for Case 2 and 7.62 mTorr for Case 7, the mean free paths for the ionization are of the order of $10 cm$ for Case 2 and $6 cm$ for Case 7, which are both comparable to or larger than the simulation length. Thus, the ionization probability approximately follows a linear profile with the distance *x*, and an approximately linear ion flux profile and a parabolic-like number density profile are formed. More details can be seen in Sec. IV. Subfigures (c) and (f) show the profiles of the electron temperature $ T e$ of cold electrons as a function of *x* for Case 2 and Case 7, respectively. As we can see, despite a variation of a factor of two in $ T e$, the ion flux to the cathode and anode remains similar (see Fig. 5). Since the mean energy of cold electrons in three dimensions is very close and approximately follows the Maxwellian distribution function (because of the isotropic scattering with the background gas and Coulomb collisions with other electrons), the effective electron temperature in 3D can be taken as $ T e real\u22481.5 T e$.

The Electron Velocity Distribution Functions (EVDFs) for the two representative cases are shown in Fig. 8. As we can see, despite some high-energy tail in the bulk plasma, the distribution functions of bulk electrons at three different positions are nearly Maxwellian for both Case 2 and Case 7, indicating that a steady state plasma is achieved. Note that in the simulations, once the beam electrons go through the ionization or elastic collisions, they are eliminated from the beam electron system and become the population of cold bulk electrons, thereby forming a fraction of the high-velocity tail part of the bulk EVDF. The electron–electron Coulomb collision frequencies for Case 2 and Case 7 are $ \nu e e= n e \sigma 90 v e b\u22487.1\xd7 10 8 s \u2212 1$ and $ \nu e e=6.7\xd7 10 6 s \u2212 1$, where $ \sigma 90$ is the cross section that makes the electrons turn a 90° angle, which is calculated by Eq. (3.3.7) in Ref. 32, whereas the electron–neutral collision frequencies are $ \nu e n=8.7\xd7 10 7 s \u2212 1$ and $ \nu e n=8.7\xd7 10 6 s \u2212 1$, respectively. The beam energy relaxation length for Case 2 can be estimated as $ l relax , e e= v e b/ \nu e e\u22484.6 mm$, which is roughly the same as the distance that is indicated by Figs. 8(a) and 8(b). The beam energy relaxation length for Case 7, on the other hand, must be calculated using the quasi-linear theory of Ref. 54. $ l relax , Q L= v e bln\u2061 \varphi 0 \u2212 1/2\gamma \u22481.5 mm$, where $\gamma =\pi \omega p e n e b/ n e ( v e b \Delta v e b ) 2$ is the growth rate of the two-stream instability and $ \varphi 0\u22481/130$. As we can see from Fig. 8(a), at $x=5 mm$, the beam energy has been relaxed and no clear peak can be observed.

As we have mentioned, Eqs. (1)–(5) can be compared with Eq. (6) only after they are integrated over the whole simulation domain. In order to check if the estimates for the ion loss to the anode agree with the simulation results, we plug the electron density and the electron temperature from our PIC simulations into Eqs. (1)–(5) for Case 2 and get the estimate of the energy loss to the anode $ Q anode theory$ by $ Q anode theory=\u222b( W e e + W ioniz + W wave$ $ \u2212 W elas \u2212 W inelas)dx=4.59\xd7 10 5 W / m 2$. This value can be compared with $ Q anode real$ that represents the simulated heat flux to the anode. We have $ Q anode real=5.08\xd7 10 5 W / m 2$, which is similar to $ Q anode theory$. Therefore, we can say that our theory achieves a good agreement with the PIC simulation. The discrepancy here mainly comes from the inaccurate calculation of $ W e e$, where the function $ \Phi ( v e b / v T e)\u2248 v T e 2/2 v e b 2$ is accurate only when $ v e b/ v T e$ is sufficiently large, which is not accurate everywhere in the plasma because the beam energy decreases when it passes through the bulk plasma [see the beam EVDF in Fig. 8].

We further plot the profiles of five terms as a function of *x* for both Case 2 [shown in Fig. 9(a)] and Case 7 [shown in Fig. 9(b)]. For Case 2, a comparison of the five terms shows that the energy exchange between electrons and neutrals, and the wave heating are much smaller compared to the energy transfer due to Coulomb collisions. The same plot is also made for Case 7 in Fig. 9(b), where we see that the wave heating dominates over the heating due to Coulomb collision. This is primarily because of a lower plasma density and a lower electron temperature in Case 7. A lower density reduces the Coulomb collision frequency, whereas a lower temperature makes the relative wave amplitude $W= E 2/4 n e T e$ larger,^{12} resulting in a stronger wave and more wave heating. Note that these observations are different from that in highly collisional plasmas because of the relatively low background gas pressure.

## IV. ESTIMATION OF *n*_{e} AND *T*_{e}

Now, we take a step further and try to further simply our theory. Despite only considering two situations here, one should keep in mind that our theory can always be simplified into other forms depending on which energy terms become dominant.

### A. Estimation of *n*_{e} and *T*_{e} in an e–e Coulomb collision dominated case

The first case is when the energy balance is dominated by the heating due to e–e Coulomb collisions and cooling due to loss to the anode, which means that we can use only these two terms to estimate the electron temperature under a relatively low pressure. Since the energy loss to the anode is sensitive to the anode sheath voltage drop, a simple model for the anode sheath voltage drop is required and is shown as follows.

*x*(thus, the ion fluxes to the cathode and anode are equal). Also, for all the thermionic electrons, no matter whether they go through the ionization and scattering or not, we assume that their energy is still higher enough to penetrate the anode sheath. Therefore, we do not need to consider the contribution of thermionic electrons to the anode sheath voltage. As a result, the electron current to the anode is balanced by the ion current to the two electrodes and then we have the current balance equation

^{32}

The anode sheath voltage drop can be calculated as $\varphi = T e/2 ln( m i / 2 \pi m e)$ in Eq. (8) so that $ v min$ can be estimated as $ v min= T e m e ln \u2061 ( m i 2 \pi m e )$. Thus, Eq. (8) can be used to approximate the anode sheath voltage drop. Please note that Eq. (8) works perfectly only when the beam electrons are sufficiently scattered and do not contribute to the current to the anode. For some cases with high beam energy and low plasma density, this scattering may not be sufficient, resulting in discrepancy between theory and simulations.

*p*and $ T g$ are the background neutral gas pressure and the neutral gas temperature. $ D a=( T i + T e) D i D e/( T i D e + T e D i)$ is the ambipolar diffusion coefficient and $ \mu i$ is the ion mobility. According to Eqs. (10) and (11), we obtain

*V*(or we say $ E e b$), the electrode distance

*L*, and the background neutral gas pressure

*p*. Please note that this equation can be invalid when $ E e b$ becomes very large because the wave heating and the e–n energy exchange may become important.

### B. Estimation of *n*_{e} and *T*_{e} in a wave heating-dominated case

### C. Comparing theoretically estimated *n*_{e} and *T*_{e} with PIC simulations

The comparison between theory and PIC simulations on electron temperature and density is summarized in Fig. 10, where $ n i$ is used for comparison since the charge quasi-neutrality condition holds for the bulk plasma. In general, we see a good match between theory and simulations. The deviation on average is about 30%. Given the assumptions we made and the numerical uncertainty of PIC simulations, we think it is a good agreement. From the figure, we see the following trend of $ T e$ and $ n i$ based on the results of our theory (blue dots): For Coulomb collision dominated cases, the electron temperature decreases with the increase of beam energy, while the plasma density increases with the increase of beam energy. For wave heating-dominated cases, the electron temperature decreases a little bit with beam energy (since the ratio $ I 0/ v e b 2$ decreases a bit with beam energy) but increases with beam current, and the plasma density increases with the increase of both beam energy and current.

As mentioned in Sec. IV A, the discrepancy may result from the fact that the beam current is not sufficiently scattered. Another reason is that the beam velocity and the beam current are not a constant across the simulation domain, but decreasing as the beam passes through the bulk plasma, whereas here when calculating the electron temperature in Fig. 10, we assumed that they are constant. Note that even though we have a relative inaccurate $ T e$ estimate in Fig. 10(c), the density estimates are not significantly affected and remain showing a good agreement. This is mainly because the density estimate is determined mostly by the first two terms on the right-hand side of the density estimate in Eqs. (13) and (15).

## V. CONCLUSIONS AND DISCUSSION

In this paper, we develop an analytical model to estimate the effective electron temperature and the plasma density for a thermionic beam-generated plasma. The model could be simplified when the energy balance is dominated by heating due to e–e Coulomb collisions or waves and cooling due to loss to the anode. In practice, for both these two cases, $ W elas$ and $ W inelas$ can be neglected because they are much smaller. Coupled with the ion diffusion model, the plasma density and the electron temperature can be obtained as functions of the current density, electrode distance, pressure, and applied voltage, giving a good match with PIC simulations. Although we only considered two simplified scenarios, one can always simplify our theory when other terms are dominant, such as when $ W ioniz$ becomes dominant at high gas pressure, which is around $100 mTorr$ for our wave heating-dominated cases.

Although the usage of quasi-linear theory for wave heating is just an approximation, we believe that this is reasonable at equilibrium because the system at equilibrium should be in the weak turbulent regime without modulational instabilities^{13,44} and the strong transit-time damping does not need to be considered.^{46,52,55–57} In this case, the random phase assumption should be valid and a good agreement between PIC simulations and theory is indeed found in our work.

We notice that our model used several other assumptions, including an approximately Maxwellian distribution function for bulk electrons and a cathode without SEE. We also realize that we need validations against experimental results, such as the beam discharge devices with chemical reactions. These extensions will be left for future investigations.

## ACKNOWLEDGMENTS

The authors express their gratitude to Dr. Igor Kaganovich at the Princeton Plasma Physics Laboratory for the invaluable suggestions and are indebted to all the developers of the EDIPIC code. This work was supported by the National Natural Science Foundation of China (Grant No. 12305223) and the National Natural Science Foundation of Guangdong Province (Grant No. 2023A1515010762). This work was supported in part by the Swiss National Science Foundation.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Haomin Sun:** Conceptualization (lead); Data curation (lead); Investigation (lead); Methodology (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). **Jian Chen:** Conceptualization (lead); Funding acquisition (lead); Investigation (equal); Methodology (lead); Resources (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Guangyu Sun:** Writing – original draft (equal); Writing – review & editing (equal). **Liang Xu:** Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

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

## REFERENCES

*Nonlinear Wave and Plasma Structures in the Auroral and Subauroral Geospace*

*Principles of Plasma Discharges and Materials Processing*

*Collision Processes in Low-Temperature Hydrogen Plasmas*(Forschungszentrum, Zentralbibliothek Jülich, 2003), Vol. 4105.