In this work, since the area of the electrode is much larger than the discharge gas gap, the drift–diffusion model can simulate the development of the discharge. Four species mentioned previously are mainly considered in this model. The Ar^* is in the 4s excited state, and the lifetime of the 4s excited state of Ar^* is 100–140 ns. In the vicinity of this excited level, there are many semblable levels, and they can transform with each other quickly; hence, we consider these particles as metastable particles. The spatial and temporal distribution of particles is described with the continuity equation: ∂ n r ∂ t + ∇ ⋅ N r = S r ⋅ (1) Here, n _r is the number density of r-th particles, including e, Ar^*, Ar⁺, and Ar₂ ⁺. S _r is the source term of r-th particles, and it is calculated by S r = ∑ i k i n r 1 n r 2 + α i N e n r 3 . (2) Here, k _i is the reaction coefficient; α _i is the Townsend coefficient; n _r1, n _r2, and n _r3 are the reaction number densities in different equations. The N _r in Eq 1 is the r-th particle flux, and it is described as N A r * = − D A r * ∇ n A r * (3) N i = − D i ∇ n i + μ i n i E (4) N e = − D e ∇ n e − μ e n e E (5) Here, i represents Ar⁺ and Ar₂ ⁺; E is the electric field; μ _e and μ _i are the mobility rates of electron and ions, respectively. The ion mobility rate is the function of the reduced electric field (E/n), and the electron mobility rate is the function of the mean electron energy, which is solved by Boltzmann equation. D _e, D _i, and D _Ar* are the diffusion coefficients of e, ion, and Ar*, respectively, which are calculated by Einstein relation: D i = k B T i μ i e (6) Here, e is the elementary charge, T _i is the temperature of the i-th particle, and k _B is the Boltzmann constant. In this work, the electric field is solved by Poisson’s equation: − ∇ 2 φ = e n A r + + n A r 2 + − n e ε (7) Here, φ is the electric potential and ε is the dielectric constant. The discharge circuit is calculated by Kirchhoff’s voltage law (KVL): V a t = V d t + R 1 i 1 t (8) Here, V _a(t) is the applied voltage, V _d(t) is the dielectric voltage and the loop circuit, and i ₁(t) is obtained by the Morrow–Sato equation.

In this model, we assume the positive ion particles are totally recombined when they impact on the cathode dielectric layer, and hence the boundary condition cathode is considered as N i ⋅ u n = 1 4 v i n i + α i 0 μ i n i E ⋅ u n ⋅ (9)

Here, v i = 8 k B T i π m i is the ion thermal velocity, m _i is the mass of ion, u _n is the normal vector which directs toward the surface of the dielectric layer, and α _i0 is the coefficient which relates to the electric field of the boundary and is expressed as α i 0 = 1 , E ⋅ u n > 0 0 , E ⋅ u n ≤ 0 (10)

The secondary electron emission induced by the positive ion impact on the cathode dielectric surface has been considered by N e ⋅ u n = 1 4 v e n e − α e 1 μ e n e E ⋅ u n − α e 2 γ N A r * + N A r 2 + ⋅ u n (11) Here, γ is the secondary electron emission coefficient; α _e1 and α _e2 are the coefficients which relate to the electric field of the boundary and are expressed as α e 1 = 0 , E ⋅ u n > 0 1 , E ⋅ u n ≤ 0 (12) α e 2 = 0 , E ⋅ u n ≤ 0 1 , E ⋅ u n > 0 (13)

The boundary of electron energy is related to the electrons, and hence they have a similar form. Since the generation mechanism of secondary electron emission is fairly complicated in dielectric barrier discharge, we simplify the mean secondary electron energy as 2.5 eV, which is approximately identical to the mean electron energy near the plasma sheath. The boundary of electron energy is simplified as N ε ⋅ u n = 1 4 v e n ε − α e 1 μ ε n ε E ⋅ u n − 2.5 α e 2 γ N A r + + N A r 2 + ⋅ u n (14) Here, v e = 8 k B T e π m e is the electron thermal velocity, N _ε is the electron energy flux, and n _ε denotes the energy density of the electron.

For the metastable atom, we assume they totally de-excitate to the ground state atoms, and the boundary condition is written as D A r * ∇ n A r * ⋅ u n = 0 (15)

