A multi-term, multi-harmonic Boltzmann equation model for kinetic behavior in intense microwave and terahertz excited low temperature plasmas Free

The goal of this study is the development and utilization of an efficient, non-stationary description of electron kinetic behavior in low temperature plasmas (LTPs) under the influence of an AC electric field from microwave to THz frequencies. For decades now, significant effort has been focused towards achieving this via Boltzmann equation (BE) based models (see Refs. 1 and 2 for a more complete historical review). Margenau and Hartmann proposed the first formal decomposition of the BE into a Fourier series in time, and a Legendre series about the anisotropy in velocity-space.³ In the high frequency limit, researchers have regularly employed the effective field approximation (EFA),^4–7 which is a special case of the Fourier-Legendre expansion, where only two Fourier terms are considered. Significant efforts have been dedicated to exploring and extending the range of applicability of the EFA and similar low frequency approximations.^8–12 However, all of the aforementioned studies considered only two terms in their Legendre expansion of the BE (i.e., the two-term approximation), which is now known to be inadequate in many cases.^13–16 Makabe and Goto extended the implementation of the Fourier-Legendre approach by considering additional terms, yet still limited to a low order truncations of both the Fourier and Legendre expansion.^17,18 This approach was employed to study a wide range of frequencies,^17–19 however, low order truncation of the Fourier-Legendre series failed to yield the desired accuracy. Owing to their comparative complexity and marginal accuracy, BE models based on the Fourier-Legendre expansion technique were abandoned following the emergence of strictly time-dependent multi-term (MT, more than two Legendre terms) BE models.^20–27 Since these initial demonstrations, time-dependent MT-BE models continue to be developed and applied.^28–34 And in general, multi-term BE models remain a standard tool in the plasma science community.^35–38 

On the other hand, modeling microwave excited LTPs and particularly microwave discharge thresholds and delays, has involved a number of strategies. In early efforts, scaling laws and empirical models were developed (e.g., Refs. 39 and 40). While BE models have been used for this purpose,^4–7,41,42 in recent years, such efforts were largely abandoned for particle based Monte Carlo collision (MCC) models,^43–46 fluid models,^47–50 and global models.^51–55 Moreover, experimental research of high power microwave breakdown in both gas and vacuum continues to be an active field of study.^50,56–61

With consideration for higher frequencies, experimental studies have exhibited a great deal of novel and unexpected phenomena. For example, the first high power microwave breakdown experiments at 110 GHz by Hidaka et al.⁶² have been the focus of extensive study owing to the unique observation of pattern and filament formation. This observation has since motivated a number of subsequent experiments^63–66 and theoretical studies^67–78 yielding unique and valuable information regarding the basic mechanics and properties of LTPs excited at this uncharacteristically high frequency. In general, there is a growing interest within the community regarding LTPs excited at higher frequencies for many applications, including reconfigurable photonic devices,^79,80 detection of concealed nuclear materials,^81–83 and even rocket propulsion.^84–87 This fact, combined with the general trend of increasing power and frequency capabilities of the high power microwave community,^88–90 necessitates improved understanding and prediction capabilities for microwave and THz plasma interactions.

Hence, highly accurate time-dependent kinetic models of electrons in a LTP exist, namely time-dependent MT-BE models and MCC models. BE models based on the Fourier-Legendre approach were also proposed, but in practice, a formal multi-harmonic (MH, more than two Fourier harmonics) implementation has never been demonstrated. The absolute capability of a MTMH-BE model is a priori unknown. However, a key advantage of a MTMH-BE model is that rather than numerically calculating electron kinetic behavior over an electric field cycle, the MTMH-BE approach seeks to reconstruct the time-periodic nature of the electron swarm via Fourier harmonics. As a result, the MTMH-BE model is not restricted to numerically resolving the electric field cycle. In fact, a MTMH-BE model may feature a numerical time step far exceeding the electric field period. This becomes increasingly advantageous at THz frequencies, where the electric field is of the order of femtoseconds, while characteristic timescales within the plasma may be microseconds or longer. Thus, with the ongoing need to model LTPs under the influence of AC electric fields, which may include THz frequencies, a reconsideration of the Fourier-Legendre approach is merited.

