A model of an electron-beam–plasma system is introduced to model the electrical breakdown physics of low-pressure nitrogen irradiated by an intense pulsed electron beam. The rapidly rising beam current induces an electric field, which drives a return current in the plasma. The *rigid-beam model* is a reduction of the problem geometry to cylindrical coordinates and simplifications to Maxwell's equations that are driven by a prescribed electron beam current density. The model is convenient for comparing various reductions of the plasma dynamics and plasma chemistry while maintaining a good approximation to the overall magnitude of the beam-created electric field. The usefulness of this model is demonstrated by coupling the rigid-beam model to a fluid plasma model and a simplified nitrogen plasma chemistry. The dynamics of this coupled system are computed for a range of background gas pressures, and the results are compared with experimental measurements. At pressures 1 Torr and above, the simulated line-integrated electron densities are within a factor of two of measurements and show the same trend with pressure as observed in experiment.

## I. INTRODUCTION

The interaction of an intense electron beam (e-beam) with partially or fully ionized gas has many applications, including plasma processing,^{1–3} interactions in the ionosphere (e.g., solar wind and aurora),^{4,5} and electromagnetic pulse events.^{6–8} In many of these natural and laboratory plasmas, air and its properties play an important role. The rapid change of the air from an insulator to an electrically conductive plasma in the presence of an intense e-beam is a complex phenomenon, which depends on the plasma chemistry and rapidly changing electromagnetic fields. In this paper, a set of approximations to the beam and electromagnetic field equations, called the rigid-beam (RB) model, are introduced. This model provides a convenient framework for examining the plasma response as system parameters are varied or when different plasma approximations are made.

At early times during an e-beam pulse, the beam's space-charge creates an electrostatic field in the gas that can cause the beam to expand radially at the front. Beam-impact ionization then produces sufficient plasma density to effectively neutralize the beam's space charge so that it can be confined radially and propagate axially.^{9} Once charge neutralization occurs, the electric field is dominated by an inductive component driven by the rapidly changing e-beam current. Electrical breakdown in this field, together with beam-impact ionization, rapidly increases the plasma density, and a plasma return current begins to flow in the increasingly conductive plasma. This plasma return current can be comparable to the beam current and contributes to the duration and magnitude of the inductive electric field inside the plasma and the electromagnetic field radiating outside of the beam–plasma system. The processes that cause the air to evolve from a low conductivity gas to a highly conductive plasma are quite complex and involve e-beam dynamics, physical and chemical changes in the gas, thermal plasma dynamics, and collisional interactions between the plasma and the gas. The beam and plasma are coupled through Maxwell's equations. To help better understand the plasma physics, it is useful to develop reduced models of the beam–plasma system that capture the essential physics. These reduced models offer relatively fast computation times, which allows us to rapidly develop plasma chemistry models and to determine pressures where a kinetic plasma model is needed and where a fluid model can suffice.

The rigid-beam (RB) model described here is a reduction of the problem geometry to cylindrical coordinates with the e-beam treated as a specified source that drives the remaining system of equations, i.e., beam dynamics are not coupled into the governing equations. The model allows the e-beam to have arbitrary time variations of current and voltage and an arbitrary current density profile. In this paper, we use spatial and time variations guided by experiments. The electron beam used in experiments has a peak current of about 4 kA, a peak voltage of about 100 kV, a pulse width of about 50 ns, and a diameter of about 4 cm.^{10} The model described in this paper is an extension of previous models developed to study the interaction of an intense e-beam with nitrogen gas.^{11–13} The e-beam is assumed to be injected parallel to the axis of a perfectly conducting cylinder of radius *R _{w}* filled with low-pressure nitrogen in the 1 to 10 Torr range. The usefulness of an earlier RB model was demonstrated by comparing two models for electron–gas collisions: one with collision rates computed from the reduced electric field (

*E*/

*p*) and a second where they were determined from the mean plasma electron energy.

^{14}The plasma chemistry model used in this paper adds dissociative recombination and beam-impact ionization, both of which were neglected in that earlier RB model. In addition to this expanded plasma chemistry, the model presented in this paper extends the zero-dimensional RB model of Ref. 13 to a one-dimensional model by allowing an arbitrary radial profile for the electron beam. The RB model employed in this paper is described more fully in Sec. II.

The complete set of collisions that affect the plasma and gas particles is referred to as the plasma chemistry. Here, the plasma chemistry will be limited to that of a weakly ionized plasma. “Weakly ionized” implies that the collision rate between electrons and molecular nitrogen in the $N2(X1\Sigma )$ ground state remains much larger than all other collisional rates. The various cross sections and rates required to model weakly ionized N_{2} discharges are readily available through tables in the literature,^{15} the LXCat database,^{16} and the Quantemol database QDB.^{17} The chemistry model is described in more detail in Sec. III. Results from fluid simulations using the RB model over a range of gas pressures are presented in Sec. IV.

The RB model described in this paper extends the 0D model of Ref. 13 to a 1D model and expands the chemistry used in Ref. 14 to include recombination and beam-impact ionization. These additions allow simulations to be performed, which maintain a higher fidelity with the experiments. The line-integrated electron densities predicted by the simulations and reported in Sec. IV are within a factor of two of the experimental measurements. The density obtained in simulations increases with increasing pressure, which is a trend that has also been observed in experiments. The important conclusions of the paper and suggestions for future work are presented in Sec. V.

## II. THE RIGID-BEAM MODEL

The RB model is a simplification of both the e-beam dynamics and Maxwell's equations that results in a set of equations that are easier and faster to solve but retain the gas breakdown physics driven by a rapidly rising, high-current-density electron-beam. The equations for the electric and magnetic fields, when coupled to dynamic equations for the plasma response and to a plasma chemistry model, provide a self-consistent set of equations for the evolution of the plasma and the fields. This model retains the essence of the gas breakdown problem and is an ideal platform to develop, compare, and contrast different plasma chemistry and response models. In the remainder of this section, the reduction of Maxwell's equations to a simpler set of equations is developed.

A major approximation used to derive the RB model is that the net particle currents are assumed to be large compared to the displacement current so that $\u03f50\u2202E/\u2202t\u226aJnet$, where $Jnet$ is the net current density obtained by summing the primary beam current density $Jb$ and the plasma current density $Jp$. This approximation is good when the plasma density is large compared to the beam density so that the beam–plasma system is nearly charge neutral. Under these conditions, the electrostatic field is negligible and the inductive electric field associated with the rapidly rising beam current is the dominant field. The RB model therefore does not apply to the early time beam–plasma behavior before space-charge neutralization is established or to subnanosecond pulses. It is assumed that this non-neutral phase is very short compared to the rise time of the beam so that the main electric field is the inductive electric field over the majority of the beam pulse. An order-of-magnitude estimate based on beam parameters produced in experiments^{10} indicates that, at 1 Torr pressure, beam-impact ionization will charge-neutralize the beam in about 1 ns.

Ignoring displacement currents in Ampère's law, the equations for the electric and magnetic fields can be written as

where SI units are used. Taking the curl of Faraday's law and substituting Ampère's law to eliminate the magnetic field allows these two equations to be combined into a single Poisson-type equation for the electric field, which is given by

