Particle-in-cell simulations of electron beam-plasma discharge in a narrow gap with varying transverse boundary conditions

1.   Introduction

Plasma techniques are extensively utilized in surface treatments to enhance material properties, such as hardness and corrosion resistance [1–5]. In specific applications, achieving good plasma uniformity is essential for improving treatment quality. However, when treating the inner surfaces of long narrow tubes, conventional plasma discharge methods, like hollow cathode discharge [6], often encounter challenges due to sheath overlap caused by the tubes’ small inner diameters. This overlap complicates both the penetration of plasma into the tube and the maintenance of uniformity [7].

Recent studies have highlighted the effectiveness of electron beam-plasma discharge as a method to overcome these challenges [8–10]. For the inner surface treatment of tubes, an intense electron beam generated by an external plasma source is injected into the tube, creating an electron beam-generated plasma (EBP). Considerable effort has been dedicated to characterizing these EBPs. For example, a platform for the study of EBP properties was developed in reference [11], which included a plasma-cathode electron source and a biased long metal tube, 200 mm in length and approximately 10 mm in inner diameter [12]. Experimental findings indicated that variations in the tube’s inner diameter and bias voltage significantly influence the EBP parameters and their distributions [13]. Due to the lack of effective diagnostic methods, measuring particle transport and beam-plasma interactions in experiments remains challenging, thus highlighting the need for kinetic simulations. Therefore, our study aims to reveal the impact of transverse boundary conditions on such transversely confined EBPs through kinetic simulations.

Despite extensive numerical studies on electron beam-plasma interaction, few of them focus on the real electron beam-plasma discharges transversely confined within a narrow gap. Under such conditions, electron beam transport behaves differently as the system length becomes comparable to or even shorter than the beam energy relaxation scale. Our previous study demonstrated the generation of oblique waves and the variation of EBP properties with changes in beam parameters [14]. However, the effects of transverse boundary conditions have not been investigated yet with a comprehensive analysis of the transport process.

Therefore, in this paper, a 2D3V particle-in-cell/Monte Carlo collision (PIC/MCC) model is employed to investigate the effects of transverse boundary conditions on the EBPs confined in a narrow gap. The study utilizes kinetic simulations, incorporating electrons, ions, and neutral atoms. The collision processes including elastic scattering, excitation, ionization and resonance charge exchange are considered [15]. The primary focus of this work is on the effects of varying gap width and bias voltage on the properties of EBPs. By conducting phase space and spectrum analyses, we aim to gain insights into the characteristics of electron beam-plasma interactions in a long narrow gap. The results will provide effective recommendations for optimizing plasma surface treatment techniques in long, narrow tubes.

The paper will be organized as follows. Section 2 describes the employed PIC/MCC model. Section 3 presents the simulation results and discussions. Section 4 provides summaries and conclusions.

2.   PIC/MCC model

In this study, the electron beam-plasma discharge in a narrow gap is simulated using the EDIPIC-2D code [16]. EDIPIC-2D is an open-source electrostatic code designed for low-temperature plasma simulations. It has been benchmarked against various codes [17–19] and utilized for simulating the EBPs [14] and waves [20, 21]. The simulations are primarily modeled off the EBP platform in reference [12]. Instead of the cylindrical shape, we adopted a slab geometry (see figure 1) as the calculation domain. Such simplification may affect some transverse processes, e.g. the transverse instabilities that are important when an axial magnetic field is applied, which can be omitted for our cases. The longitudinal length of the simulation domain is fixed as L = 200 mm, while the half of transverse gap width D is varied from 5.0 mm to 8.7 mm. In the subsequent sections of the article, “the gap” will be used as a shorthand for “half of the gap”. A symmetric boundary condition is adopted for the left boundary (x = 0 mm) to reduce the computation time. Namely, only half of the real configuration is simulated and the other half is obtained from the symmetry along the plane of symmetry (x = 0 mm). Two metal walls biased with the negative voltage −U varying from 0 V to 500 V are located at the lower and right boundaries. Particles that hit these two walls are assumed to be fully absorbed. Secondary electron emission is turned off for the present simulations to simplify our simulation and focus on the process of beam-plasma interaction, as the contribution of SEE is complicated and it is worthy to be discussed separately in future work. An electron beam is uniformly injected from a grounded left injection electrode with a length of 3.3 mm and ionizes the background argon gas within the gap. The incident electron beam energy ε_b = 4 keV, beam current density J_b = 300 A/m² and working gas pressure p = 3 Pa are chosen from the experiments. Electrons will go through elastic scattering, excitation, and ionization collisions with a neutral gas, which are modeled using a null collision algorithm based on Monte Carlo methods [15]. The ion-neutral collisions are considered as resonance charge exchange in the simulations. The cross-section data are taken from references [22, 23].