For the electric field, the dielectric surface can simulate the charge, and the electric flux in the dielectric surface is discontinuous. D g − D a ⋅ u n = σ (16) Here, D _g is the electric flux in the ground electrode, D _a is the electric flux in the applied excitation electrode, and σ is the surface charge density obtained by σ = ∫ t 1 t 2 e N A r + + N A r 2 + − N e ⋅ u n d t (17)

The related species reactions and their associated transport coefficients in this work are shown in Table 1. The sketch of the DBD model in this work is shown in Figure 1.

We first verify the validity of the model by comparing the numerical and experimental results for similar conditions in single LF or RF DBD. The discharge plasma is generated and sustained in a geometrically symmetric parallel-plate DBD configuration. The dielectric material we adopt is polyethylene terephthalate (PET), whose relative permittivity (dielectric constant) and thickness are ε = 3.4 and d _B = 1 mm, respectively. The gas gap between the two layers is fixed as d _g = 2 mm. The secondary electron emission coefficient of the dielectric material is 0.07. The gas temperature at the dielectric layer boundary is fixed at T _g = 400 K ( ≈ 0.04 eV). The initial ions and electrons are uniformly distributed in the discharge region with the same density of 10⁷ cm^-3. Additionally, 600 equidistant grids are used to resolve the 1-mm inter-electrode gap, and 2×10⁵ time-steps are applied within an applied voltage cycle to resolve the temporal dynamics of all plasma species.

Figure 2 presents a typical discharge waveform, including applied and discharge current density, driven by a 50 kHz@1.2 kV LF (the upper row) as well as a 5 MHz@0.7 kV RF (the lower row), respectively, and the corresponding spatio-temporal distributions of the electric field and electron density during two discharge periods. The discharge current density is characterized by a short pulse per half cycle in the LF driving case and approximately sinusoidal in the RF case. As the driving frequency increases from ∼kHz to ∼ MHz, both the electric field and the electron density increase gradually, due to the higher electrical oscillation compressing the sheath electron density. Both numerical results of the individual LF or RF discharge presented here are clearly in good agreements with previous numerical and experimental results under the similar condition [24].

The governance pattern of DF discharges can be simply described by the ratio of LF to RF voltage amplitudes with fixed dual frequencies. In this work, with f _L = 50 kHz/f _H = 5 MHz, we gradually vary the LF (50 kHz) voltage from 0 to 1500 V and alter the RF (5 MHz) voltage from 0 to 400 V. Figure 3 presents typical discharge waveforms under various voltage amplitude ratios of LF/RF. For the LF-dominant discharge mode with a ratio of V _LF = 850 V/V _RF = 150 V, as shown in Figure 3A, the discharge current density (red line) oscillates in the high frequency of the RF and is enveloped by a pulse-like waveform during the low-frequency period. When reducing the ratio to V _LF = 700 V/V _RF = 300 V, although the voltage curve (blue line) is still LF-dominant, the discharge current density waveform is in an LF-modulated quasi-sinusoidal oscillation, as shown in Figure 3B. We then call the discharge mode in this case as a DF-comparable mode. The discharge current pattern clearly distinguishes various discharge mechanisms and plasma parameter regimes. Furthermore, lowering the ratio to V _LF = 100 V/V _RF = 400 V, the discharge is then clearly dominated by the RF component, and the discharge current density exhibits a completely high-frequency oscillation, with the discharge characteristic not essentially different from the RF-only driving. Hence, its feature is similar to that of the HF-only discharges.

To better understand and distinguish the aforementioned governance patterns, we show in Figure 4 the spatial distributions of feature parameters, including electron/argon ions/metastable ion densities and electron temperature, in the discharge gas gap, which are temporally averaged over an LF period. Obviously, one can see a sheath region near the boundary, characterized by the significant variations of particle densities (n _e, n _Ar+, and n _Ar*). For the LF-dominant discharge, as shown in Figure 4A, the electron density decreases substantially near the boundary to form a sheath. In addition, the LF component dominates the discharge and induces the cathode sheath expansion and decreases the plasma bulk region and the density compared with the DF-comparable mode and the HF-dominant mode, as shown in Figures 4B, C, respectively. This can be explained by the increase in the ion loss due to the drift to the cathode by the strong LF electric field [24]. The electron temperature in the bulk plasma decreases to its bottom level at x = 0.5–1.5 mm, due to the higher electron power absorption in this region. For the single-frequency dielectric barrier discharge, the electron temperature has a strong coupling with electron power absorption, i.e., if the power absorption increases so does the electron temperature. For the dual-frequency discharge, the higher frequency oscillation within a period of low frequency can effectively capture the electron and trap electrons within plasma sheaths. The sheath thickness is similar to that in the RF-only DBD discharge, but the electric field in the sheath region is much stronger. The peaks of particle densities are located at the central region where the electron temperature reaches its minimum.

