# Mathematical Modelling of RF Plasma Flow at Low Pressures with 3D Electromagnetic Field

Authors: A. Yu Shemakhin; V. S. Zheltukhin
Journal: Advances in Materials Science and Engineering
 (2019)

---

## Abstract

In this study, a hybrid mathematical model of a low-pressure RF plasma jet in transition mode between continuum and free molecular flow at a Knudsen number of 8·10−3 ≤ Kn ≤ 7·10−2 for a carrying gas is described. The model takes electrons, ions, metastable atoms, and potential and curl electromagnetic fields into account. The model is based on both a statistical approach for the atoms in the ground state and a continuum model for other components. The results of plasma flow calculations in an undisturbed jet are described. The distributions of the electrodynamic and electrostatic parts of the electric field are given. It has been observed that the plasma jet has a layered structure along the stream: a positive charge region is formed at the beginning of the jet, followed by a negative charge region, and then a positive one again. The reason for the formation of a layered structure is the fast flow expansion when the plasma inflows into the vacuum and the difference in electron and ion pulse.


---

## Body

## 1. Introduction

RF plasma at low pressures (0.15–150 Pa) is used for the modification of various materials such as metals, polymers, and composites [1–5]. Plasma has the following properties: electron density (ne) ranging from 1015 to 1018 m−3, ionization degree ranging from 10−6 to 10−3, electron temperature (Te) ranging from 1 to 4 eV, and atom and ion temperature (Ta) ranging within (3-4)·103 K in the plasma torch and within (3.2–10)·102 K in the plasma stream [6].Low-pressure RF plasma has specific features. Firstly, experimental results [6] showed that the plasma jet is both an inductive and capacitive coupled RF discharge because the electron density in the plasma jet is multiple orders higher than that in decaying plasma, and both axial and azimuthal components of the magnetic strength and plasma current are detected. Hence, all 6 components Ex,Ey,Ez,Hx,Hy,Hz of the electromagnetic field should be considered in common cases. Secondly, the neutral atom flux is in the transient mode between continuous and free molecular flow because the Knudsen parameter (Kn) ranges between 8·10−3 and 7·10−2 [1] while the charged components and metastable atoms satisfy the continuous media flow. The Navier–Stokes equations for carrier gas are not applicable in this case.Therefore, particle-in-cell (PIC) method is used for simulation of the inductive coupled radio-frequency (ICRF) discharge plasma to analyze the mechanism of power absorption in the pressure range from 1.33·10−2 to 13.3 Pa [7, 8].Maxwell’s equations are usually reduced to vector and scalar potential equations [9] or reduced to a wave equation [10] or analyzed in conditions of neglecting the radial components of the magnetic field [11]. Poisson’s equation is used for calculation of the electric field strength in [7, 8]. 2D Maxwell’s equations for ICRF discharge in the pressure range from 13.3 to 133 Pa and are reduced to a system of elliptical equations for the modulus, phase and angle functions of the magnetic field, and electric strength vectors in [6].Different models of jet inductive coupled RF discharge at low pressure have been considered in our previous studies [12–16]. Direct Monte Carlo simulation (DSMC) [17] is used for simulation of carrier gas flow in these studies. The dynamics of a carrier gas was considered in [12, 13], taking into account the internal distributed heat source due to the transfer of energy from the electron gas. The pressure, the temperature, and the velocity of carrier gas both in the undisturbed flow and in the flow with a sample are calculated. As a result of numerical experiments, the effect of upwarming of the peripheral region of the flow near the inlet was discovered and experimentally confirmed [14]. The cause of gas upwarming is jet dumping due to a collision with a stationary gas in the vacuum chamber.It is known that excited atoms (metastable) are essential for the balance of particles and energies in RF discharges in argon [18–22]. The spatial distribution of metastable particles in the jet of an inductive coupled RF at low pressure was studied in [15]. It was found out that the metastable density nm is two orders higher than the density of electrons, whereas their spatial distributions correlate with each other everywhere except in the vicinity of the inlet about 0.15 m along the stream and 0.1 m by diameter. Maximum nm is reached at 0.09 m from the inlet, whereas maximum ne is at the inlet.The mathematical models considered in [12–15] are not self-consistent because the distribution of the electromagnetic field was approximated by the experimental data from [6]. Distribution of the RF electromagnetic field in the low-pressure inductive coupled RF plasma jet is calculated in [16] using the Biot–Savart law.In addition, studying [12–16] is performed under the assumption of quasi-neutral plasma flow. It is known that a positive sheath is formed near the electrodes in capacitive coupled RF discharges. The reason for the sheath formation is fast oscillations of electrons referred to slow ions [23]. Because the jet of a low-pressure RF plasma contains a factor of capacitively coupled discharge, it is interesting to investigate the effect of the capacitive component on the distribution of charged particles and the electromagnetic field in the plasma. In this regard, the present paper aims in constructing a model of low-pressure ICPRF plasma to analyze the role of electrons and ions in sustaining the discharge, taking into consideration the interaction with electromagnetic field.


