Stochastic and Nonlinear Dynamics in Low-Temperature Plasmas

Abstract

Low-temperature (LT) plasmas have a substantial role in diverse scientific areas and modern technologies. Their stochastic and nonlinear dynamics strongly determine the efficiency and effectiveness of LT plasma-based procedures involved in applications such as etching, spectrochemical analysis, deposition of thin films on substrates, and others. Understanding and controlling complex behaviors in LT plasmas have become a serious research problem. Modeling their behavior is also a major problem. However, models based on hydrodynamic equations have proven to be useful in their study. In this chapter, we expose the use of fluid models taking into account relevant kinetic processes to describe out from equilibrium LT plasma behavior. Selected topics on the stability, stochastic, and nonlinear dynamics of LT plasmas are discussed. These include the coexistence of diffusive and wave-like particle transport and delayed feedback control of oscillatory regime with relaxation.

1. Introduction

A plasma is a complex system consisting of mutually interacting particles. Some or all of the following can coexist in a plasma: electrons, photons, positive ions, negative ions, metastables, free radicals, and neutral species. In the laboratory, plasmas can be classified into two categories, namely, high-temperature and low-temperature plasmas. The former is obtained in fusion devices, and the plasma is generally in local thermal equilibrium (LTE), which means that each of the existing species is in thermal equilibrium with each other. Typical LTE plasma temperatures range from 5000 K to 20000 K. Low-temperature plasmas, produced in well-known electrode configurations, are not in thermal equilibrium (non-LTE) as each of the species has its own temperature. In these systems, electrons are characterized by their higher temperatures than those of other heavier species (ions, atoms, etc.). The thermodynamic state (thermal plasma or nonthermal plasma) of the plasma mainly depends on its pressure $$p$$. This is understood if one considers that in a high-pressure plasma, many collisions occur leading to a more frequent interchange of energy between species than in nonthermal plasmas. The pressure determines not only the thermodynamic state of the plasmas but also the distance $$d$$ between electrodes. In fact, the product $$pd$$ can be used to distinguish the thermal and nonthermal states. Besides pressure and electrode separation, the working gas, the applied electric (or electromagnetic) field, and the configuration of the containing device also determine the thermodynamic state of the plasma. Among the wide variety of gas discharges, in this chapter, we will discuss nonthermal plasmas. Among them, we find the following: direct current (DC) glow discharges, capacitively coupled radio-frequency (RF) discharges, pulsed glow discharges, atmospheric pressure glow discharges, corona discharges, dielectric barrier discharges, and magnetron discharges. For an overview of the physical description and applications of both thermal and nonthermal plasmas, see Ref. [1]. A DC glow discharge is a nonthermal plasma wherein the main role is played by electrons. Normally, they have high temperatures which lead to inelastic collisions determining the ionization rate in the system. In fact, these collision processes sustain the plasma. The temperature of heavier particles is lower than that of electrons. Their kinetics has importance in applications such as etching, spectrochemical analysis, and deposition of thin films on substrates through processes such as sputtering, ion implantation, deposition, etc. When one or both of the electrodes are insulating materials, positive or negative charges accumulate and get charged causing plasma annihilation. This process is avoided by introducing an alternating voltage to the electrodes giving rise to the so-called capacitively coupled discharges which are obtained in radio frequencies. This kind of discharge is especially used for the spectrochemical analysis of nonconductive materials or deposition of thin films on dielectric electrodes [2]. If the applied alternating voltage is in the form of pulses, a very similar discharge to DC glow discharge is obtained in a short period of time. In this way, higher peak voltages and currents can be applied, and more efficient sputtering, ionization, and excitation processes are obtained. Because of the existence of lapses of time where the discharge burns out, the pulsed discharges do not excessively heat up in such a way that the temperature of the gas is lower than the electron temperature. This is conducive to the nonthermal state of the plasma. Atmospheric pressure glow discharges occur at higher (e.g., atmospheric) pressures. Increasing the operating pressure makes it necessary to reduce the distance between electrodes in such a way that discharge conditions are maintained. In this kind of discharge, one (or two) electrode is made of a dielectric material, and the device must be operated with an alternating voltage. However, under these operating conditions, a stable glow discharge can be obtained when some working gases (e.g., helium) are used. In other cases (nitrogen, oxygen, argon), it takes the form of a filamentary glow discharge, however, in a stable state. Atmospheric pressure plasmas are especially suitable for the deposition of both organic and inorganic coatings [3] and the treatment of industrial and biological surfaces [4]. We end the description of nonthermal plasma discharges by considering the dielectric barrier discharge. The device used to produce this kind of discharge is quite similar to that of the atmospheric pressure glow discharge. The main difference is that the latter is homogeneous along the electrodes, while the dielectric barrier discharge produces filamentary discharge of a very short duration. Note, however, as mentioned earlier, that the atmospheric pressure discharge may also consist of those filamentary discharges for certain gases. In each of the filamentary discharges, the high electron energy can excite the gas-provoking chemical reactions and emission of radiation which are of interest in various applications [5].

In this brief overview, it may be appreciated that nonthermal plasmas have a substantial role in diverse scientific areas and modern technologies. This has attracted the interest of both theorists and experimentalists giving rise to the discovery of very complex behaviors which were unknown until 20 years ago. It has also attracted the attention to nonlinear dynamics of nonthermal plasmas, since it strongly determines the efficiency and effectiveness of nonthermal plasma-based procedures. Because of this, understanding and controlling complex behaviors in nonthermal plasmas have become a serious research problem. We find a great variety of complex dynamics in nonthermal plasmas, namely, deterministic chaos, mixed-mode oscillations, and homoclinic chaos as well as order–chaos transitions, stochastic resonance, coherence resonance, etc. In this chapter, we will describe some of them and the control techniques developed to mitigate or even enhance those kinds of complex behaviors in low-temperature discharges. In Section 2, the basis of fluid modeling of nonthermal plasmas is exposed. The dynamic equations of particle densities and electron energy as well as constitutive expressions for the fluxes can be found in the same section. In Section 3, a brief discussion on stability aspects of nonthermal plasmas is given. In Section 4, the plasma dynamics is analyzed as a stochastic process. The stochastic equations are solved by means of a high-order numerical method. We also discuss the coexistence of diffusive and wave-like particle transport regimes in the stationary state. Finally, Section 5 is devoted to the nonlinear behavior of nonthermal plasmas, and several control methods are described. Then the delayed feedback control (DFC) is applied to the oscillation with relaxation regime under specific values of physical parameters of the plasma, resulting in an efficient control of the oscillations.

2. Fluid modeling of nonthermal plasmas