In this study, the framework of a MTMH-BE model, based on the Fourier-Legendre approach is developed and validated via model gases. The MTMH-BE model is utilized as a basic kinetic tool to explore electron kinetic behavior in atmospheric pressure air, and specifically explore the limits of the EFA. In the second part of this study, the MTMH-BE model is utilized to predict the temporal evolution of electron kinetic behavior from a cold initial condition, as would be relevant to the application of this model within a LTP global model. To highlight the advantage of this work, the MTMH-BE model is utilized for the prediction of the formative delay time of 95 GHz high power microwave breakdown. In this demonstration, the numerical time step of the MTMH-BE model is approximately six orders of magnitude larger than is possible with a conventional, strictly time-dependent kinetic model.

The focus of this work is the development and utilization of a MTMH-BE model of a LTP excited by a microwave or THz electric field. The electron phase space is expanded into a Fourier series in time. However, the velocity distribution functions (VDFs) retain a time dependence, meaning that while the phase space is assumed to be periodic with the electric field, it is also permitted to vary arbitrarily in time. This enables the calculation of kinetic behavior from an arbitrary initial condition, to a steady-state periodic solution without the need to resolve the electric field cycle. In this section, the basic mathematical framework of the MTMH-BE model is presented. Secondly, the proposed MTMH-BE model is validated using standard model gas benchmarks.

To begin, consider the BE for electrons in a neutral background

Here, $F(r→,v→,t)$ is the phase space of the electrons, for the electron position, velocity, and time, $r→, v→$⁠, and t, respectively. The collision operator, which accounts for the effect of electron-neutral interactions on the electron phase space, is given as C[F] (given in Ref. 38). For this study, the electrons of mass, m, and charge, e, are accelerated (i.e., $a→$⁠) by a z-directed sinusoidal electric field, as $E(t)=E cos (ωt)$⁠, where E is the magnitude of the field and ω is the angular frequency of the electric field. Second, only spatially uniform conditions are considered. The BE may then be rewritten as follows:³ 

Here, the cosine of the anisotropic angle is given as $x=cos θ=vz/v$⁠, where v_z and v are the z-velocity and absolute velocity of the electrons, respectively. As in previous studies,³ the phase space may be expanded into a Fourier series in time, and a series of isotropic and anisotropic VDFs. The electron density, n_e(t), is assumed to be spatially uniform and separable from the VDFs

Zero-order Legendre polynomials about the cosine of the anisotropic angle are given as P_ℓ(x). In previous studies, the number of considered Legendre and Fourier terms (N_ℓ and N_k, respectively) was limited to low order truncations.^17–19 Here, an arbitrary number of Legendre and Fourier terms may be used.

Velocity distribution functions of the even time-harmonics of the ℓ-th Legendre term and k-th Fourier harmonic are given as f_ℓ,k, whereas VDFs of the odd time-harmonics are given as g_ℓ,k. Each VDF retains an arbitrary temporal dependence, while the time-harmonic content is accounted for via the Fourier expansion. In order to account for the effect on non-conservative collisions on the VDFs, the time-harmonic ionization frequency, ν_i(t) = (∂n_e/∂t)/n_e, must be considered. This effect appears in the time-derivative of the electron phase space given below, normalized by the electron density (i.e., via applying a time-derivative to Eq. (3) and dividing by the electron density)

The ionization frequency is defined as the integral of the isotropic portion of the VDF (i.e., VDFs with ℓ = 0) with the reactive part of the collision operator. The reactive cross-section can be defined as σ_R = σ_i − σ_a, where σ_i and σ_a are the ionization and attachment cross-sections, respectively. The isotropic VDF is harmonic in time, meaning that the ionization frequency is also harmonic in time, and may be written as follows:

The ionization frequency for the even and odd time-harmonic VDFs is as follows:

The typical procedure for exploiting the orthogonality of the Legendre and Fourier content may be used to yield a coupled system of ordinary differential equations for both the even and odd time-harmonic VDFs. Note that the impact of the time-dependence of the VDFs, f_ℓ,k(v, t) and g_ℓ,k(v, t), on the Fourier orthogonality was neglected

The Kronecker delta is given as δ_p,q, where δ_p,q = 1 for p = q and δ_p,q = 0 for p ≠ q. The Ξ coefficients are conditional parameters which describe coupling between the time harmonic content through the harmonic ionization frequency, given by the following integrals:

The time dependence of the VDFs is discretized in a conventional implicit finite difference scheme, $∂fℓ,k(v,t)/∂t=(fℓ,k(v,t+Δt)−fℓ,k(v,t))/Δt$⁠, where Δt is the numerical time step. The same approach is taken for g_ℓ,k(v, t). As in prior work,³⁸ Eqs. (8) and (9) are solved numerically using a finite difference scheme with respect to energy. The maximum considered energy in the energy grid was selected such that the magnitude of the distribution at the peak energy was negligible (i.e., below ∼10⁻¹⁰). Typically, 1000 energy grid points were utilized, as the use of additional energy grid points yielded no observable difference in the solutions.

The MTMH-BE model utilized N_ℓ = 4 Legendre terms (i.e., ℓ = 0, 1, 2, 3). As in previous studies,³⁸ additional Legendre terms affected the predicted quantities by less than one percent. Physically relevant VDFs have an even sum of their Legendre and Fourier indices (i.e., ℓ + k = even).^3,17,18 A total of four Fourier harmonics were considered for each Legendre term. For instance, for (ℓ, k), for ℓ = 0, the (0, 0), (0, 2), (0, 4), and (0, 6) terms were considered. Similarly, for ℓ = 1, the (1, 1), (1, 3), (1, 5), and (1, 7) terms were considered. This is equivalent to N_ℓ = 4 and N_k = 8.

This combined harmonic temporal approach circumvents the restriction that the numerical time step must finely resolve the electric field period. In this study, it was observed that the time step needed to be sufficiently short to resolve the temporal development of the VDFs. Once steady-state periodic conditions were achieved, no restriction to the time step was observed.

In Sec. II B, a benchmark of the MTMH-BE model is performed using the conservative Reid ramp model gas⁹¹ and the non-conservative ionization model gas of Lucas and Salee.⁹² To check calculations in air, the independently developed MCC model, METHES, is used,⁹³ along with a second, strictly time-dependent MT-BE model based on the framework of MultiBolt.³⁸ In all cases, the time-dependent MT-BE model utilized N_ℓ = 4 Legendre terms. Although not shown here, the strictly time-dependent MT-BE model was also verified using AC benchmark data for the Reid ramp and Lucas-Salee model gases. The MTMH-BE model calculations are consistent with each of these models.

To validate the MTMH-BE model, benchmark calculations are made using the conservative Reid ramp model gas,⁹¹ and the non-conservative ionization model gas of Lucas and Salee.⁹² The Reid ramp model gas is defined by a constant momentum transfer cross-section, σ_m = 6 × 10⁻²⁰ m², and a linearly increasing excitation cross-section, σ_ex = 10(ε–0.2)10⁻²⁰ m² with a 0.2 eV threshold energy (ε is the electron energy, in eV). The neutral mass is 4 AMU, and neutrals are assumed to be stationary (i.e., gas temperature, T = 0). For all calculations with the Reid ramp model gas, the reduced electric field is $E/N=102 cos (ωt) Td$⁠. (N is the neutral density and 1 Td = 10⁻²¹ V × m²).

The MTMH-BE calculated electron mean energy and drift velocity are given in Fig. 1, along with the strictly time-dependent MT-BE results of White et al.,²² for a wide range of reduced angular frequencies, ω/N. Note that in this study, the electron phase space is not expanded about gradients in the electron density. As a result, transport coefficients associated with density gradients (e.g., the diffusion tensor, bulk drift velocity, etc.) cannot be calculated (see Refs. 94 and 95). It is emphasized that the goal of this work is kinetic modeling of microwave and THz discharges, and not swarm and transport theory; hence these coefficients are not needed.

The MTMH-BE is in reasonable agreement with the results of White et al.,²² with outstanding agreement at the higher reduced angular frequencies. At lower frequency, the periodic behavior becomes strongly non-sinusoidal. As one might expect, representing this phenomenon using a truncated Fourier series yields some limitation. Still, the MTMH-BE model reasonably predicts swarm behavior within the desired accuracy. On the other hand, at higher reduced angular frequency, the solution is notably more sinusoidal and very well reproduced by the MTMH-BE model.