where it has been assumed that $\u2207\xb7E\u22430$. Equation (3) shows that the initial source of the inductive electric field is the rapidly rising, e-beam current density $Jb$. Another important term that determines the overall magnitude and duration of the electric field is the plasma current density $Jp$. During the rising part of the beam pulse, the plasma current density opposes the beam current, which tends to reduce both the net current and the magnitude of the electric field. During the falling part of the beam pulse, the plasma current is in the same direction as the beam current and can enhance the magnitude of the electric field.

It is often assumed that nearly complete current neutralization occurs when modeling the transport of high current-density electron beams in low pressure gas.^{18} In this work, however, no assumptions are made about the magnitude of the induced plasma current $Jp$. Rather, both the plasma current and the induced electric field are modeled using a set of coupled equations.

To simplify Eq. (3) even further, it is assumed that the beam is cylindrically symmetric and flows in the positive z-direction. It is also assumed that spatial gradients in the direction of the beam propagation are small compared to radial gradients. This also implies that the pulse is long enough so that rapid temporal variations do not induce steep axial gradients. With these assumptions, both the beam and plasma current densities can be written as

With these assumptions, Eq. (3) simplifies to

Once the beam and plasma current densities are known (*J _{b}* specified and

*J*computed) and boundary conditions are specified, this equation can be solved for the electric field. The boundary condition at

_{p}*r*= 0 is that $\u2202Ez/\u2202r|r=0=0$ and, assuming a perfectly conducting wall at

*R*, the other boundary condition is $Ez(Rw)=0$. The plasma model described in Sec. III is used to compute the plasma current density required to solve Eq. (6). After solving Eq. (6) for

_{w}*E*, Eq. (1) can be used to calculate $B\theta $.

_{z}## III. PLASMA RESPONSE MODELS

The simplifications described in Sec. II deal with how the e-beam and the electromagnetic fields are treated in the RB model. In addition to these approximations, a plasma response model is needed so that the plasma current can be computed for the source term in the RB model.

As the plasma density grows, the plasma current can become a significant fraction of the beam current. Since the plasma electrons are much more mobile than plasma ions, and the beam timescale is much shorter than the plasma hydrodynamic response time, ions will be treated as an immobile species that provides overall charge neutrality. Under this assumption, the plasma current density will be determined entirely from the motion of the plasma electrons. The plasma current density can then be written as

where *e* is the magnitude of the electron charge, *n _{e}* is the plasma electron density, and

*V*is the electron drift velocity that occurs as a result of acceleration in the electric field and scattering collisions. The electron density and the drift velocity are defined by moments of the electron energy distribution function, which itself depends on the electric field and the details of the plasma chemistry. In the remainder of this section, the details of the plasma chemistry and the fluid model used in the results section are presented.

_{e}### A. Plasma chemistry

The plasma chemistry used in this paper includes any physical changes to the individual molecules, such as excitation or ionization, that result from energetic collisions between electrons and gas particles. This includes elastic scattering as well as inelastic collisions where a fraction of the electron energy is transferred to the molecule to produce rotational, vibrational, and electronic excitations, and ionization. Since the purpose of the present paper is to introduce the rigid-beam model as a way to study electron-beam–produced plasmas, for simplicity it will be assumed that the e-beam is injected into a gas consisting entirely of ground state N_{2} molecules. The role of oxygen and other trace constituents of air (such as water vapor) in these plasmas will be the subject of future work.

A complete plasma chemistry model describes processes that can take place between the plasma electrons and all heavy species, neutral and ionized, that may be present. In addition to collisions between thermal plasma electrons and the gas, collisions that result from the direct impact of the beam electrons and the gas can also be significant. In the remainder of this section, the electron-impact driven interactions that will be used in this paper are described.

### B. The weakly ionized plasma model

The weakly ionized plasma model is the simplest plasma chemistry model for an electrical discharge. In the weakly ionized model, only reactions between electrons and N_{2} gas need to be followed. The possible end products of these reactions are either rotationally excited, vibrationally excited, or electronically excited molecules, or molecular ions. It is assumed that the densities of these excited or ionized species are small compared to the ground state neutral density so that the reactions between electrons and excited species can be neglected. For weakly ionized N_{2}, the complete set of collisions can all be written in the form

where $N2*$ represents an excited or ionized state of nitrogen and *a* is an integer that indicates how many electrons are produced in these reactions. For molecular nitrogen, *a* = 1 for excitation processes (rotational, vibrational, and electronic excitation) and *a* = 2 for ionization. For elastic collisions, where there is no change in internal energy of the molecule, *a* = 1 and $N2*=N2$. The kinetic energy transfer between electrons and N_{2} molecules during an elastic collision is small because of the small electron-to-molecule mass ratio. While elastic collisions are much more frequent than inelastic collisions, the loss of electron energy due to these collisions is less than 7% of the inelastic energy loss for the rate tables used in this paper. When the mean electron energy is greater than 0.5 eV, the loss is even smaller, less than 0.1%. The energy loss due to elastic collisions is therefore neglected in this paper. Inelastic collisions result in quantized energy transfers between the plasma electrons and molecules, increasing the internal energy of the molecules by the amount $\epsilon *$. In addition to the reactions given in Eq. (8), electrons with energies above the 9.75 eV dissociation threshold can drive reactions, which result in the dissociation of the nitrogen molecule into two nitrogen atoms or ions. The atomic products of these dissociations are not tracked in this model.

In the weakly ionized plasma model, it assumed that, in general, the densities of all the excited species that result from collisions with electrons are small so that reactions with these excited products can be neglected. However, for diagnostic purposes and to prepare for future versions of the model where the weakly ionized approximation breaks down, the density of the ground-state neutrals as well as the densities of all excited species is tracked using rate equations. The time-dependent rate equations for the type of reactions described by Eq. (8) are given by

where $\nu N2,j*=nN2kN2,j*$ is the collision frequency for the $jth$ excited state, and $kN2,j*$ is the rate coefficient obtained by averaging the energy-dependent cross sections over the electron energy distribution function. The density of N_{2} at *t* = 0 is set based on the pressure of the background gas, which is one of the main input parameters for this model. The densities of all other excited states at *t* = 0 are set to zero. Beam impact ionization as well as dissociative recombination will be treated in Sec. III C 2 and will result in additional terms being added to Eqs. (9) and (10).

The rate coefficients are defined by

where $\sigma N2,j*(\epsilon )$ is the cross section for the reaction, $v\epsilon =2\epsilon /m$ is the speed of an electron at energy *ε*, and $f0(\epsilon )$ is the electron energy distribution function (eedf). The eedf is normalized so that

The time dynamics of the eedf is described by the Boltzmann equation. However, in this paper, we will use a simplification where the eedf takes the form of a Maxwellian distribution with a time-dependent temperature. Further details will be described in Sec. III C 3.

Each type of collision is represented by a cross section and a threshold energy. The N_{2} collisional cross section data used in this paper are the Phelps data set downloaded from the LXCat website.^{16} The individual cross sections in this database come from a number of sources, which can be found in Ref. 15. In this data set, there are a total of 24 different collision types with each identified by the excited product that results from the collision. The complete set of reactions and thresholds is summarized in Table I. Molecular term notation is used to describe both the ground state ($X1\Sigma $) of the N_{2} molecule and the electronic excited states. There is one momentum transfer collisional cross section and one for the aggregate rotational mode of $X1\Sigma $ labeled as *rot*. There are cross sections for eight vibrational modes of $X1\Sigma $ with vibrational quantum numbers $v=1\u20268$, for thirteen molecular electronic excited states, and there is one total ionization cross section. Even though the ionization reaction includes dissociative ionization events as well as ionization into excited states of $N2+$ such as $A2\Pi $ and $B2\Sigma $, the ionization product state is labeled in this paper by $N2+$. The electronic triplet state $A3\Sigma $ is divided into three lumped vibrational bands as indicated in the table. No vibrational resolution is included for any of the other electronic states.

