Numerical solution of the model problem of CCRF-discharge at atmospheric pressure

This work describes a 1D mathematical model of capacitive coupled RF discharge between symmetrical electrodes in argon at atmospheric pressure in a local approximation. Electrons, atomic and molecular ions, metastable atoms and argon dimmers as well as ground-state atoms are considered. A simplified diagram of argon excited states when two metastable and two resonance states are replaced with the uniform level. Such diagram is frequently used to simulate argon plasma due to efficient mixing of these layers at electron impacts. Velocity factors of electron impact processes were calculated using Boltzmann equation with a glance to electron-electron collisions. This work describes numerical algorithm of mathematical model implementation, which is based on finite-dimensional approximation of the problem using difference schemes together with iteration process. The software was developed to implement iterative processes using MatLab. Characteristics of atmospheric pressure capacitive coupled RF discharge at interelectrod distance 20 mm are calculated.


Introduction
Low-temperature plasma is widely used for various spheres of science and technologies. Plasma is created by various types of gas discharges. An important role is assigned to high-frequency discharges, and particularly, capacitive coupled radio-frequency (CCRF) discharges. Argon is often used as a plasma gas [1][2].
Experimental and calculation methods, which complement each other, are used to tie internal and external parameters of the discharge and allow solving a lot of problems of physics and chemistry of lowtemperature plasmas.
In our earlier works we found that, as a rule, all investigations, which are dedicated to simulation of CCRF-discharge in argon, considered that plasma consists of electrons, atomic ions, metastable atoms and ground state atoms. Dependencies of molecular and atomic ions densities in CCRF argon plasma at atmospheric pressure on gas temperature is researched in the article [3].
It is found that molecular ions prevail at the temperatures near to 500 К, and density of atomic ions increases at the temperatures equal to 1500 К and above. This fact is the reason for including molecular ions and dimmers to kinetic diagram of argon in this study. Moreover, earlier the authors tested different types of boundary conditions as for balance equations of charged particles and for balance equation of gas temperature. [4,5] Self-consistent 1D model of CCRFdischarge in argon at atmospheric pressure is offered because of this investigation.
The authors have developed approximate algorithm of numeric implementation of the model; this algorithm is based on finite-dimensional approximation together with using of iterative method. It should be noted that different iteration methods intended to solve non-linear problems including partial derivatives, were offered earlier (see, for example, [6][7][8][9][10] at al.). But the problem studied herein has a range of peculiar features, and particularly, availability of different time scales of change of the main characteristics of the steady state of CCRF-discharge at low pressure. Moreover, gradients of charged particles density and electric field intensity, electron temperature in near-electrode layers at the boundary of computational domain is the specific feature of the problem too.

Setting the problem
Here we study CCRF-discharge between two parallel plate electrodes, one of which is grounded and another is connected to CCRF-generator, and distance between electrodes is less then size of the electrodes. Under such conditions electric field is close to potential one and discharge is uniform along electrodes thereby we may use 1D model, which allows considering kinetics of plasma-chemical processes in diffusive-drift performance.
Assessment of time and distance, required for electrons to lose energy obtained from field, shows that to simulate CCRF-discharge at atmospheric pressure, we can use local approximation, under which the parameters of plasma electronic component (diffusion, mobility as well as moderate energy, velocity of plasma and chemical reactions, etc.) depend on local value of the reduced electrical field E/N, where E is electric field intensity, N is density of ground states atoms [11].
The following features of CCRF-discharge is considered in the model: presence of domains without quasi-neutrality, fast time change of applied voltage, as well as availability of processes involving metastable atoms, molecular ions and dimmers. Efficient mix of four lowest electronically excited states allows replacing them with a unified layer [12,13]. Let us assume that b is the distance between electrodes, the grounded electrode locates at the point x = 0 and the charged one locates at x=b (axis Ox is directed perpendicular to the surface of the electrodes).
Processes, occurring inside of CCRF-discharge at atmospheric pressure, can be described by the following initial-boundary value problems and Cauchy problem.

Convection-diffusion equation for atomic ions
where ne, n+, n2+ are densities of electrons, atomic and molecular positive ions, respectively, nm is density of metastable atoms, G+=-D+dn+/dx + μ+n+E is density of atomic ion flow, μ+ is mobility of atomic ions, D+ is diffusion coefficient of atomic ions, E=-dφ/dx is intensity of electric field, φ is electric field potential.
Boundary conditions at x=0 shall be chosen as: and at x=b -shall be chosen as: and where m+ is atomic ion mass, k is a Boltzmann constant. We assume that temperatures of atoms, atomic ions and excited atoms are equal Ta.