Next, the non-conservative ionization model gas of Lucas and Salee⁹² is employed. In this model, the momentum transfer cross-section is given as σ_m = 4ε^−1∕2 10⁻²⁰ m², while the excitation and ionization cross-sections are defined as σ_ex = 0.1(1 − F)(ε–15.6) 10⁻²⁰ m² and σ_i = 0.1F(ε–15.6) 10⁻²⁰ m², respectively. The excitation and ionization cross-sections have a threshold energy of 15.6 eV. F is a parameter of the model gas that defines the fraction of the inelastic cross-section that is made up by the ionization cross-section. For F = 1, the ionization cross-section accounts for the entire inelastic cross-section. For F = 0, the ionization cross-section is zero, and the inelastic cross-section is entirely defined by the excitation cross-section. The neutral mass is taken to be 1000m. Again, neutrals are assumed to be at rest (T = 0). The reduced electric field for the Lucas-Salee calculations is $E/N=10 cos (ωt) Td$⁠.

The steady-state periodic mean electron energy, drift velocity, and ionization rate (defined as k_i = ν_i/N), are given in Fig. 2. Similar to the Reid ramp model gas, the MTMH-BE model is able to predict desired swarm coefficients fairly well, with excellent agreement at higher frequency. Interestingly, at lower frequency, Gibbs phenomenon-type artifacts in the solution are apparent, with the MTMH-BE model oscillating about the target solution [e.g., see mean energy in Fig. 2 (left)]. Additional Fourier terms improve the prediction capability of the MTMH-BE model, however, Fourier artifacts remain.

Overall, the model gases serve to highlight the general performance of the MTMH-BE model and validate its accuracy. The MTMH-BE model is observed to very accurately predict electron swarm behavior, especially at higher frequency where solutions become increasingly sinusoidal. At very low frequency, solutions become highly non-sinusoidal and the MTMH-BE model features some limitation, but is still able to reasonably predict electron behavior.

In this section, a basic kinetic phenomenon in microwave and THz excited LTPs in air is explored using the MTMH-BE model. The gas is composed of N₂:O₂ in a 0.78:0.22 fractional composition. The cross-sections of Biagi,^96,97 which were developed via MCC calculations, are utilized.

Similar to the model gases, the MTMH-BE model is able to accurately predict basic plasma properties, such as the mean electron energy and ionization rate over a wide range of frequencies, although, some limitation is observed when solutions become highly non-sinusoidal. For atmospheric pressure air, this behavior is observed for microwave frequencies in the range of f = 1 GHz (ω = 2πf), and relatively high reduced electric field, see Fig. 3.

For f = 1 GHz and high reduced electric field, the sharp decay in mean electron energy is not entirely reproduced by the MTMH-BE model, although the general behavior is well-represented. At this same frequency and lower reduced electric field, the MTMH-BE better reproduces the results of the strictly time-dependent MT-BE model. This behavior is still observable for f = 10 GHz, although the deviation is less pronounced. Examples of the MTMH-BE model predictions at higher frequency are given in Sec. III B.

The strong modulation in macroscopic properties of the LTP over a microwave cycle suggests that the EFA is not well justified. Recently, Kourtzanidis and Raja have shown that the assumption of a stationary distribution over a microwave cycle is satisfied only in a limited regime for atmospheric pressure air.⁹⁸ Still, many researchers have observed that even when kinetic behavior is highly non-stationary over a microwave cycle, time-averaged quantities are only weakly influenced by the modulation when the microwave frequency is on the order of the energy relaxation frequency or higher.^8,19,99 This may be mathematically attributed to the weak coupling between the isotropic VDF and higher order harmonic VDFs.⁹ Using the MTMH-BE model, quantitative calculations may be made regarding the absolute accuracy of the EFA and associated time-averaged properties of the swarm.

Since the EFA is simply a truncation of the Fourier expansion of the BE to only two time-harmonics, benchmarking the EFA is equivalent to benchmarking the Fourier expansion to only two terms. Thus, exploring the limits of the EFA is analogous to early multi-term BE model investigations, where the error in low-order truncation of the Legendre (or spherical harmonic) expansion of the BE was first quantified.^100,101

