Particle-in-cell simulations of the ionization process in microwave argon microplasmas

In the last few decades, microwave plasmas have played an important role in numerous areas, such as semiconductor fabrication,^1,2 combustion,³ and biomedical treatment,⁴ where breakdown is created intentionally. On the other hand, in electronics, such as transmission lines,⁵ magnetrons,⁶ and antennas,⁷ it is crucial to prevent microwave breakdown. In the late 1940s, MacDonald and Brown performed microwave experiments in air to understand microwave breakdown phenomenon.⁸ Since then, many experimental,^9–11 theoretical,^12–15 and simulation studies^16–18 have investigated the threshold voltage/power for microwave breakdown and interpreted the possible mechanisms. For high-frequency fields, breakdown is typically assumed to occur when the production and loss rates of electrons balance.¹⁹ 

Recently, interest in microwave microplasmas has grown due to the rapid development of micro- and nanotechnology. For DC microscale breakdown, the required strong electric field strips electrons from the cathode by field emission, which increases gas ionization and generates more secondary electrons, resulting in a lower breakdown voltage than predicted by Paschen's law.^20–25 In microwave microgaps, one important parameter is the critical frequency, where the amplitude of oscillating electrons becomes the same as the electrode gap.^26–28 Below the critical frequency, one must incorporate field emission and secondary emission into the microwave breakdown condition in microscale gaps because of the large electron-drift wall loss driven by the electric field.^28,29 However, electron wall loss is minimized above the critical frequency, and breakdown rarely depends on secondary emission.²⁹ Regardless of the difference between the breakdown mechanisms, it is critical to accurately predict the ionization rate and its dependence on the electric field, frequency, and gas density or pressure when assessing breakdown.

Although some studies have measured and modeled the ionization rate, they typically focus on gaps ranging from millimeters to centimeters,^30–33 making the data inapplicable to microscale gaps. While recent studies have investigated microwave breakdown in microscale gaps,²⁹ they rely on the standard empirical DC formula to calculate the ionization coefficient,³⁴ which is neither optimum for microscale gaps³⁵ nor microwave frequencies. Thus, a comprehensive assessment of the ionization rate with gap distance, electric field, and field frequency is required to accurately calculate the breakdown voltages. This paper fills this gap by assessing the ionization rate in microwave microplasmas using one-dimension in space, three-dimensions in velocity (1D-3v) particle-in-cell/Monte Carlo collision (PIC/MCC) simulations. Section II provides background on the theory. Section III summarizes calculations of $ν i$ in argon at atmospheric pressure for gap sizes ranging from 2 to 10 μm and field frequencies ranging from 1 to 1000 GHz using PIC/MCC simulations. We make concluding remarks in Sec. IV.

We point out that Eq. (2) is based on empirically fitting experimental data from previous studies assuming constant electron mobility.^19,39 This approach gives an average collision frequency and only works for a certain electron energy range. The complete equation for collision frequency should be $ν m=N ⟨ σ m v ⟩=p∫ σ m(v) v f(v) d v / k T$⁠, where k is Boltzmann's constant, T is gas temperature, f(v) is electron probability distribution function, and v is electron velocity, making $C=∫ σ m(v) v f(v) d v / v$⁠, which depends on the electron temperature. Other studies found $C=f( E / p)$ by assuming the electron energy follows the Maxwellian distribution, the system is at room temperature, and all the electron-neutral elastic scattering is isotropic.¹³ However, those formulas have more complicated forms, and the variation of C is usually within an order of magnitude for an electron energy from 0.01 to 1 eV (10²–10⁴ K).⁴⁰ Thus, the collision frequency calculated from Eq. (2) serves as a reasonable estimate for further calculation. For more precise values of the collision frequency, one can perform the integration of $σ mv$ based on the electron energy distribution function. The precision obtained here suffices for our current goal of predicting $ν i$ for use in microwave, microscale gas breakdown theories.

Note that while $ν i$ is sinusoidal due to its dependence on E, the electric field oscillates so rapidly since $f> f c r$ that $ν i$ cannot keep up with the electric field oscillation. This makes $ν i$ essentially constant, allowing us to assume that $ν n e t$ is constant with respect to time.