The importance of nonthermal gas discharges in various industrial and technological applications makes understanding and controlling such plasmas necessary. This requires to model a dissipative system out from thermodynamic equilibrium where several nonlinear kinetic processes take place giving rise to complex behaviors. Among others, nonlinear kinetic processes include electron attachment and ionization, ion–ion recombination, etc. Models describing gas discharges can be classified as fluid methods [6], particle methods, and hybrid methods [7] which are combinations of the former and the latter. Hybrid methods are particularly useful when there exist several time scales in the system. In this kind of model, the fast and energetic electrons are treated as particles while the slow components (heavier species) as coexisting fluids. Fluid models (FMs) are used for describing collective movement of particles, such as wave effects, which cannot be obtained from a kinetic equation of the electron behavior. In fluid models, all the species are considered as fluids, and they make use of continuity equations and the energy balance equation (normally written only for electrons). Some additional equations are necessary in order to complete these sets of equations, namely, the Poisson equation for the self-consistent (self-generated) electric field and constitutive equations relating the dissipative fluxes with the thermodynamic forces. These equations resort to well-known phenomenological laws such as Ohm’s law, Fourier law, etc. [8, 9]. To complete the framework, expressions for the transport and rate coefficients are necessary. An often-used assumption for these coefficients is the local field approximation. In this approximation, the coefficients are supposed to depend on the reduced local electric field $$E/p$$ with $$p$$ the pressure of the working gas . This assumption has a limited validity since it implies that ionization can occur in the cathode sheath solely. Nevertheless, it is known that bright luminous regions of high charge density in the negative glow can be produced by nonlocal ionization processes [6, 10, 11]. Instead, sometimes the local mean energy approximation is used. In this approximation, expressions for the transport and rate coefficients come from the Boltzmann equation as functions of the electron temperature which is obtained by solving the energy balance equation of electrons.

The continuity equations can be written as

$$\frac{\partial n_j}{\partial t}+\nabla \cdot \overrightarrow{J_j}=S_j$$E1

where $$n_j$$ denotes the densities of the species in the plasma, with $$j=\text{p,}\text{e,}\text{m,}\text{n,}\text{g},\mathrm{...}$$ for positive ions, electrons, metastable molecules, negative ions, working gas, etc., respectively. $$S_j$$ is the volume gain or loss of particles. The particle fluxes are given by the constitutive equations which, if only drift and diffusion are considered, become [12]

$$\begin{array}{l}{\overrightarrow{J}}_{\text{p}}={\mu }_{\text{p}}{n}_{\text{p}}\overrightarrow{E}-D_{\text{p}}\nabla n_{\text{p}}\\ {\overrightarrow{J}}_{\text{n}}=-{\mu }_{\text{n}}{n}_{\text{n}}\overrightarrow{E}-D_{\text{n}}\nabla n_{\text{n}}\\ {\overrightarrow{J}}_{\text{e}}=-{\mu }_{\text{e}}{n}_{\text{e}}\overrightarrow{E}-\nabla \left(D_{\text{e}}{n}_{\text{e}}\right)\\ {\overrightarrow{J}}_{\text{m}}=-D_{\text{m}}\nabla n_{\text{m}}\end{array}$$E2

where $${\mu }_{\text{p}},$$ $${\mu }_{\text{n}},$$ $${\mu }_{\text{e}}$$ are the mobilities of ions and electrons; $$D_{\text{p}},$$ $$D_{\text{n}},$$ $$D_{\text{e}},$$ $$D_{\text{m}}$$ the corresponding diffusion coefficients (see Section 3); and $$\overrightarrow{E}$$ the electric field. Note that $$D_{\text{e}}$$ is included in the gradient operator, since it depends on the mean electronic energy $$U_{\text{e}}$$ which is a function of the spatial position. The same must be said for the electronic mobility $${\mu }_{\text{e}}.$$ The balance equation for the mean electronic energy is given as

$$\frac{\partial (n_{\text{e}}{U}_{\text{e}})}{\partial t}+\nabla \cdot {\overrightarrow{J}}_{U_{\text{e}}}=-{\overrightarrow{J}}_{\text{e}}\cdot \overrightarrow{E}-n_{\text{e}}{S}_{U_{\text{e}}}$$E3

The work done by the electric field on the electrons and the electronic energy changes due to elastic and inelastic collisions (the term with $$S_{U_{\text{e}}}$$). In the same equation, the energy flux is given by

$${\overrightarrow{J}}_{U_{\text{e}}}=-{\mu }_{\text{e}}{}^{*}{n}_{\text{e}}{U}_{\text{e}}\overrightarrow{E}-\nabla \left(D_{\text{e}}{}^{\ast }{n}_{\text{e}}{U}_{\text{e}}\right)$$E4

with $${\mu }_{\text{e}}{}^{*}$$ and $$D_{\text{e}}{}^{*}$$ being functions of the mean electronic energy.

The source terms in the continuity equations depend on several kinetic processes. Electrons are mainly generated by ionization (collisions of electrons with the working gas molecules). Negative ions are generated by the attachment of electrons by the working gas molecules and lost by ion–ion recombination. Metastable molecules are generated by electron collisions with the working gas molecules and lost by superelastic collisions of electrons with metastables. In the last process, associative detachment processes of metastables with the working gas molecules are also of importance. The total electric flux (electric current density) depends on particle fluxes as shown in the following:

$$\overrightarrow{j}=q\left({\overrightarrow{J}}_{\text{p}}-{\overrightarrow{J}}_{\text{n}}-{\overrightarrow{J}}_{\text{e}}\right)$$E5

which obeys the condition

$$\nabla \cdot \overrightarrow{j}=0$$E6

The total electric flux can be used to eliminate $$\nabla \cdot {\overrightarrow{J}}_{\text{e}}$$ in the continuity equation for the electrons. In this way, the electronic transport coefficients, which depend on the energy $$U_{\text{e}}$$, are no longer necessary.

Respecting boundary conditions, particle (electrons, ions) fluxes must be specified at the electrodes of the device, often considered as perfect absorbers of particles [11]. The energy flux (due to electrons) must also be specified at electrodes. In the case of conducting electrodes, the boundary condition for Poisson equation should fix the level of the electric potential (including the imposed potential) giving rise to Dirichlet-type and Neumann-type conditions when they are dielectrics (perpendicular electric field imposed on them).

Much is known on low-temperature plasmas by fluid model-based investigations. Fluid models (FM) have been used to study the fundamental properties and dynamic behavior of several gases: argon, oxygen, hydrogen, and helium in a diversity of configurations (see Section 1) and its results confirmed by experiments. The following is not an exhaustive list. In RF discharges, the FM of the plasma has reported the experimentally observed dependence of the harmonic part of the discharge current on the geometry of the device [11]. Coupled with external circuit equations, FMs have been used to reproduce the commonly observed mixed mode oscillations and chaos in glow discharges [13, 14]. The influence of impurities constituted by small amounts of nitrogen on the dielectric controlled discharge at atmospheric pressure in helium has been investigated by means of FM in one-dimensional models [15]. The temporal nonlinear behavior in atmospheric glow and dielectric barrier discharges such as period multiplication and secondary bifurcations has been studied by using a FM approach [16, 17, 18]. Another interesting problem analyzed with FM is the stability of homogeneous state in the positive column of a glow discharge [12]. The results agree with the observation that a hysteresis occurs at the transition from the H-mode to the T-mode.

3. Fluid modeling and instabilities in low-temperature plasmas

3.1. Ionization instability in the positive column of a glow discharge

