Influence of the Dufour effect on striations formation in radio-frequency discharges

In recent years, interest in striation phenomena in radio-frequency (rf) discharges has risen due to the availability of new experimental data and implementation of new computational models. Depending on the conditions, different mechanisms of discharge striations are realized. These are the ionization instability, the instability due to the electron attachment to electronegative gases or the instability due to thermoelectric transport. Although the first two mechanisms were modeled quite extensively in recent years, the understanding of the influence of Dufour effect originating from plasma density gradients on stability of radio-frequency discharges in long tubes remains poor. In this paper, the influence of this mechanism on the longitudinal striations of radiofrequency discharge is presented using a one-dimensional model of argon discharge driven with rf excitation under intermediate pressure conditions of 0.5 Torr. It is found that striation formation is sensitive to the value of the thermoelectric heat transport coefficient in the low electron temperature range. The critical value of this coefficient necessary for the instability onset is derived using the linear stability analysis.


I. Introduction
Striations in glow discharges is one of the most fascinating phenomena in discharge physics. 1,2They have been observed in rare and molecular gases, at high and low gas pressures, in magnetized and unmagnetized plasmas of direct current (dc) and radio frequency (rf) discharges (see Refs. 3,4,5,6,7 and references therein).Due to such disparity of physical scales, different physical mechanisms govern the formation of striations. 2 These can be kinetic effects, non-linear dependence of the ionization frequency on the electron density, electron attachment to electronegative gases, thermoelectric heat transfer etc.
Striations of positive column of glow discharges was studied in numerous publications using fluid, kinetic and hybrid approaches. 8,9,10,11,12Less attention was devoted to the computational analysis of rf discharge striations.The striations of electronegative gas CF4 was studied in Ref. 3. It was found that the striations are generated due to the resonance between the rf frequency and the eigenfrequency of the ion-ion plasma.This frequency establishes a modulation of the plasma parameters and the electron transport.
Pattern formation in low-pressure inductively coupled plasma was studied in Ref. 6. Linear stability analysis has shown that striations appear at low and intermediate pressures due to electron energy transport.This mechanism was used in Ref. 13 to explain the formation of striations in rf discharges of argon and nitrogen forming perpendicular to the applied electric field.The authors used a two-dimensional fluid model coupled with the Poisson's equation.In this model, the discharge instability was enforced by the random perturbation of the electron density by 1%.The longitudinal striations of rf discharge due to the ionization instability was modeled in Ref. 7 using the one-dimensional fluid drift-diffusion model.The nonlinear dependence of the ionization frequency on the electron density was used in these studies.Standing striations were observed for a wide range of rf frequencies and gas pressures.The electron kinetics for similar conditions was studied in Ref. 12.
Thermocurrent instability of direct current discharge in various gases was studied in Refs.14,15,16.It was shown that this instability arises due to the competition between the diffusion and thermodiffusion electron fluxes.The electron density and temperature gradients are directed in opposite directions.Then, if the thermodiffusion flux exceeds the diffusion flux, any plasma perturbation will grow with time resulting in thermocurrent instability.
To our knowledge, the self-consistent modeling of striations formation in rf discharges in long tubes was attempted only in Refs.7 and 12.In Ref. 7, the mechanism of discharge striations was prescribed by introducing the special form of the ionization frequency, while the mechanism of discharge striations in Ref. 12 was not clear.In the present paper, the influence of Dufour effect originating from plasma density gradients within the plasma on striations of rf discharge at the gas pressure of 0.5 Torr is studied by the drift-diffusion fluid model.At such a pressure, the fluid approach to the description of plasma is justified due to the frequent electron-neutral collisions.
Also, it is important to note that for the conditions of our studies, the angular rf frequency is much higher than the electron energy relaxation frequency (see discussion in Section III.A).Therefore, the plasma density remains constant while the electric field oscillates with the rf frequency. 17e paper is organized as follows.Section II discusses the fluid model used in the present studies.Section III presents the analysis of the influence of the thermoelectric heat transfer on the discharge stratification.Section IV summarizes the results of the present studies.