Collision type . | Products . | Threshold (eV) . | Beam-impact production efficiency, $\xi N2,j*$ . |
---|---|---|---|

Elastic | $X1\Sigma $ | 0 | ⋯ |

Rotational^{a} | rot | 0.02 | 20.25 |

Vibrational^{a} | v = 1, 2, 3, 4 | 0.29, 0.59, 0.88, 1.17 | 5.237, 0.391, 0.147, 0.104 |

v = 5, 6, 7, 8 | 1.47, 1.76, 2.06, 2.35 | 0.100, 0.075, 0.070, 0.049 | |

Electronic | $A3\Sigma $ ($v=0\u22124$), $A3\Sigma $ ($v=5\u22129$) | 6.17, 7.0 | 0.167/3,^{b} 0.167/3^{b} |

$B3\Pi $, $W\u20093\Delta $, $A3\Sigma $ ($v\u226510$), $B\u20323\Sigma $ | 7.35, 7.36, 7.8, 8.16 | 0.148, 0.112, 0.167/3,^{b} 0.030 | |

$a\u20321\Sigma $, a $1\Pi $, w $1\Delta $, C $3\Pi $ | 8.4, 8.55, 8.89, 11.03 | 0.028, 0.081, 0.031, 0.057 | |

E $3\Sigma $, $a\u20331\Sigma $, sum | 11.87, 12.25, 13.0 | 0.001, 0.010, 0.664 | |

Ionization | $N2+$ | 15.6 | 1 |

Dissociative recombination^{c} | Atomic N | 0 | ⋯ |

Collision type . | Products . | Threshold (eV) . | Beam-impact production efficiency, $\xi N2,j*$ . |
---|---|---|---|

Elastic | $X1\Sigma $ | 0 | ⋯ |

Rotational^{a} | rot | 0.02 | 20.25 |

Vibrational^{a} | v = 1, 2, 3, 4 | 0.29, 0.59, 0.88, 1.17 | 5.237, 0.391, 0.147, 0.104 |

v = 5, 6, 7, 8 | 1.47, 1.76, 2.06, 2.35 | 0.100, 0.075, 0.070, 0.049 | |

Electronic | $A3\Sigma $ ($v=0\u22124$), $A3\Sigma $ ($v=5\u22129$) | 6.17, 7.0 | 0.167/3,^{b} 0.167/3^{b} |

$B3\Pi $, $W\u20093\Delta $, $A3\Sigma $ ($v\u226510$), $B\u20323\Sigma $ | 7.35, 7.36, 7.8, 8.16 | 0.148, 0.112, 0.167/3,^{b} 0.030 | |

$a\u20321\Sigma $, a $1\Pi $, w $1\Delta $, C $3\Pi $ | 8.4, 8.55, 8.89, 11.03 | 0.028, 0.081, 0.031, 0.057 | |

E $3\Sigma $, $a\u20331\Sigma $, sum | 11.87, 12.25, 13.0 | 0.001, 0.010, 0.664 | |

Ionization | $N2+$ | 15.6 | 1 |

Dissociative recombination^{c} | Atomic N | 0 | ⋯ |

^{a}

These are rotational and vibrational excitations of the ground state, N_{2}(X $1\Sigma $).

^{b}

The beam-impact excitation of the A$3\Sigma $ state is not vibrationally resolved in Ref. 12. In this model, 1/3 of the excitations are attributed to each of the three vibrational bands.

^{c}

Unlike the other reactions in this table, the target species for this reaction is $N2+$.

In Refs. 15 and 16, the authors include a cross section for excitation to a sum of singlet states (labeled *sum* in Table I). The singlet states have threshold energies between 12.5 and 14.8 eV, which are well above the 9.75 eV dissociation energy for the nitrogen molecule. As a result, the excited singlet states are highly unstable and often result in the dissociation into two nitrogen atoms. In Ref. 19, the fraction of the excited singlet sum states that dissociates is estimated to be 0.73. This dissociation pathway as well as the pathway through high lying triplet states is ignored. Future work will consider these pathways for regimes where this is important.

Collisions between electrons and excited molecular states of nitrogen (momentum transfer, stepwise excitation and stepwise ionization) as well as collisions between electrons and atomic nitrogen are beyond the scope of the weakly ionized plasma assumption. Some of these processes were considered in Ref. 13, and they will also be the subject of future papers. There are many additional processes which affect the excited state populations but are not included in the plasma chemistry model used in this paper. These neglected processes include collisional de-excitation, spontaneous radiative de-excitation, and heavy–heavy reactions such as charge exchange. Future work should include the physics of these processes in order to more accurately track the densities of the excited states.

There is one additional reaction beyond those described above that is tracked in this paper and that is dissociative recombination of molecular nitrogen ions,

where the electron temperature *T _{e}* is in eV. Including this reaction provides an estimate for the production of atomic nitrogen via this pathway, and it also captures the reduction in plasma density due to recombination. This is the only reaction in this plasma chemistry model which allows for the decay of electron density, a phenomenon which is experimentally observed. With this reaction, we also have a rate equation for atomic nitrogen,

The factor of two on the right hand side comes from the fact that the dissociation of one molecular ion produces two atoms of nitrogen. Loss terms due to dissociative recombination are also added to the electron and ion rate equations. Note that since the ion density $nN2+$ is equal to the electron density in this quasi-neutral model, the source term for atomic nitrogen and the associated loss terms in the electron and ion rate equations are second order in *n _{e}*. All other processes considered in this paper are first order in

*n*.

_{e}### C. Plasma dynamics

The way that the plasma chemistry is used depends largely on the model used for the plasma dynamics. For example, in a kinetic model, the plasma chemistry is represented by a set of energy-dependent collisional cross sections. In a fluid model, a set of collisional rates is determined by Eq. (11) with assumptions about the electron energy distribution function. The fluid model described in this section will be used to describe the plasma dynamics in the remainder of this paper.

#### 1. The fluid equations for the thermal plasma

The plasma electrons can be treated as a fluid provided the collision rate is sufficiently large such that the mean-free time between collisions is small compared to the rise time of the beam, and the mean-free path is short compared to the length scales in the problem. If $V\u2032$ and $U\u2032$ are the mean velocity and energy of the secondary electrons created from ionization events by thermal electrons, and $Vs$ and *U _{s}* are the mean velocity and energy of secondary electrons from beam-impact ionization, then the thermal electron fluid equations for the beam–plasma system in conservative form are given by

where *p _{e}* is the plasma electron pressure, $\nu N2+$ is the total thermal electron ionization rate,

*ν*is the momentum transfer frequency, $\nu \epsilon $ is inelastic energy exchange frequency, and $dnedt|beam$ is the beam-impact ionization rate. The last two terms on the right-hand sides of each equation represent the volumetric rates for adding or removing electron density, momentum, and energy during ionization and recombination events. To close the equations, the pressure must be given in terms of the fluid variables or it can be neglected if appropriate for the parameters of interest (Sec. III C 3).

_{m}#### 2. Beam-impact ionization