The ionization instability occurs in the positive column of a glow discharge. It is a consequence of the nonlinear dependence of the ionization rate on the electron temperature. The main dynamic processes behind the instability are the fluctuations of the ionization rate provoked by fluctuations of the electron temperature. Electrons diffuse and are drifted in the electric field being heated and cooled by collisions with other species in the plasma. In this way, the electron temperature becomes a fluctuating quantity and thus the ionization rate. The self-consistent electric field also contributes to this process. The ionization instability can be described through a properly selected set of dynamic variables and the associated fluid equations [19]. The relevant properties are the electron and positive-ion densities, the ionization rate, the electric field, and the electron temperature. Fluid modeling description is based on i) the equation describing the diffusion, production, and recombination of charged particles Eq. (1); ii) an expression of the ionization rate as a function of the electron temperature and the ionization potential; iii) an equation for the self-consistent electric field in terms of the inhomogeneities of the electron density, namely, Poisson equation; and iv) the energy balance equation for electrons, Eq. (3) [19], describing their heating and cooling by the electric field and collisions respectively. Under specific conditions, the solutions of this set of equations show that a local fluctuation of $$n_{\text{e}}$$ results in a fluctuation of the temperature $$T_{\text{e}}$$ on the anode side, and then an ionization perturbation propagates toward the anode with a certain group velocity.

3.2. Transition from a uniform nonthermal RF discharge to an arc

The power density in a uniform RF discharge is limited by instabilities which destroy the homogeneity of the plasma by arcing. Several kinds of instabilities have been identified, the thermal instability and the $$\alpha$$ to $$\beta$$ state transition, being the two most encountered. The thermal instability occurs in a sequence of processes [20] leading to a (positive) fluctuation of the electron density $$\delta n_{\text{e}}$$ which increases the electric current density $$\delta j_{\text{e}}$$. It results in an increment of power dissipation $$\delta (j_{\text{e}}E)$$. This results in an increase in gas temperature $$\delta T$$ and leads to a decrease in the gas density $$\delta n_{\text{g}}$$. In turn, the reduced electric field fluctuates positively by $$\delta (E/n_{\text{g}})$$ causing an increase in the electron temperature $$\delta T_{\text{e}}$$. Hence, an increment in ionization occurs and a further increasing in the electron density $$\delta n_{\text{e}}$$. It has been assumed that the initial fluctuation of the electron density occurs rapid enough to consider that the electric field remains constant. It must be mentioned that, in fact, the described amplification process can be initiated with a fluctuation of any of the involved physical properties. The discharge is unstable if the rate of heating of the working gas, $${\nu }_{\text{h}}$$, is higher than the cooling rate, $${\nu }_{\text{c}}$$, divided by a factor, $$\nu *$$, known as the logarithmic sensitivity, given by [20]:

$$\nu *=\frac{d\left(\ln \left({\nu }_{\text{i}}\left(E/n_{\text{g}}\right)-{\nu }_{\text{a}}\left(E/n_{\text{g}}\right)\right)\right)}{d\left(\ln \left(E/n_{\text{g}}\right)\right)}$$E7

where $${\nu }_{\text{i}}$$ is the ionization frequency and $${\nu }_{\text{a}}$$ the attachment frequency. The cooling rate $${\nu }_{\text{c}}$$ enters the energy balance equation of the working gas which is written as follows:

$$c_{\text{p}}\rho \frac{\partial T_{\text{g}}}{\partial t}={\overrightarrow{j}}_{\text{e}}\cdot \overrightarrow{E}-c_{\text{p}}\rho \left(T_{\text{g}}-T_0\right){\nu }_{\text{c}}$$E8

The two terms on the right-hand side describe Joule heating and cooling of the gas through conduction and convection. At this point, it is necessary to introduce the growth rate of thermal fluctuations, $$\mathrm{\Omega }$$, in such a way that $$\delta T_{\text{g}}\approx e^{\mathrm{\Omega }t}$$. Using the two previous equations, it can be proved that the growth rate of thermal fluctuations depends on the heating and cooling frequencies as follows:

$$\mathrm{\Omega }={\nu }_{\text{h}}(\nu *+1)-{\nu }_{\text{c}}$$E9

Clearly, if $$\text{Ω}<0$$, the discharge is stable, and if $$\mathrm{\Omega }>0$$, thermal instability occurs, and the discharge is unstable. Calculated values for $$\nu *$$ are in the range of 5–10 [21]. It has been shown that if there is an increment of 10%–20% in the gas temperature, thermal instability occurs.

The $$\alpha$$ to $$\beta$$ state transition is another known instability of the RF gas discharge. This instability occurs between two modes of sustaining the discharge. In the $$\alpha$$ mode, the discharge is sustained through volumetric ionization, while in the $$\beta$$ mode, sustaining takes place through secondary electron emission. The transition appears when the sheath suffers a breakdown after the electric field reaches a critical value. The instability associated with the $$\alpha$$ to $$\beta$$ state transition takes place in a smaller time scale than that of the thermal instability described above.

3.3. Linear stability of a DC-driven electronegative discharge

Consider the homogeneous and stationary state of an electronegative discharge with axial symmetry [12], where the homogeneity of the system is in the axial direction. The equilibrium condition can be expressed in terms of the averaged source terms as $${\overline{P}}_{\text{e}}={\overline{P}}_{\text{n}}={\overline{P}}_{\text{m}}=0$$ (these averages are written in terms of the cross section averaged density of the working gas, the details in [22]). Except for equilibrium, the generation of metastables by super-elastic collisions (see Section 2) has practically no influence on the linear properties of the discharge. In order to describe the properties of stability in a cylindrical geometry, the radial and axial dynamics are separated by introducing a new set of fields, $$N_{\text{j}}$$, in Eq. (1) (we follow the work of Bruhn et al. [12]). They are defined as

$$\begin{array}{l}{n}_{\text{e}}(r,z,t)=N_{\text{e}}(z,t)\frac{g(r)}{\overline{g}}\\ n_{\text{n}}(r,z,t)=N_{\text{n}}(z,t)\frac{h(r)}{\overline{h}}\\ n_{\text{m}}(r,z,t)=N_{\text{m}}(z,t)\frac{l(r)}{\overline{l}}\end{array}$$E10

where $$g(r)$$, $$h(r)$$, and $$l(r)$$ are the corresponding radial profile functions [23, 24, 25]. The bar over each one denotes a cross-section average. Then, it is assumed that the current density $$\overrightarrow{j}$$ (Eq. (5)) has the radial electron profile since electrons are the dominant components of the axial electric current. In this way

$$\overrightarrow{j}=j\frac{g(r)}{\overline{g}}{\overrightarrow{e}}_{\text{z}}=\frac{I(t)}{\pi r_0^2}\frac{g(r)}{\overline{g}}{\overrightarrow{e}}_{\text{z}}$$E11

where $$r_0$$ is the radius of the discharge and $$I(t)$$ is the electric current. The same assumption is also made for the mean electron energy $$U_{\text{e}}=U_{\text{e}}(z,t)$$, while the electric field is written as $$\overrightarrow{E}(r,z,t)=E(z,t){\overrightarrow{e}}_{\text{z}}+E(r,t){\overrightarrow{e}}_{\text{r}}$$ since its radial component cannot be neglected. Substituting Eq. (10) into Eq. (1), the cross-section averaged balance equations for the axial part of the particle densities are found