## 2. Mathematical Model of Low-Pressure RF Plasmas with Electromagnetic Field

As mentioned above, a carrier gas flows in low-pressure RF plasma in a transition mode between continuum and free-molecule flux. Navier–Stokes equation is not correct at Knudsen’s parameter Kn > 0.01 [24, 25]. On the contrary, charged particles satisfy continuity because their movements are controlled not only by atom flow but also by long-range Coulomb force that prevents charge separation [26]. Electrical field does not have an effect on excited atoms, but atom exciting and metastable quenching results are due to collision with electrons so that their movements are controlled by both ground-state atoms and electron collisions. Hence, hypothesis of continuity is also valid for metastable atoms.We neglect Hall’s effect, electron pressure gradient, radiation energy loss, electron attachment, multiply charged ions, and ions slipping. Direct electron impact and step ionization are assumed as the basic mechanisms of charged particle formation. Penning’s ionization, triple recombination, and photo recombination are also considered.We assume that plasma is generated in a discharge tube and then is passed to a vacuum chamber. A model of the plasma flow in the vacuum chamber is considered.Let the radius and length of the cylindrical vacuum chamber be denoted byRvk and Lvk, respectively, the radius of the vacuum chamber inlet denoted by Rrk, and the subscripts inlet, outlet, and walls be used for parameter values of inlet, outlet, and walls of the vacuum chamber, respectively.The model involves the following initial and boundary value problems:(1)
Boltzmann’s transport equation for neutral atoms:(1)∂f∂t+c⋅∂f∂r+F˜⋅∂f∂c=Sf,t>0,r∈Ω,c∈R3,fc,r,tr∈Ω,t=0=fc,r,tr∈∂Ω,t=0=f0c.(2)
The equation of electron continuity:(2)−∇⋅De∇ne+μeEcap−vane=R1nena+R2nm2+R3nenm−R4neni−R5ne2+ni,neinlet=ne,in,neoutlet=0,newalls=0.(3)
The equation of ion continuity:(3)−∇⋅Di∇ni−μiEcap+vani=R1nena+R2nm2+R3nenm−R4neni−R5ne2+ni,niinlet=ne,in,nioutlet=0,niwalls=0.(4)
The equation of metastable continuity:(4)−∇⋅Dm∇nm−vanm=R6nena−R2nm2+R3nenm−R7nm−R8nmna−R9nm,nm∂Ω=0.(5)
Maxwell’s equations for curl component of electromagnetic field in inhomogeneous wave view [23]:(5)Δ−μ0ε0∂2∂t2Er,t=μ0∂∂tjr,t,Er,tr∈∂Ω=−iωμ04π∫jselr−sd3s.(6)
Irrotational electrical field:(6)Ecap=−∇φ,Δφ=eε0ne−ni,φwalls=φoutlet=0,φinlet=φin,here, va and r are vectors of the velocities and coordinates of ground-state atoms, respectively, c is the vector of velocities in 3d Euclidean space R3, Ω is the domain of the vacuum chamber, fc,r,t is the velocity distribution function of ground-state atoms, f0c is Maxwell’s velocity distribution function, Sf is the collision integral, F˜=−1/ma∇WT is the reduced force field which has effects on the ground-state atoms at elastic collisions with electrons, t is time, ∂Ω is the Ω boundary ∂Ω=∂Ωinlet∪​∂Ωoutlet∪​∂Ωwalls, ne,ni,nm are electron, ion, and metastable density, respectively, na is the ground-state atom density, ne,in is a the prescribed value of the electron density on the inlet, De,Di,Dm are the electron, ion and metastable diffusion, νc is the elastic collision frequency of electrons and atoms, Ecap is the strength of the irrotational electrical field, φ is the electrical potential, e is the electron charge, ε0 is the electric constant, E is the curl electric field strength, Efull = Ecap+E,E=Efull, d3s is the elementary volume, R1 is the impact ionization rate, R2 is the Penning ionization rate, R3 is the step ionization rate, R4 is the photorecombination rate, R5 is the triple recombination rate, R6 is the metastable excitation rate, R7 is the radiative recombination rate, R8 is the collisional quenching rate, R9 is the electron de-excitation rate, jel=σEfull is the plasma electron current, and σ is the plasma conductivity,(7)σ=nee2νcmeνc2+ω2,var,t=∫−∞∞cfc,r,tdc,Tar,t=ma3kBc¯2−va2,c¯2r,t=∫−∞∞c2fc,r,tdc,pr,t=kBNTa,where me is the electron mass, ω=2πf is the cyclic frequency, and f is the generator frequency.During the elastic collisions, electrons impart energy to atoms:(8)Ec=32kBδνcneTe−Ta,where(9)δ=me2ma,Te=e2E23kBmeδνc2+ω2.The specific power of the distributed heat source can be written as(10)WT=∫VcellEcdVdt,where dV is a volume element of the mesh cell.It is assumed that transport coefficients and kinetic rates are functions of the reduced electrical fieldE/N [15, 23, 27–35]. The boundary conditions for equation (6) are set by the Biot–Savart law. The mobility of electrons and ions are calculated as(11)μe=emeνc,μi=emaνia,where νia is the collision rate between ions and neutral atoms [23].The direct simulation Monte Carlo (DSMC) modelling [17, 36–38] is used for the numerical solution of the kinetic equation (1). The method is based on the splitting of Boltzmann’s equation (1) on the movement and collision processes that allow us to describe the gas-dynamic processes in the transition mode for a neutral environment [17]. Therefore, both modification of Bird’s method to the plasma flows and concordance results with the continuous models of other particles are required. The software package using libraries of OpenFoam open source CFD software [39] was created for calculating the plasma parameters.