Ionization of neutral gas by energetic beam particles plays a significant role in the overall plasma evolution, especially at higher gas pressure. Early in time, the energetic beam electrons collisionally ionize the background gas and start to generate a plasma. This plasma produced by the electron beam also provides seed electrons for thermal excitation and ionization processes.

The fluid model in this paper specifies the beam-impact ionization rate as well as the average momentum and energy of the secondary electrons from beam impact events. These quantities are needed for the source terms in the electron fluid equations. One common model is to represent the beam impact rate by

where $Jb$ is the beam current density, *ρ _{G}* is the gas density in g/cm

^{3},

*S*is the beam-electron-energy-dependent stopping power in eV-cm

^{2}/g, and

*W*is the average energy expended by the beam per electron–ion pair.

^{12}The stopping power approach assumes that energetic secondary electrons produced by individual beam-impact events both excite the gas and produce additional electron–ion pairs as they slow down and join the thermal population. A detailed analysis for N

_{2}gas shows that the value of $W\u224336\u2009eV$ is nearly constant for beam energies above 1 keV.

^{12}The value of

*W*is quite a bit higher than the ionization energy since there are many competing inelastic excitation processes which act as energy sinks as the energetic secondary electrons produced by the beam slow down. Since Eq. (19) assumes that secondary electrons created by beam impact rapidly join the thermal population, it is reasonable to take $Vs=Ve$ and

*U*=

_{s}*U*in Eqs. (17) and (18).

_{e}To solve Eq. (19), the stopping power of electrons in nitrogen is needed. The total stopping power in N_{2} (collisional plus radiative bremsstrahlung losses) can be calculated by the NIST electron stopping power program ESTAR.^{21} The Bethe–Bloch formula^{22} is an approximate expression for the stopping power, which agrees well with the collisional part of the ESTAR stopping power below 100 MeV in nitrogen. Since this paper is concerned with e-beam energies below 100 keV, radiative stopping can safely be ignored in our analyses, and the Bethe–Bloch formula for stopping power from Ref. 22 is used.

The beam-impact ionization rate in Eq. (19) is valid for the case where the beam electrons are monoenergetic and have no spread in velocity (which could be due, for example, to scattering of the beam electrons in the foil when they enter the gas chamber). A spread in either energy or velocity could change both the stopping power and the effective path length of beam electrons. These effects are beyond the scope of this work.

In addition to ionization events, the beam can directly excite the gas. For N_{2} gas, the number of excited species produced by the e-beam per electron ion pair is relatively constant for beam energies above 1 keV.^{12} The number of excited species of type *j* created per electron–ion pair is called the production efficiency, which is denoted $\xi N2,j*$. A term is added to the rate equations for the neutral species to account for these beam-impact reactions,

The beam excitation production efficiencies are taken from Ref. 12 for a 1 keV beam and modified slightly to be compatible with the Phelps cross section set that is used to compute thermal rates.

The single Phelps ionization cross section is a lumped cross section, which includes all final states of the $N2+$ ion as well as dissociative ionization to N^{+}, whereas each of these ionization channels is tracked separately in Ref. 12. Instead of tracking these final states separately, we sum all of the ionization channels from Ref. 12 to obtain a single beam-impact ionization production efficiency, with the final state being a generic “ion.” In our model, this product species is added to the $N2+$ species, which now represents a generic ion final state.

Two other modifications are also made. First, the triplet state A$3\Sigma $ in the Phelps dataset is vibrationally resolved in three bands, while the A$3\Sigma $ production efficiency in Ref. 12 is not vibrationally resolved. In this paper, 1/3 of the beam-driven excitation of the A$3\Sigma $ state is assigned to each of the three vibrational bands. Second, the dissociation channels in Ref. 12 must be reconciled with the sum-of-singlets state that appears in the Phelps dataset. Both of these processes are consistent with the dissociation cross section in Ref. 19, with the Phelps sum-of-singlets-state cross section being larger than that dissociation cross section by (1/0.73), where 0.73 is the dissociation fraction in Ref. 19. Therefore, in our model, we take the total beam-impact dissociation efficiency from Ref. 12, divide by the dissociation fraction of 0.73 from Ref. 19, and attribute the result to the Phelps sum-of-singlets-state. The production efficiencies used in our model are summarized in Table I.

#### 3. Collisional frequencies

As mentioned above, the fluid equations (16)–(18) need frequencies that describe the interaction of the thermal electrons with the background gas. There are four frequencies for thermal interactions that need to be specified for the electron fluid: the ionization frequency, the momentum transfer frequency, the total inelastic energy transfer frequency, and the dissociation frequency. The dissociation frequency is $\nu dr=nN2+kdr$, where the dissociative rate constant *k _{dr}* is given in (14).

The collisional ionization frequency by the thermal electrons is given by

where $kN2+$ is the rate coefficient defined in Eq. (11).

The electron momentum equation requires the momentum transfer frequency for collisions between thermal electrons and neutral gas molecules. This frequency is defined by

where *σ _{m}* is the total momentum transfer cross section. This cross section is from all processes including elastic collisions, ionizing collisions, and all other inelastic collision processes.

In the energy equation, an expression for the inelastic energy transfer between thermal electrons and gas molecules is needed. The inelastic energy-transfer frequency is defined by

where $\epsilon N2,j*$ is the threshold energy and $\nu N2,j*$ is the inelastic collision frequency for the process with product $N2,j*$. The sum is over all inelastic collision processes. The inelastic collision frequency is given by

where $kN2,j*$ is the rate coefficient for producing the excited state $N2,j*$.

These frequencies cannot be evaluated without an estimate of the eedf, *f*_{0}. If the ionization fraction is low and Coulomb collisions are infrequent, the dynamics of the eedf is dominated by the energy transfer between the electrons and the neutral gas molecules via inelastic collision processes. In this case, a steady-state two-term approximation for the eedf $f0(\epsilon )$ can be used to calculate the collision frequencies. A code like BOLSIG+ can be used to calculate the eedf in this approximation.^{23} If the ionization fraction is high and Coulomb collisions between electrons and other free electrons in the plasma dominate the total inelastic collision frequency, then the Coulomb collisions in the plasma will drive the eedf to a Maxwellian distribution. In this case, the eedf is

where $Te=23Ue$ is the electron temperature. Figure 1 shows the differences between the steady state eedf computed by BOLSIG+ and a Maxwellian eedf for two different values of the mean electron energy *U _{e}*. At high mean energy (

*U*= 10 eV, shown in red), there is not much difference between the two eedfs, and it is valid to approximate the eedf by a Maxwellian. At low mean energy (

_{e}*U*= 1 eV, shown in blue), the steady-state eedf differs significantly from Maxwellian. Energy losses during inelastic collisions cause much of this difference. If there are no electron–electron Coulomb collisions to cause $f0(\epsilon )$ to spread in energy and relax to a Maxwellian, then the shape of $f0(\epsilon )$ is set by a balance between acceleration of electrons in the electric field and energy loss due to inelastic collisions.

_{e}In order to decide which approximate eedf to use in this work, a comparison can be made between the inelastic collision frequency and an estimate of the electron–electron collision frequency. In cgs units (but with temperature in eV), this collision frequency is^{24}