$$\begin{array}{l}\frac{\partial N_{\text{m}}}{\partial t}=D_{\text{m}}\frac{{\partial }^2{N}_{\text{m}}}{\partial z^2}+{\overline{S}}_{\text{m}}\\ \frac{\partial N_{\text{n}}}{\partial t}=D_{\text{n}}\frac{{\partial }^2{N}_{\text{n}}}{\partial z^2}+{\mu }_{\text{n}}\frac{\partial }{\partial z}(N_{\text{n}}{E}_{\text{z}})+{\overline{S}}_{\text{n}}\\ \frac{\partial N_e}{\partial t}=D_e\frac{{\partial }^2{N}_{\text{e}}}{\partial z^2}+(D_{\text{e}}-D_{\text{n}})\frac{{\partial }^2{N}_{\text{n}}}{\partial z^2}+{\overline{S}}_{\text{e}}-{\mu }_{\text{e}}\frac{\partial }{\partial z}(N_{\text{e}}{E}_{\text{z}})-({\mu }_{\text{e}}+{\mu }_{\text{n}})\frac{\partial }{\partial z}(N_{\text{n}}{E}_{\text{z}})\end{array}$$E12

The coefficients in Eq. (12) have been defined in Section 2. Since the explicit expressions, in their most general form, for the averaged source terms in Eq. (12) shed light on the kinetic reactions involved in the dynamics of species existing in the plasma, they are written here:

$$\begin{array}{l}{\overline{S}}_{\text{m}}=k_{\text{m}}\overline{n}{N}_{\text{e}}-k_{\text{s}}{c}_1{N}_{\text{e}}{N}_{\text{m}}-k_{\det }{c}_2{N}_{\text{n}}{N}_{\text{m}}-{\nu }_{\text{w}}^{\text{m}}{N}_{\text{m}}\\ {\overline{S}}_{\text{n}}=k_{\text{at}}\overline{n}{N}_{\text{e}}-k_{\det }{c}_2{N}_{\text{n}}{N}_{\text{m}}+({\nu }_{\text{w}}^{\text{n}}-{\underset{\_}{\nu }}_{\text{w}}^{\text{n}})N_{\text{n}}-k_{\text{rec}}{N}_{\text{n}}(c_3{N}_{\text{e}}+c_4{N}_{\text{n}})\\ {\overline{S}}_{\text{e}}=k_{\text{io}}\overline{n}{N}_{\text{e}}-k_{\text{att}}\overline{n}{N}_{\text{e}}+k_{\det }{c}_2{N}_{\text{n}}{N}_{\text{m}}-({\nu }_{\text{w}}^{\text{e}}-{\underset{\_}{\nu }}_{\text{w}}^{\text{e}})N_{\text{e}}-\left(1+\frac{{\mu }_{\text{p}}}{{\mu }_{\text{n}}}\right){\nu }_{\text{w}}^{\text{n}}{N}_{\text{n}}+\left(1-\frac{D_{\text{p}}}{D_{\text{n}}}\right){\underset{\_}{\nu }}_{\text{w}}^{\text{n}}{N}_{\text{n}}\end{array}$$E13

where $$c_1\mathrm{...}{c}_4$$ are the so-called shape factors and $${\nu }_{\text{w}}^{\text{m}}$$, $${\nu }_{\text{w}}^{\text{n}}$$, ν $${}_{\text{w}}^{\text{n}}$$, $${\nu }_{\text{w}}^{\text{e}}$$, ν $${}_{\text{w}}^{\text{e}}$$ denote the losses of different species in the walls. The metastable working gas molecules are generated by electron impact at the rate $$k_{\text{m}}$$ and they are lost by superelastic collisions at $$k_{\text{s}}$$, as well as by associative detachment at $$k_{\det }$$. The negative ions are generated by attachment at $$k_{\text{att}}$$, and they are lost by ion–ion recombination at $$k_{\text{rec}}$$, where positive and negative ions recombine to form molecular and atomic oxygen. The main generation process of electrons is ionization at the rate $$k_{\text{io}}$$. After averaging, the balance equation for the mean electronic energy can be written as

$$\frac{\partial (N_{\text{e}}{U}_{\text{e}})}{\partial t}=\frac{{\partial }^2\left(D_{\text{e}}^{*}{N}_{\text{e}}{U}_{\text{e}}\right)}{\partial z^2}+\frac{\partial \left({\mu }_{\text{e}}^{*}{N}_{\text{e}}{U}_{\text{e}}{E}_{\text{z}}\right)}{\partial z}+E_{\text{z}}\frac{\partial \left(D_{\text{e}}^{*}{N}_{\text{e}}\right)}{\partial z}+{\mu }_{\text{e}}{N}_{\text{e}}{E}_{\text{z}}^2-N_{\text{e}}{\overline{S}}_{\text{u}}^{\text{e}}-N_{\text{m}}{\overline{S}}_{\text{u}}^{\text{m}}$$E14

where $${\overline{P}}_{\text{u}}^{\text{e}}$$ and $${\overline{P}}_{\text{u}}^{\text{m}}$$ are averaged loss terms given by

$$\begin{array}{l}{\overline{S}}_{\text{u}}^{\text{e}}=\frac{{\mu }_{\text{e}}^{*}}{{\mu }_{\text{e}}}\left(1+\frac{{\mu }_{\text{p}}}{{\mu }_{\text{n}}}\right){\nu }_{\text{w}}^{\text{n}}{U}_{\text{e}}+\frac{{\mu }_{\text{e}}^{*}}{{\mu }_{\text{e}}}\left(\frac{D_{\text{p}}}{D_{\text{n}}}-1\right){\underset{\_}{\nu }}_{\text{w}}^{\text{n}}{U}_{\text{e}}+c_5\left(D_{\text{p}}-D_{\text{n}}\right)\frac{D_{\text{e}}}{{\mu }_{\text{e}}}+c_6\left({\mu }_{\text{p}}+{\mu }_{\text{n}}\right){\left(\frac{D_{\text{e}}}{{\mu }_{\text{e}}}\right)}^2\\ {\overline{S}}_{\text{u}}^{\text{n}}=P_{\text{u}}\left(U_{\text{e}}\right)+{\underset{\_}{\nu }}_{\text{w}}^{\text{u}}{U}_{\text{e}}\end{array}$$E15

The quantities $${\mu }_j$$, $${\mu }_j^{*}$$, and $$D_j$$ were introduced in Section 2.

Equations (12) and (13) can be used to study the stability of the discharge by introducing nondimensional variables which measure how much the system is far from equilibrium. The particle density deviation is nondimensionalized through the corresponding particle density at equilibrium, while the electric field and the mean electronic energy are nondimensionalized by the values at equilibrium. If the new set of variables are denoted by $$n_j^{*}$$, with $$j=\text{m,}\text{n,}\text{e,}{\text{U}}_{\text{e}}\text{,}\text{E}$$, the continuity equations, Eq. (12), can be written up to first order as

$$T_{jk}\frac{\partial n_k^{*}}{\partial \tau }=\left(A_{jk}+B_{jk}\frac{\partial }{\partial \xi }+C_{jk}\frac{{\partial }^2}{\partial {\xi }^2}\right)n_k^{*}+P_{jk}\frac{1}{l}{\displaystyle \underset{0}{\overset{l}{\int }}{n}_k^{*}d\xi }$$E16