II. Numerical model
In the present studies, the one-dimensional fluid model for two-component plasma consisting of electrons and argon ions (Ar + ) is used.The dominant background neutral (Ar) "gas" density was kept constant and homogeneous throughout the entire domain.Also, the generation of the electronically excited argon species and consequent stepwise ionization were neglected in the present model.These processes may be important for the conditions of our studies, but do not affect the main objective of our study in elucidating the role of the thermoelectric heat transfer coefficient on striation formation. 17e model solves fluid equations for the species continuity with the drift-diffusion approximation: .
Here,  , and Γ , are the number density and number flux of electrons and ions, respectively,   is the ionization rate coefficient tabulated as the function of the electron temperature   ,   is the electron number density and   is the background gas density.The last term on the right-hand side of Eq. ( 1) simulates the charged species ambipolar diffusion in the radial direction.This diffusion time is defined as , where Λ =   /2.4 is an effective transverse discharge dimension (  is the discharge tube radius) and   is the ambipolar diffusion coefficient.The density of the background gas argon,   , is assumed to be homogeneous and constant.
The species flux terms for both charged and neutral species are calculated using the driftdiffusion approximation: Here,  , and  , are the species mobility and diffusion coefficients,   = 1 and   = -1.Electron transport properties are computed as the functions of electron temperature using the Boltzmann equation solver Bolsig+, 18 which solves the Boltzmann equation in the two-term approximation.
The Ar + transport properties were obtained from Ref. [19].The electrostatic potential is obtained by solving the Poisson's equation: where   is the elementary charge.The electron energy density   is defined as where   is in the units of eV.It is calculated by solving the electron energy conservation equation: 13 Here, Γ  is the flux of electron energy and   is the source term.They are defined by 18 In Eq. (7),  is the electric field obtained as  = −/,   and   are the thresholds for Ar ionization and excitation to the first electronic level (in eV),   is the rate coefficient of the excitation reaction of the first electronic level of Ar,   is the electron-argon momentum transfer collision frequency,   is the electron mass,   is the argon mass, and   is the background gas temperature.Also,   and   are, respectively, the electron energy mobility and diffusion coefficient calculated by Bolsig+.
The first term on the right-hand side of Eq. ( 7) describes the electron heating by the electric field (further ohmic heating), the second term describes the contribution of inelastic collisions, and the last term describes the contribution of elastic collisions.For further discussion, it is also convenient to introduce the thermoelectric energy transport coefficient defined as Here,   and   are, respectively, the coefficients of heat diffusion and coefficient of thermoelectric transport for electrons.According to Ref. [18], they are defined by solving the Boltzmann equation as   = ̅   and   = ̅   , where ̅ is the average electron energy.
Here, it is important to note that Eq. ( 5) in which   and   are obtained from the EEDF is consistent with the two-term approximation of the Boltzmann equation.Another form was used in Ref. [13] which explicitly includes the thermoelectric flux defined as the coefficient (8)   multiplied by the electron density gradient.This form is inconsistent with the two-term approximation. 18netic Maxwellian flux condition combined with secondary electron emission flux is used to specify the boundary flux of the electrons to electrode surfaces: Here,  ̂ is the unit normal vector pointing towards the electrode surface from the plasma.In the present studies the secondary electron emission from the walls was neglected.
Mobility limited flux condition is imposed at the walls for ions using For the electron energy equation, the energy flux incident at the walls is given by Here, Γ   is the electron number flux at the walls given by Eq. (10).
Equations ( 1) and (5) were discretized using the Scharfetter-Gummel spatial discretization scheme. 20,21A uniform mesh with 2000 grid cells is used for the simulations and the time step of 10 -10 s was used for all cases considered.
In the present studies, the interelectrode distance was set equal to 30 cm.The background gas pressure was 0.5 Torr and the background gas temperature was 300 K.The right boundary was kept grounded while the rf potential was applied to the left boundary,   () =  0 sin (2) with  0 = -400 V and  = 13.56MHz.The tube radius was assumed equal to   = 1.1 cm.The secondary electron emission from the boundaries was neglected in the present studies.These coefficients were obtained by using the solver Bolsig+. 18Three cases were considered.
First, it was assumed that the electron energy distribution function (EEDF) is Maxwellian.In this case, one obtains   =   , 18 and, consequently,   = 0.In case 0, the nonequilibrium conditions were considered without accounting for the electron-electron Coulomb collisions.In case 1, to understand the influence of the electron-electron collisions on the Dufour effect the plasma ionization degree of 10 -4 was set in the solver Bolsig+.Here, it is important to note the rate coefficients of the electron-neutral reactions considered in our studies also depended on the EEDF's choice.