Figure 2 shows a log –log plot of the inelastic collision frequency computed using a Maxwellian eedf (red), the inelastic collision frequency computed using a steady-state eedf (blue), and estimates of $\nu ee$ made using ionization fractions of $ne/Ng$ = 0.01%, 0.1%, and 1%, and $ln\u2009\Lambda =10$ (dashed black lines), where *N _{g}* is the number density of the background neutral gas. As this plot shows, electron–electron collisions occur much more frequently than inelastic collisions for $Ue\u22721\u2009eV$, even for plasmas with low ionization fraction. The electron–electron collisions drive the eedf toward a Maxwellian distribution, while the inelastic collisions drive it away from Maxwellian. In the low energy region where the Maxwellian and two-term collision frequencies differ more, the electron–electron collisions occur more often than the inelastic collisions, and we would expect them to drive the eedf toward a Maxwellian. Because of this, we will use collision rates based on a Maxwellian eedf assumption. More accurate collision frequencies could be computed by coupling this model to a time-dependent solver for the eedf; however, this is beyond the scope of this paper.

Figure 3 shows a summary of the collision frequencies used in this work computed using a Maxwellian eedf. The recombination frequency in Fig. 3(a) is shown for an ionization fraction of $ne/Ng=0.1$%. The total inelastic energy transfer frequency shown in Fig. 3(b) as a solid red line is the sum of contributions from all of the inelastic processes. The dashed lines show how the various types of process contribute to the inelastic energy transfer frequency. Since energy transfer due to elastic collisions between electrons and the neutral molecules is reduced by their mass ratio, it is very small compared to inelastic energy transfer and is neglected in this work. Note that the ionization frequency is nonzero even for mean energies below the ionization threshold. This is due to the tail of the eedf extending above the threshold, even when the mean energy is below threshold.

#### 4. Simplified fluid equations

To proceed with the fluid equations, it is also necessary to make some assumptions about the average velocity and energy for secondary electrons produced during ionization events by thermal electrons. It is common to assume that secondary electrons are produced with an isotropic distribution. With this assumption, the average secondary velocity is given by $V\u2032=0$. Since momentum is conserved during the ionizing collision, any momentum lost by the primary electron must show up in the ion. Since *ν _{m}* represents the momentum transfer for both inelastic and elastic collisions, the momentum transfer to ions is already included. Ignoring any kinetic energy in the ion, the average secondary kinetic energy is approximated by $U\u2032=Ue$ in the electron energy equation. Energy lost to ionization is already one of the terms in the inelastic energy loss. It is also assumed that $Vs=Ve$ and

*U*=

_{s}*U*.

_{e}Because of the symmetry assumptions in the rigid-beam model, Eqs. (4) and (5), the fluid equations simplify because the divergence terms on the left-hand sides of Eqs. (16)–(18) are zero. Another simplification comes from comparing the $Ve\xd7B$ term with the term due to momentum transfer with the neutrals. The ratio of these terms is proportional to the Hall coefficient, $CH=\omega ce/\nu m$, where the electron cyclotron frequency is $\omega ce=eB/me$. For the parameters being considered, the Hall coefficient is fairly small, and so the $Ve\xd7B$ term will be neglected. However, at lower pressures, the Hall coefficient can reach a maximum value greater than one, as will be discussed in Sec. IV C. This is an indication that, at lower pressures, the magnetic field and the radial component of the electric field may need to be included.

The pressure gradient term in Eq. (17) will also be neglected in this model. This approximation is justified because the pressure gradient gives rise to a radial electric field, which holds the electrons back and accelerates the ions outward radially. The resulting expansion of the plasma happens at the ion sound speed, which is small for the parameters of this plasma.

With these simplifications, the fluid equations (16)–(18) for the plasma electrons become

where we have defined the effective momentum transfer frequency as

Note that only the *z* component of the momentum equation (29) remains and that it can be used to obtain an equation for the plasma current,

In this equation, we have used the definition of the conductivity,

The second two terms on the right-hand side of Eq. (31) are orders of magnitude smaller than the momentum transfer frequency *ν _{m}* for the parameters of our beam. Early in the beam pulse, the beam impact ionization term can be larger than the other terms but it quickly becomes small as the electron density rises. Because these terms are small, they will be neglected in the remainder of this paper, and we will assume that $\nu eff\u2243\nu m$.

## IV. NUMERICAL SIMULATIONS

The equations for the rigid-beam model, together with the fluid response model for the plasma, have been solved numerically for parameters selected to match experiments that were performed at NRL. In these experiments, a pulsed electron beam is injected into a cylindrical chamber filled with low-pressure gas. Experimental diagnostics measure various aspects of the plasma response, including the net current, the line-integrated plasma density at several locations, and framing camera images and spectra of the light emitted by the plasma. The pulse length of the electron beam is about 100 ns, the peak current in the pulse is about 4 kA, and the beam diameter is about 4 cm full width at half maximum. The pressure in the gas chamber ranges from about 100 mTorr to 10 Torr. Details of the experiments and analysis of the experimental results were reported elsewhere.^{10}

The simulations reported in this section are driven with beam current and voltage histories derived from experimental measurements. The time variations of the diode current are shown in Fig. 4 (blue curve), with a blue shaded region that represents experimental measurement uncertainty. Shot-to-shot variations in the diode current give a smaller uncertainty and are thus neglected in this work. Electrons are extracted from the diode and injected into a gas-fill transport chamber through a pair of 8.52-*μ*m-thick aluminized Kapton foils. The beam is attenuated and scattered by the foils, which results in a reduction of both the energy and current that the beam has as it enters the transport section. Monte Carlo simulations of the current loss provided an estimate of the beam current entering the gas cell. This current is shown in red in Fig. 4, along with a shaded region that represents the uncertainty in the diode current measurement. The normalized radial profile of the beam is inferred from experimental measurements and shown in Fig. 5.

The energy of the beam electrons is estimated from the measured diode voltage, which is shown in Fig. 6 (blue curve) along with an estimate of the uncertainty in the voltage measurement (blue shaded region). The Monte Carlo calculations described above also give an estimate of the scattering and slowing of the beam electrons in the foils. The statistical nature of the energy loss in the foils introduces a spread in the energy of the electrons that enter the gas chamber. The red curve in Fig. 6 shows the mean energy of these beam electrons. The red shaded region comes from propagating the uncertainty in the measured diode voltage through the Monte Carlo simulations. Since it is calculated using nonlinear uncertainty propagation in this way, it accounts for both the energy spread induced by electron scattering as well as the experimental uncertainty in the voltage measurements. When the diode voltage is low, the electrons are completely stopped in the foil. This leads to a narrower beam pulse as compared to the diode current, which can be seen in Fig. 4.

In the simulation results that are reported in this section, estimates of uncertainty in the simulation outputs are obtained by running a series of simulations with different combinations of beam histories. For beam current (and beam energy), there are three options: the mean value, the upper uncertainty bound, and the lower uncertainty bound. Combining these options gives a total of nine simulations that are run in order to get a simulation result with an estimated uncertainty. In the plots of the simulation results shown in this section, the mean of these nine results will be shown with a solid line, and the standard deviation will be shown using shaded uncertainty bounds.

Note that there are additional parameters of the simulations (such as tabulated collision rates, gas pressure, initial electron density, or beam radial profile) that have experimental uncertainties associated with them. These uncertainties are not accounted for in the simulations presented here. Only the uncertainties in the beam current and energy have been propagated through the simulation. A more complete and systematic investigation of these other uncertainties is beyond the scope of this work.

### A. Custom turboPy modules for the rigid-beam model