where the coefficients $$A_{jk}$$, $$B_{jk}$$, $$C_{jk}$$, $$P_{jk}$$, and $$T_{jk}$$ are known quantities (see Appendix C in [12]), and $$\xi$$ and $$\tau$$ are dimensionless variables of position and time, respectively. The perturbations $$n_k^{*}$$ are expressed as a linear combination of plane waves as follows:

$$n_k^{*}(\xi ,\tau )={\displaystyle \sum _{N=-\infty }^{\infty }{n}_k^{(N)}{e}^{i(k_N\xi -{\omega }_N\tau )}}$$E17

where $$k_N$$ is the real number and $${\omega }_N$$ is the complex number, and $$n_k^{(N)}$$ is the amplitude of the $$k$$ th Fourier component of each specific property. The dispersion relations are obtained from Eq. (16) with the expansion of Eq. (17). The complex roots can each be identified with a potential unstable plasma behavior. The specific value of the root $${\omega }_N$$ from the dispersion relation determines the time evolution of the perturbation. In case $$\mathrm{Re}\left(i{\omega }_N\right)>0$$, the perturbation grows in time which implies an instability. In this way, $$Re\left(i{\omega }_N\right)$$ represents the rate of growth of the perturbation. The characteristic $$\mathrm{Im}\left(i{\omega }_N\right)$$ determines if the wave behavior of perturbations is propagative or dispersive in nature. The calculation of the roots $${\omega }_N$$ may be difficult since it may involve high-order algebraic equations. The analysis of time scales of different relaxation processes can reduce the order of them. For instance, characteristic time of electron energy relaxation is in the range $${10}^{-10}-{10}^{-8}\sec$$, which is very short when compared with the time required for ionization, attachment, etc. [26]. Thus, the time response of the electron energy density to a perturbation is almost instantaneous compared with the time required for a change that occurs in the electric charge density. This allows to omit the term in the time derivative in Eq. (14), which reduces the order of the dispersion relation. The stationary state is stable if all $$\mathrm{Im}\left({\omega }_N\right)$$ are non-positive and if the instability boundary is given by one vanishing $$\mathrm{Im}\left({\omega }_N\right)$$. In Figure 1 (from [12]), it is shown the boundary of the mode $$N=1$$ for $$k>0$$.

The only mode which has influence on the coupling of the plasma with the external circuit is the mode $$N=0$$, which has purely imaginary eigenvalues $${\omega }_0$$ in the considered range of plasma parameters. In Figure 2, the instability curves $${\omega }_0=0$$ of the mode for different values of the external resistor can be seen. Above the curves, all the sets of parameters $$(I_0,k)$$ represent stable mode $$N=0$$. The line marked as lower critical point corresponds to the lower critical current in Figure 1. In this way, all the states of the mode $$N=0$$ inside the instability loop of the mode $$N=1$$ are stable. The stability of the mode $$N=0$$ determines the nonlinear properties of the system since only in that case a Hopf bifurcation arises at the critical points of the mode $$N=1$$ [12].

The linear stability analysis has been used to study the stability of the low-pressure positive column in plasmas of electronegative gases [27, 28]. Testrich et al. [28] considered as relevant components of the plasma: electrons, negative ions (O⁻), positive ions (O₂⁺), and metastable oxygen molecules. Other excited oxygen states were considered to have negligible concentrations. The kinetics included the following: electron attachment and detachment, oxygen ionization, metastable molecules production, and recombination O⁻–O₂⁺. In Figure 3, the comparison of theoretical with experimental results at 0.5 Torr can be seen. As mentioned earlier, the theoretical loop (right) defines the stable and unstable discharge regimes. Waves with dimensionless wave number k within the loop are linearly unstable. On the left, the detected frequency spectra depending on the discharge electric current in different investigated operating modes are shown. The instability windows were in good agreement with experimental data.

Figure 1.

Instability curves for the mode $$N=1$$ for (1) vanishing diffusion coefficients $$D_{\text{p,n}}=0$$ and (2) $$D_{\text{p,n}}\ne 0$$. Critical points are indicated by arrows (from Bruhn et al. [12]).

Figure 2.

Instability curves for the mode $$N=0$$, for three different values of the external resistor (from Bruhn et al. [12]).

Figure 3.

Comparison of experimental results with theoretical investigation on the stability of a DC oxygen discharge. On the left, detected frequency spectra depending on the discharge electric current in different investigated operating modes are shown. The colored plot represents the averaged relative intensity of electric field, 0.1 (blue), 0.25 (cyan), 5 (green), and 100 (orange), in arbitrary units. On the right, the theoretical loop defines the stable and unstable discharge regimes. The experiments were carried out in a pressure range from 0.5 to 0.9 Torr within a discharge current interval from 0.5 to 90 mA. From Testrich et al. [28].

4. Stochastic dynamics in glow discharges

4.1. Plasma dynamics as a stochastic process

It is well known that the random fluctuations of certain observables give information about the transport processes occurring in a given system. This fact is actually a direct consequence of the well-known fluctuation–dissipation theorem which relates the random fluctuations of a given variable to the transport coefficient of the corresponding conjugated thermodynamic variable [29]. It has been used to obtain information on the transport coefficients in glow discharges from the time and space autocorrelation function and their corresponding auto-power spectrum of random fluctuations in plasmas [30]. Moreover, it has been shown experimentally that the random current fluctuations exhibit, under certain conditions, the Markovian property in the fully developed turbulence regime [31].

The above-mentioned works have in common that the random fluctuations are assumed to be homogeneous and isotropic in space. Actually, it is known that glow discharges are highly inhomogeneous, and therefore the fluctuations might vary in space which might be due to nonlocal effects, finite size confinement, instabilities, etc. The study of the influence of inhomogeneities in glow discharges has revealed the presence of two particle transport regimes, namely, a purely diffusive transport regime and a wave-diffusive transport regime [32, 33]. The non-locality, memory effects, and other phenomena inducing inhomogeneities in glow discharges have been introduced in the referred works by considering a space-dependent diffusion coefficient. Within this treatment, it has been found that every type of transport regime induces a characteristic type of decay in the correlation function. Indeed, the purely diffusive regime induces a decaying mode, while the wave-like transport regime induces an oscillatory decay mode. These two types of behavior have actually been observed in fusion plasmas in tokamak and non-tokamak devices ([32] and references therein). It is important to point out that it is possible to transit from one behavior to the other by tuning a parameter value [32]. Moreover, it was also found theoretically that it is actually possible to observe the two transport regimes (diffusive transport and wave-like transport) in a single device at different locations where the plasma density has substantial differences.

4.2. FM model for fluctuations in glow discharges