Figure  1.  Schematic diagram of the simulation domain. D is the half of transverse gap width varying from 5.0 mm to 8.7 mm and L = 200 mm is the longitudinal length.

DownLoad: Full-Size Img PowerPoint

Simulations start from an initial plasma with the parameters: the initial plasma density n_p,0 = 10¹⁶ m⁻³, the initial electron temperature T_e,0 = 4 eV, and ion temperature T_i,0 = 0.5 eV. Cases with different initial parameters are performed and the results show that they have little influence on the profiles at steady state. The cell size ∆x = 55 μm for U = −300 V and −500 V cases or ∆x = 83 μm for others is used to resolve the Debye length (≈ 70 or 86 μm respectively) and the smallest collision-free path (≈ 0.047 m), and the time step ∆t = 7.8×10⁻¹² s or 1.2×10⁻¹¹ s is used to resolve the electron plasma frequency (ω_e∆t ≈ 0.1 < 1) and the smallest collision period (≈ 1.8 ns). The particle weight is fixed to be 1.4×10⁵. Initially, 50 macro-particles per cell are created for each species. The total macro-particle number increases until the steady state is attained. Namely, the balance between the produced ions and lost ions is achieved. Each simulation takes 7–14 days for 45 μs simulation time with 240 cores in Sugon high-performance cluster.

3.   Results and discussions

In this section, we provide the discussions on the profiles of beam and plasma parameters, electron phase space, and wave spectrums varying with the bias voltage U on the transverse wall (x = 5 mm) and the gap width D. All the presented data are taken from the steady state, and the profiles at the negative x range are obtained using the plane of symmetry at x = 0 mm.

3.1   Typical density and potential profiles

Profiles of beam density, ionization-produced electron density, and electrostatic potential for the typical case with J_b = 300 A/m², ε_b = 4 keV, p = 3 Pa, U = −100 V, and D = 5 mm are depicted in figure 2. In figure 2(a), it can be observed that at y < 100 mm, beam electrons are focused towards the center due to the pinch electric field created by the bias voltage. Consequently, the produced plasma also exhibits a converged shape at y < 100 mm, as shown in figure 2(b). At y > 100 mm, the coherent wave becomes oblique, eliminating the transverse pinch electric field, which weakens the further focusing of the beam electrons and instead induces the oscillation of the beam density (figure 2(a)). The produced plasma impinges on the end wall at y = 200 mm and forms a bell structure, similar to that observed in free-burning arc discharge. In this region, the bulk plasma density increases by ionization, likely because the beam electrons undergo cycles of scattering and convergence by the oblique wave (see figure 2(a)), thereby experiencing more collisions with neutral gas [14]. The effects of the wave excitation and beam deflection on the plasma parameters will be further discussed later in this paper.

Figure  2.  (a) 2D beam electron density profile, (b) 2D electron density profile, and (c) 2D plasma potential profile. The operating conditions are J_b = 300 A/m², ε_b = 4 keV, p = 3 Pa, U = −100 V, and D = 5 mm.

DownLoad: Full-Size Img PowerPoint

3.2   Plasma profiles with varying bias voltages U and gap width D

The profiles of electron and ion densities for cases with varying bias voltage U and gap width D are shown in figure 3. Due to the existence of the bias voltage on the transverse walls, the Child-Langmuir high voltage sheath forms. As shown in figure 3(a), the Child-Langmuir sheath thickness s increases with −U, qualitatively consistent with the analytical prediction [24]. In the bulk region (x < 1 mm) shown in figures 3(a) and (b), density oscillations are clearly observed, indicating that some transverse waves are excited. These waves might be ion acoustic waves generated by the sheath-induced instability [25, 26] and require further investigation.