Convection-diffusion equation for ionized dimers
where n2* is density of argon dimers, μ2+ is mobility of molecular ions, D2+ is diffusion coefficient of molecular ions, G2+=μ2+n2+E-D2+dn2+/dx is density of molecular ion flow. By the same way as above, let us set the following boundary conditions: Boundary conditions at x=0 shall be chosen as: and at x=b -shall be chosen as: and Here m2+ is mass of atomic ion.

Balance equation for metastable density
Boundary conditions for this equation are: Let us assume that electrode surface is catalytic, when metastable atoms decay completely on it, then nm(x,0)=0.

Kinetic equation for argon dimer density
Let us take zero initial conditions as well as for metastable atoms: n2*(x,0)=0.

Kinetic equation for neutral atoms
Since fluctuations of atomic temperature in scope of average values taken within the period of field change are negligibly small, the balance equation for atomic temperature can be considered as average temperature within the period of variation of electric field, identified with frequency generated by a high-frequency generator. Moreover, we consider that coefficient of energy transfer from ions to atoms equals to 1.

d(-λadTa/dx)/dx=[jiE]+[QNne],
where ji=e(G++G2+) is ion current, λa is heat transfer ration of atom-ion gas, Q is energy obtained by the heavy particles because of elastic collisions with electrons. Let us assume that electrode is cooled with water. Therefore, the boundary conditions for this equation are the following: -λadTa/dx=-χ(Ta -Tw), x=b , where χ is complete heat transfer factor, Tw is the temperature of water cooling electrode.

Simulation of the discharge
Numerical algorithm of simulating of CCRF-discharge at atmospheric pressure is the same as the ones used for simulation of the same problem at low pressure which is described in the paper [18,19]. This algorithm is based on the method of transferring of non-linearity in equation coefficients to the lower time layer. Non-linear square terms in the right hand of the equations are linearized due to Newton scheme. Solution of equations for gas temperature was realized once per T period using Jacobi's iteration process.
Calculating of charged particle flow is needed for determine of conductivity current. Density of averaged ion flow is also included into the first part of equation for atomic and ion temperature. The problem is that numerical differentiation results significant inaccuracy in calculation of flow density due to high gradients of solution and coefficients of balance equations. So, we used Scharfetter-Gummel's approach to calculate flow density [20]. We determine charged particles densities by using implicit difference scheme, then we determine charged particle flow per the algorithm. This calculations sequence allows us rejecting a restriction on a scheme time step, which is associated with Courant condition.

Calculation data
Results of simulation of CCRF discharge at atmospheric pressure, interelectrode distance 2 mm, and frequency f=13.56 MHz, came near results of researches [21,22] and gave good fit with them. We perform simulation of the problem at atmospheric pressure, interelectrode distance 20 mm, and frequency f=13.56 MHz (Fig.1)   Fig. 1. Distribution of averaged density of charged particles at interelectrode distance of 20 mm.
It is found that maximum of gas temperature about 2000 К is achieved in the middle of the discharge, and the minimum about 930 К is achieved near the electrodes. Difference between maximum and minimum of gas temperature becomes significantly affecting on charged particles distribution. The density of charged particle distributions in case of interelectrode distance d=20 mm becomes more observable then in case of d=2 mm.
The most significant changes occur in distribution of Ar2 + . The density in the middle of the charge is almost constant in case of d=2 mm, but local maximum emerges close to the near-electrode layer when interelectrode distance is in case of d=20 mm, and local minimum is reached in the middle of charge area.
Ratio of minimum to maximum density of molecular positive ion of argon value is reached about 0.43, in other words, density of molecular ions in the middle of discharge gap drops twice. However, quasi-neutrality of plasma in the ambipolar area survives.
Ratio of dimmers to excited atoms also changes in case of d=2 mm, the maximum is equal approximately 0.78, but in case of d=20 mm, maximum value is about 0.38. Note that dimmer density can be considered permanent in case of d=2 mm, but in case of interelectrode distance is 20 mm, the drop of density occurs at the same point, where gas temperature is reached maximum.