In this section, we explore the occurrence of the above-described phenomenon, namely, the coexistence of the diffusive transport regime and the wave-like transport regime in glow discharges. For such purposes, we implement numerical experiments with a closely related model to that used in Ref. [32]. The model we consider here corresponds to the case of Eqs. (1)– (3) where stochastic properties are introduced through additive random terms in the constitutive equations as a stochasticity source of the particle fluxes (Eq. (2)). The stochastic terms are intended to model spontaneous rapid variations in the particle fluxes which might be due to some phenomena which are not present in the constitutive equations but which would be relevant for the dynamics of the plasma fluctuations. It is clear that Eq. (2) can be considered as an approximation of a larger set of constitutive equations [9, 32]. Of course, we introduce the stochastic terms as a simple way to model (in a rather heuristic way) such phenomena here. Actually, there is a formal way to do this, but our intention in this section is just to study the phenomena induced by the introduction of those stochastic “corrections.” We leave, as a future work, the introduction, in a more rigorous way, of those stochastic terms. For the sake of completeness, we write in the simplest case the stochastic fluid model obtained from the above-mentioned considerations [34]. We consider only two species, namely, positive ions and electrons in an external electric field generated by parallel plates. First, we have the continuity equations given by

$$\begin{array}{l}\frac{\partial n_{\text{p}}}{\partial t}+\nabla \cdot {\overrightarrow{J}}_{\text{p}}=S_{\text{p}}\\ \frac{\partial n_{\text{e}}}{\partial t}+\nabla \cdot {\overrightarrow{J}}_{\text{e}}=S_{\text{e}}\end{array}$$E18

where the source terms are written to include only the ionization and the electron–ion recombination rates.

Next, we introduce the stochastic version of the constitutive equations as

$$\begin{array}{l}{\overrightarrow{J}}_{\text{p}}={\mu }_{\text{p}}{n}_{\text{p}}\overrightarrow{E}-D_{\text{p}}\nabla n_{\text{p}}+{\widehat{f}}_{\text{p}}\\ {\overrightarrow{J}}_{\text{e}}=-{\mu }_{\text{e}}{n}_{\text{e}}\overrightarrow{E}-\nabla \left(D_{\text{e}}{n}_{\text{e}}\right)+{\widehat{f}}_{\text{e}}\end{array}$$E19

where $${\widehat{f}}_{\text{p}}$$ and $${\widehat{f}}_{\text{e}}$$ are assumed to be $$\delta$$ -correlated white noise terms with zero mean and Gaussian distribution. Formal aspects of the inclusion of stochastic terms in the constitutive equations can be found in the context of the general equation of the reversible–irreversible coupling, GENERIC [35, 36]. Another approach based on the Markovian method can be seen in Ref. [31]. Within the first formalism, it is shown that fluctuating phenomena which are described by dynamical equations having the GENERIC structure are stationary, Gaussian, and Markovian processes. In general, as it is the case here, with the use of this kind of stochastic terms in Eq. (19), one can incorporate the statistical properties of stationary, Gaussian, and Markovian processes to the mean Eqs. (1)– (3). The set of equations is completed with the Poisson equation describing the self-consistent electric field and the electronic energy conservation equation (Eq. (3)). The former is given by

$${\nabla }^2\phi =\frac{e}{\epsilon }\left(n_{\text{e}}-n_{\text{p}}\right)$$E20

where $$\phi$$ is the electric potential function.

We close this subsection with the reduced expressions for the source terms in Eq. (18). In fact, we have

$$S_{\text{p}}=S_{\text{e}}=k_{\text{io}}\overline{n}{n}_{\text{e}}-k_{\text{att}}\overline{n}{n}_{\text{e}}$$E21

In the last equation, $$k_{\text{io}}$$ and $$k_{\text{att}}$$ are the ionization and recombination ($$e-p$$) coefficients, respectively, and $$\overline{n}$$ the bulk gas density.

4.3. Numerical method of solution of FM model

Our main goal is to study the phenomenon of coexistence of the two regimes of transport previously reported in plasmas [32]. To this end, we numerically simulate the dynamics of Eqs. (18)– (21) in a one-dimensional “box” of length L = 1 (see Figure 4) with a noise level of D = 0.10. The boundary conditions needed by Eqs. (18)– (20) are specified at each electrode. Accordingly with Figure 4, the grounded anode is at x = 0, and the powered cathode is at x = L. At the anode, the electron density is taken as zero, since it is assumed here that the recombination rate $$k_{\text{att}}$$ is infinite as it corresponds to a grounded conducting surface. It is also assumed that the ion flux at position x = 0 is due only to drift. $$\phi =0$$ is the given value of the plasma potential at x = 0. At the cathode, the electronic flux is proportional to the ion flux which is also assumed to be due only to drift.

The numerical method used here is based on the spectral Chebyshev collocation (SCC) method which assumes that an unknown partial differential equation (PDE) solution can be represented by a global, interpolating, Chebyshev partial sum. This form is often convenient, particularly if the solution is required for some other calculation where its rapid and accurate evaluation is important. The global character of spectral methods is beneficial for accuracy. Once the finite Chebyshev series representing the solution is substituted into the differential equation, the coefficients become determined so that the differential equation is satisfied at certain points within the range under consideration. In this spectral method, the PDE equation is required to be satisfied exactly at the interior points, namely, the Gauss–Lobatto collocation points given by

$$x_i=\cos \left(\frac{i\pi }{N}\right)$$E22

with $$i=1,\mathrm{...},N-1.$$

Figure 4.

The device with a parallel-plate configuration. The anode and the cathode are conducting surfaces. The anode is grounded and the cathode powered by the shown power supply. The anode is situated at x = 0 in the reference frame and the cathode at x = L. The figure also shows the narrow box within which the simulations are made with the use of Eqs. (18)– (21). The separation between the plates is L.

Figure 5.

Spatial distribution of electron and ion densities. The anode is at x = 0 and the cathode at x = 1 (nondimensional position). Both the distributions show maxima at x = 0.14. The plasma is not in the neutral condition along the spatial domain.

In equation (22), N denotes the size of the grid. The number of points is chosen so that, along with the initial or boundary conditions, there are enough equations to find the unknown coefficients. The positions of the points in the range are chosen to make the residual obtained small when the approximate solution is substituted into the differential equation. The range in which the solution is required is assumed finite, and, for convenience, a linear transformation of the independent variable is made to make the new range (−1, 1).

Finally, an additional property of the spectral methods is the easiness with which the accuracy of the computed solution can be estimated. This can be done by simply checking the decrease of the spectral coefficients. There is no need to perform several calculations by modifying the resolution, as is usually done in finite difference and similar methods for estimating the “grid convergence.” A further explanation of this spectral method can be found in Refs. [37] and [38]. The spatial derivative terms in Eqs. (18) and (19) were expressed on derivative matrices expanded on Chebyshev polynomials. The matrix-diagonalization method was used to solve the coefficient equation system in physical space directly. A coordinate transformation was necessary to map the computational interval to 0 <x <1.

Figure 6.

Spatial distribution of the plasma potential and the self-consistent electric field. The electric field has the expected behavior, namely, it decays almost linearly from the cathode, and then it becomes negative in the anode sheath. The distributions seen in Figures 5 and 6 indicate that the negative glow does not exist for the plasma conditions assumed here.