Figure  3.  Plasma density profiles for the cases with (a) different bias voltages U and fixed D = 5 mm and (b) different gap widths D and fixed bias voltage U = −100 V along the x-axis at y = 100 mm, and (c) different bias voltages and fixed D = 5 mm and (d) different gap widths D and fixed U = −100 V along plane of symmetry at x = 0 mm. The operating conditions are J_b = 300 A/m², ε_b = 4 keV and p = 3 Pa.

DownLoad: Full-Size Img PowerPoint

As shown in figure 3(b), as the gap width D increases from 6.2 mm to 8.7 mm, the sheath edge gradually moves outward while the thickness remains unchanged (s ≈ 1.5 mm). With the fixed sheath thickness and bias voltage, the Child-Langmuir sheath flux remains almost the same (Γ_i ~ ⎜U⎜^3/2s⁻²). Therefore, varying the gap width shows limited effects on the bulk plasma density, as seen in figure 3(b).

Figure 3(c) shows that for the case of no bias voltage (U = 0 V), plasma density monotonically increases along the plane of symmetry. When the bias voltage is applied, n_e exhibits a “V-shaped” profile with lower values in the middle, which aligns with the experimental measurements [12]. The density hump at small y occurs because beam electrons are focused and concentrated, leading to enhanced ionization. Therefore, as −U increases, the density gradient increases due to the stronger beam deflection. In contrast, as depicted in figure 3(d), increasing the gap width can smooth out the density gradient and reduce the density hump at small y because the pinch field is weakened.

The potential profiles along the plane of symmetry for the cases with varying U and D are shown in figures 4(a) and (b), where coherent wave structures are clearly observed. As shown in figure 4(a), the increase of −U results in lower plasma floating potential and lower wave amplitude. As the bias voltage increases to above 300 V shown in figure 4(a), the resulting beam deflection suppresses the potential oscillations. Since the total density changes little between 0 V and 100 V, the wavelengths remain almost unchanged. As illustrated in reference [27], the wave number of the streaming waves generated in a finite beam-plasma system can be calculated by

Figure  4.  Plasma potential profiles for (a) different bias voltages U and fixed D = 5 mm and (b) different gap widths D and fixed bias voltage U = −100 V along the plane of symmetry at x = 0 mm. The operating conditions are J_b = 300 A/m², ε_b = 4 keV and p = 3 Pa.

DownLoad: Full-Size Img PowerPoint

Here, ω_e is the electron frequency, v_b is the beam velocity and $L_{\text{n}}=L\omega_{\text{e}}\text{/}v_{\text{b}}$ is a normalized length. With the fixed beam energy and longitudinal length, Re(k) ~ ω_e ~ n_e^1/2. Figure 4(b) shows that the ratio of wavelengths for the cases with D = 5.0 mm and D = 8.7 mm is about $\lambda_{\text{5.0 mm}}/\lambda_{\text{8.7 mm}}\approx \text{15/12.5}\approx\text{1.2}$, consistent with the ratio of their densities ${\left({{n}}_{\text{e}\text{,}\text{5.0 mm}}\text{/}{{n}}_{\text{e}\text{,}\text{8.7 mm}}\right)}^{{-1/2}}\approx \sqrt{\text{3/2}\text{}}\text{}\approx \text{}\text{1.22}$. Moreover, one can see that the wave amplitude shows a non-monotonic variation with the gap width (increases from D = 5.0 mm to 6.2 mm but decreases from D = 6.2 mm to 8.7 mm), suggesting the complex saturation mechanism of the streaming waves in such a system. Such a phenomenon is still not clear and needs to be investigated in future studies.

The electron temperature profiles along the plane of symmetry for the cases with varying U and D are shown in figures 5(a) and (b). The increase of −U leads to the significant increase in electron temperature at y < 100 mm, indicating that the generation of the pinch electric field leads to strong electron heating. The underpinning mechanism remains unclear and requires further investigation. Moreover, it can be observed that the electron temperature is about 5 eV, which is lower than the ionization energy for argon. Consequently, the major contribution to the ionization processes comes from beam electrons rather than plasma electrons. Given their significant role, it is important to investigate the porperties of beam electrons in the tube.

Figure  5.  Electron temperature profiles for (a) different bias voltages U and fixed D = 5 mm and (b) different gap widths D and fixed bias voltage U = −100 V along the plane of symmetry at x = 0 mm. The operating conditions are J_b = 300 A/m², ε_b = 4 keV and p = 3 Pa.