For the DF-comparable mode, as shown in Figure 4B, as the proportion of the HF component increases, all particle densities increase while the electron temperature steeply decreases in the plasma bulk region. The minimum value of electron temperature occurs at the sheath edge due to the enhanced direct ionization and strong RF electric field. Furthermore, the discharge in this mode is characterized by a cathode reduction on both sides of the plasma bulk region, meaning the density of ions is always large enough to cause the RF electron oscillation, which induces a space–charge electric field. For the HF-dominant mode, as shown in Figure 4C, the spatial distributions of time-averaged particle densities and electron temperature are similar to those in the RF-only DBD discharges [41, 42]. Hence, the discussion for features of the HF-dominant mode is also similar to that in the previous RF-only DBD studies.

In Figure 5, we further show the time-averaged particle densities varied with LF and RF voltage amplitudes. In general, it is observed that all time-averaged particle densities are positively proportional to V _RF. For LF voltage, however, although n _e, n _Ar*, and n _Ar+ increase with the increase in V _LF when V _RF = 350 V, they decrease first and then increase with a turning point approximately at V _LF = 200 V when V _RF = 250 V or 300V. The possible reason is due to the uncompensated ion density loss at a low LF component voltage, as shown in Figure 5D. The ratio of electrons and Ar ions shows the Ar ion density has a rapid increase when V _LF>300 V. As V _RF reaches a high level, e.g., at V _RF = 350 V, the RF component is large enough to compensate the ion loss, and hence the particle densities keep increasing with the LF voltage amplitude, and also the ratio of electrons and Ar ions presents a relatively steady increase. The electron density is enhanced by a factor of 4.6 (the ratio of V _RF = 350 V, V _LF = 1400 V, and V _RF = 250 V/V _LF = 100 V) by changing the LF/RF voltage amplitude ratio. Correspondingly, the densities of Ar^* and Ar⁺ densities can be enhanced by a factor of 2.7 and 2.5, respectively.

To gain a clear picture of the density evolution, the spatial distributions of particle densities developed in two LF periods are shown in Figure 6. In Figure 6A, one can find that for discharges dominantly controlled by LF, the maximum electron density is located near the sheath region. Modulated by the RF component, the electron density in the bulk plasma region maintains a high level, different from the electron density distribution in only LF-driven DBDs. The similar spatial distribution tendency can be observed for n _Ar+. For the distribution of n _Ar+, the maximal value is at ∼400 μm from the instantaneous cathode, near the sheath region. In the LF-dominant mode, n _Ar+ has a similar distribution to only LF-driven DBDs, with one peak in each LF period, but of a higher intensity. The electron-depleted region oscillates with the RF driving. Moreover, even a weak RF component can notably expand the bulk plasma region, causing further ion loss at the instantaneous cathode.

For the DF-comparable mode, the electron is mostly distributed in the plasma bulk region, with a significant electron density decrease near the sheath region, as shown in Figure 6D. In this case, LF-only DBD discharge characteristics completely disappear due to the considerable RF modulation. The Ar* species majorly distributes near the sheath region and is continuously created in the LF period. It is worth noting that all particle densities have a reduction in comparison with the LF-dominant mode. However, electron and argon ion densities only decrease slightly, while the argon metastable density decreases obviously by almost an order of magnitude. This suggests that the LF component mainly affects the metastable particle generation.

Figures 7A, B shows the spatio-temporal profiles of the power absorption (P _e,abs = J _e· E ) as a function of the RF voltage amplitude in two LF periods with V _RF = 250 V and 350 V, respectively. It is observed that the major power absorption occurs near the plasma bulk region, which mainly results from RF ohmic heating since the power absorption intensity increases with the RF voltage. Thus, electrons can be heated throughout the LF period. To further investigate electron dynamics, the power absorption in five continuous RF periods is shown in Figures 7C, D. One can find that the discharge remains in an α-γ transition mode. In the α-γ transition mode, the plasma is driven by the electron current density and SEE. The power absorption is majorly observed in both the expansion and collapse areas of the sheath and extended into the sheath. As the RF component is enhanced, the plasma sheath is compressed and the power absorption inside the sheath increases obviously. By further increasing the RF component, the discharge then completely enters into the γ-mode (not shown here).

