Analytical model for estimating the equilibrium plasma quantities in an electron beam–plasma system

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.³⁰ 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.³² 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,²⁷ 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.

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:

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

2. Heating due to the ionization by the electron beam.

3. Wave heating through the wave–particle interaction.

Three cooling sources can be important:

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

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

3. 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.

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

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

The obtained energy balance of the considered beam–plasma system will be validated against the numerical simulation results in Secs. III and IV.

We perform one-dimensional PIC simulations of an electron beam–plasma interaction system using the PIC code, EDIPIC,⁴⁷ 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.⁴⁸ EDIPIC incorporates the Boris scheme for particle pushing, and the field is solved using the Thomas algorithm.⁴⁹ 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 − 3$⁠, $T e 0=1 eV$ for Coulomb collision heating-dominated cases and $n e 0= 10 17 m − 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 $Δt=1× 10 − 3 ns≈1.12× 10 − 3 ns=0.2 ω p e 0 − 1$⁠, $Δx=1.13× 10 − 3 mm$ $=1/2 λ D e 0$⁠, while the time steps and cell sizes of wave heating-dominated cases are $Δt=1× 10 − 2 ns≈1.12× 10 − 2 ns=0.2 ω p e 0 − 1$⁠, $Δx=7.84× 10 − 3 mm=1/2 λ 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× 10 22 s − 1 m − 2$ (see Fig. 5). Assuming the SEE coefficient is $γ S E E=0.1$ (reasonable assumption, see, for instance, Ref. 50) and suppose the SEE flux $F S E E=1.5× 10 18 s − 1 m − 2$⁠. Comparing to the electron beam flux $F e b= n e b v e b=1.559× 10 23 s − 1 m − 2$⁠, the SEE flux is much smaller.

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.

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., $∂ Γ i ∂ 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.

An important assumption of our model for wave heating is the $I∼ I 0 ( k / k 0 ) − 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≈1.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 $ν e e= n e σ 90 v e b≈7.1× 10 8 s − 1$ and $ν e e=6.7× 10 6 s − 1$⁠, where $σ 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 $ν e n=8.7× 10 7 s − 1$ and $ν e n=8.7× 10 6 s − 1$⁠, respectively. The beam energy relaxation length for Case 2 can be estimated as $l relax , e e= v e b/ ν e e≈4.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⁡ ϕ 0 − 1/2γ≈1.5 mm$⁠, where $γ=π ω p e n e b/ n e ( v e b Δ v e b ) 2$ is the growth rate of the two-stream instability and $ϕ 0≈1/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=∫( W e e + W ioniz + W wave$ $− W elas − W inelas)dx=4.59× 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× 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 $Φ( v e b / v T e)≈ 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,¹² 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.

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.

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.

The anode sheath voltage drop can be calculated as $ϕ= T e/2 ln( m i / 2 π m e)$ in Eq. (8) so that $v min$ can be estimated as $v min= T e m e ln ⁡ ( m i 2 π 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.

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).

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.

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.

The authors have no conflicts to disclose.

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).

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