To quickly implement and test the various parts of the rigid-beam model, we use the Python physics simulation framework turboPy.^{25,26} This framework manages boilerplate aspects of the code such as file I/O and problem configuration, while also providing basic computational physics functionality such as defining the main simulation loop, a clock, and a grid.^{27} For this rigid-beam problem, we have created custom turboPy PhysicsModules which solve the electric field equation, the rate equations, and the electron fluid equations. Additional PhysicsModules read in tabulated rate coefficients and perform the rate table lookups, while others calculate the stopping power and beam-impact ionization and excitation rates. In addition to these PhysicsModules, several custom Diagnostics were written to compute quantities, such as the line-integrated electron density (for comparison with experiment) and the electron–electron collision frequency (for checking the validity of the weakly ionized model). While most of these PhysicsModules use standard numerical techniques, there is one whose details are worth describing.

The equations for the electric field and the plasma current are closely coupled in this model. In particular, the electric field equation (6) does not have the form of a temporal evolution equation, but rather is a Poisson-like second-order differential equation that must be inverted to find *E _{z}*. In the code, this is done by combining it with the equation for the plasma current to obtain

The electric field and current densities are discretized onto a radial grid that extends from 0 to 8.65 cm with 100 grid points. The operator acting on *E _{z}* on the left-hand side of this equation is written as a finite difference matrix, and the electric field is updated at each time step using the solve function from the SciPy package scipy.linalg.

^{28}This function uses the industry-standard LAPACK library to solve the matrix equation using LU decomposition. The values of the electric field and plasma current at time step

*n*are used when calculating the plasma current at time step

*n*+ 1,

The rate equations for the density of the electrons, ions, ground state neutral N_{2}, and all excited states are solved numerically using reaction rates that are tabulated as functions of the mean electron energy. The mean electron energy is obtained by solving Eq. (30) numerically. The forward Euler method is used for solving these density and energy equations.

The simulations shown in the remainder of this section were produced using various open-source packages, most notably turboPy, NumPy, Matplotlib, and xarray.^{25,26,29–31}

### B. Simulation results

To compare the results from the rigid-beam model with experiments, the rigid-beam turboPy program described above was run with background gas pressures of 0.7, 1, 3, 5, 7, and 10 Torr. The grid spacing and size of the time step were chosen so that the simulation is numerically converged (see the Appendix for details). While the experiments were performed in dry air, and the model only includes nitrogen, the model results are expected to at least qualitatively reproduce the experimental measurements.

Figure 7 shows the radial profiles of several quantities during the beam rise, at 40 ns, from the simulation with 1 Torr background gas pressure. The beam current [shown in blue in Fig. 7(a)] is confined within about 3 cm radius, as are the plasma electrons [Fig. 7(c)]. The axial electric field [Fig. 7(b)] drives an axial current in the plasma [shown in red in Fig. 7(a)], which cancels some of the beam current, resulting in a net axial current [shown in black in Fig. 7(a)] which peaks near a radius of 1.7 cm.

Figure 8 shows a comparison of the measured net current^{32} (green), the beam current that drives this simulation (red), and the simulated net current (black), at 1, 5, and 10 Torr. The net current initially rises with the beam current, but then the plasma current starts to cancel out a portion of the beam current, giving a lower net current. Once the beam pulse has ended, the inductive electric field keeps the current flowing in the plasma and the plasma current slowly decays away because of the finite plasma resistivity. The differences between the simulated and measured net currents, especially at early times, could be due to approximations in the model or in the current pulse used to drive the simulations. These differences will be examined more closely in future work.

The line-integrated electron density is measured in the experiments using an interferometer with a line-of-sight along a chord through the plasma. A synthetic diagnostic was used in the simulation to obtain a similar output, where the chord is taken through the center of the plasma. Figure 9 shows this line-integrated plasma electron density as a function of time for both simulation (black) and experiment (blue) at 1, 5, and 10 Torr. The generation of plasma through beam-impact and thermal ionization gives densities that rise during the beam pulse, and then dissociative recombination causes the densities to decrease later in time. At lower gas pressures, beam-impact ionization has a smaller effect, and the density peaks at a smaller value. These trends are seen in both the experimental data and the simulation outputs. The simulations reproduce the experimental trends with gas pressure, and the simulated peak value of line-integrated density is about a factor of two lower than the experiment. There are several factors that could contribute to these differences, such as uncertainties in the rates used in the plasma model, or the fact that the model only includes nitrogen, but the experiments were in air. Further work is planned to study these differences.

In addition to experimentally motivated synthetic diagnostics, additional diagnostics can be added to the simulation in order to probe the dynamics of the model. For example, the value of the induced electric field on axis as a function time is shown in Fig. 10 for the 1, 5, and 10 Torr simulations. The rapidly rising beam current density initially induces a very large electric field. As the plasma becomes more conductive, the plasma current tends to respond quickly enough to cancel out rapid changes in the beam current, and thus the magnitude of the induced electric field is smaller later in time. Once the beam pulse has ended, the stored magnetic energy in the chamber slowly dissipates through the resistive plasma. The electric field associated with this dissipating energy can be seen in the long, late time tail in Fig. 10, and is generally much smaller than the peak field early in the pulse.

The plasma response model can also be probed through the use of custom diagnostics. Figure 11 shows the mean energy of the plasma electrons on axis for the 1, 5, and 10 Torr simulations, with peaks in energy that correlate with the induced electric fields, showing Ohmic heating of the plasma electrons. The peak at early time, around 20 ns, corresponds to the large electric field that is induced during the rise of the beam current. The smaller peak at late time, around 120 ns, corresponds to the end of the beam pulse, which also induces an electric field which heats the electrons.

In Figs. 12 and 13, the time histories of the on-axis densities of the N_{2}(*v* = 1) vibrational state and the N_{2}(B$3\Pi $) electronic state of nitrogen are shown for simulations at 1, 5, and 10 Torr. These states exhibit changes in density due to both beam-impact and thermal excitation. As can be seen by the increasing densities late in time, the thermal plasma electrons continue to excite these states even after the pulse ends. At 10 Torr, excitation of the N_{2}(B$3\Pi $) happens almost exclusively through beam-impact events, which can be understood by examining the energy of the thermal electrons in Fig. 11. This state has a threshold energy of 7.35 eV, but the electrons in the 10 Torr simulation are significantly lower energy than this for nearly the entire simulation. Note that the model used for these simulations does not include any means by which the densities of the excited neutral states could decrease. Future work should include, for example, the physics of collisional de-excitation in order to more accurately track the densities of these states and estimate optical emission spectra. The densities of all the excited states at 125 ns (near the peak of the line-integrated densities) are summarized in Table II for the 10 Torr simulation.

Species . | Density at 125 ns (%)^{a}
. |
---|---|

X $\u20091\Sigma $ | 95.15 |

rot | 2.51 |

v = 1, 2, 3, 4 | 1.12, 0.36, 0.24, 0.14 |

v = 5, 6, 7, 8 | 0.11, 0.08, 0.04, 0.02 |

A$\u20093\Sigma $ ($v=0\u22124$), A$\u20093\Sigma $ ($v=5\u22129$) | 0.01, 0.01 |

B $\u20093\Pi $, W $\u20093\Delta $, A$\u20093\Sigma $ ($v\u226510$), $B\u2032\u20093\Sigma $ | 0.02, 0.01, 0.01, 0.00 |