As a figure merit, consider the modulation depth of the mean electron energy, defined as $1−ε¯peak/ε¯avg$⁠, where, $ε¯peak$⁠, is the peak of the mean electron energy over a microwave cycle, and, $ε¯avg$⁠, is the mean electron energy, time-averaged over a microwave cycle. The electron energy modulation depth for LTPs in air at 760 Torr is given in Fig. 4, for reduced electric fields up to 1000 Td and electric field frequencies from 1 GHz to 10 THz. Similar to the conclusions of Kourtzanidis and Raja,⁹⁸ the assumption of a stationary mean electron energy is poorly justified at almost any reduced electric field level for frequencies below 100 GHz. At 100 GHz, the error of the modulation depth exceeds 10%, even for modest reduced electric field levels below 200 Td. Also at 100 GHz, for reduced electric fields exceeding ∼600 Td, the modulation depth exceeds 30%. The modulation depth is below 10% for electric field frequencies of 0.8 THz and higher and a reduced electric field up to 1000 Td.

The time-averaged ionization rate versus reduced electric field for various frequencies for atmospheric pressure air is given in Fig. 5. Note that the MTMH-BE model with N_ℓ = 4 and N_k = 2, is a multi-term BE model utilizing the EFA. The most complete models, that is the MTMH-BE model with N_k = 8 and the strictly time-dependent MT-BE model, are in very good agreement. On the other hand, both BOLSIG+¹⁰² and the MTMH-BE model utilizing the EFA (N_k = 2) are in a very good agreement. This suggests that the error associated with both BOLSIG+ is not a consequence of the two-term approximation, but rather a consequence of the EFA. The error of the EFA is apparent at 1 GHz, where models based on the EFA may under-predict the time-averaged ionization rate by more than a factor of two. Even though the electron swarm demonstrates notably stronger modulation at higher reduced electric fields, the error associated with the EFA decreases with increasing reduced electric field.

To summarize this subsection, the general kinetic behavior of microwave and THz excited LTPs has been explored for atmospheric pressure air using the MTMH-BE model. These calculations were also verified using a strictly time-dependent MT-BE model. Over a wide-range of frequencies, well into the THz regime, the assumption stationary kinetic behavior over a microwave cycle is not well justified. Still, the non-stationary nature of the electron swarm has overall little impact on the time-averaged quantities, as has been demonstrated in earlier studies. However, at the low frequency limit, rates associated with the EFA may be incorrect by more than a factor of two. Hence, utilizing EFA based rates in a fluid model, or a similar model, in this range, is not recommended. In such instances, time-averaged rates from a MTMH-BE model, or a comparable kinetic model should be utilized.

As previously discussed, a key application of the MTMH-BE model is its use for the efficient kinetic modeling of microwave and THz excited LTPs. In this section, the ability of the MTMH-BE model to predict electron heating in air is explored. The initial condition for these calculations is equivalent to a Maxwellian distribution with an electron temperature of zero (T_e = 0). Utilizing a finite T_e (i.e., room temperature) has a negligible effect on the calculation results.

The large vibrational cross-sections of molecular nitrogen and oxygen result in electrons very rapidly achieving equilibrium. Unipolar discharges in air at atmospheric pressure, achieve equilibrium of the electrons within ∼10 ps,^103,104 which is on the order of an electric field cycle for microwave and THz frequencies. Hence, the temporal development of an electron swarm responding to an instantaneous microwave or THz electric field is a novel regime for investigation.

The development of the electron mean energy at 760 Torr for moderate-high RMS reduced electric field levels is given in Fig. 6. Figure 6 also shows MCC calculations using METHES,⁹³ and time-dependent MT-BE model calculations. Excluding obvious statistical artifacts in the MCC calculations, the MTMH-BE model demonstrates excellent agreement with both the time-dependent MT-BE model, and MCC calculated mean electron energy. Consistent with the findings reported in Sec. III A, even at a microwave frequency of 110 GHz, strong modulation of the mean electron energy is observed in atmospheric pressure air.

The temporal development of the electron mean energy for a 110 GHz microwave frequency and 100 Torr pressure is given in Fig. 7. Here, electron heating occurs over a longer, tens of picoseconds timescale, still on the order of the microwave period. In general, the MTMH-BE model agrees reasonably well with the strictly temporal BE model, although some deviation is observed.