We used the PIC/MCC code XPDP1, which is a bounded electrostatic code for simulating 1D planar plasma devices,⁴⁸ to determine $ν i$⁠. The PIC/MCC technique⁴⁹ is a kinetic approach that is commonly used to simulate plasmas over a wide range of operating conditions and parameters. The 1D computational domain was discretized into computational cells containing computational particles representing larger numbers of either electrons or ions with the ratio of real to computational particles varied to maintain the number of computational particles below 10⁶. The background gas was distributed uniformly in the computational domain and is not treated as a particle. All simulations presented here included three electron-neutral collisions (elastic scattering, excitation, and ionization) and two ion-neutral collisions (elastic scattering and charge exchange), as summarized in Table I. Although excitation collisions were included, the densities of excited species were not tracked since they did not participate in any subsequent collisions. We initialized the simulation with a uniform number density of ions and electrons with the computational particles randomly distributed in the domain. At each time step, we used a leap-frog algorithm to update the computational particle locations and velocities, followed by collisions with the background gas. Poisson's equation was then solved to update the electric potential and electric field distributions in the domain. Finally, new particles that are injected due to boundary processes are introduced in the computational domain.

For the simulations reported here, the electrode at $x=0$ was connected to the AC voltage source $V= V 0 sin( ω t)$⁠, and the electrode at $x=d$ was connected to a capacitor of 1 F to filter any DC signal, as shown in Fig. 2. We divided the computational domain into 100 cells based on the desired spatial resolution for the pre-breakdown simulations. We used time steps between 10⁻¹⁷ and 10⁻¹⁵ s to satisfy the standard Courant condition for numerical convergence of the PIC/MCC technique.⁵⁰ The electrons were injected into the gap from the electrode at $x=d$ with a given current density $j 0$⁠. Once we obtain the instantaneous value for electron number, we fit the results and calculate the ionization rate based on Eq. (12). We used argon gas at atmospheric pressure and consider $2≤d≤10μ m$ with $1≤f≤1000 GHz$⁠.

Previous studies calculated the ionization coefficient for DC fields by using the ratio of electron to ion current density at steady state.^22,24,35 However, this formula becomes invalid here because the current always oscillates with the electric field at microwave frequencies. Thus, we report below a new method to evaluate the ionization process in microwave fields.

Figures 3(a) and 3(b) show the variation of electron number with time in 8 μm gaps for $f< f c r$ and $f> f c r$⁠, respectively, where $V 0=100 V$ and $f c r=15.3 GHz$⁠. For $f< f c r$⁠, the electron number oscillates with the electric field, while for $f> f c r$⁠, the electron number exhibits predominantly exponential growth (although relatively minor oscillations remain) with time. Figures 4(a) and 4(b) show phase-space plots of velocity across the gap as a function of position in the gap to demonstrate the frequency dependence of this behavior for $f< f c r$ and $f> f c r$⁠, respectively. Figure 4(a) demonstrates the collection of the electrons by the electrode at a certain time [when $V ≈ − V 0 , t ≈ ( 4 n − 1 ) π / ( 2 ω )$] during each cycle when $f< f c r$⁠, whereas Fig. 4(b) indicates that most of the electrons always remain confined within the gap when $f> f c r$⁠. Since f plays a key role in the loss mechanism, we divide our simulation results into two regimes based on $f c r$⁠.

For $f< f c r$⁠, the electrodes collect the electrons during each cycle and the electron number oscillates with the electric field, as shown in Fig. 3(a). Next, one can obtain $ν n e t , 0 / ω$ by fitting the simulation results to Eq. (22), as mentioned in Sec. II. Figure 5 shows the variation in electron number with time for a gap with $d=8μ m$ and $f=10 GHz$ for various applied voltages compared to a curve fitted to Eq. (22). Although the fits represent the main characteristics of simulations, smaller peaks that deviate from the fitted curves appear for higher $V 0$⁠. This difference may be attributed to numerical errors or limitations of the simplified model used in the study (such as potentially not fully accounting for changes in gap capacitance). This could also arise due to beating caused by the mismatch between the frequencies of $ν i$ and E. These hypotheses could be further analyzed in future studies. For now, this simple model reproduces much of the dominant behavior.