We then obtain the density of electrons and ions as well as the plasma potential and the self-consistent electric field resulting in the spatial distributions shown in Figures 5 and 6. The simulation in time at three different positions, namely, at x = 0.14, x = 0.5, and x = 0.86, of the positive ion density results in time series representing the density variations in time. Representative numerical time series can be observed in Figure 7a. Each time series is different from each other due to the spatial inhomogeneities of the plasma and the presence of fluctuating phenomena. These behaviors can be qualitatively compared with experimental results that we have obtained in a DC discharge shown in Figure 7b. The measurements were made in a hollow cathode device (90 mm long, 29.5 mm radius) which is depicted in Figure 8. Then we can calculate the autocorrelation function of each time series and compare the resulting autocorrelation functions in the different locations we used to obtain the density variations (Figure 9). The smoothed curves shown in Figure 9 have been averaged over a large number of numerical realizations corresponding to each of the considered positions.

4.4. Discussion

As it can be seen in Figure 5, the electron and ion densities have a maximum at x = 0.14, while in the other two points (at x = 0.5 and x = 0.86) the densities show smaller values. From the graph shown, it can be concluded that the plasma is in a non-neutrality condition in most of the domain [0,1]. In the sheath regions, the electron density is nearly zero. The plasma potential and the self-consistent electric field are shown in Figure 6. The electric field has the expected behavior, namely, it decays almost linearly from the cathode, and then it becomes negative in the anode sheath. The distributions seen in Figures 5 and 6 indicate that the negative glow does not exist for the plasma conditions assumed here.

Figure 7.

Time series of the positive ion density at three positions in the profile. The series are generated by a) numerical simulation of the time evolution of the density through stochastic version of the model, Eqs. (18)– (21), and b) experimental measuring of the floating potential at each position. The insets show the rapid variations occurring during a longer lapse of simulation in each case. Time is in a.u.

In Figure 9, we observe that the density autocorrelation function has different behaviors at different positions. Though the shapes of the shown curves are very similar, the comparison reveals that in the plasma two different particle transport processes are present, the so-called diffusive and the wave-like transport regimes. These two regimes are characterized by its autocorrelation function as anticipated in Ref. [32]. The autocorrelation function for the pure diffusive transport would be a strictly monotonous decreasing function of the delay time, while the correlation function for the wave-like transport would be an oscillating decaying function. These pure modes are not indeed present in the density variations we have obtained from the numerical simulation we performed, but a transport regime resulting from the coexistence of both of the mentioned regimes is observed. We name this the wave-diffusive regime. The tendency shown by the three correlation functions in Figure 9 reveals that the “less oscillatory” curve (solid line) corresponds to the position at x = 0.14 where the electronic density is maximum. At this position, one should indeed expect a “more diffusive” behavior. We note that the oscillating decaying autocorrelation functions correspond to measurements where the plasma has small diffusivity gradients [33], while in the diffusive transport regime where the autocorrelation function has a purely decaying mode, the plasma has high diffusivity gradients.

Figure 8.

Schematic diagram of the hollow cathode DC discharge device.

Figure 9.

Numerical autocorrelation function for the plasma density fluctuations. The autocorrelation function is calculated from the measurements of the density as a function of time at three different locations in the box. The autocorrelation with a more oscillating decaying mode (solid line) corresponds to measurements near the anode at x = 0.14 where the electronic density is maximum. Other measurements (at x = 0.50 and x = 0.86) result in a less oscillating decaying autocorrelation function. This means that the diffusive transport regime is present in regions where the plasma has a high density of electrons, while the wave-diffusive transport regime is present in regions where the electron density is low. Time is in a.u.

Finally, we should mention that this work can be extended in several directions. For example, if we think of the plasma as a system out of equilibrium, we could explore the phenomenon of resonant response reported in Ref. [39]. This phenomenon occurs when the system relaxes to its stationary state in an oscillatory decaying mode. It has been proved that such a frequency actually can be used to stimulate the system in the linear response regime, thus leading to a response in the density fluctuations of the system.

5. Control of nonlinear dynamics in discharge plasmas

5.1. Overview of nonlinear dynamics in discharge plasmas

Among the very interesting nonlinear systems are gaseous plasmas. As mentioned, plasmas are governed by dissipative effects and are far from equilibrium. They exhibit complex or regular behavior depending on the controlling parameters such as discharge current, voltage, and pressure. The transition from one state to another takes place for a change in any one of the parameters [40]. An important feature is the understanding of the complex chaotic behaviors commonly observed in glow discharges. The transition to chaos and constructive effects of noise, such as stochastic resonance, and coherence resonance have been studied in many glow discharge plasma systems [41, 42].

Nonlinear dynamical systems present several routes to chaos such as period doubling, intermittency, quasiperiodicity, etc. [43]. Very important and different transition sequences from periodic to chaotic patterns have also been observed. An example is the so-called alternating periodic chaotic (APC) sequence. Much of the behavior concerning APC sequences is understood in terms of a homoclinic chaotic scenario [44]. Homoclinic orbits are associated with erratic behavior (chaos) in a dynamical system [45]. Homoclinic chaos in discharge tubes dates back to more than twenty years, as well as the chaotic behavior and period doubling in pulsed plasma discharge [46]. Currently, much attention has been given to the generic problem of phase synchronization using these devices. In particular, phase synchronization has been demonstrated under different conditions [47]. Given a chaotic oscillator for which an angle coordinate can be suitably introduced as a state space variable, it is often the case that the phase of this oscillator synchronizes with the phase of an external periodic perturbation, or pacer [48, 49].

Chaos may be undesirable for industrial applications where cycle-to-cycle reproducibility is important, yet for the treatment of cell-containing materials including living tissues, it may offer a novel route to combat some of the major challenges in medicine such as drug resistance. Chaos in low-temperature atmospheric plasmas and its effective control are likely to open up new vistas for medical technologies. In the general framework of nonlinear dynamical systems, a number of strategies have been developed to achieve active control over complex temporal or spatio-temporal behavior [50]. Many of these techniques apply to plasma instabilities. There are some methods of controlling chaos where the main idea is to convert chaotic behavior to periodic behavior by inducing a small perturbation in the system. Depending on the physical mechanism of the specific instability in each case, an appropriate control strategy is chosen out of a variety of different approaches; in particular, discrete feedback, continuous feedback, or spatio-temporal open-loop synchronization [51, 52]. In particular, time-delayed feedback plays a prominent role in controlling chaos; this is one of the successful applications that knowledge acquired in nonlinear science has provided to plasma physics. Pyragas [53] proposed the time-delayed feedback technique, which is based on feedback perturbation in the form of the difference between a delayed output signal and the output signal itself; this is appropriate for laboratory experiments conducted in real time. This method is robust to noise and does not require real-time computer processing to calculate a target unstable periodic orbit (UPO); therefore, it can act on the experimental system continuously over time. The feedback perturbation signal that is applied to the nonlinear system is proportionally adjusted to the difference between the two successive values of an arbitrary dynamic variable. Since plasma is a typical nonlinear dynamical system with a large number of degrees of spatiotemporal freedom, various unexpected phenomena are observed when time-delayed feedback is applied in plasma.

Therefore, investigations into the behavior of nonlinear systems with regard to time delay are currently required. From a theoretical point of view, discharges have been traditionally described taking into account the complex processes involved in the plasma recombination and electric conductivity. Such descriptions require, as mentioned in Section 2, the use of coupled partial differential equations involving spatial and time variables, the transport of momentum and energy of plasma components, the continuity equation, diffusion equation, Poisson equation, etc.