Well into the THz regime, the electron heating profiles for a 0.67 THz electric field and atmospheric pressure are given in Fig. 8. Again, at this higher pressure, electron heating occurs on a picosecond timescale, on the order of the electric field period. Overall, the MTMH-BE model is again in reasonable agreement with the temporal BE model. Again as reported in Sec. III A, even for a 0.67 THz electric field, strong modulation of the mean electron energy is still observed in atmospheric pressure air.

Finally, electron heating profiles for a 0.67 THz electric field and 100 Torr pressure are shown in Fig. 9. Though not shown, the MTMH-BE model features roughly the same degree of agreement as before during the electron heating process, with very good agreement in the periodic steady-state solution. Again, for a 0.67 THz electric field, and a reduced air pressure of 100 Torr, some modulation of the swarm over an electric field cycle is observed.

Finally, to relate this work to practical experimental application, theoretical formative delay times, τ_f, are compared with the recent W-band (95 GHz) experimental breakdown data of Kim et al.,¹⁰⁵ see Fig. 10. Note that these experiments were conducted in argon. Again, the cross-sections of Biagi are utilized.^96,97 An initial electron density of n_e = 10³ cm⁻³, is assumed and the temporal development of the electron density is numerically calculated from the time-dependent ionization frequency, ν_i(t). The formative delay time corresponds to the time necessary for an initial electron density to avalanche into a fully developed plasma which strongly reflects the incident microwave (i.e., when the plasma frequency exceeds the microwave frequency).

The MTMH-BE model predicted formative delay times are shown in Fig. 10, along with the experimental data of Kim et al.¹⁰⁵ Overall, the MTMH-BE model demonstrates fair agreement with the experimental data, within the apparent measurement accuracy. Note that for an electric field frequency of 95 GHz, the electric field period is approximately 10 ps. Though not shown, the strictly time-dependent MT-BE model utilized in this study required Δt ∼ 20 fs to achieve convergence (approximately 500 samples per field cycle). Again, by making use of the a priori known periodic nature of the electron kinetic behavior, the MTMH-BE model does not feature such restrictions on the time step. The MTMH-BE calculations in Fig. 10 utilized a numerical time step of 20 ns, a factor of 10⁶ larger than a strictly time-dependent kinetic model.

In total, this study details the first demonstration and validation of a formal, multi-term, multi-harmonic BE model based on a Fourier-Legendre expansion of the BE. In the first portion of this study, the MTMH-BE model was utilized to elucidate the fundamental kinetic behavior in microwave and THz excited LTPs in atmospheric pressure air. It was demonstrated that the assumption of stationary kinetic behavior over a microwave cycle (i.e., the EFA) is not well justified in a wide parameter range. However, consistent with earlier studies, the non-stationary behavior over a microwave cycle was found to have limited influence over the time-averaged properties of the swarm. Still, in atmospheric pressure air, for a microwave frequency of 1 GHz, the EFA was observed to under-predict the time-averaged ionization rate by more than a factor of two. However, these deviations diminished at moderate to high reduced electric fields. Overall, the MTMH-BE model was demonstrated to very accurately predict time-averaged swarm properties, with notable improvement over conventional EFA based models. Consequently, the use of a MTMH-BE model for the accurate prediction of rate and transport coefficients for fluid models should be considered in ranges where EFA models are known to fail.

A second focus of this study was the demonstration of the MTMH-BE model as a fundamental kinetic modeling tool for microwave and THz excited LTPs. The MTMH-BE model was demonstrated to accurately predict the temporal evolution of kinetic behavior in microwave and THz induced electron heating. Further, the MTMH-BE model was demonstrated to accurately and efficiently predict high power microwave induced plasma formation delays. At the present time, global models continue to be a fundamental tool in the plasma science community.¹⁰⁶ With consideration for ongoing efforts to develop kinetic global models,^107,108 the integration of the MTMH-BE model (or even the strictly temporal BE model) within such a global model offers the opportunity to bring true kinetic fidelity to current global modeling efforts.

The author thanks Ron White for helpful discussion regarding his benchmark data, Mohamed Rabie for his help with METHES, and Dongsung Kim and EunMi Choi for the useful discussion regarding their experimental data. The author also thanks Guy Rosenzweig and Michael Shapiro for their careful review of this manuscript.

This work was supported by the Air Force Office of Scientific Research Program on Plasma and ElectroEnergetic Physics under Grant No. FA9550-15-1-0058.