III. Results and discussion
A.

Maxwellian vs. non-Maxwellian electron energy distribution function
In this subsection, we compare case 0, i.e. the case of non-Maxwellian EEDF without accounting for the electron-electron collisions, with the case of the Maxwellian EEDF.In the steady state, the electron energy flux defined by Eq. ( 6) may be simplified. 6,9Let us consider the quasi-neutral plasma of positive column.Then, the transport of both electrons and ions is described by the ambipolar flux defined as while the ambipolar electric field within the plasma is Then, the electron energy flux becomes: Here, the thermoelectric energy transport coefficient is defined by Eq. ( 8) and the thermal conductivity is defined as Equation (14) shows that in the case of non-equilibrium EEDF the flux consists of two competing components proportional to the electron density and temperature gradients.The term with   leads to the flux similar to the Dufour effect discussed, for instance, in Ref. 22.For the Maxwellian EEDF one obtains   = 0, while the thermal conductivity is defined as   = − i.e. it is defined by the electron temperature gradient only.Namely, it is directed opposite to  The spatial distribution of the period-averaged plasma parameters for the Maxwellian and non-Maxwellian EEDFs is shown in Figure 2 and Figure 3, respectively.One can see that for   = 0 the discharge structure resembles the structure of dc glow discharge in long tubes. 23mely, there are the sheaths of non-neutral plasma in the vicinity of the electrodes, and the long positive column of quasi-neutral plasma.For the observation of such a structure, the last term on the right-hand side of Eq. ( 1) describing the plasma ambipolar radial diffusion to the side wall is crucial.Without accounting for this term, the plasma density reaches unrealistically large values, ⁓10 20 -10 21 m -3 due to the slow plasma diffusion from the discharge center to the left and right electrodes.
Unlike the dc discharge, rf discharge does not have the anode layer. 23Here, both electrodes act as the anode each half of the rf period.As in the case of dc discharge, the sheaths are crucial for the discharge maintenance since the external electric field is shielded by quasi-neutral plasma.
Therefore, the main potential drop is obtained near the electrodes [see Figure 3   The switching to the non-equilibrium EEDF changes the discharge structure drastically (see Figure 2).One can see the stratification of the positive column into 14 standing peaks in the electron density.Unlike in Ref. 7 where the generation of striations was enforced by the nonlinear dependence of the ionization rate coefficient on the electron density, in the present studies they are initiated self-consistently.As will be discussed in more detail in Subsection III.B, they originate at the sheath edges and facilitate the electron heat transfer from the sheaths towards the discharge center.Here, it is also important to note that in the simulations presented in Ref. 13, the transverse striation formation was initiated by randomly perturbing the electron density, i.e. by introducing the numerical noise.In our simulations, the striations are formed self-consistently without the explicit introduction of any noise although there still can be some noise associated with the numerical discretization of plasma equations.
Figure 3(b) shows the period-averaged electron ohmic heating.One can see that for the stratified discharge, the dominant heating is still obtained within the plasma sheaths where the electric field is the highest.However, now there is also significant electron heating within the plasma bulk where the electric field ~400 V/m is obtained [Figure 3  Figure 4 shows the electron density, electric field, electron temperature and power deposition from the electric field to the electrons at two different times of the rf period.One can see that the electron density does not change with time, while the electron temperature exhibits modulations around its average value.  ≈ 1.4×10 6 s -1 .Thus, the collision frequency is   ≈ 1.3×10 7 s -1 , which is much larger than   .
One can conclude that for the conditions of the present studies, discharge operates in the so-called high-frequency regime 24 since  >   in which the plasma density is defined by the "effective" electric field and do not depend on its instantaneous values.This is seen in Figure 4(a).
The discharge striations change the mechanism of electron mass and heat transfer between the electrodes.Figure 5(a) shows the comparison between the period-averaged drift and diffusion components of the electron flux in Eq. ( 2) and Figure 5(b) shows the difference between their absolute values.One can see that unlike in the case of   = 0 where the electrons are carried by the small electric field presented within the plasma, in the case of non-Maxwellian EEDF the main mechanism of electron transport is their diffusion due to the density gradient within one striation.
Note that the diffusion component does not oscillate since the electron density is independent of the oscillating electric field (see Figure 4).The drift component exhibits oscillations due to the modulations of the electric field like those shown in Figure 4(d) for the ohmic heating term.

B. Influence of the electron-electron collisions on discharge structure
In order to clarify the nature of the obtained striations, this subsection analyzes the influence of the electron energy transport coefficient on the rf discharge structure.For this, two cases discussed in Section II were considered.In case 0 (  =  ,0 ), the influence of the electronelectron collisions on the electron energy distribution function was ignored.In case 1 (  =  ,1 ), these collisions were considered which resulted in a small value of the thermoelectric coefficient for   < 3 eV [see Figure 1(c)].Figure 6 shows the comparison between the electron density and the electron temperature obtained in these two cases.One can see that there is no discharge stratification for   =  ,1 .
Equation (14) shows that the electron energy flux is defined by two competing components proportional to the electron density and temperature gradients.The comparison between these two components for two cases considered in this subsection is shown in Figure 7.The sum of these two terms is balanced by electron ohmic heating and the electron energy dissipation in inelastic collisions.The momentum transfer term is much smaller than other terms due to the small ratio between the electron and neutral masses and can be neglected in our analysis.does not change with time in the plasma bulk because the electron density is constant here.
At the same time, the thermal flux follows the oscillations of the electron temperature gradient.
Therefore, there are the times during one rf period when both fluxes add up which facilitates the energy transfer.This distance agrees with the mean free path of electrons being accelerated within the sheaths and entering the quasi-neutral plasma bulk.For energetic electrons obtained within the sheaths, this is also the electron energy relaxation length.For these electrons, the dominant collision is the ionization one which is ~2×10 -20 m 2 .Then, for the gas pressure of 0.5 Torr the mean free path is estimated as  ≈ 0.3 cm.
After propagating the distance of 0.25 cm, hot electrons enter the plasma with a weak electric field within it and start dissipating their energy.This leads to the local decrease of the electron temperature [see Figure 6(b)] and the electron energy density [see Figure 7(a)].Since the electron density at the sheath edge is still growing [see Figure 7(a)], the thermoelectric flux does not change its direction.Also, at the sheath edge the gradient of the electron density starts growing faster than within the sheath.Figure 1(c) shows that for   =  ,0 the thermoelectric coefficient remains small in the temperature range 1 -2 eV.Thus, the thermoelectric flux obtains the local maxima at the sheath edge.This energy is carried by the thermoelectric flux from the sheath to the striations closest to the electrodes.Similar explanation is valid for the flux      which is directed opposite to the electron temperature gradient.

C. Linear analysis of the instability growth
For the linear analysis of the instability growth, let us consider the steady state quasi-neutral plasma of the positive column.Then, the ambipolar flux and electric field are defined, respectively, by Eqs. ( 12) and ( 13), and the electron energy flux is defined by Eq. ( 14).The solutions of plasma equations are looked in the form: Here,  0 and  0 are, respectively, the unperturbed values of the electron density and temperature, Ω =  +  with  is the angular frequency and  is the instability growth rate, and  is the wave number.
Let us assume that the EEDF is close to the Maxwellian.Then, the electron thermal .Also, let us assume the Einstein relation between the electron diffusion coefficient and mobility, i.e.   =     .Then, using the observation that  0   ≫  0   9 and neglecting the term  2   compared to   ′ , one obtains from the particle balance equation ( 1): Thus, the system of equations ( 17) and ( 18) defines   and   .The determinant of this system must be zero which results in the system of equations for the angular frequency and the instability growth rate.
The equation for the imaginary part gives  = 0 which is expected for the rf discharge where the standing striations are obtained (see Section III.A).The equation for the real part allows obtaining the following dispersion relation: Analogous relation was obtained in Ref. 9 for dc discharges.One can see that the instability growth rate does not depend on the electron density but is defined by the ionization and the electron heat and mass transfer.The maximum value of the wave number is: , this equation can be written as Then, the maximum value of the instability growth rate is which is analogous to the one derived in Ref. 9.
For plasma to be unstable, the growth rate should be positive.This condition defines the critical value of the thermoelectric transport coefficient necessary for the instability onset: The expression on the right-hand side is positive.Therefore, the first observation from Eq. ( 23) is that the thermoelectric transport coefficient should also be positive [see Figure 1 -2)×10 2 m 2 /eV-s (see Figure 1) and   ~ 1.6 m 2 /eV-s.Then, Eq. ( 24) requires that  0 > 660 eVm 2 /s. Figure 1(c) shows that this condition is satisfied for the case 0 and is not satisfied for the case 1.This explains the results discussed in Section III.B.
Based on the results presented in this section and Section III.B, we can conclude that the mechanism of discharge stratification at the given conditions is the following.If the thermoelectric transport coefficient exceeds the critical value defined by Eq. ( 23), then the thermoelectric flux starts dominating the flux of thermal conductivity.This leads to the energy flow towards the regions with the largest electron density and temperature from the regions with the smaller density and temperature.This leads to the increase of the electron temperature in these regions which, in turn, leads to the increase of the ionization frequency.The latter results in further increase of the electron temperature in these regions etc.The instability stabilizes when both fluxes are balanced by the electron heating by the electric field and the energy dissipation in collisions with neutrals.

D. Influence of the non-linear dependence of the ionization rate coefficient on the electron density
This subsection briefly discusses the influence of the ionization rate coefficient on the electron density.Such dependence is obtained due to the influence of electron-electron collisions on the EEDF 25 Here,   (  ) is the ionization rate coefficient obtained from Bolsig+,   controls the rate of the non-linear dependence and  1 defines the saturation value.In the present studies, we have used the values   = 8×10 15 m -3 and  1 = 10  as in Ref. 7. For the gas pressure of 0.5 Torr, the chosen value of   implies the plasma ionization degree exceeding 10 -6 .Therefore, here we compare the case 1 (see Section II) obtained under the assumption of non-Maxwellian EEDF with the accounting for the Coulomb collisions and the Maxwellian EEDF (for the transport parameters, see Figure 1).In both cases, the ionization rate coefficient was determined by Eq. ( 25) with   (  ) obtained for a given assumption on the EEDF.
Here, it is important to note that the striations formation driven by the Dufour effect is only obtained for small discharge currents when the electron-electron collisions are negligible.Hence, the thermoelectric transport coefficient is large, and the electron diffusion flux destabilizes the electron thermal diffusion flux in the electron energy balance equation.
Figure 8 shows the comparison between these two cases for the rf frequency of 13.56 MHz, background gas pressure of 0.5 Torr and three values of the rf voltage amplitude.Our simulation results have shown that for both cases, striations exist only in the narrow window of applied voltages (discharge currents).Namely, for the Maxwellian EEDF it is in the range  0 ≈ 1200 -2500 V, while for the non-Maxwellian EEDF it is in the range  0 ≈ 1400 -2200 V.In both cases, the discharge current is approximately in the range 100 -200 mA.This observation is in qualitative agreement with the experimental results reported in Ref.This is obtained because for this case the electron density within the positive column is ⁓10 16 m -3 which exceeds  1 in Eq. (25).Therefore,   saturates at   >  1 and the ionization rate coefficient does not depend on the electron density.For the Maxwellian EEDF, the plasma density is below  1 which results in the striations formation.

IV. Summary
The influence of Dufour effect on the striation of radio-frequency discharge was studied using the drift-diffusion fluid model for electrons and singly charged argon ions.
It is obtained that the striations formation is sensitive to the value of the thermoelectric heat transport coefficient in the range of electron temperatures 1 -4 eV typical for the positive column of radio-frequency discharge in long tubes.The critical value of the thermoelectric heat transfer coefficient was obtained using the linear stability analysis.The simulation results have shown that for the non-equilibrium electron energy distribution function the thermoelectric heat coefficient exceeds this critical value.Therefore, the standing striations originated at the sheaths edges and filled in the entire interelectrode space.
The accounting for the electron-electron collisions led to the thermalization of electrons which resulted in the decrease of the thermoelectric heat transfer coefficient.Consequently, its value dropped below the critical value for the instability onset.Therefore, the striations were not obtained.
Thus, the Dufour effect is responsible for the discharge stratification only for small discharge currents when the electron-electron collisions are negligible.For higher discharge currents, these collisions are non-negligible which results in the non-linear dependence of the

Figure 1 .
Figure 1.(a) Electron mobility, (b) electron diffusion coefficient and (c) thermoelectric energy transport coefficient obtained for three different assumptions on the electron energy distribution function.

Figure 1
Figure 1 shows three sets of the electron transport coefficients used in the present studies.

5 2𝑘
.Then, the energy flux for the Maxwellian EEDF is reduced to . leads to the energy flow from the regions of hotter plasma to the colder regions.In such plasmas, the Dufour effect cannot be obtained.In the case of non-Maxwellian EEDF, flux      may have an opposite direction and, as will be discussed in Section III.B, destabilize discharge.16

Figure 2 .
Figure 2. (a) Period-averaged electron density and (b) period-averaged electron temperature obtained with and without accounting for the thermoelectric effect.
(a)].Hence, the highest electron ohmic heating is also obtained within the sheaths [see Figure 3(b)], where the oscillating electric fields heat up the plasma electrons reaching the electrodes during the sheaths collapse.This results in the highest electron temperature [see Figure 2(b)] within the sheaths.

Figure 2 (
Figure 2(a) shows that the peaks of the plasma density are obtained near the electrodes.

Figure 3 .
Figure 3. (a) Period-averaged electric field and (b) period-averaged ohmic heating of electrons obtained with and without accounting for thermoelectric effect.
Figure 3(b) that the total energy gain by the electrons within one striation is comparable to that

Figure 4 .
Figure 4. Spatial distributions of (a) the electron density, (b) electric field, (c) electron temperature and (d) power deposition from the electric field obtained near the center of the simulation domain at two different times of rf period.Dash lines show the period-averaged values.

Figure 4 (
b) shows the oscillations of electric field during one rf period which results in the oscillations of the electron ohmic heating term [Figure4(d)].Such dependence of the plasma parameters can be understood by comparing the typical frequencies with the rf frequency.For the conditions of our studies, the frequency of electron loss to the radial wall is estimated as   =     Λ 2 ~ 0.3×10 6 s -1 , where the ion mobility is   ≈ 1.6 m 2 /Vs.The collision frequency is estimated as   =   +   +   , where for argon gas  = √2  /  ~ 5.2×10 -3 ,   is the ionization frequency, and   is the electronic excitation frequency.For   = 4 eV [see Figure2(b)], one has   ≈ 2.4×10 9 s -1 ,   ≈ 3.2×10 4 s -1 and

Figure 5 .
Figure 5. (a) Comparison between the period-averaged drift and diffusion components of the electron flux and (b) the difference between the absolute values of the drift and diffusion components of the electron flux for case 0 in the vicinity of one striation.

Figure 6 .
Figure 6.The comparison between (a) the electron density and (b) the electron temperature obtained for two values of the thermoelectric coefficient.

Figure 7 (
Figure 7(a) shows that for   =  ,0 both fluxes are directed in the opposite direction.The

Figure 7 (
Figure 7(b) shows that for   =  ,1 both fluxes are non-zero only near the electrode, at x

Here 1 𝜏.
)   is the derivative of the ionization frequency at   =  0 .We also took into account the fact that in the steady state the local electron density in the positive column is defined by the balance between the ionization and wall losses, i.e.     ( 0 ) = By neglecting the elastic collision term in the electron energy balance equation, we find: 5 2    0  0 +    0   ′ )   = 0. (18) (c)].In many cases, the Arrhenius form of the ionization frequency   =  exp (−     ) approximates well the ionization frequency obtained by solving the Boltzmann equation.Then, one can obtain that .(23) can be re-written as  0 > 2   0 (  +  0 ( of the present studies,  0 ~ 3 -4 eV [see Figure 2(b) and Figure 6(b)],   ~ (1 7.

Figure 8 (
Figure 8(b) shows the comparison between the period-averaged ion densities obtained for

Figure 8 .
Figure 8. Period-averaged ion densities obtained for the rf frequency of 13.56 MHz for three values of the voltage amplitude: (a) 1000 V, (b) 1600 V and (c) 2400 V.