5.2. Time-delayed feedback control

Several algorithms to control chaotic dynamics have been proposed in the last 20 years ([54, 55] and references therein). They include OGY method [56], linear feedback [57, 58], nonfeedback method [59, 60], adaptive control [61, 62], backsteeping method, and sliding mode control. Delayed feedback control (DFC) has been applied to a large variety of systems in physics, biology, medicine, and engineering [63, 64], in purely temporal dynamics as well as in spatially extended systems ([65] and references therein). DFC has also been used to control purely noise-induced oscillations and pattern formation in steady-state regimes in chemical systems [65, 66], electrochemical cells [67], neural systems [68], laser diodes [69], and semiconductor nanostructures [70].

DFC control comprises various aspects such as the stabilization of unstable periodic orbits embedded in a deterministic chaotic attractor, stabilization of unstable fixed points (steady states), or control of the coherence, and timescales of stochastic motion. Moreover, it has been shown to be applicable also to noise-induced oscillations and patterns. DFC has been widely used to control chaos in discharge plasmas. Several interesting reports can be found in the literature [51, 52, 65, 67]. Among them, for instance, Fukuyama et al. [52] applied the method to study the nonlinear periodic regime associated with the current-driven ion acoustic instability (IAI) in a double-discharge plasma system. The discharge chamber was divided at the center into a driver region and a target region by a separating grid that was kept at a floating potential. A DC potential was applied to the other grid in order to excite the instability. The current-driven IAI appeared when the DC potential reached a threshold. Once the instability is established (the threshold was 23 V), the dynamics of the floating potential shows a variety of behaviors. First, a limit cycle appears and persists for values of the DC potential between 23 and 35 V. The amplitude of the limit cycle switches stochastically between two values for the DC potential between 36 and 39 V, and a larger limit cycle appears for 40 V. After this, as the DC potential increases, the system gradually falls into disorder. Then the chaotic state appears at 54 V when the system is completely disturbed. When the potential is 67 V, the instability disappears with a decrease in the noise level. DFC is applied to the nonlinear periodic regime for values of the DC potential between 40 and 45 V. The feedback perturbation signal $$F(t)$$ applied to the nonlinear system was proportionally adjusted to the difference between two successive values of a dynamic variable $$x(t)$$, i.e.,

$$F(t)=k\left[x(t-\tau )-x(t)\right]$$E23

where $$\tau$$ is the time delay and $$k$$ a proportionality constant. The behavior of the system was investigated by changing the values of $$\tau$$ and $$k$$ when time-delayed feedback is applied to the nonlinear periodic regime. Intermittency and chaos were observed.

In the remaining of this section, we examine the effects of DFC method on a special periodic regime in discharge plasmas which are often observed in experiments. This can be seen in Figure 10, where the autonomous dynamics as reflected by the time evolution of the positive ion density when the discharge voltage ($$V_{\text{D}}$$) has been applied is shown. The gas is supposed to be in a parallel electrode device (Figure 4). One of the electrodes is grounded, and the other is at $$V_{\text{D}}$$. Here, we test a DFC method on periodic oscillations with relaxation predicted by the FM for the discharge plasma given by Eqs. (18)– (21). The graph in Figure 10 is obtained by using the deterministic form of Eqs. (18)– (21) with an additional term in continuity equations through which the control is applied as suggested by Pyragas [53]. The continuity equation for the positive ion density is then modified as follows:

$$\frac{\partial n_{\text{p}}}{\partial t}+\nabla \cdot {\overrightarrow{J}}_{\text{p}}=S_{\text{p}}+F(x,t)$$E24

where $$F(x,t)$$ is the parameter control given by the formula

$$F(x,t)=k\left[n_{\text{p}}(x,t-\tau )-n_{\text{p}}(t)\right]$$E25

with $$k$$ a constant. As it is well known, Eq. (25) for the control term ensures that the system will not drive to a regime where the target dynamics were naturally stable. Equations (24) and (19)– (21) are solved by means of the spectral Chebyshev collocation numerical method described in Section 4.3. The value of the constant $$k$$ is chosen to optimize the effect of the feedback on the system’s oscillatory regime. This process is made by trial and error. The value of the delay time $$\tau$$ is directly related with the period of the rapid oscillations observed between the oscillations in Figure 10.

Figure 10.

Time evolution of the positive ion density. The dynamics, obtained by numerical solution of Eqs. (18)– (21), shows a periodic behavior with relaxation. The discharge voltage is fixed in 460 V.

The observed periodic oscillation behavior with relaxation in Figure 10 was experimentally found by, for instance, Nurujjaman and Sekar Iyengar [71] in the cathode sheath in a cylindrical DC discharge device by measuring the floating voltage in the plasma which is directly related with the charge density (see also [42] and [72]). It can be seen in Figure 11 that the power spectrum and the phase space plot were included.

Figure 11.

a) Measured floating potential in a discharge plasma at 292 V, b) power spectrum, and c) phase space plot. PS, power spectrum. From Nurujjaman and Sekar Iyengar [71].

Figure 12 shows the result of the application of the control method sketched in Eqs. (24) and (25) to the oscillatory regime between $$t=10000$$ and $$t=20000$$. As observed, the stabilization of the dynamics is easily achieved with great efficiency. The optimum value of the constant $$k$$ used to obtain Figure 12 is 0.1, while other constants are $$\tau =5$$. We leave, for future work, implementing the control method to an experimental plasma. In the laboratory, the density $$n_{\text{p}}(x,t)$$ should be replaced by an accessible parameter of the system, e.g., the floating potential and the control should be externally applied on the system.

Figure 12.

Delayed feedback control is applied to the time series of positive ion density from Figure 10. The control is turned on at t = 10000 and turned off at t = 20000. Time is in a.u.

6. Conclusion

Low-temperature (LT) plasmas have a substantial role in diverse scientific areas and modern technologies. Industrial applications of LT plasmas form a very important part of the productive infrastructure of advanced economies all over the world. In LT plasmas, a large number and diversity of reactive species are generated that activate physical and chemical processes hard to obtain in ordinary chemical environments. Their stochastic, nonlinear dynamics, and chemical kinetics strongly determine the efficiency and effectiveness of LT plasma-based procedures. Typical applications are etching, spectrochemical analysis, deposition of thin films on substrates, and others. Much is known on low-temperature plasmas by fluid model-based investigations (Sections 2 and 3). Fluid models have been used to study the fundamental properties and dynamic behavior of several gases: argon, oxygen, hydrogen, and helium in a diversity of configurations and its results confirmed by experiments. In Section 4, we analyzed the diffusive and wave-like transport of particles in LT plasmas (glow discharges). To this end, we used the well-known relationship between the random fluctuations of physical observables and the transport processes occurring in the system. We concluded that both the diffusive transport and the wave-like transport of particles coexist in the numerical simulations and experiments that were carried out. Section 5 was devoted to analyzing the state of the art of the nonlinear dynamics of LT plasmas; in this section, besides discussing the wide variety of nonlinear behaviors shown by LT plasmas, we illustrated the case of particle density oscillation control through a time-delayed feedback technique. The understanding of the many complex chaotic behaviors commonly observed in glow discharges still offers an exciting field for scientific research.