DownLoad: Full-Size Img PowerPoint

In summary, the simulation results indicate that increasing −U leads to V-shaped profiles which mean larger non-uniformities along the plane of symmetry, while decreasing −U or increasing D results in a stronger excitation of oblique waves. Based on these results, in order to realize a surface treatment of high uniformity, we recommend using a relatively narrower gap width tube and applying an appropriate bias voltage. These operations can prevent the formation of a V-shaped profile while also suppressing waves.

3.3   Electron phase space

To highlight the effects of the pinch electric field on the beam dynamics, beam electron density distributions for cases with different bias voltages are displayed in figure 6. As shown in figures 6(a)–(d), electrons near the edge of the beam are strongly deflected, thereby forming the envelope oscillation. Owing to the shield of the electrons near the edge, those near the beam axis are less affected. These central electrons maintain relatively good coherence and can further interact with plasma, exciting the waves downstream (y > 100 mm). With the increase of the bias voltage, more and more beam electrons are deflected, and therefore, waves are suppressed.

Figure  6.  Beam electron density profiles for the cases with different bias voltages. (a) U = −0 V, (b) U = −100 V, (c) U = −300 V and (d) U = −500 V. The operating conditions are J_b = 300 A/m², ε_b = 4 keV and p = 3 Pa.

DownLoad: Full-Size Img PowerPoint

Figure 7 presents the phase space y-v_x and y-v_y for beam electrons under different bias voltages, where v_x and v_y are the transverse and longitudinal velocities, respectively. In the absence of bias voltage, both v_x and v_y primarily oscillate downstream (see figures 7(a) and (d)), indicating that beam electrons are strongly scattered by the oblique wave. For the case with U = −300 V, as shown in figures 7(b) and (e), the coherent oscillations of v_x and v_y disappear. Instead, three-band structures characterized by the sudden rise of v_x and fall of v_y become pronounced, corresponding to the deflection of beam electrons by the pinch electric field. When the bias voltage rises to U = −500 V, the pinch electric field becomes stronger, causing the decrease of effective betatron length (i.e. the length of the single beam envelope). Consequently, two more bands are observed.

Figure  7.  Space phases y-v_x ((a)–(c)) and y-v_y ((d)–(f)) under different bias voltages. (a) and (d) U = 0 V, (b) and (e) U = −300 V, (c) and (f) U = −500 V. The operating conditions are D = 5.0 mm, J_b = 300 A/m², ε_b = 4 keV and p = 3 Pa.

DownLoad: Full-Size Img PowerPoint

Figure 8 shows the phase space of the beam electrons under different gap widths. Given the relatively small bias voltage U = −100 V, the band structures are not pronounced. At the same time, as shown in figure 8(a), when y > 0.15 m, the band also oscillates, indicating the complex relationship between coherent wave excitation and the beam deflection process, and that coherent wave excitation also impacts the beam deflection. Additionally, it can be seen that when the gap width increases from 5.0 mm to 6.2 mm, the band structure in the v_x direction disappears, and the v_x oscillations recover first. The oscillations in the v_y direction do not appear until the gap width reaches 8.7 mm, as shown in figure 8(f). This suggests that the oscillations in the v_y direction are more sensitive to the beam deflection compared with v_x direction.

Figure  8.  Space phases y-v_x ((a)–(c)) and y-v_y ((d)–(f)) under different bias voltages. (a) and (d) D = 5.0 mm, (b) and (e) D = 6.2 mm, (c) and (f) D = 8.7 mm. The operating conditions are U = −100 V, J_b = 300 A/m², ε_b = 4 keV and p = 3 Pa.

DownLoad: Full-Size Img PowerPoint

3.4   Wave spectra

In this subsection, wave spectra are analyzed by applying the fast Fourier transform (FFT) to the potential profiles within the space [x = 0 mm, 2 cm < y < 17 cm] and the 50 ns time interval. The typical wave spectrum and its schematic diagram are plotted in figure 9. As shown in figure 9(a), three types of wave branches can be categorized:

Figure  9.  (a) Wave spectrum under 1 keV beam energy and (b) schematic diagram of the spectrum.

DownLoad: Full-Size Img PowerPoint

(1) Electron beam mode whose frequency is $\omega\mathrm{_b}$ (red line in figure 9(b)) with a linear dispersion relation