## 3. Methods and Numerical Experiments

The calculations of charged particles, metastable atoms, and electromagnetic field was performed in the vacuum chamber atRvk = 0.1 m, Lvk = 0.3 m, and Rrk = 0.012 m. Flow input parameters are the following: the plasma forming gas is argon, gas flow rate G is 0.12–0.24 g/s, pressure Pinlet is 35–85 Pa, the temperature Tinlet is 400–600 K, and the degree of ionization is 10−4. The initial pressure in the vacuum chamber P0 is 3.5–8.5 Pa. The model parameters correspond to the experimental setup [1, 6]. As shown by calculations, steady state flow is established for t = 10−2 s [8].Velocity, temperature, and pressure of the carrier gas are closed to earlier works [12–15] because the ionization ratio is too small for the charged and metastable particles to have a significant effect on the atoms in the ground state. The vector plot on the left side of Figure 1 shows the real part of electrodynamical part E of the electrical field at 0.02 m from the inlet. One can see that vectors of electrical tense are winded around the chamber axis. On the right side of Figure 1, curves of ReEz,ImEz, and potential φ of electrostatic part of the electrical field along the flow axis is showed. Based on exponential law, the electrodynamical part of the electrical field is decreased approximately just as the electrostatic potential is decreased to a weaker level with distance from the inlet.Figure 1
Transverse section of curl electric field in ICRF plasma jet at z = 0.02 m from inlet (left side) andReEzblue curve,ImEzred curve,andφgreen curve along the flow axis (right side).
(a)
(b)Spatial distribution of electron and ion densities is shown in Figure2. On the left side, the longitudinal cross section of the electron density slice is presented. It is apparent that the plasma jet is almost uniformly diffused in the width direction. On the right hand of Figure 2, electron and ion densities along the flow axis are showed. Both positive and negative regions in the plasma flow are observed. At the beginning and ending of the flux, ion density is greater than electron density, whereas in the middle, the electron density exceeds the ion ones. The reasons of these phenomena are the fast flow expansion after inflow into the vacuum chamber [40], bounded ionization and recombination rates, and a high flow rate. Ion pulse is more than 7·104 times greater the electron ones, whereas transit time of a particle through the vacuum chamber is less than 10−3 s. Therefore, an ion fails to adapt to electron density during flowing through a vacuum chamber.Figure 2
The distribution of the electron density in the ICRF plasma flow in longitudinal cross section (left side) and distributions of electron (red curve) and ion (black curve) densities along the flow axis (right side).
(a)
(b)


## 4. Conclusion

A hybrid model of the low-pressure RF plasma flow in transition mode between continuum and free molecular flow is constructed. Electrons, ions, metastable atoms, and potential and curl electromagnetic fields are taken into consideration. The model is based on both a statistical approach for the atoms in the ground state and a continuum model for other plasma characteristics. Calculations of plasma parameters in an undisturbed jet at local approximation of transportation coefficients and kinetic rates are performed. The results show that hard ionization nonequilibrium is observed in the plasma stream. The reasons of these phenomena are the fast flow expansion in the vacuum chamber, bounded ionization and recombination rates, and a high flow rate. Based on exponential law, the electrodynamical part of the electrical field is decreased approximately along the jet just as the electrostatic potential is decreased to a weaker level with distance from the inlet.



---

*Source: 7120217-2019-01-16.xml*