When the discharge is dominantly controlled by the LF driving, the power absorption presents a different spatio-temporal profile, as shown in Figures 8A, B. In the figures, we can see that major power absorption occurs near the sheath region, and the major excitation peaks around each LF power maximum, which is quite different from the DF-comparable mode. Increasing the RF component to 350 V, the RF oscillation begins to be apparent, as shown in Figure 8B. To show the process more clearly, the power absorption in five continuous RF periods is plotted in Figures 8C, D. In the figures, one can find an asymmetric structure of the electron power absorption. For the instantaneous cathode side, the electron power absorption occurs near the sheath region; for the instantaneous anode side, a weaker power absorption occurs at both expansion and collapse regions of the sheath. Furthermore, by increasing the RF component, the power absorption is enhanced at both sides and gets extended to the bulk plasma region from the instantaneous anode, which suggests that more secondary electrons are trapped by the enhanced RF electric field.

Furthermore, we investigate the influence of the DF discharge on the plasma parameters. Fast Fourier transform (FFT) of the discharge current and their associated power spectra for different discharge modes is carried out and shown in Figure 9. As a rapid and effective method to analyze non-linear phenomena, we and other groups have also applied FFT in previous studies [23, 26]. From the frequency spectrum (0–25 MHz) shown in Figure 9A, one can see that, for the discharge dominantly controlled by the LF component, apart from the fundamental frequencies of ω _L = 50 kHz and ω _H = 5 MHz, the third RF harmonic (3ω _H = 15 MHz) and sum + beat of the LF and the second RF harmonic (2ω _H + ω _L and 2ω _H - ω _L), as well as side-bands, are also excited. Nevertheless, the second RF harmonic (2ω _H) itself is too weak to be observed. These components are related to the RF amplitude modulated by the LF and reveal the frequency interaction in the discharge process by DF driving [23]. Clearly, the maximum of the power spectrum density (PSD) occurs at the fundamental RF frequency. Nevertheless, PSD at high harmonics is also substantial, as shown in Figure 9C. These results suggest that the excited high harmonics also significantly contribute to the high electron density generation.

By weakening the LF component, the discharge achieves the DF-comparable mode. In this regime, the second RF harmonic structure is quite different from the LF-dominant mode. As shown in Figure 9B, not only the associated sum and beat frequencies, 2ω _H + ω _L and 2ω _H - ω _L, respectively, but also the second RF harmonic, 2ω _H, are all clearly excited. Moreover, the third RF harmonic, 3ω _H, is also strengthened. In the corresponding PSD shown in Figure 9D, one can find that the second RF harmonic and the LF fundamental are comparable, and other high harmonics are also more intensive than those in the LF-dominant case.

The FFT amplitudes of the discharge current as a function of LF and RF voltage amplitudes are presented in Figure 10, including some major excitation spectral lines, such as the fundamentals (ω _L and ω _H), the second RF harmonic (2ω _H), associated beat (2ω _H - ω _L) and sum (2ω _H + ω _L), and the third RF harmonic (3ω _H). As shown in Figure 10A, V _LF varies from 500 V to 1500 V, with a fixed V _RF = 350 V, and one can observe that as the LF component increases, the amplitudes of the LF-related components, namely, ω _L, 2ω _H - ω _L, and 2ω _H + ω _L, also have an obvious increase. However, the RF component and its frequency doubling are substantially suppressed by enhancing the LF component. On the other hand, as shown in Figure 10B, V _RF varies from 0 to 400 V, with a fixed V _LF = 150 V. In the figure, one can see that, for V _RF = 0, the discharge does not occur due to the excessively low input power. As the RF component, i.e., the fundamental, increases, its second harmonic and associated sum and beat components with the fundamental LF also increase, while the fundamental LF component is suppressed substantially. Meanwhile, the third RF harmonic is also suppressed. Combining these results with those in Figure 5, we find that, in the electron density modulation, restraining the high harmonic and increasing the sum and beat components can be an effective mechanism to obtain the desired electron density.