where $k$ is the wave vector and ${v}_{\mathrm{b}}$ is the beam velocity calculated as $v_{\mathrm{b}}=4.19\times10^5\times\sqrt{E_{\mathrm{b}}(\mathrm{eV})}\mathrm{\ m}/\mathrm{s}=1.32\times10^7\ \mathrm{m}/\mathrm{s}$ for the beam whose energy is $E_{\mathrm{b}}\ =\ 1\ \mathrm{keV}$ beam, matching the slope of the first wave in figure 9(a).

(2) Langmuir wave whose frequency is $\omega$ marked by the blue line in figure 9(b) with the dispersion relation

where $v_{\text{th}}\approx4.19\times\text{10}^{\text{5}}\ \text{m/s}$ is the thermal speed of the electron and $\omega_{\mathrm{p}}\approx\text{10}^{\text{10}}\ \text{rad/s}$ is the electron plasma frequency. Since the electron density and electron temperature are not uniform, broad spectra are seen in figure 9(a) for this branch.

(3) Electron acoustic wave whose frequency is $\omega\mathrm{_{ea}}$ marked by the green line with the dispersion relation [28]:

Here, ${{ \lambda }}_{\text{Db}}$ is the Debye length for beam electron. Electron acoustic modes are known to be excited in plasma with cold and hot electrons. As shown in figure 7, beam electrons show a broadened velocity range, indicating that they are heated and serve as the hot component that supports the electron acoustic modes. In the long wavelength limit ($k\mathit{\lambda}_{\text{Db}}\ll$ 1), ${{ \omega }}_{\text{ea}}\approx{{ \omega }}_{\text{p}}{{ \lambda }}_{\text{Db}}{{k}}\propto{{k}}$, matching the Fourier spectrum of the third branch in figure 9(a). It is worth noting that typically the electron acoustic mode with high k has a negative growth rate, making them unobservable [29].

In addition, one can see from figure 9(a) that the intersections of the spectra of different modes show very high intensity (deep red regions), suggesting the reactive couplings between different waves. Specifically, as a negative energy wave, electron beam mode can react with the Langmuir wave (a positive energy wave), leading to the growth of both these two waves [30]. A similar mechanism applies to the coupling between electron beam mode and electron acoustic mode [31].

Figure 10 presents the wave spectra for the cases with different bias voltages. As can be observed, varying the bias voltage has a minimal effect on the beam mode and electron acoustic mode. Since the range of electron density expands and electron temperature increases with the bias voltage (see figures 3 and 5), a higher frequency of Langmuir wave can be reached under a large bias voltage. A broader frequency range is also observed in the cases with high bias voltages. Wave spectra are scarcely changed with varying gap widths and thus not plotted here.

Figure  10.  Wave spectra along the plane of symmetry for different bias voltages. (a) U = −100 V, (b) U = −300 V, (c) U = −500 V. The operating conditions are J_b = 300 A/m², ε_b = 4 keV and p = 3 Pa.

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

In this study, we investigated the effects of transverse boundary conditions, specifically the bias voltage on the transverse wall and the gap width, on the properties of electron beam-generated plasmas confined in a narrow gap using the PIC/MCC method. Such plasmas are used for the inner surface treatment of tubes in practical applications. The key findings are summarized as follows:

(1) Beam deflections are observed when a bias voltage is applied to the transverse wall, leading to the formation of a beam oscillation envelope and the band structures in beam electron velocity space.

(2) Beam electrons near the plane of symmetry are less affected by the pinch electric field due to the shield of the edge beam electrons. These central electrons continue to interact with the produced plasma, leading to the excitation of oblique waves. Three branches of electrostatic waves, i.e. the electron beam mode, the Langmuir wave, and the electron acoustic mode are identified. Increasing the bias voltage would expand the range of the electron density, therefore broadening the frequency range of the Langmuir wave.

(3) Both the beam deflections and oblique wave excitation significantly modify beam electron transport, thus causing the plasma non-uniformity. Increasing the bias voltage and decreasing the gap width enhance the deflections of the beam by the electric field, thereby lifting the electron temperature. Additionally, as the bias voltage increases, more beam electrons are deflected, resulting in the suppression of waves.

It is worth noting that secondary electron emission and the cylindrical effects are neglected in this study, which may cause some changes in the beam deflection process and will be investigated in our future work.