$a\u2032\u20091\Sigma $, a $\u20091\Pi $, w $\u20091\Delta $, $C3\Pi $ | 0.00, 0.01, 0.00, 0.01 |

E $\u20093\Sigma $, $a\u2033\u20091\Sigma $, sum | 0.00, 0.00, 0.07 |

$N2+$ | 0.05 |

Atomic N | 0.10 |

Species . | Density at 125 ns (%)^{a}
. |
---|---|

X $\u20091\Sigma $ | 95.15 |

rot | 2.51 |

v = 1, 2, 3, 4 | 1.12, 0.36, 0.24, 0.14 |

v = 5, 6, 7, 8 | 0.11, 0.08, 0.04, 0.02 |

A$\u20093\Sigma $ ($v=0\u22124$), A$\u20093\Sigma $ ($v=5\u22129$) | 0.01, 0.01 |

B $\u20093\Pi $, W $\u20093\Delta $, A$\u20093\Sigma $ ($v\u226510$), $B\u2032\u20093\Sigma $ | 0.02, 0.01, 0.01, 0.00 |

$a\u2032\u20091\Sigma $, a $\u20091\Pi $, w $\u20091\Delta $, $C3\Pi $ | 0.00, 0.01, 0.00, 0.01 |

E $\u20093\Sigma $, $a\u2033\u20091\Sigma $, sum | 0.00, 0.00, 0.07 |

$N2+$ | 0.05 |

Atomic N | 0.10 |

^{a}

Density is expressed as a percent of the initial neutral density, $3.53\xd71017\u2009cm\u22123$.

### C. Simulation results as functions of pressure

Because of the simplifications and reductions in the rigid-beam model, the numerical simulations described above run quickly, and parameter sweeps can be performed without significant computational cost. This allows for the rapid evaluation of chemistry models, and pressure scans can be easily performed. To illustrate this, a set of simulations were performed for gas pressures ranging from 0.7 to 10 Torr. For each simulation, the peak values of several quantities are recorded and are shown as functions of pressure in Figs. 14–19. At each pressure, the dots show the mean of the nine simulations with different combinations of beam current and energy histories, and the error bars show standard deviation. The dashed lines connect the mean values to illustrate the variation with pressure. Note that simulations at pressures lower than 0.7 Torr are not shown because at those low pressures, the mean plasma electron energy reaches values greater than 100 eV. This is the maximum energy value for the reaction rates tables used in these simulations, and reaching this value indicates that the chemistry model used in this paper is likely insufficient at these low pressures.

Figure 14 shows the peak value of the line-integrated electron density. Since the beam-impact ionization rate is proportional to the number density of the background gas, the density of plasma produced by beam-impact ionization should increase with pressure. As seen in the figure, this is what happens for pressures above about 3 Torr. At lower pressures, the effect of thermal ionization starts to become more important, and the density is higher than would be expected from beam-impact ionization alone. The trend of density increasing with gas pressure in the experiments is captured in the simulations. However, the peak line-density from the simulations is lower than the measured value at all pressures. This discrepancy is slightly larger than the ±10% uncertainty in the measurements. The overall comparison of the model with the data is within a factor of about 1.5 at all pressures and the experimental trends are reproduced. This is an indication that the model captures the essence of the physics and provides an excellent framework to do systematic studies of the beam–plasma system. An effort is under way to further investigate the source of the remaining discrepancy. Also, note that the experimental data at pressures lower than 0.7 Torr are outside the range of model validity and are thus not shown in Fig. 14.

Figure 15 shows the peak value of the on-axis electric field. While this value does vary somewhat with pressure, the values do not vary by more than about 20%. This is not the case for the mean electron energy, however. Figure 16 shows the on-axis value of the mean electron energy. This quantity does strongly depend on pressure, varying by nearly a factor of 10, from a minimum of about 6 eV at 10 Torr to a maximum of above 100 eV at 0.7 Torr. For the simplified plasma chemistry that is being used in this work, there are important electron scattering processes at higher energy that have been neglected. At mean energies above about 10 eV, there are likely to be processes such as stepwise ionization or ionization to excited states of $N2+$, which can significantly affect the electron energy. The results in Fig. 16 indicate that this simplified chemistry should be used with care below about 5 Torr, and below about 2 Torr, the model could be giving incorrect results because of the missing plasma chemistry. The results at 0.7 Torr give peak energies above 100 eV. Since the rate tables used in the work only extend to 100 eV, these results indicate that the model should not be used at any lower pressures.

The magnetic field in each simulation was calculated from the net current using Ampère's law Eq. (2). This was used to calculate the electron cyclotron frequency so that the Hall coefficient C_{H} for each simulation could be calculated. As shown in Fig. 17, the maximum value (for any radius and time) of C_{H} is less than one at higher pressures but above one for simulations with pressures of 3 Torr and below. This is an indication that the assumptions made in this model may not be valid at these lower pressures and that the $Ve\xd7B$ term in the momentum equation should be retained.

There is another quantity that can be extracted from the simulations in order to obtain insights into the physics of formation of the plasma. Since thermally driven ionizations are computed separately from beam-impact ionizations, the code can track how much of the plasma was created by thermal processes and how much by beam impact. The fraction of ionization events due to thermal processes is calculated by taking the ratio of electron density produced thermally to the total electron density. The on-axis value of this ratio at 200 ns is show in Fig. 18 as a function of pressure. Above about 2 Torr, the fraction of plasma produced thermally drops below 50%, indicating that at higher pressure beam-impact ionization is the dominant process by which the plasma is produced. This is consistent with the approximately linear scaling of peak plasma density with pressure shown in Fig. 14.

The last simulation output that is useful to examine is a summary of the dynamics of the excited states of N_{2}. Plots such as those in Figs. 12 and 13 can be helpful for predicting experimental outputs that depend on the densities of excited states, such as spectral measurements of the optical emission once the necessary processes are included. In order to condense this information in a way that would be helpful for discerning trends, a series of sparkline plots were created and are shown in Fig. 19. Each row is a different state, and the states are sorted in order of increasing threshold energy (the state name and threshold energy are listed to the right of each corresponding row). The left column of this plot shows the density vs time sparkline from the 1 Torr pressure simulation. The second column shows density vs time from the 10 Torr simulation. This column shows some differences from the 1 Torr case. Most noticeably, the fact that beam-impact excitation is more important for the states with threshold above ∼6 eV at 10 Torr than it is for 1 Torr. This is seen in the difference in the shapes of the density curves; in the 1 Torr case, the density of these states continues to increase after the pulse is over, which for the 10 Torr case the density is constant after the pulse. The third column of sparklines shows how the peak value of density varies with pressure. These are sparkline versions of plots that are analogous to the plot in Fig. 14 for electron density. Several of the species, such as the $N2+$ ion and the N_{2} (*rot*) state, increase with pressure. Others, however, such as the B$3\Pi $ and C$3\Pi $ electronic states, actually decrease with pressure. It would be interesting if these trends could be observed experimentally through spectroscopic or other means. In order to use this rigid-beam code to predict optical spectra, a more complete model of the plasma chemistry is required, including a model of the radiation dynamics.

Plots like Fig. 19 can be valuable for identifying areas where the plasma chemistry model should be improved. For example, in the *p* = 10 Torr case, densities for all excited states with thresholds above 6.17 eV are unphysically constant after the beam is off, showing clearly the impact of neglecting de-excitation processes.

## V. CONCLUSIONS