Based on the fitted parameters, we subsequently calculated the ionization frequency. Figure 6 shows that $ν n e t , 0 / ω$ does not change much with f as a function of the reduced electric field $E e f f / p$ for 2–10 μm gaps. For larger gaps, the curves do not exactly overlap. This may occur because the critical frequency is lower for larger gaps, causing the field frequency (1–10 GHz) to fall in the transition from the boundary-controlled region to the diffusion-controlled region. Also, as we mentioned earlier, the scattering of the data could also result from limitations of the simplified model used in this study or slight numerical errors that are amplified in the simplified model. Figure 7 shows that $ν n e t , 0$ does not exhibit universal behavior with $E e f f / p$ since it varies monotonically with f over this same range of gap distances. However, for all cases, the ionization rate plateaus or even decreases with increasing $E e f f / p$ due to the change in the ionization cross section or the appearance of runaway electrons at high fields,⁵¹ which resembles the results in DC breakdown.³⁹ 

The results suggest that $ν i / p$ is roughly proportional to $E / p$⁠. When the electron energy reaches steady state, the energy gained from the electric field between each collision balances the energy loss due to collisions. This steady state energy $ε$ is proportional to the electric field and also the mean free path, which is inversely proportional to the pressure or gas density; thus, $ε∝ E / p$⁠. Since the energy loss is proportional to $ε$⁠, assuming most of the energy loss produces ionization, then the number of ionization event per molecule is proportional to $E / p$⁠. Strictly speaking, by assuming a Maxwellian distribution, the ionization coefficient may be calculated analytically as $ν i=[ c 1 + c 2 ( E / p ) c 3] exp[ c 4 / ( E / p )]$⁠, where $c 1$⁠, $c 2$⁠, $c 3$⁠, and $c 4$ are gas-dependent constants.¹⁵ Here, we used a simpler format similar to previous studies to fit the results, which captures the main characteristic and is handy to use. Equation (17) is a general form to fit the ionization coefficient for different gas species in macroscale gaps, while Eq. (25) is specific for the microwave discharge of argon in microscale gaps based on the simulations in this work. Figure 9(b) compares our results with argon data reported by others^30,33 and also shows that the ionization frequency of argon exceeds that of air,¹³ primarily due to argon's larger ionization cross section. For both cases at either $f< f c r$ or $f> f c r$⁠, the results validate the similarity law, which is usually used to map how physical laws remain the same at systems of different discharge conditions.^29,52

In summary, this study uses PIC/MCC simulations to determine the ionization rate in argon for microscale gaps at microwave frequencies at atmospheric pressure. We have developed a method to evaluate the ionization rate in AC fields that differs from the approach used at DC. The ionization rate differed for frequencies below and above the critical frequency $f c r$⁠. For $f< f c r$⁠, the electrodes collect the electrons during each cycle, causing the electron number to oscillate with the AC field and $ν i / ω$ to scale approximately with $E e f f / p$⁠. For $f> f c r$⁠, the electrons remain confined within the gap, causing the electron number to grow exponentially and $v i / p$ to become a function of $E e f f / p$⁠. In the future, we will investigate the ionization rate for more gas species at different pressures in microscale gaps. This more accurate representation of ionization could then be incorporated into microwave microscale breakdown theories^14,53 to better guide device design and plasma formation.

This material was based upon the work performed at Purdue University that was supported by the Office of Naval Research under Grant No. N00014-21-1-2441. We also thank Dr. Abbas Semnani for helpful discussions.

The authors have no conflicts to disclose.

Haoxuan Wang: Formal analysis (lead); Investigation (lead); Methodology (equal); Writing – original draft (lead). Ayyaswamy Venkattraman: Formal analysis (supporting); Investigation (supporting); Methodology (equal); Writing – review & editing (equal). Amanda M. Loveless: Investigation (supporting); Supervision (supporting); Writing – review & editing (equal). Allen L. Garner: Conceptualization (lead); Funding acquisition (lead); Project administration (lead); Supervision (lead); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.