In this paper, we define a reduced model of the electrical breakdown of a low-pressure gas caused by the injection of an intense electron beam. Several simplifying approximations are made to the equations for the beam dynamics and the electromagnetic fields. Much of the physics is driven by the electric field that is induced by the changing beam current, and the *rigid-beam model* expresses the dynamics of this field through a one-dimensional second-order Poisson-like differential equation.

The rigid-beam model, when combined with a model for plasma that is created by beam-impact reactions and reactions driven by the electric field, forms a complete set of equations for the dynamics of the beam/plasma system. A fluid model of the plasma electrons is presented, along with a weakly ionized plasma chemistry which describes the electron-driven scattering processes in low pressure molecular nitrogen. The prescribed properties of the injected beam, provided by experimental measurements, appear in source terms to the electron fluid equations.

The coupled rigid-beam equations and the plasma response are solved with a code written using the turboPy framework. Numerical results are presented for experimentally relevant parameters, showing that this model is useful for computing the beam-driven breakdown of a low-pressure gas. Because of the simplifications inherent to the rigid-beam model, the code runs quickly. This enables rapid numerical exploration of the parameter space, which includes variation with pressure, electron beam current density, and the plasma chemistry model. In addition to running simulations for a range of model parameters, it is also straightforward to expand the code through additional PhysicsModules which implement different reduced plasma models.

There are several ways in which the model described in this paper expands on previous work. Previous papers which used a rigid-beam model to simulate e-beam driven breakdown had several limitations which prevented them from accurately modeling the experiments described in this paper. Reference 13 used a 0D model for the beam, which has been expanded to a 1D radial model in this paper. Reference 14 used a simplified plasma chemistry, which has been extended here to include beam-impact ionization and dissociative recombination. Additionally, a more realistic radial profile and temporal pulse shape have been used for the beam. By extending the rigid-beam model in these ways, higher-fidelity simulations are able to be performed. Even with all the simplifications in this model, the simulation results agree favorably with experimental measurements. The simulated line-integrated electron density increases with increasing pressure, which is a trend observed in experiments. While the simulated densities are consistently lower than measured in experiments, the results are within a factor of two of the measurements.

Future work is planned where this program will be used to map out regions of validity in parameter space for the reduced plasma response models used in this work. Extensions to the program are also planned, which will implement models that account for physics which has been neglected. For example, electron–ion collisions will be added to the conductivity, de-excitation processes will be added for tracking excited state populations, and advanced chemistries will be added which include state-to-state transitions. Additionally, this program can be used to perform sensitivity studies, for example, to determine the effect of uncertainties in the electron-scattering cross sections that are used in the weakly ionized plasma chemistry.

The RB model described in this paper employs a set of assumptions about the electromagnetic fields and the electron-beam dynamics which are useful for modeling the breakdown of low-pressure gas by an intense electron beam. This model also lays a foundation for future work studying approximations related to the plasma electron dynamics or related to the plasma chemistry. Systematic investigation of advanced plasma chemistries—which examine the effects of state-to-state transitions for example—can now be performed and compared *ceteris paribus* to results obtained with simpler plasma chemistry approximations.

## ACKNOWLEDGMENTS

The authors wish to thank Dr. P. Ottinger and Dr. T. Haines for useful discussions.

This work was funded by the Defense Threat Reduction Agency and the U.S. Naval Research Laboratory Base Program.

## DATA AVAILABILITY

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

### APPENDIX: CONVERGENCE

The grid and the time step are tested at various pressures in order to ensure that the numerical solution is acceptably converged. The asymptotic approach in Ref. 33 allows us to relate the approximate result to the expected result as we refine the grids in space and time. The method is a uniform refinement method generalized from Richardson extrapolation theory. Here, the grid convergence index (GCI) gives some statistical significance to the percent error away from the converged solution. The advantages of using this method for the rigid-beam problem are that it allows for quick and confident testing of space and time resolution used to produce the results in the paper.

The electric field on axis was examined in this grid spacing convergence study. The electric field in the rigid-beam model is computed using a second order discretization in space; thus, this is a natural quantity to use for this convergence study. The grid extends from *r* = 0 to *r* = 8.65 cm. The grid spacing is refined from 10 grid points, to 30 grid points, to 90 grid points, and the results at different pressures are reported in Table III. Here, the theoretical order of convergence is $O(2)$ and what we found as the actual order of convergence is $O(1.6)\u223cO(1.7)$. The extrapolated value *h*_{0} is what the approximate converged value of the electric field is on axis at the end of the simulation, and the GCI is the error band from the fine grid (90 points) to the medium grid (30 points). The results show that as the pressure increases, the error gap improves and the solution is closer to the converged value, but all pressures are acceptably converged.

Pressure . | Order of convergence ($O$) . | Extrapolated value, h_{0}
. | GCI (%) . |
---|---|---|---|

1 | 1.677 | 475 | 8.71 |

5 | 1.579 | 2063 | 1.57 |

7 | 1.675 | 2159 | 0.857 |

10 | 1.704 | 1991 | 0.614 |

Pressure . | Order of convergence ($O$) . | Extrapolated value, h_{0}
. | GCI (%) . |
---|---|---|---|

1 | 1.677 | 475 | 8.71 |

5 | 1.579 | 2063 | 1.57 |

7 | 1.675 | 2159 | 0.857 |

10 | 1.704 | 1991 | 0.614 |

Because of the optimized solver, vectorized operations, and relatively small number of grid points used in these simulations, increasing the number of spatial grid points did not affect the runtime of the simulations very much. For the 1 Torr simulation with time step of $1\xd710\u22123$ ns, the simulations with 10, 30, and 90 grid points took 441.6, 442.9, and 438.9 s to run, respectively.

Similarly, for time step refinement, we choose the net current as the variable to study for this experiment. In this experiment, we refine the time step from $2\xd710\u22123$ to $1\xd710\u22123$ ns, and to $5\xd710\u22124$ ns. Here, we expect the theoretical order of converge to be $O(1)$ because of the forward Euler method used for the time update. The results at different pressures are reported in Table IV. The results show the actual order of convergence is almost exactly the predicted value, and the error bands of this experiment are all under 1% with the best case being the 5 Torr result.

Pressure . | Order of convergence ($O$) . | Extrapolated value, h_{0}
. | GCI (%) . |
---|---|---|---|

1 | 0.9424 | 2354 | 0.42 |

5 | 1.008 | 3313 | 0.111 |

7 | 0.9855 | 3307 | 0.208 |

10 | 0.9824 | 3240 | 0.351 |

Pressure . | Order of convergence ($O$) . | Extrapolated value, h_{0}
. | GCI (%) . |
---|---|---|---|

1 | 0.9424 | 2354 | 0.42 |

5 | 1.008 | 3313 | 0.111 |

7 | 0.9855 | 3307 | 0.208 |

10 | 0.9824 | 3240 | 0.351 |

Changing the size of the time step used in these simulations had a much larger effect on the run time. Simulations with 100 grid points at 1 Torr pressure and time steps of $2\xd710\u22123$, $1\xd710\u22123$, and $5\xd710\u22124$ ns took 240.7, 485.6, and 949.4 s to run, respectively.

For the rigid-beam model and the results shown in this paper, we find that the use of 90 grid points and time step of $5\xd710\u22124$ ns is acceptable. For the results shown in Sec. IV, we used $5\xd710\u22124$ ns as the time step and rounded up to 100 grid points for the spatial resolution.