A GLOBAL MODELING FRAMEWORK FOR PLASMA KINETICS: DEVELOPMENT AND
APPLICATIONS
By
Guy Morland Parsey

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Physics – Doctor of Philosophy
Electrical Computer Engineering – Doctor of Philosophy
2017

ABSTRACT
A GLOBAL MODELING FRAMEWORK FOR PLASMA KINETICS: DEVELOPMENT AND
APPLICATIONS
By
Guy Morland Parsey
The modern study of plasmas, and applications thereof, has developed synchronously with computer capabilities since the mid-1950s. Complexities inherent to these charged-particle, manybody, systems have resulted in the development of multiple simulation methods (particle-in-cell,
fluid, global modeling, etc.) in order to both explain observed phenomena and predict outcomes
of plasma applications. Recognizing that different algorithms are chosen to best address specific
topics of interest, this thesis centers around the development of an open-source global model framework for the focused study of non-equilibrium plasma kinetics. After verification and validation
of the framework, it was used to study two physical phenomena: plasma-assisted combustion and
the recently proposed optically-pumped rare gas metastable laser.
Global models permeate chemistry and plasma science, relying on spatial averaging to focus
attention on the dynamics of reaction networks. Defined by a set of species continuity and energy
conservation equations, the required data and constructed systems are conceptually similar across
most applications, providing a light platform for exploratory and result-search parameter scanning. Unfortunately, it is common practice for custom code to be developed for each application—
an enormous duplication of effort which negatively affects the quality of the software produced.
Presented herein, the Python-based Kinetic Global Modeling framework (KGMf) was designed to
support all modeling phases: collection and analysis of reaction data, construction of an exportable
system of model ODEs, and a platform for interactive evaluation and post-processing analysis. A
symbolic ODE system is constructed for interactive manipulation and generation of a Jacobian,
both of which are compiled as operation-optimized C-code.
Plasma-assisted combustion and ignition (PAC/PAI) embody the modernization of burning fuel
by opening up new avenues of control and optimization. With applications ranging from engine

efficiency and pollution control to stabilized operation of scramjet technology in hypersonic flows,
developing an understanding of the underlying plasma chemistry is of the utmost importance.
While the use of equilibrium (thermal) plasmas in the combustion process extends back to the advent of the spark-ignition engine, works from the last few decades have demonstrated fundamental
differences between PAC and classical combustion theory. The KGMf is applied to nanoseconddischarge systems in order to analyze the effects of electron energy distribution assumptions on
reaction kinetics and highlight the usefulness of 0D modeling in systems defined by coupled and
complex physics.
With fundamentally different principles involved, the concept of optically-pumped rare gas
metastable lasing (RGL) presents a novel opportunity for scalable high-powered lasers by taking
advantage of similarities in the electronic structure of elements while traversing the periodic table. Building from the proven concept of diode-pumped alkali vapor lasers (DPAL), RGL systems
demonstrate remarkably similar spectral characteristics without problems associated with heated
caustic vapors. First introduced in 2012, numerical studies on the latent kinetics remain immature.
This work couples an analytic model developed for DPAL with KGMf plasma chemistry to better understand the interaction of a non-equilibrium plasma with the induced laser processes and
determine if optical pumping could be avoided through careful discharge selection.

Copyright by
GUY MORLAND PARSEY
2017

Cette thèse est dédié à ma grand-mère, Françoise Bratt

v

ACKNOWLEDGEMENTS

The journey I have taken through graduate school was only made a reality with the support and
kindness of many people along the way. Professors and fellow graduate students at Michigan State
University have taught and changed me into the person I have happily become.
Firstly, I wish to express my deepest gratitude to my adviser Prof. John Verboncoeur for the
years of patience and support you have given me. It is through your invaluable guidance and encouragement that I have been able to complete my journey through graduate school, including
through moments of my self-induced frustration. Your tolerance of my shortcomings and willingness to allow for creativity, and the time it may require, has allowed me to grow academically,
professionally, and as a person. The information and skills you taught me, or motivated me to
learn, has advanced me beyond what I ever expected my time at MSU to provide me with. I cannot
thank you enough for your mentorship and friendship.
I want to thank Prof. Andrew Christlieb, for the knowledge, advice, and realism you shared
with me in class or conversation. Together with my other committee members, Profs. Martin Berz,
Philip Duxbury, and Brian O’Shea, I want to thank you all for your time and patience with my
dissertation and writing process.
Amongst the other professors, advisers, and students who I have had the fortune of learning
from, I want to thank Dr. Steve Lund and Dr. Francesc Roig for helping me along the academic
path that brought me to graduate school. I wish to extend this gratitude to Dr. Yaman Güçü and
(soon to be Dr.) Connor Glosser for the knowledge and philosophy shared with regards to work,
code, and life.
Finally, I want to thank my family for their unrelenting support that has pushed me, calmed
me, and made all of this possible. To my girlfriend, Melissa, thank you for holding me in hard
times and pushing me through tired ones. To my parents and sister, thank you for your endless
love, support, and motivation for me to try my hardest.

vi

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xi

CHAPTER 1 INTRODUCTION . . . . . . . . . . . . .
1.1 Plasma Fundamentals . . . . . . . . . . . . . . .
1.2 Plasma Simulation . . . . . . . . . . . . . . . . .
1.2.1 Kinetic . . . . . . . . . . . . . . . . . .
1.2.2 Fluid . . . . . . . . . . . . . . . . . . . .
1.2.3 Hybrid . . . . . . . . . . . . . . . . . . .
1.2.4 Global/Chemistry . . . . . . . . . . . . .
1.3 Kinetic Global Modeling Framework Motivation .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

. 1
. 3
. 5
. 5
. 8
. 12
. 12
. 13

CHAPTER 2 KINETIC GLOBAL MODELING FRAMEWORK . . . . . . .
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1.1 Overview of available software for global modeling of plasmas
2.1.2 Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Global Model Equations . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 Continuity and Energy Conservation Equations . . . . . . . .
2.2.2 Reaction Rates . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2.1 Rates as Symbolic Expressions . . . . . . . . . . .
2.2.2.2 Reaction Rates Computed from Cross Section Data .
2.2.3 Additional Physics Effects . . . . . . . . . . . . . . . . . . .
2.2.3.1 Pseudo-spatial Effects . . . . . . . . . . . . . . . .
2.2.3.2 Non-standard Processes . . . . . . . . . . . . . . .
2.3 Process of Operation . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Data Acquisition and local database . . . . . . . . . . . . . .
2.3.2 Model Source Generation . . . . . . . . . . . . . . . . . . . .
2.3.3 Preparing for Evaluation . . . . . . . . . . . . . . . . . . . .
2.3.4 Evaluation and Usage . . . . . . . . . . . . . . . . . . . . . .
2.3.4.1 Interactive Usage . . . . . . . . . . . . . . . . . . .
2.3.4.2 Batch Operation . . . . . . . . . . . . . . . . . . .
2.4 Example Usage with Verification . . . . . . . . . . . . . . . . . . . .
2.4.1 Plasma Verification . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 Gas Phase Verification . . . . . . . . . . . . . . . . . . . . .
2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

15
15
16
17
18
19
21
21
22
24
24
25
25
26
27
29
31
32
32
33
33
35
39
40

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

41
41
42
43

CHAPTER 3 PLASMA-ASSISTED COMBUSTION WITH KGMF
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Classical Combustion . . . . . . . . . . . . . . . . . . . . .
3.2.1 Chain Reactions . . . . . . . . . . . . . . . . . . . .
vii

.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

44
46
49
52
54
56
56
58
59
60
62
63
65
66
70
72

CHAPTER 4 RARE GAS METASTABLE LASERS WITH KGMF
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 General Laser History and Concepts . . . . . . . . . . . . .
4.2.1 Laser Background . . . . . . . . . . . . . . . . . . .
4.2.2 Laser Fundamentals . . . . . . . . . . . . . . . . . .
4.3 Motivation of this Work . . . . . . . . . . . . . . . . . . . .
4.3.1 DPAL . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 OPRGL . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Derivation of Lasing Model . . . . . . . . . . . . . . . . . .
4.4.1 Lasing Model . . . . . . . . . . . . . . . . . . . . .
4.4.2 Pumping Model . . . . . . . . . . . . . . . . . . . .
4.5 Modeling Results . . . . . . . . . . . . . . . . . . . . . . .
4.5.1 Pulsed Discharge . . . . . . . . . . . . . . . . . . .
4.5.2 Continuous Discharge . . . . . . . . . . . . . . . .
4.5.3 EEDF driven . . . . . . . . . . . . . . . . . . . . .
4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

75
75
76
76
77
81
81
83
89
91
93
94
95
97
98
99

CHAPTER 5 CONCLUDING REMARKS AND FUTURE WORK
5.1 Future KGMf . . . . . . . . . . . . . . . . . . . . . . . . .
5.1.1 Self-consistent EEDF . . . . . . . . . . . . . . . . .
5.1.2 Uncertainty Quantification . . . . . . . . . . . . . .
5.1.3 Spatial - Fluid Code . . . . . . . . . . . . . . . . . .
5.2 Future PAC/RGL . . . . . . . . . . . . . . . . . . . . . . .
5.2.1 PAC . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2.2 RGL . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

102
103
103
104
104
105
105
106

3.3

3.4

3.5

3.6

3.2.2 Explosion Limits . . . . . . . . . . . .
Plasma Assisted Combustion . . . . . . . . . .
3.3.1 PAC Fundamentals . . . . . . . . . . .
3.3.2 Physical system for study of kinetics . .
3.3.3 Nanosecond Discharge System . . . . .
KGMf overview with PAC . . . . . . . . . . .
3.4.1 KGMf analytic additions . . . . . . . .
3.4.2 Mechanisms . . . . . . . . . . . . . . .
3.4.2.1 Argon buffer gas . . . . . . .
3.4.2.2 H2 - O2 . . . . . . . . . . . .
3.4.2.3 Adding N2 . . . . . . . . . .
Modeling . . . . . . . . . . . . . . . . . . . .
3.5.1 Rate Variation From EEDF dependence
3.5.2 O2 (a1 ∆g ) Injection and Ignition Time .
3.5.3 Nanosecond discharge . . . . . . . . .
Conclusion . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
APPENDIX A
Appendix to KGMf . . . . . . . . . . . . . . . . . . . . . . . . . . 109
viii

APPENDIX B

Appendix to PAC . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

ix

LIST OF TABLES

Table A.1 Species and reaction counts for a range of chemistry models demonstrating a
huge increase in complexity with included species. . . . . . . . . . . . . . . . . 114
Table A.2 GCC-5.4 Compile and evaluation timing. Runs marked with X failed to compile due to a lack of memory on a machine with 64GB of RAM available. . . . . 115
Table A.3 Clang-3.8 Compile and evaluation timing. Runs marked with X failed due to
a segmentation fault during function loading. . . . . . . . . . . . . . . . . . . . 116

x

LIST OF FIGURES

Figure 1.1

Types of plasma relative to electron density and temperature [118]. . . . . . . .

2

Figure 2.1

Cross section data for Ar excitation from the LXCat IST Lisbon [71] database.
Horizontal lines are species energy levels, vertical lines are electron excitation reactions, dashed lines are radiative reactions, and red boxes highlight
incomplete orbitals or orbitals constructed from grouped excitation. . . . . . . . 28

Figure 2.2

Single data mapping for Ar with Biagi data . . . . . . . . . . . . . . . . . . . . 28

Figure 2.3

Model generation for Ar with Biagi database [16]. . . . . . . . . . . . . . . . . 29

Figure 2.4

KGMf flow chart of the latter three user stages. Standard usage involves
modifying control text files for extra reactions and system parameters, while
every section is accessible for modification in a Python instance. Red dashed
lines mark the future work of incorporating a Boltzmann equation solver. . . . . 31

Figure 2.5

Basic KGMf integrator usage. . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

Figure 2.6

Transient square wave behavior for a period of τ = 100µs, duty ratio of 25%,
and 500W of total absorbed power from: the original paper (x markers), faithful model with guessed unreported parameters (dotted) and tuned parameters
(solid), and recreation with cross section dependent rates (dashed). . . . . . . . 35

Figure 2.7

Argon ionization rates defined with: Arrhenius rate used by Ashida et al.
(solid black) and KGMf spline implementations for convolutions of an assumed Maxwellian EEDF with cross sections from: the LXcat Biagi database
(solid blue) and estimates of the cross section extracted from a figure in the
Arrhenius rate source [127]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

Figure 2.8

Comparison between results from different assumed EEDFs. With duty ratios of 10% ≤ α ≤ 50%, the faithful model is plotted in black while the cross
section dependent model for a Maxwellian (x = 1) and Druyvestyn (x = 2)
EEDFs are in blue and green respectively. . . . . . . . . . . . . . . . . . . . . 36

Figure 2.9

Example output a CH4 -O2 -Ar mechanism with 90% argon and a stoichiometric fuel/oxidizer mixture initiated at a pressure of 1 atm and temperature
of 1400 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

Figure 2.10 Density-scaled ignition delay times for a mixture of CH4 -O2 -Ar combustion
with ζ = 1/3. Experimental data points were extracted from figure 5 in [29],
the black line is the reported fitting function, and the blue line represents
scaled estimates from the KGMf recreation model. . . . . . . . . . . . . . . . . 38

xi

Figure 3.1

Experimental explosion limits of stoichiometric H2 -O2 combustion in a spherical KCl-coated vessel from [30] with original data credited to [103]. Dotted
line corresponds to where 2k2 = k5 nM . . . . . . . . . . . . . . . . . . . . . . . 47

Figure 3.2

Schematic of plasma assisted combustion enhancement pathways taken from [80].
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

Figure 3.3

Comparison of chain propagation reactions with electron impact dissociation
rates from [149] using an EEDF found using a two term Boltzmann solver
assuming different field strengths, E/N, and a fixed gas composition. . . . . . . 50

Figure 3.4

Experimental extension of explosion limits of H2 -O2 combustion under electrical discharge from [149] with original reference to [56]. . . . . . . . . . . . . 53

Figure 3.5

Types of discharges with regards to electron density and electric field values,
taken from [152]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

Figure 3.6

Experimental demonstration of the bypassing of extinction limits. Plotted
against fuel mass fraction, XF , with different cases of oxidizer mass fraction, XO , left demonstrates the classical S-curve behavior, while the case
with higher oxidizer levels, right, shows monotonic behavior; taken from [159].

55

Figure 3.7

Electron impact cross sections for Argon from the Morgan LXcat database [119]. 60

Figure 3.8

Electron impact cross sections for Argon from the Morgan LXcat database [119]. 61

Figure 3.9

O generation from gas phase reactions. Reactions involving M collision
species are scaled assuming 100 Torr for matching units. . . . . . . . . . . . . 62

Figure 3.10 Plot of EEDF, eq. 2.9, with different values of X and values of Te. . . . . . . . . 64
Figure 3.11 Effect of EEDF selection on different electron impact rates: excitation of
molecular oxygen on the left, O2 → O2 (a1 ∆g ), and dissociation through
quenching on the right, O2 (a1 ∆g ) → O + O(1 D). . . . . . . . . . . . . . . . . . 65
Figure 3.12 Effect of EEDF selection on different electron impact rates, shown through
relative error from rates calculated with a Maxwell-Boltzmann distribution
(X=1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
Figure 3.13 Ignition time in 10 Torr of H2 :O2 ::2:1 with equilibrium and non-equilibrium
EEDFs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
Figure 3.14 Rates of radical producing reactions with (a) and without (b) singlet delta
oxygen supplementation. Without external power, the rates appear quite similar but differ in scale. The total system number density is used for the M
collision species. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

xii

Figure 3.15 Changes to ignition time due to injection of 6% O2 (a1 ∆g ) at 10 Torr. . . . . . . 69
Figure 3.16 Effect of O2 (a1 ∆g ) injection on the rates of reactions involved. . . . . . . . . . 70
Figure 3.17 Baseline results at 100 Torr with different EEDFs. . . . . . . . . . . . . . . . . 71
Figure 3.18 Pulsed power 60 kW 100 Torr. Electron impact reaction terms jump in magnitude due to temporary increases of Te during pulses. . . . . . . . . . . . . . . 71
Figure 3.19 Pulsed Ignition time at 100 Torr with a 60 kW 25 ns pulse. Contains modeling results for stoichiometric combustion, φ = 1, with varied EEDF, and
examples of lean and rich behavior. . . . . . . . . . . . . . . . . . . . . . . . 73
Figure 3.20 Species densities in 1000K NSD ignition. . . . . . . . . . . . . . . . . . . . . 73
Figure 3.21 Pulsed 1000 K Temperatures. . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
Figure 4.1

Two and three level laser schemes. . . . . . . . . . . . . . . . . . . . . . . . . 78

Figure 4.2

Electronically excited levels of a 3-level DPAL system with ground state
n2 S1/2 and further excited levels n2 P1/2 and n2 P3/2 [94] . . . . . . . . . . . . 82

Figure 4.3

Excited level diagram for an argon oRGL system. . . . . . . . . . . . . . . . . 84

Figure 4.4

Experimental time-resolved traces of 5p[1/2]1 ↔ 5s[3/2]2 gain in Kr, with
colors denoting total pressures of 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6,
1.8, and 2 bar [64]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

Figure 4.5

Experimental time-resolved traces of 5p[1/2]1 ↔ 5s[3/2]2 intensity in Ar,
with colors denoting total discharge voltages of covering [-1 kV,-2 kV] in
steps of 100V [63]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

Figure 4.6

Optical gain along microplasma vs pumping intensity; 2% Ar - 98% He at
769 Torr with flow of 3.7 mmoles/s and total microwave power of 9W, data
adapted from [135]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

Figure 4.7

Assumed two-way passing schemes for induced lasing emission, adapted
from [62]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

Figure 4.8

Excited levels of argon considered in the model. Levels involved in the laser
process are in red. Data comes from the LXcat Biagi database [16]. . . . . . . . 95

Figure 4.9

Optimized delay time between discharge and optical pulses depends on species
equilibration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

xiii

Figure 4.10 Effect of adding helium to the system over a range of ratios. The line type
in the right plot transitions from solid to dotted with the same labels as the
colors in the left plot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
Figure 4.11 Species densities and involved intensities as a result of repetitive optical pumping. 98
Figure 4.12 Equilibration of species and intensities with 50% helium buffer gas. . . . . . . . 98
Figure 4.13 Electron impact cross sections involved in the laser process from LXcat Biagi [16]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
Figure 4.14 Same cross sections fig. 4.13, zoomed in to look at relevant lasing states. . . . . 100
Figure A.1 Example simulation input file for an argon plasma. . . . . . . . . . . . . . . . . 112
Figure A.2 Example reaction input file for an argon plasma. . . . . . . . . . . . . . . . . . 113

xiv

CHAPTER 1
INTRODUCTION

Plasma, described as the fourth state of matter or a charge-balanced gas of ions and electrons,
makes up the majority of the known visible universe. This profusion of matter has a broad range
of defining characteristics due to the combined effects of the underlying mechanical, electrical and
quantum forces and atomic and molecular collisions. Existing both naturally (e.g. astrophysical
plasma) and in man-made forms, plasmas have been a focus of study since before the term was
coined by Irving Langmuir in 1928 [98]. The complex nature of the driving forces limits analytical
solutions to most phenomena of interest, which restricted this field to experimental research until
the advent of computational methods of numerical simulation. The first numerical studies were
performed using simplified models by the plasma simulation pioneers in the 1950s [162, 23] and
60s [17, 37]. The evolution of computational hardware since has led to an exponential growth in
the field, resulting in a wide range of applications and simulation methods.
Though the forces involved in plasma are relatively well defined at a fundamental level, the
sheer number of particles leads to an impossibility of simulating all forces involved in an ab initio
manner. Each method of simulation is reliant on a certain set of imposed assumptions in order
to balance accuracy and computational evaluation time. The wide range in the types of plasmas
results in overlaps between physical behavior and assumptions built into computational methods.
These overlaps are used to choose simulation methods given a specific system of interest and,
with advancing computational hardware, motivated the development of hybrid methods combining beneficial aspects. The majority of development has focused on understanding the physical
dynamics of plasma systems, commonly resulting in simplifying assumptions placed on the underlying chemistry involved in gas-plasma interactions. To be fair, modeling this chemistry is
still unfortunately restricted by the availability of experimental and numerical fundamental data as
opposed to operating on fundamental principles.
This chapter provides an overview of plasma fundamentals and the main types of simulation

1

Figure 1.1 Types of plasma relative to electron density and temperature [118].
methods, followed by motivations to develop a plasma-chemistry simulation platform with features
differing from available tools. The key capabilities and limitations of each method are described
along with examples of available software implementations. With the realization of plasma applicability to a wide range of scenarios, commercial and academic software exist for each method,
designed to tailor towards specific regimes. Emerging interest in regimes requiring integration of
involved plasma-chemistry motivates a simulation platform with a wide acceptance of data and
an intrinsic flexibility to quickly incorporate different physics. If designed in a manner to rapidly
evaluate different parameters, a quick chemistry engine could be included in a hybrid modeling
method for higher fidelity simulation.

2

1.1

Plasma Fundamentals

Langmuir described plasma as the “region containing balanced charges of ions and electrons,”
found in between the electrodes of an RF discharge encased by “sheaths” containing predominantly
ions [98]. This description has been expanded to define plasma as a combination of free electrons
with ionized and excited molecular, or atomic, species interacting with each other and external
forces [19]. The set of charged species in a given plasma exhibits collective behavior due to longrange electrical forces and yet remains macroscopically neutral in a large number of cases. In the
same manner that a solid or liquid undergoes a phase transition with the introduction of enough
thermal kinetic energy, a plasma is formed through a system absorbing enough kinetic energy to
exceed the binding energies of both gaseous molecules and orbital-bound electrons. The kinetic
energy transfered by inelastic collisions results in the production of electrons and ions along with
rotational, vibrational, and electrical excitations of atomic and molecular species. As the reaction
dynamics can differ widely between ground and excited species, the introduction of electron impact
processes can dramatically alter chemical kinetics with high electron energies and/or densities.
Aside from the plethora of possible excited species depending on the reaction network of a
given system, plasmas can be further classified by their ionization fraction, whether or not they
exhibit thermal behavior, and if the electrons are in equilibrium with their own ensemble. With
the ionization fraction defined as the ratio of heavy particles that are ionized, the classification of
thermal vs. non-thermal plasmas corresponds to the difference between electron and heavy species
temperatures. Thermal plasmas are characterized as having approximately equal electron, heavy
excited/ionized and ground species temperatures, Te

Ti

Tg , and obeying the laws of thermody-

namic equilibrium. In these cases, gas temperatures have increased to the order of 1 eV ≈ 10,000
K or greater from joule heating through electron-neutral collisions [50]. Lightning is a natural
example of a thermal spark discharge, with local temperatures increasing to above 24,000 K [165].
Non-thermal plasma are characterized by electron and excited species temperatures differing from
the ground state temperatures, Te

Ti

Tg , resulting in a dominance of excited species and elec-

3

tron impact reactions as opposed to equilibrium gas-phase reactions. The expression for energy
exchange under elastic collisions between two gas species is written as, [6],
dU
m1 m2
k(T1 − T2 )n2 < (v1 − v2 )σm >
= −3n1
dt
(m1 + m2 )2

(1.1)

with cross section σ , <> denoting averaging with the velocity distribution, and mass m, density n,
temperature T and velocity v for each species. The mass difference between electrons and heavy
species results in an inefficient momentum transfer to heavy species, meaning a large portion of
electron energy resides in excited states and electrons do not thermalize well with heavy species.
An increasingly common example is a medical atmospheric plasma jets, where, in the discharge
region, gas temperatures can approach 2000 K and electron temperatures are estimated to be in excess of 10 eV [88]. Paired with a difference in temperature, electron velocity distributions may not
follow a Maxwell-Boltzmann, or equilibrium, distribution. This is primarily driven by the lack of
thermodynamic equilibrium between forward and reverse inelastic collisions [54]. Consider a gas
discharge system where fast electrons are generated in order to ionize and excite the gas in the bulk
volume and the primary loss of electrons and radiation are to the wall. With reasonable gas pressures, the electron-electron collision rate is too small to thermalize low and high-energy electrons.
Though electron-electron collisions support a Maxwellian distribution in general, a comparably
high inelastic collision frequency will perpetuate a non-equilibrium, or kinetic, electron energy
distribution.
The utilization of plasma has spread across a wide range of technologies; to name a few:
medical applications of sterilization and tissue regrowth [174], semiconductor manufacturing [96],
hypersonic jet engines [74], thrusters [25], fusion [40, 51] and material modification [68]. Alterations in pressure, temperature, gas mixtures, and the type of electrical discharge all contribute to
physical behavior of the plasma due to the underlying dominant forces, fig. 1.1. Selection of the
best simulation method for a numerical study depends on the primary physics of interest and the
computational resources and time available.

4

1.2

Plasma Simulation

A complete numerical description of a plasma requires resolving a number of underlying physical
forces and reaction processes. As a system of neutral and charged species, the particle motion
must be represented in a manner which responds to electromagnetic and thermal effects. This
requires a solution of Maxwell’s or Poisson’s equations for the forces acting on species described
as either particles or fluids. A description of the different species’ kinetic and potential energies
then dictates collisional processes when paired with reaction probabilities. Reaction probabilities
are used as the addition of quantum processes for electron impact reactions would make simulation
intractable regardless of the method used. Simulation methods are categorized into one of two
methods: kinetic or fluid descriptions. However, interest in different physics, while factoring in
computational run time, has also lead to the development of hybrid and reduced methods.

1.2.1

Kinetic

A kinetic description of the electrons and heavy particle motion inherently implies moving individual, or clustered, particles around a domain as a response to electromagnetic, collisional, and
boundary forces. Using a predominantly fundamental description retains very good physics; numerical solutions to the coupled Newton-Lorentz and Maxwell’s equations (or Poisson’s equation
in the electrostatic case) are found with very few approximations [86], resulting in a self-consistent
evaluation of the particle velocity distributions, f (x, v,t). This high-fidelity method comes with
a cost of computational intensity with regards to the total number of simulated particles. Due to
memory limitations, it is not possible to set up a 1:1 relationship between plasma and computer
particles for most plasma pressure regimes. To resolve this issue, the majority of current kinetic
methods applied to plasmas rely on the particle-in-cell (PIC) scheme. Formalized in the 1970s by
Birdsall and Langdon [18] and Hockney and Eastwood [43], the PIC scheme takes advantage of
the collective behavior existing between charged particle species to simulate macro-particles as a
statistical representation of 105−7 real particles.

5

The typical solution scheme in PIC is reliant on the particles existing in continuum velocity
and position space. With the macro-particles initialized to some set of positions and velocities,
fields from both external forces and charged particles are found at discrete locations in space, at
the nodes of some imposed grid, and then interpolated into collective forces acting on simulated
macro-particles. Since both fields and particle data are known at discrete time intervals, their
respective solutions are temporally separated by half of a time step ∆t/2 allowing for leapfrog,
center-difference, integration of the equations of motion [169],
ut+∆t/2 − ut−∆t/2
q
ut+∆t/2 + ut−∆t/2
=
Et +
× Bt
∆t
m
2γ t
ut+∆t/2
xt+∆t − xt
=
∆t
γ t+∆t/2

(1.2)

for each macro-particle. Once the particles have been moved to their next configuration in time,
boundary conditions and collisional process schemes are used to transfer energy and mass into
different species states. The new set of particles at distinct positions and velocities then starts the
cycle again until the total time of interest has been resolved. Selection of the weighting scheme
used to transmit forces back to particles requires a balance of computational efficiency, physical
accuracy, and simulation noise in the system [28] with the latter often suppressed actively or passively. Although short-range Coulomb forces between macro-particles in the same grid cell are
underestimated [169], the simulations are fully kinetic and and take significantly less time than a
numerical evaluation of direct Coulomb’s law between all particles as the operations performed are
O(N log N) and O(N 2 ) respectively. Collisions are typically handled using a Monte Carlo collision
scheme which determines particle interaction results from statistical probabilities. Further speeding up simulation time, a fictitious null-collision is used such that the total collision frequency of
each species is independent of incident particle energy and uniform in space. This allows for collisional processes to be calculated on a random subset of macro-particles as opposed to performing
multiple operations on each one. Considerations can also be made to fix certain species such that
their full dynamics are ignored; heavy species can be fixed to the background accounting for relative immobility compared to electron motion, or, a background fluid of electrons can be assumed
6

when ion and heavy particle behavior is of interest.
Though representing the highest-fidelity modeling of plasmas, PIC methods are not without
stability requirements and evaluation limitations. The discrete nature of particle kinetics brings
√ −1
numerical fluctuations which converge proportionally to N for N macro-particles, requiring
large numbers of macro-particles for a given simulations. Grid spacing ∆x, is limited to being
smaller than the smallest characteristic length scale of the system to be resolved (e.g. Debye
length, sheath length, or Larmor radius). Time steps ∆t, must be chosen to satisfy a number of
conditions: resolving the shortest characteristic time scale (e.g. inverses of plasma frequency,
driving frequencies, or maximum collision frequency), satisfying the particle Courant condition,
v∆t < ∆x, with particle velocity, v, and restricting particles to a single collision per time step . The
major limitation of PIC simulations lies in the computational time required for large simulations,
controlled predominantly by the number of macro particles and selection of ∆t and ∆x. Attempts
to reduce the computational intensity relative to a a fully-kinetic molecular dynamics simulation
induces numerical issues. Noise in the system from finite-particle numerical fluctuations is only
reduced through increasing particle count, higher-order algorithms, or filtering high-frequency effects in plasma, the latter of which can also filter out certain physics of interest. The noise in the
system further causes artificial cumulative heating which can dramatically effect the behavior of
macro-particles confined for long periods of time. Finally, the number of particles must be large
enough to adequately describe a statistical representation of particle velocity distributions which
can range in several orders of relative magnitude.
In the context of plasma simulation, PIC is commonly used when high-fidelity modeling is
required or the system is either at low-pressure or demonstrating non-equilibrium behavior. Lowpressure systems have the advantage of needing fewer particles to adequately describe the kinetics. The self-consistent description of the particle distribution is capable of demonstrating nonequilibrium effects accurately. Finally, the minimal number of assumptions made results in PIC
methods being the most accurate description if particle count, grid size, and time steps are chosen
correctly. The largest constraint on PIC simulations is simply the computational resources and

7

clock time required for attaining a solution.
Several simulation packages exist for PIC simulation, ranging from open-source academic
codes to commercial applications. Some of the first examples came from the Plasma Theory Simulation Group (PTSG) (e.g. XES1, XPD[P,C][1,2] and XOOPIC) [18, 170, 166, 171], laying the
foundation for almost any PIC implementation. The desire for application-specific considerations
led to the development of a number of simulation platforms for the study of plasma physics, such
as: WARP for modeling particle beams with high-space charge intensities [58] and PICCANTE,
as an effort towards massively parallel fully-relativistic PIC [137] extending from ALaDyn [14].
There also exist a number of commercial options, such as TechX’s VSim (formally VORPAL) [33]
and Orbital ATK’s MAGIC Tool Suite [5], originally designed for microwave engineering applications.

1.2.2

Fluid

Reducing computational intensity relative to kinetic methods, a fluid description describes the evolution of macroscopic quantities coming from velocity moments of the Boltzmann equation. The
first three moments become the species continuity, momentum transfer and heat transfer equations. As with PIC, these equations are coupled with Maxwell’s equations (or Poisson’s equation)
for self-consistent electric and magnetic fields [86]. The Boltzmann equation is defined as
B=

∂
+ v · ∇ + a · ∇v f (x, v,t) =
∂t

∂f
∂t col

(1.3)

where f (x, v,t) is the particle distribution, (∂ f /∂t) is the collision term, a and v are particle accelerations and velocities, and ∇ and ∇v are del operators acting in position and velocity space
respectively [19]. The zeroth velocity moment, (B) d3 v, defines the continuity equation for each
species
∂n
Sm
+ ∇ · (nu) =
=
∂t
m

with species density n =

∂f
d3 v
(1.4)
∂t col
f (v) d3 v, average velocity u = 1n v f (v) d3 v, and collisional mass

source term Sm . In a similar fashion, the momentum and heat transport equations are found from
8

the first, m (vB)d 3 v, and second, 12 m (v2 B) d3 v, velocity moments respectively,
mn

∂ (u)
+ mn(u · ∇)u + ∇P − nq(E + u × B) = S p = m
∂t

1
γ −1

∂p
1
+ ∇ · (pu) + (P · ∇)u + ∇ · L = Sε = m
∂t
2

(v − u)

∂f
d3 v
∂t col

(1.5)

∂f
d3 v
(1.6)
∂t col
after substituting the equations into each other. Analogous to Sm , the collisional momentum trans(v − u)2

fer S p and collisional heat transfer Sε are found from velocity moments of the distribution time
derivative with respect to collisions. Completing the equation, the pressure tensor P and heat flux
L are defined as
P=m

1
L= m
2

(v − u)(v − u) f (v) d3 v

(v − u)(v − u)2 f (v) d3 v

(1.7)

with scalar pressure defined through pδi j = Pi j and γ is the ratio of specific heats. Energy is related
to other macroscopic quantities through 2ε = mnu2 + 2/(γ − 1)p. It should be noted that the use
of macroscopic quantities averaged over the distribution implicitly assumes local thermodynamic
equilibrium [139].
Unlike the kinetic formulation, the general fluid equations are not closed and the distribution is
not found self-consistently. Demonstrating the closure problem, determining L requires the third
velocity moment of the Boltzmann equation. This problem lies in the inability to describe both
momentum P and heat L fluxes in general and only in terms of macroscopic quantities; instead
they require a prescription dependent on the physical regime and effects of interest [38]. Though
higher-order closure methods exist, and are required for descriptions of turbulence [109], solution
methods commonly rely on the first two or three moments, a declaration of distribution function,
and approximative methods. Closure relations come from taking a hydrodynamic limit of the
Boltzmann equation, commonly in the form of the Navier-Stokes equation by using Fourier’s law
for treating heat flux. These approximations and imposed distribution result in a loss of information
relative to kinetic methods.
Depending on the regime of interest, approximations to terms in eqs. 1.4, 1.5 and 1.6 can drastically simplify simulation evaluation. The selection of which approximations, and hence forces
9

to exclude, directly relates to both the accuracy and computational intensity of the simulation. As
with PIC methods, species do not need to be declared in the same manner, neglecting small terms
relative to differences in energy, mass, and reactive frequencies. A good set of relevant approximations can result in high-accuracy solutions at a considerably lower computational cost than
kinetic methods. Approximations range from treatments of dimensionally-coupled force terms to
reduction of the total number of PDEs required for simulating the system. Depending on the regime
considered, the dynamics of some species can be neglected altogether. While the general equations
above apply to all species, a two-fluid plasma model can be made when only considering electrons
and ions [140]. This assumption forms the basis for the magnetohydrodynamic formulation [48].
Except in certain plasma regimes, such as for strongly magnetized systems such as in tokamaks [124], shear and viscous forces can be neglected [19]. The pressure tensor, P, is then reduced to purely diagonal elements representing an anisotropy in the random velocities of species.
Further neglecting any dependence on thermal conductivity, this combination results in an ideal
compressible Eulerian system for closure relations. If the momentum transfer collision frequency
is larger than the electric driving frequency, the drift-diffusion approximation can be made. This
requires that inertia of the particles are neglected and assumes they respond instantaneously to
external fields. Reducing the time derivative of the distribution to a relaxation time from an equilibrium, (∂ f /∂t) ≈ ( feq − f )/τ, introducing species mobility, µ = qτ/m, neglecting the magnetic
kB T
field, and replacing the square velocity average with the equilibrium value, n1 v2 f (v) d3 v = e m ,

results in the simplified drift-diffusion formulation of eqs. 1.4 and 1.5 such that species and flux
are defined as
∂ nk
Sm
+ ∇ · Γk =
∂t
m

(1.8)

Γk = µk nk E − ∇(Dk nk )

(1.9)

for each species nk , where the diffusion coefficient, Dk = µk kB T /q, comes from the Einstein relation [55]. These equations can provide analytic solutions and insight into complicated behavior
when considered in reduced dimensions (e.g. ambipolar-diffusion from the drift diffusion formu-

10

lation). Though these derivations do not assume an equilibrium distribution, the use of equilibrium
velocity limits validity to small perturbations of the distribution from equilibrium.
The limits of fluid methods come from the approximations used to close the governing equations along with the methods used to solve them. Solutions to these PDEs are commonly found
through the use of finite difference, finite element, or finite volume methods. Finite difference
methods rely on representing derivatives through truncated Taylor expansions and allow the simplest implementation [67]. A variety of implicit (e.g. Crank-Nicholson) and explicit (e.g. RungeKutta) numerical methods can be used, each carrying their own limitations/restrictions [172]. Solving a major limitation of finite difference methods, finite element and volume methods are capable
of complicated geometries from using unstructured grids. Finite element relies on the assumption
of an approximate solution to the discretized problem with the method looking for a solution to the
integral form of the PDEs through residual calculations [145]. Finite volume methods instead discretizes the integral forms of the equations, which becomes advantageous when the conservation
laws can only be described in integral form [67]. While advantageous in situations with discontinuities, finite volume methods run into trouble when describing derivatives, limiting use to when
viscosity terms can be neglected.
The wide range of approximations that can be made results in a surplus of fluid simulation
packages for plasma applications. The general applicability of computational fluid dynamics has
resulted in a number of general platforms which have been adapted to include plasma physics.
Though a good number of codes remain in-house, a growing number have become available as
open-source packages: plasmaFOAM [168] built from the openFOAM fluid dynamics framework
for low temperature plasmas, BOUT++ [42] designed for curvilinear systems using computer architecture parallelization and ITM-TF [89], built on the open workflow system Kepler [36] and intended for analysis of ITER. Commercial applications include the COMSOL Plasma Toolkit [70],
Tech-X’s USim [32], and ESI’s CFD-ACE+ [59].

11

1.2.3

Hybrid

Expanding upon both methods above, hybrid simulation methods are made by combining kinetic
and fluid methods. The intention is to get methods which have both accurate particle kinetics and
the evaluation speed of fluid methods. As with fluid methods, the physics of interest dictates the
pairing of particle species with solution method, with completely different overall descriptions.
Electrons can be considered using PIC while considering a background ion fluid [146]. By comparison, studies on the ion energy distribution can treat ions with a PIC scheme and electrons as a
background Maxwellian fluid [131]. Aside from the particle and fluid method limitations of each
hybrid component, communication between the two methods must be given consideration, as reducing the time step size to the smallest requirement results in excessively long simulation times.
The built in assumptions of hybrid models limits the ability for general-application software, however a number of codes are available such as Pegasus [95], designed for multi-fluid quasi-neutral
astronomical plasmas, HPEM [97], designed for low-pressure plasma processing reactors, and
HPHall-2 [125], designed for Hall-effect thrusters.

1.2.4

Global/Chemistry

Volume-averaged global models represent the simplest forms of the plasma fluid description. Neglecting spatial variation requires a sizable number of approximation models and physical data
to account for relevant physics while greatly reducing the computational intensity of evaluation.
The simple description allows for emphasis on the scalar mass and energy flux terms, Sm and Sε
respectively. In order to speed up PIC and fluid methods, limited reaction chemistry is typically
considered. This is not a problem when the dynamics of particle motion are the main interest,
however applications in higher pressure gas mixtures depend on a non-linear collection of reactions which are not easily simplified a priori. With a minimal increase of computational intensity,
continuity equations can be written for each considered species, coupled with one or more energy
equations, to describe a complete picture of the plasma-chemistry interaction. Due to their inherent

12

simplicity relative to other methods, global models are used for: fundamental studies of collisional
effects on species populations, determining dominant reaction pathways to use in higher-fidelity
solution methods, and analysis of reaction data sensitivity.
It is the emphasis on plasma-chemistry which makes global models ideal for exploratory analysis of low-temperature, high-pressure, nearly-homogeneous plasma systems. In these regimes, reactions between heavy species and electron-impact reactions both play an important role in system
evolution. The number of distinct species required limits the applicability of a multi-fluid approach.
Similarly, the number of reactions required for a complete picture would dramatically slow down
kinetic methods. Not without their own caveats, global models require a large amount of reaction
chemistry data paired with declaration of the energy distribution function. Along with continuity
equations, energy equations can be derived for both electrons and heavy species. The energy distribution of heavy particles is commonly assumed to be in equilibrium, while the implementation
of the electron energy distribution tends to separate global model platforms. Explicit selection of
the electron energy distribution allows for algebraic electron impact rate coefficients, allowing for
plasma global models using the chemical kinetics package CHEMKIN-III [85]. Codes such as ZDPlasKin [123] and GlobalKin [41] represent some of the higher-fidelity global modeling softwares,
coupling a Boltzmann equation solver for a self-consistent description of the electron distribution
along with databases for electron impact reactions. Larger packages, such as PLASIMO [76] and
COMSOL [70], are able to generate global models by reducing higher-dimensional descriptions.
In many cases, the data used is restricted to pre-confirmed sets compiled by the code developers.

1.3

Kinetic Global Modeling Framework Motivation

In the context of modeling systems with complicated plasma chemistry, global models are the
best platform for developing an understanding of the species and reaction chain details of the
underlying chemistry. The simplicity of the implementation further allows for the addition of new
physics modules with relative ease. These modules are combined with the set of approximate

13

spatial contributions which close the system of equations. The variable nature of approximations
leads to a requirement of understanding implemented terms and the effect they impose on final
solutions. What was found to be lacking in available software was the ability to modify reaction
data, analytic approximation terms, and electron energy distributions along with easy avenues for
dissecting contributing terms.
The data dependence of global models motivates the study of reaction data validity and system
sensitivity. Multiple databases exist with different data describing the same process, differing by
acquisition method or regime considered. The speed of global models should allow for quick
replacement of data to quickly compare relative effects.
Analytic approximations have the capability of drastically altering simulation results. The
assumptions built into these approximations tend to be extremely regime specific, requiring modification or replacement under different scenarios. Extending to the context of academic research,
the incorporated models should be made transparent through open-source code. Even with the
correct selection of approximations, the inclusion of new physics through additional terms should
be made possible with minimal code alteration.
Finally, treatment of the electron energy distribution can dramatically change the underlying
chemistry. While examples exist of global model software coupled to Boltzmann equation solvers,
the approximations made to solve for the distribution are not always valid in considered regimes.
For example, the common two-term approximation for the Boltzmann equation breaks down in
systems with a large collisional mass flux (e.g. combustion). Even without a self-consistent description of the EEDF, sensitivity to energy distributions should be regularly considered, especially
in systems with non-equilibrium chemistry.
Bringing all of these considerations to fruition required the development of a new generalpurpose global model platform with an emphasis on variation of assumptions and interactive usage.

14

CHAPTER 2
KINETIC GLOBAL MODELING FRAMEWORK

2.1

Introduction

In the last few decades, advances in computer architecture and scaling have drastically pushed
the boundaries of what can be simulated of known physics. The study of plasma in particular
requires the combination of electromagnetics, fluid dynamics, atomic and molecular collisions,
and statistical mechanics in order to begin describing a physical model. There are three main
components to a computational description of a plasma system: the environment and external
forces, dynamics of the particles in response to their surroundings, and composition of the particle
ensemble. Due to the high number of particles found even in a low pressure system, plasma
simulations are performed using some type of approximating method to reduce the computational
complexity as direct solve molecular dynamics simulations are usually not feasible. Fluid and
Particle-in-Cell (PIC) methods remain the most common, for macroscopic plasmas, by simplifying
the dynamics of individual particles into ‘groups’ while maintaining the complexity of transport
and spatial variation.
Taking the extreme of low spatial resolution, global (volume-averaged) models reduce all spatial dependence to analytic approximations and provide a focused platform for studying reaction
kinetics. Though the lack of spatial variations can lead to inaccurate results with regard to parameters of a specific simulated system, the intention is to study physical trends across variations of
inputs. With the greatly reduced computational intensity of these simulations, global models can
be used to map out a system’s parameter space either in an exploratory manner or retroactively
matching experimental results by iterating through possible input changes. An example of the latter is a “collisional radiative model” [141], wherein parameters of the global model are adjusted
until a simulated spectrum matches that of an experiment. In the context of interest in high fidelity simulations, such as PIC, global models can be used to reduce the included reaction set to a

15

minimized configuration maintaining relevant behavior.
The reaction kinetics of a given system describe the overall changes in particle densities and
temperatures as a result of various forms of energy transfer. In the context of classical chemistry,
this typically involves decomposition, displacement, and synthesis reactions involving molecular and atomic species. Plasma chemistry then expands the number of tracked species to include
excited species created through processes such as electron impact ionization and excitation. Reactions involving electrons are herein referred to as ‘plasma phase’ whereas reactions between only
heavy species are referred to as ‘gas phase.’ Further effects, such as diffusion or charge species
absorption by sheaths, are treated through the use of approximate models. A reaction mechanism
is then created by compiling a collection of all relevant reactions and effects involving species of
interest, represented by a global model in the form of a coupled set of ordinary differential equations (ODE). While the simplest of models can include a handful of species and reactions, certain
systems require hundreds of species and thousands of reactions for an accurate description.

2.1.1

Overview of available software for global modeling of plasmas

Defined by a set of species continuity and energy conservation equations, required data and constructed global models are conceptually similar across most applications [78]. Due to the relative simplicity of smaller reaction mechanisms, it is common practice for custom code to be
developed for each application; this results in a duplication of efforts which can negatively affect
quality control and code maintenance. There otherwise exists a few published softwares for generating and evaluating plasma-chemistry global models; research codes such as PLASIMO [76],
ZDPlasKin [123], and PLASMAKIN [128] and commercial options, such as COMSOL [70] and
Quantemol-P [132] (based on GLOBALKIN [41]).
Global models have a rich history of use for studying the chemical kinetics of applications
ranging from plasma-reactor etching/deposition systems employing different reactive discharges
(e.g. radio frequency [100], inductively coupled plasmas [161], high power impulse magnetron
sputtering [133]) to plasma assisted combustion [177] and Hall-effect thrusters [27]. With growing
16

interest in applying plasmas in biomedical, technological, and environmental applications, global
models present an excellent opportunity for better understanding complex plasma chemistry interactions. While used for direct simulation of systems, the relative computational simplicity, as
compared to spatially resolving codes, allows for quick sensitivity analysis [164] and experimentmatching iterative operation. Not only useful for research purposes, global models can be an
invaluable tool for teaching chemical kinetics and building intuition about reaction network behavior. With a generality of application, access must be provided to both examine sub-components of
underlying equations and to experiment with features and physics beyond the original design considerations. Closed-source codes, no matter how complex, do not offer enough flexibility due to an
inability to modify methods. The use of global models as both a research tool and an educational
platform calls for an open-source code written in an interactive language.

2.1.2

Design

It should be evident that global models obfuscate quickly as the number of species and reactions
increases. In order to handle these complexities, the Kinetic Global Modeling framework (KGMf)
was developed with two primary goals in mind: 1) minimize confusion when trying to gather and
compare data for a given reaction network and 2) provide a symbolic description of the global
model allowing for data exchange, tuning, and the ability to incorporate custom effects into user
modules. These two goals call for the creation of a local and self-checking database of reaction and
species data, to be paired with multi-level symbolic equation generation. Secondary considerations
during KGMf development were to: maintain as much operational flexibility as possible, allow for
on-the-fly user interaction, and avoid redundant calculations.
The framework is written in Python, using Numerical Python (NumPy) [167], Scientific Python
(SciPy) [79], Symbolic Python (SymPy) [117], Cython [12], and H5py [31]. The usage can be
broken up into four stages:
1. importing of data into the local KGMf database for reference

17

2. gathering the set of relevant reactions and processes
3. generating symbolic and numerical forms of the global model equations
4. evaluation of the global model and post-processing.
In principle the first stage is only done with initial installation or with the addition of new data
for the database. The second stage is initially performed automatically with regards to internal
databases, generating two humanly-readable input files for model setup (herein referred to as reaction and simulation input files) and allowing for user truncation of the rate mechanism, addition
of non-imported reaction data, selection of non-rate dependent options. The third stage is largely
automatic; input files define everything from physics assumptions and pre-processing steps to the
logic of the reaction network, allowing the KGMf to generate all symbolic equations and terms
within an interactive Python instance followed by compilation into a callable form from generated
C-code. The fourth stage involves integrating the system of ODEs in time, likely being performed
multiple times with variations on inputs (from integration time to data used), and post-processing
analysis by decomposing the ODE system into evolved components or collecting results from multiple integrations.
Taking advantage of the nature of Python, each user stage can be performed as a script or interactively via command line. While scripted usage has obvious benefits with regards to repetitive
tasks, interactive use exposes KGMf subroutines, data, and results for plotting, symbolic manipulation and extraction (post-processing), and user additions beyond those natively supported by the
KGMf. To bridge the gap between the two methods, each stage can be completed with a single
function call or split into sections for analysis or tweaking.

2.2

Global Model Equations

Global models are composed of a set of coupled non-linear ODEs representing species continuity
and energy conservation equations. Using the notation of the system state vector, Y , containing all

18

evolving species densities, temperatures, and additional terms for which a differential equation is
included, the global model system of ODEs, F, and corresponding Jacobian, J, are defined as
F(t,Y ) =

dY
dt

and

J(t,Y ) =

dF
dY

(2.1)

By using a symbolic representation of F and all functions therein, the Jacobian is found using the
standard derivative functions included in the SymPy library. Of increasing importance with the
stiffness of F, J greatly improves integrator performance and is found by default in the KGMf
routines.

2.2.1

Continuity and Energy Conservation Equations

From a given species set and reaction network, a governing set of equations can be defined to
describe the propagation of energy in a system. Evolution of each tracked species density, nk , is
described by a continuity equation of the form:
n˙ k =

dnk
= ∑ νik Ki ∏ n j + (G − L)k
dt
j
i

(2.2)

where i sums over all reactions including species nk , Ki is the reaction rate for a given process,
j creates a product of reactant species for reaction i, νik is the integer growth or loss factor for
species nk from reaction Ki , and (G−L)k represents gain and loss terms independent of the reaction
network. If relevant to given system, the imposition of quasi-neutrality allows for electrons to be
defined as the charged weighted sum of all ionized species, removing the need for an additional
ODE.
Depending on the system being modeled, energy conservation equations are implemented to
maintain physicality. From the assumption of thermodynamic equilibrium, all massive species
are assumed to have the same rotational and translational temperatures. The energy conservation
equation for this shared gas temperature is defined as:
dTg −Tg ∑k=1 Hˆk (Tg ) − 1 n˙ k + Q˙ + P˙proc
=
dt
ˆ
∑N
k=1 C pk (Tg ) − 1 nk
19

(2.3)

where Tg is the gas temperature in kelvin, k sums over all tracked species, Hˆk = Hk /RTg and
Cˆpk = C pk /R are the species specific enthalpy and specific heat in normalized NASA polynomial
format [24], Q˙ represents external heat flux, and P˙proc represents heat flux from thermodynamic
process (e.g adiabatic) energy transfer. It should be noted that n˙ k , in eqs. 2.2 and 2.3, represents
the full continuity equation for species nk , resulting in a large expression as the number of species
increases. This equation is derived in appendix A.1. Though appropriate for massive species,
electrons usually have a different temperature from the background gas and can be far from having
a Maxwellian distribution. Completing the energy picture is a description of the effective electron
temperature, Te , from the electron energy equation,
P
d 3
ne Te = abs − ne ∑ ni Ki ∆Ei j + Q˙ e
dt 2
V
i

(2.4)

defined as the sum of absorbed power, Pabs (in power density units with volume, V ), the kinetic
energy lost during electron impact reactions, and kinetic energy lost to the system from electron
absorption, Q˙ e . The second term is defined as the sum over all electron impact reactions of the
product between larger species’ potential energy flux, ∆Ei j , and average reaction frequency, from
the product of reaction rate, Ki , and involved reactant electron, ne , and large species, ni , densities.
In the case of momentum transfer reactions, a modified form of Ki ∆Ei j is constructed from reaction
cross section and involved electron energy. This equation is solved for dTe /dt using the product
rule, resulting in the electron continuity equation, dne /dt, appearing on the RHS as in eq. 2.3. The
inclusive repetition of ODE equations is used for function optimization prior to compilation.
It should be noted that terms such as Q˙ in eq. 2.3 and Pabs in eq. 2.4 are not restricted to a
constant value and can instead be defined as functions of time and/or the state vector Y . While
there are included methods to provide time dependent pulsing terms (square, gaussian, etc.) and
other known forms of external effects (e.g. microwave power absorption dependent on both electron density and electron-neutral collision frequencies [120]), any analytic expression can be used
through symbolic substitution.

20

2.2.2

Reaction Rates

By their construction, global models are entirely data driven with regards to the reaction rates making up the majority of the continuity equations. These rates represent a wide range of processes,
with defining data found through both empirical fitting and theoretical calculation, resulting in multiple forms of mathematical expressions. Data for these expressions ranges from constant values
and parameterization to numerical function convolution with an excited species energy distribution function. Relying on data availability, the KGMf was designed to keep as much flexibility and
variability as possible for both physical accuracy and integration error handling.
The range of both gas and plasma phase chemistry covers a large number of different uni-, bi-,
and termolecular processes. Gas phase chemistry is largely responsible for mass transfer between
species, either through synthesis or decomposition reactions, potentially affecting the gas temperature. Plasma phase reactions diversify gas species into ionized, electronic, and vibrationally
excited levels resulting in spontaneous emission and charged species accumulation. Though data
for reactions can be stored in a multitude of forms, they can all be thought of as either expressions
of system parameters/terms or cross section convolution with an energy distribution function.

2.2.2.1

Rates as Symbolic Expressions

The more directly handled rates are those that follow a specific analytic form with the ability to
differentiate symbolically. The majority of empirically found rates, both in gas and plasma phase,
are fitted with the modified Arrhenius temperature dependence equation:
β

Ki (Ta ) = Ai Ta i exp

−Ei
Rc Ta

(2.5)

where Ai , βi , and Ei are fitted numerical values, Ta is the temperature (either Tg or Te ), and Rc is
either the universal gas constant, R, or Boltzmann constant, kB , depending on the unit balancing
required. The simplest usage of eq. 2.5 is for constant rates, setting all but Ai to zero, while more
complex terms can be used depending on data availability. If changes in enthalpy, ∆H and entropy,

21

∆S, are known for involved ground state species, the reverse reaction rate, Kri , is defined as:

Kri =

∆Sl0

Ki
= Ki exp
Kcl

R

∆Hl0
−
RT

Patm
RT

−1
νk
∑N
k=1 i


(2.6)

using the change in the Gibbs free energy for the reaction, where Kcl is the equilibrium rate, R
is the universal gas constant, Patm is the pressure at one atmosphere, and ∆Sl0 and ∆Hl0 represent
changes in entropy and enthalpy respectively [52]. Recombination reactions are considered as a
collision of two reactants, the high-pressure limit k∞ , or as a three-body reaction with another
species participating as a collisional partner, the low-pressure limit k0 . The arbitrary pressure
reaction rate is defined using the Lindemann approach [105]:
K = k∞

Pr
F
1 + Pr

k [M]
with Pr = 0
k∞

(2.7)

where the Arrhenius rates for each limit are related to the reduced pressure Pr , and equivalently the
gas mixture density. In the Lindemann form, F is set to unity, however increased flexibility can be
achieved with the addition of Troe [163] and SRI International [156] forms of F as data dependent
parameterizations of gas temperature and reduced pressure.
To account for other parameterizations existing in the literature, the KGMf supports expression parsing from the reaction input file or the addition of symbolic rate expressions composed
of system variables prior to compilation. All expressions used within the KGMf retain their variable dependence in order to generate a system Jacobian and to be accessible for post-processing
analysis as functions.

2.2.2.2

Reaction Rates Computed from Cross Section Data

When cross section data is available for a given process and the impact particle energy distribution
is known, the rate is defined as a convolution of the impact velocity v(ε), reaction cross section
σ (ε), and energy distribution function (EDF) f (ε, Y˜ ):
Ki (Y˜ ) =

∞
0

dε vα (ε)σi (ε) fα (ε, Y˜ )
22

(2.8)

where ε is the relative impact kinetic energy approximated as the absolute impact species energy
and Y˜ can represent any expression composed of state vector components (e.g. Te or the ionization
fraction, ne /ntotal ). Though most commonly used with electrons as the energetic impact species,
the KGMf supports cross section defined ion impact reactions through EDF declaration for any
charged impact species. Since electron impact reactions are prevalent in global model reaction
networks and the effective electron energy evolves with a conservation equation, the rest of this
work will discuss convolved rates as being functions of the electron energy distribution function
(EEDF), however it all applies to any rate defined through a cross section and an associated EDF.
Cross sections can be defined either numerically or analytically and, analogously to eq. 2.6, cross
sections for reverse inelastic excitation processes can be calculated using the principle of detailed
balance [112].
The largest present caveat with the KGMf is the lack of a Boltzmann equation solver for a
self-consistent evaluation of the EEDF. Though presently fixed with regards to a given system, the
KGMf supports dynamic EEDFs through parameterization of one evolving system term. Without
extra parameterization of the EEDF beyond ε, the convolved rate reduces to a constant value.
An external Boltzmann solver, such as BOLOS [110] or BOLSIG+ [60], can be used to generate
numerical EEDFs. Multiple generated EEDFs can be a function of some variable, such as Te or the
ionization fraction, or analytic expressions can be used, such as the generic form [120]:
3

1

5 )) 2
xε 2 (Γ( 2x
exp −
fe (ε, Te , x) =
3
5
3
3
2
2
( 2 Te ) (Γ( 2x ))

5)
Γ( 2x
ε
3 ) 3T
Γ( 2x
2 e

x

(2.9)

where Te is used as the evolving parameter, x is a fixed shape parameter (x = 1 corresponds to
a Maxwell-Boltzmann distribution, x = 2 yields a Druyvesteyn distribution, etc.), and Γ is the
gamma function. In either case, reaction rate convolutions are performed for a range of values of
the parametrized variable as a pre-processing step. These convolution results are then converted
into either third order b-splines or PCHIP-splines, as functions of the parameterization. These
splines are used to maintain continuous derivatives, using the Cox de-Boor recursion formula [20],
for use within the system Jacobian. Though not entirely self-consistent, this method is accurate
23

assuming the EEDF does not evolve beyond the scope of the chosen parameterization.

2.2.3

Additional Physics Effects

Extending beyond reaction rates, the KGMf contains a few different modules for including terms
for pseudo-spatial effects and non-standard processes. In the simplest form these contribute to
the (G − L) terms in eq. 2.2, however, they can extend with the addition of new species in the
state vector, Y , and corresponding continuity equations into the ODE system, F and J. Aside
from whatever methods are used to acquire reaction data, these modules encapsulate the majority
of assumptions found within the global model reaction network. Though the modules discussed
below are incorporated into the KGMf source code, methods exist for a user to use custom modules
through specification in the simulation input file without requiring source code modification.

2.2.3.1

Pseudo-spatial Effects

Since global models are inherently independent of spatial characterization, effects due to geometrical configuration, gradients, and gas flow must be approximated. Two such native examples are
the flow and sheath/diffusion modules. The flow module adds rates corresponding to inflow of
sourced species and outflow of all species from the system with assumed flow rates and a linear
flow progression across the global model system. Excited neutral species diffusion is treated as
a function of mass, gas temperature, neutral collision cross section, and a characteristic diffusion
length of the assumed geometry. If a plasma sheath profile can be derived analytically or fit to
an experiment or spatial simulation, losses of charged particles nc due to a plasma sheath can be
described as:
Lsheath = uB nc g(hcS , A,V )
where uB =

(2.10)

(eTe /Mc ) is the Bohm velocity, hcS is the edge-to-center density ratio of species nc ,

g is a function of area A, volume V , and sheath density ratio [53].

24

The form of G largely depends on the system and the method of determining the sheath profile.
Analytical expressions are typically found assuming constant ionization of a single species without electron-ion recombination. Though this assumption breaks down with higher gas pressures,
higher pressures also result in the charged species flux terms becoming negligible, countering the
inaccuracy. For example, Godyak and Maximov [53] derived analytical expressions for hcS and
g as functions of length, radius, and ion-neutral mean free path in the context of a low-pressure
discharge and cylindrical geometry.

2.2.3.2

Non-standard Processes

Unlike the relatively simple spatial terms, the inclusion of some effect modules requires the treatment of new evolved species along with their corresponding equations. A simple example is the
specification of plasma voltage and current, dVp /dt and dI p /dt respectively, when the system is
powered with a circuit model. A more complicated example is the module for optical pumping
and stimulated emission. Following the derivation of a two-way averaged plane-wave model [179]
for a three-level laser system, the intracavity intensity, Ψ, is tracked as an evolved species with its
corresponding continuity equation:
Ψ
σ21 c2 hν21
dΨ(t)
= (tL Rm exp{[2β ]} − 1) + n2
dt
τrt
lg

(2.11)

where τrt is the residence time, tL is the transmission coefficient, Rm is the mirror reflectivity, lg
is the gain length, σ21 is the optical emission cross section between states n2 and n1 , and β =
g
1

lg σ21 (n2 − g2 n1 ) is the laser gain. In this notation, n1 is the bottom level, n3 is the pumped level,
and n2 is the lasing level with tracked output for the n2 → n1 transition.

2.3

Process of Operation

As mentioned before, usage of the KGMf can be broken into two pre-processing stages, a system
setup stage, and an evaluation stage. The creation of the local KGMf database and generation of

25

model input files are the two large pre-processing steps for using the framework. Both attempt to
simplify the process of building confidence in a reaction set by using automated selection methods,
human confirmation, and computer sanity checking. The local KGMf database provides a foundation of trusted species and reaction data which can be accumulated from a wide variety of data
sources. Due to potential differences in formatting, syntax, methods, and data type provided from
each data source, the KGMf database uses (automated) methods to map related data for quickly
gathering everything of relevance or detailing similarities and differences between overlapping
datasets. Not only does this reduce the number of data-mismatch mistakes, the mapping is used to
generate reaction and simulation input files with self-consistent checking. After user modification
of the model input files, the third stage of setting up the system of equations for the global model
is largely automatic relative to the user. Finally, the last stage incorporates evaluating the system
of equations and post-processing of terms and results.

2.3.1

Data Acquisition and local database

At the core level, the KGMf is designed for simulating a reaction network between a set of species.
Though extra modules are considered for approximating external effects and spatial variance (e.g.
flow diffusion, wall/sheath losses, etc), almost all data usage corresponds to the rates, or frequency
of effect, for reactions involving one or more species as reactants and products. Regardless of
project scale, physically realistic simulations depend on multiple types of data usually found from
different sources. Adding the prospect of comparing data from multiple sources and varying chemical species makes for complicated bookkeeping. A benefit of the KGMf is the ability to maintain
a local database of species and reaction data by combining generic datafile parsing routines and
consistency checking into a file-and-forget reference of trusted species and reaction data.
With regards to plasma phase reaction data, generic parsing routines are used to import and convert data regardless of source database format. Routines are selected based on expected data type,
known keywords, and formatting similarities to industry standards (NIST, LXcat, and CHEMKIN)
and output is saved into a known internal format in standard units. This is predominantly com26

posed of reference data from the NIST ASD/MSD database [92] being compared against cross
section database files. This process eliminates manual entry error and guarantees unit consistency
when using data from multiple sources. Data is organized and stored by major species, element or
molecule, allowing for overlapping/redundant data from different sources and compressed dictionary accessing. For each major species, a reference dataset is created from the energy structure of
excited levels to which all other datasets will be mapped (First half of listing 2.2). The reference
dataset is used to check parsed-in rates when values are missing (e.g. threshold reaction energy for
an inelastic collision) and to determine group state reactions (e.g. inelastic excitation to an entire
orbital level, requiring degeneracy splitting between individual excited species).
Though accumulation of data into the KGMf database is largely automated and only needs
to be performed once for a given database file, user confirmation that all data has been loaded
correctly is recommended before first use. Aside from checking saved data and mapping files,
plotting routines are used to show data truncation, grouped states, and possible gaps in the overall
description (fig. 2.1 generated by the later half of listing 2.2).

2.3.2

Model Source Generation

In the context of the KGMf, a given model is a combination of two model input files (simulation
and reaction files) and relevant data. The reaction input file is composed of three listings: species
involved in reactions, species variable mapping (e.g. pseudo species representing a combination
of ODE species), and reactions involving species with explicit or referenced unique rate data.
Defining only the reaction network, the reaction input file is completely independent of geometry
and physical approximations to be used. The simulation input file provides configuration for KGMf
modules, such as incorporating sheath losses or lasing, ODE parameters (e.g. fixed or free species
declaration, initial conditions, etc.), and reference data overwriting (e.g. substituting values for
energy levels). Both of these files are human readable and possible to create from scratch, however
automatic generation methods are included for convenience and to reduce human error. Examples
of input files are explained in Appendix A.2.
27

Figure 2.1 Cross section data for Ar excitation from the LXCat IST Lisbon [71] database. Horizontal lines are species energy levels, vertical lines are electron excitation reactions, dashed lines
are radiative reactions, and red boxes highlight incomplete orbitals or orbitals constructed from
grouped excitation.

import src
# Species mapping
mapper = src. KGM_exe . mapping_species . SpeciesMapping . MappingPureKGMdata ()
savename = ’Cross_LXcat__LX_IST_Lisbon_all . txt_Ar .npz ’
mapper . MapSpecies (’Ar’,savename )
# Plotting of DB mapping
data_plot = src. KGM_visualize . PlotMapping . PlotPureSpeciesMapping ()
layout_dat = data_plot . All2PrepPlot (’Ar’,savename )
flat_data = [ subelem for elem in layout_data for subelem in elem]
handle = data_plot . PlotLayout (flat_data ,tit=’LXcat IST Lisbon Ar Data ’)

Figure 2.2 Single data mapping for Ar with Biagi data

28

import src
# Model generation
generator = src. KGM_exe . model_gen . SpeciesModeling . GeneratePlasModelFiles ()
generator . AllPLASProcessing ([’Ar’], DBstyle =’Biagi ’,Save=True)
# Inner access for debugging
# maplist = generator . LocateDataSources ([’Ar ’], DBstyle =’ Biagi89 ’)
# gatdat = generator . GatherSpeciesData ( maplist )
#out = generator . GeneratePlasModelFile2 (gatdat ,maplist ,save=False)

Figure 2.3 Model generation for Ar with Biagi database [16].
Taking advantage of the KGMf database of previously mapped species data, generic model
input files can be generated automatically. In the purely plasma phase case, a user specifies which
major species are used and from which mapped source to pull reaction data, listing 2.3. This
results in input files incorporating every relevant reaction from the KGMf database, which takes
into account species decomposition and the expansion of the excited states. The inclusion of gas
phase reactions is done by parsing CHEMKIN-like mechanism and thermo files [85], which are
then mapped on the fly to the excited species in the local KGMf database. More advanced methods
can be used to: truncate energy levels of interest, use multiple data sources to fill reaction gaps,
and alter mapping between gas-phase and plasma-phase species. These automatic methods extract
the relevant data from the KGMf database and compress them into a reaction data file, which no
longer needs reference to the main database and can be packaged with the input files for portability
across multiple machines.
The input files can then be easily edited to alter species of interest, ignore certain reactions,
or add reactions which were not included in the KGMf database. In the same manner as would
be done when not using the KGMf database, extra reaction data files can be created and linked to
reactions in the reaction input file.

2.3.3

Preparing for Evaluation

As was alluded to previously, the system of equations defining the global model is first set up
symbolically and then rendered into a callable function. Resulting in a compiled version of the
29

ODE system, this setup only needs to be done once for a given reaction mechanism and can otherwise be referenced during later evaluation. First a generic system of equations is created which is
dependent solely on the reaction network being used, equivalently expanding continuity equations
with symbolic place holders for each term. This generic system creates the symbolic skeleton from
which a data dependent symbolic system is created. Rates which are defined via expressions are
simply substituted for their placeholder symbol while rates defined by convolution are wrapped
as SymPy generic symbolic functions. These symbolic forms are created to: generate a symbolic
Jacobian for the extremely stiff system of ODEs, take advantage of existing functions in SymPy
to generate Cython wrapped c-code and to use as reference during post-processing. Since any
collection of symbols and terms is accessible as a function, the user is able to follow the symbol
substitutions to look at how any term evolved within ODE integration.
In order to speed up both compilation and evaluation time, the KGMf takes advantage of the
known symbolic structure to optimize the generated C-code. Since nearly every reaction will
appear in multiple continuity equations, and some continuity equations contain others in their
entirety, the symbolic system is recursively processed for common-sub-expressions to be evaluated
separately. This implies that no collection of operations is ever performed multiple times in a
single function call to either the ODE system or the Jacobian. Reducing the operation count by
over an order of magnitude, paired with optimization flags for the compiled C-code, results in a
significantly faster evaluation time than remaining in pure Python.
Finally, the ODE system and Jacobian are each wrapped into Python functions broken into
two steps that occur with every call. For a given set of parameter values at some time step, all of
the spline rate coefficients are evaluated in a vectorized form. These rate values, paired with the
parameter values, are then passed to the compiled code to yield the differential change for each
tracked species.

30

Custom
User
Control

KGMf
Model
Generation

Custom
User Data

Auto Data
Collection

Reaction
Data

Control
Files

EEDF
Recovery

Integrate
Rates

EEDF
Solve

Compiled
ODE *

Batch
Methods

KGMf
Interactive

Symbolic
ODE *

External
Linking

Figure 2.4 KGMf flow chart of the latter three user stages. Standard usage involves modifying
control text files for extra reactions and system parameters, while every section is accessible for
modification in a Python instance. Red dashed lines mark the future work of incorporating a
Boltzmann equation solver.
2.3.4

Evaluation and Usage

Acknowledging the potential stiffness of the equations, the Python wrapped ODE system and
1

Jacobian are integrated in time to describe the overall evolution of species densities and energies.
Presently native KGMf functions exist for using SciPy’s various integrators under odeint and ode
(see listing 2.5), however the functions can be passed by reference to any integrator library. While
the default is to use odeint for convenience, usage of the ode class allows for on-the-fly analysis
(e.g. termination upon reaching steady state) and dynamic time-stepping for handling extremely
stiff problems.

31

mfode = src. KGM_exe . model_eval . ModelEval . ModelFile2ODE ( NumCPUs = 12)
modnam = ’BASIC_PLAS__Ar ’ # Directory of model files
init , aux , sims = mfode. All2PrepSpecific ( ModelName =modnam , \
SimFile =’SI.txt ’,ReacFile =’RF.txt ’)
td = [0.0 ,1.e −5,1000]
ode_out ,tvals ,sfo ,tfo , pltmod = mfode. KGMode . IntegrateODESystem (init , \
aux , timedat = td , pltresult =False)

Figure 2.5 Basic KGMf integrator usage.
2.3.4.1

Interactive Usage

Since the symbolic versions of the species continuity and energy equations are created from the reaction input file, references to symbols, and symbolic expressions, can be used for post-processing.
In a similar manner as to how the system of symbolic ODEs can be converted into a callable numeric function, any expression composed using the same symbols can be rendered callable. This
allows for visualization of numeric and analytic rates, along with rate terms (including the product
of reactant species), being compared to the overall system evolution. This capability is invaluable to both developing an understanding of dominant reactions and determining if reaction data is
correct.

2.3.4.2

Batch Operation

Aside from being able to access the symbolic equations of the global model, at various points
during the middle two stages, hash values are generated corresponding to macroscopic (e.g. full
reaction input file) and microscopic (e.g. selection of EEDF) quantities in order to determine
recovery points. During cases when a simulation is running multiple times with varying inputs, the
integer hash values are used to quickly determine which of the previously saved data can be reused.
This could be as simple as a change in initial condition, which would result in the same final Python
wrapper being evaluated, or as involved as redefining the EEDF, which results in new numerical
reaction rates being found and wrapped with the previously compiled C-code. Due to the reduced
evaluation time from pre- processed reaction rates and recovery points, batch parameter scanning

32

and target value search methods can be used over relatively large search spaces.

2.4

Example Usage with Verification

In order to instill confidence in the KGMf and to demonstrate functional capabilities, a few example models are included with the KGMf for tutorial and verification purposes. In the simplest
configuration, these models provide a sandbox for familiarization with input files and interactive
usage, while packaged input files attempt to cover the range of possible inputs and called modules.
Two systems in particular are used as verification of the plasma and gas phase chemistry methods:
a recreation of the argon plasma global model presented by Ashida, Lee, and Lieberman [4] and an
example of classical combustion using the validated GRI-Mech 2.11 reaction set [143] recreating
an original validation case [29].
A script within the main KGMf directory, “quick_rebuild.py” is used to populate a local
database, with packaged data files from NIST [92], the Biagi v8.97 dataset hosted on LXCat [16],
and to generate models for demonstrating code capabilities and verification. Example scripts and
supplemented control files are located within the examples sub-directory “data/Example_Source/”.

2.4.1

Plasma Verification

The work by Ashida et al. focused on studying “the behavior of argon plasmas being driven by
time modulated power” using a global model. With an emphasis on plasma reactors, a cylindrical
geometry is assumed with a microwave-driven plasma composed of argon in the ground, ionized,
and electronically excited states truncated to include the 4p level. With a constant gas temperature,
this system is used to verify implementations of plasma chemistry components and eq. 2.4. Published with two modeled species representing the grouped four 4s and ten 4p energy levels, KGMf
comparison is performed with a faithful-recreation model and a model having distinct species for
each energy level. In the first case, Arrhenius fitted reaction rates are used from the original work,
while in the latter case, cross sections from the LXcat-hosted Biagi v8.97 database [16] are used

33

in conjunction with an assumed EEDF for each electron-impact reaction.
Example files are included in the KGMf to completely recreate the original work using both
models; control files for each model along with scripts demonstrating generation of a base model,
interactive usage, batch running configuration, and access of batch saved data compressed with
H5Py. With a pressure of 5mTorr, gas temperature of Tg =600 K, and reactor with dimensions of
R=15.25cm and L=7.5cm, four different absorbed power assumptions are tested: continuous wave
to steady state, square wave, smooth change (square with coupled sinusoid segments to soften
edges), and square with nonzero“off”. In the latter three cases, multiple runs are performed using
differing duty ratios and periods while maintaining the same total absorbed power. First looking at
the transient behavior of the plasma, fig. 2.6 shows the originally published results (fig. 1.b) along
with the faithful recreation and model with cross section data and an assumed Maxwellian EEDF.
Though in very good agreement, differences between the original and faithful recreation can be
attributed to rounding differences on reaction threshold energies and initial conditions passed to
the integrator.
Exactly the same reaction pathways are allowed between the faithful recreation and the model
using cross section data, the reaction rates differ due to both cross section data and functional representation. Ashida et al. used rates from the work by Kannari et al. [83] which relied on modified
Arrhenius rate fitting of semi-empirical and experimental cross sections, from the 1960s and 70s,
being convolved with a Maxwellian EEDF. By comparison, the second model relies on cross section data updated in 2002. The difference between rate implementations is highlighted in fig. 2.7
which plots the reaction rate for argon ionization from ground versus the effective electron temperature. For the pressure and temperature considered in the original work, a Maxwellian EEDF is
a valid assumption, but for the purposes of better understanding the effects of the assumed EEDF,
examples include batch runs with both Maxwellian and Druyvestyn distributions, with square wave
time averaged results compared in fig. 2.8. As opposed to needing to find new rates in the literature
which assume a specific EEDF, rates are convolved using the same cross section data, guaranteeing
consistency.

34

Figure 2.6 Transient square wave behavior for a period of τ = 100µs, duty ratio of 25%, and
500W of total absorbed power from: the original paper (x markers), faithful model with guessed
unreported parameters (dotted) and tuned parameters (solid), and recreation with cross section
dependent rates (dashed).
Example scripts relating to recreating the work by Ashida et al. are denoted with “*VV_ALL.py”.
Control files are located within the “BASIC_PLAS__Ar” generated model directory with “*control*” and “*reALL*” denoting the faithful model and KGMf recreation respectively.

2.4.2

Gas Phase Verification

Verifying the implementation of gas-phase chemistry components and eq. 2.3, example models for
classical CH4 -H2 -O2 -Ar combustion are generated from the GRI-Mech 2.11 mechanism. Though
more modern mechanisms exist, such as GRI-Mech 3.0 and NUI Galway [126], GRI-Mech 2.11
is small enough for quick validation while demonstrating caveats of an incomplete mechanism
relative to pressure and temperature regimes. Due to the mass exchanging reactions (between
major species), the defining system of ODEs for a given gas-phase model can be extremely stiff
and require special care to integrate in a stable manner. Especially when the gas temperature
equation is coupled with species continuity equations, values of the ODE state vector can range by
35

Figure 2.7 Argon ionization rates defined with: Arrhenius rate used by Ashida et al. (solid black)
and KGMf spline implementations for convolutions of an assumed Maxwellian EEDF with cross
sections from: the LXcat Biagi database (solid blue) and estimates of the cross section extracted
from a figure in the Arrhenius rate source [127].

Figure 2.8 Comparison between results from different assumed EEDFs. With duty ratios of 10% ≤
α ≤ 50%, the faithful model is plotted in black while the cross section dependent model for a
Maxwellian (x = 1) and Druyvestyn (x = 2) EEDFs are in blue and green respectively.

36

Figure 2.9 Example output a CH4 -O2 -Ar mechanism with 90% argon and a stoichiometric fuel/oxidizer mixture initiated at a pressure of 1 atm and temperature of 1400 K.
over twenty orders of magnitude, as can be seen from fig. 2.9, resulting in a sensitivity to integrator
parameters, time steps, and initial conditions.
The first gas phase model to be considered is a replication of the experimental work by Cheng
et. al [29] which was used by the GRI-Mech group for mechanism validation. Using a shock tube,
autoignition induction times (ignition times) were found for a variety of gas mixtures covering
pressures of 1 to 3 atm and temperatures of 800 to 2400 K. An induction time correlation formula
was fitted with dependence on fuel and oxidizer densities, ratio of H2 and CH4 , ζ , and equivalence
ratio, φ
τ[CH4 ]0.48(ζ −1) [H2 ]−0.14ζ [O2 ]1.94−1.38ζ

= (1.19e − 12)(1−ζ ) (1.54e − 4)ζ e

46.4−29.2ζ
RTg

(2.12)

where R is the universal gas constant and [∗] denotes species density in units of mole/cc. While
example scripts generate models which compare against each of the eleven different gas ratios used
by Cheng et al. and recreate the four published induction time plots, herein we present a recreation
of the case with fuel defined as 66%CH4 -33%H2 (ζ = 1/3, denoted as mixture number 3 in the
original work). Fig. 2.10 plots the original experimental data, fitting function eq. 2.12, and scaled
37

Figure 2.10 Density-scaled ignition delay times for a mixture of CH4 -O2 -Ar combustion with
ζ = 1/3. Experimental data points were extracted from figure 5 in [29], the black line is the
reported fitting function, and the blue line represents scaled estimates from the KGMf recreation
model.
ignition times from the KGMf model for a range of initial temperatures and initial pressure of 1
atm. Since pressures were not defined for each point in the original work, the y-axis is scaled with
species densities matching fig. 5 presented by Cheng et al. The deviation of the fitting function
from the experimental data points can be explained by both their extraction from the original work
and the fitting originally done with two different gas mixtures.
Similarly to the plasma model, a script is provided to integrate multiple conditions and generate
curves for ignition delay times to compare with original mechanism validation work [29]. Example scripts for the GRI-Mech 2.11 validation models are denoted with “*VV_GRI.py”. Control
files are located within the “GP_GRI_Ar_CO2_CH4” generated model directory with “*cheng*”
denoting original validation tests.

38

2.5

Conclusion

Due to interest in using global models for the study of kinetic plasma systems, the KGMf was
developed for analysis of mixed-phase chemistry as both a tool for research and an educational
platform. Implementation in Python allows for both scripted and interactive usage while remaining
open-source allows for user contribution to development and guaranteed reproducibility. In using
a symbolic description, the KGMf is able to: import user-defined analytic terms, automatically
calculate global model Jacobians, minimize the operation count of evaluation, and provide access
to sub-components of the numerical system. With the largest caveat of global models being data
dependence, the KGMf was designed to quickly compare the effects of different data sources and
assumptions. The design of KGMf allows for the addition of new physics, as demonstrated by
the laser model. While operable in nearly a fully automatic manner, treatment of global models
as initial value problems requires system specific tuning with regards to integrator parameters and
initial conditions.
Examples are included to demonstrate KGMf features and both verification and validation of
generated models in two contexts. Verifying species continuity and electron temperature equations, the modeling work by Ashida et. al [4] was recreated with originally published data as well
as with cross section and EEDF declaration. This tested the implementation of conservation equations and time-dependent power absorption envelopes (constant, square, and smooth square) using
an argon gas at 5mTorr. Requiring tuning for unpublished parameters, KGMf matched the published global model well aside from integration equilibration from incorrect species initialization.
Using different data sources and explicit representation of states demonstrated variation in results
from disagreement between experimental and analytical data. Agreeing with mechanism validation efforts by GRI Mech [143], the KGMf was used to model methane combustion at atmospheric
pressure as reported by Cheng et. al [29]. Modeling results agreed with experimental results in
systems with high initial hydrogen concentrations; methane systems agreed qualitatively however,
demonstrated longer ignition times than originally reported. This is attributed to low chain branch-

39

ing rates associated with hydrogen radical generation.

2.6

Acknowledgments

This work was supported by AFOSR-BRI grant FA9550-12-1-0343, DOE Plasma Science Center
grant DE-SC0001939, and a Michigan State University Strategic Partnership Grant.

40

CHAPTER 3
PLASMA-ASSISTED COMBUSTION WITH KGMF

3.1

Introduction

The study of plasma-assisted combustion (PAC) and plasma-assisted ignition (PAI) is motivated
by the goal of optimizing combustion characteristics for novel engine concepts. The generality
of this statement lies with the different characteristics of interest for different applications, while
developing a better understanding of the underlying chemistry applies to every scenario. For over
a century, spark-ignition internal combustion engines have use thermal-equilibrium plasmas in the
form of arc discharges; however, these discharges rely on localized gas heating to initiate and propagate combustion, and provide limited options for control of ignition process [115]. Starting in the
20th century, the use of non-equilibrium plasma for ignition and combustion has received interest
due to possibilities for modified ignition pathways, flame stabilization, and byproduct control. The
seminal work of Haselfoot and Kirkby [66] introduced the concept of electrical interaction with the
combustion process, resulting in an explosion of work to study the physics involved and potential
applications. More recently in the 1980s, the effect of non-equilibrium plasma discharges on combustion kinetics was demonstrated by the pioneering experimental works by Kimura et. al [87] and
Behbahani et. al [11] with plasma jet applications to high speed propulsion and NOx destruction
respectively.
A plethora of experimental studies have demonstrated characteristics of different discharges
with different gas mixtures in different regimes, however, even with a growing list of methods that
can be shown to benefit the combustion process, it remains unclear which plasmas are optimal
for specific applications. While the literature contains multiple review papers on PAC/PAI studies
and applications [152, 129, 153, 149, 1, 81], the largest gap in PAC/PAI theory that remains is an
understanding of the intricacies of plasma-combustion chemistry. Numerical modeling presents an
opportunity to dissect components of the PAC/PAI processes while neglecting or approximating

41

complexities which arise from completely coupled systems.
This work provides an overview of PAC fundamentals in the context of volume-averaged global
modeling followed by an analysis of the kinetic reaction chemistry effects resulting from nonequilibrium plasma discharges. Starting with an introduction to classical combustion theory, this
chapter then continues into the general motivation and fundamentals of PAC/PAI. Nanosecond discharges are chosen for modeling following the progression of experimental works highlighting the
focused study of non-equilibrium plasma chemistry in the ignition process. Numerical modeling
results follow a description of the specifics of the global model formulation and systems used.

3.2

Classical Combustion

Prior to the end of the 18th century, the concept of combustion relied heavily on the legacy of
Greek philosophers to whom fire was a classical element. Developed formally by Georg Ernest
Stahl (following the works of Johann Joachim Becher) phlogiston theory attempted to explain
oxidation as the release of a fire-like element through the process of dephlogistication [175]. The
lack of an explanation for reaction products having a larger mass than the original reactants, such
as in the creation of metal oxides from heating metals in air, seeded doubt and eventually led to
the theory’s abandonment. While supporters attempted to define a negative mass of phlogiston in
their defense, the root of the problem was in not knowing the composition of the air.
Starting with experiments in the 1770s attempting to explain the results of phlogiston-supporter
Joseph Priestly, Antoine-Laurent Lavoisier determined that the composition of air contained a mixture of at least two gases: an inert gas and one on which combustion was reliant. The latter, found
in most acids and supporting both respiration and combustion, he named “oxygène,” from the
Greek words for acid generator. With his public attack on phlogiston theory beginning in 1783
based in a demand for strict involvement of the scientific method, Lavoisier presented the modern
foundations of combustion in his book Traité élémentaire de Chemie [99]. Though Lavoisier proposed empirical arguments, the development of a kinetic theory of gases, by James Clerk Maxwell

42

and Ludwig Boltzmann in the 19th century, led to an explanation of the heat and energy involved
in combustion.
Combustion is an exothermic chemical reaction incorporating both a fuel and oxidizer, complicated by a dependence on both mechanical and chemical processes. While a complete picture of
the combustion process requires mechanical processes , such as mass and heat transfer, the scope of
this work aims to develop an understanding of the chemical processes and demonstrate how plasmas might provide augmentation of combustion characteristics. In the context of volume-averaged
modeling, these mechanical processes are handled with approximations (e.g. loss of particle density due to gas flow) as opposed to having fully-resolved dynamics (e.g. turbulent mixing of fuel
and oxidizer). With the implementation of mechanical approximations discussed later, this section
will provide a brief introduction to chain branching and explosion limits.

3.2.1

Chain Reactions

Prior to delving into the methods by which plasma can affect a combustion process, it is important
to first understand how these processes evolve in a classical context. Consider, for example, the
stoichiometric global reaction describing the interaction between molecular hydrogen and oxygen:
2 H2 + O2 −−→ 2 H2 O
While descriptive of the overall process, this reaction does not exist as a single step. Instead, the
combustion process relies on radical species, such as H, O, and OH, to propagate energy release.
Describing the system adequately, independent of pressure regimes, requires eight species and
sixteen elementary reactions. The inclusion of radical species, or chain carriers, in the reaction
mechanism provides the dynamics of the combustion process. These elementary reactions support
an equilibrium composition of reactants and products at a given temperature and pressure, allowing
for a reverse rate to be calculated from the change in Gibb’s free energy. In the limiting case
of low pressure, a simplified mechanism of seven reactions adequately describes the system and

43

introduces the concept of chain reactions:
k

H2 −−0→ 2 H
k

O2 −−0→ 2 O
k1

O2 + H −−→ OH + O
k2

H2 + O −−→ OH + H
k3

H2 + OH −−→ H2 O + H
k4

H + wall −−→ 0.5 H2
k5

H + O2 + M −−→ HO2 + M

(R.1)
(R.2)
(R.3)
(R.4)
(R.5)
(R.6)
(R.7)

where M is a collision partner representing any major species, therefore directly related to the
overall system pressure, and k∗ represent reaction rates. Without some amount of chain carrier
species, the reaction mechanism would not progress through a combustion process due to a lack
of used elementary reactions. These reactions fall into one of four categories with regards to
progression through the chain reaction. Triggering the mechanism, chain initiation reactions, R.1
and R.2, providing an initial concentration of radical species. Reactions R.3 and R.4, classified
as chain branching, further produce radical species for accelerating the mechanism. Classified
as chain propagation, reaction R.5 produces the final product, H2 O. Finally, chain termination
reactions, R.6 and R.7, remove carrier species from the system.
When comparing the rates of the chain branching and propagation reactions, R.3 progresses
significantly slower than both R.4 and R.5, resulting in O and OH radicals being consumed faster
than H species. This means that H production determines the majority of the combustion mechanism propagation in the low-pressure regime. Other species become limiting factors as pressure
increases due to different sets of reactions which become dominant.

3.2.2

Explosion Limits

The asymptotic behavior of radical species production directly defines the transition between deflagration and detonation, or the explosion limits, of a given system. While difficult to perform with
a complete reaction mechanism, explosion limits determined analytically come from considering

44

regime specific simplified mechanisms and assumptions of steady state behavior. First, writing the
species continuity equation for each species, using nα to represent concentrations of species α,
dnH
= 2k0 nH2 − k1 nO2 nH + k2 nH2 nO + k3 nH2 nOH − k4 nH − k5 nH nO2 nM
dt
dnO
= 2k0 nO2 + k1 nO2 nH − k2 nH2 nO
dt
dnOH
= k1 nO2 nH + k2 nH2 nO − k3 nH2 nOH
dt
Imposing the assumption of steady-state behavior of O and OH species relative to H, by setting
their continuity equations to zero, results in steady-state values of
nO |ss =

2k0 nO2 + k1 nO2 nH
k2 nH2

and

nOH |ss =

2k1 nO2 nH + 2k0 nO2
k3 nH2

where nO |ss has been substituted into nOH |ss in order to remove all dependence on O and OH
species. Plugging the two steady-state forms into the continuity equation for H and simplifying the
result
dnH
= 2k0 (nH2 + 2nO2 ) +(2k1 nO2 − (k4 + k5 nO2 nM ))nH
dt
s

b

t

where s represents a general source term, b denotes chain branching and t chain termination. The
behavior of H production reduces to two cases, depending on the relationship between b and t,
when analyzing the analytic solution of
nH =

s
1 − e−(t−b)t
t −b

In the case of t > s, the H concentration grows to an asymptotic value of s/(t − b) limiting the
overall speed of the mechanism. While in the case of s > t, H species experiences exponential
growth which in turn further increases the mechanism velocity resulting in an explosion as opposed
to a slow and incomplete burn.
Not only does this analysis use a simplified mechanism neglecting reactions dominant at higher
pressures, but the relationship between s and t depends on both total system pressure, through M
collision partners, and O2 concentration, resulting in extremely different behavior under different
45

pressure regimes. Repeated asymptotic analysis in each of the pressure regimes, using subsets of
the 16 reactions (Appendix B.1) corresponding to their dominance at different pressures, results
in the explosion limits shown in fig. 3.1. Here, the first limit occurs as chain carrier production
exceeds wall termination due to pressure increasing. As pressure increases, the relatively-inert
HO2 production, R.5, consumes carrier species and impedes run-off behavior, resulting in the
second limit. Further increasing the pressure, heat release from the exothermic reactions cannot
be dissipated, resulting in higher rates of endothermic reactions dissociating HO2 and H2 O2 to
produce more carrier species and an explosion.
While asymptotic analysis provides analytical solutions, the effect of simplifying reactions
mechanisms alters the exact locations of the explosion limits. The dotted line in fig. 3.1 corresponds to the second explosion limit in a stoichiometric mixture, where 2k2 = k5 nM , if all other
reactions considering dissociation of HO2 are neglected. Without analytical solutions to complete
combustion mechanisms, the determination of the explosion limits requires the use of numerical
integration or experimental studies.

3.3

Plasma Assisted Combustion

While understood for the last two centuries, humanity’s methods of implementing classical combustion processes has not drastically changed since the patents of Nikolaus Otto and Rudolf Diesel
in the late 19th century. Until recently, relying on thermal effects to drive the combustion process
has not limited the ability to optimize designs of modern combustion platforms with regards to
environmental, economic, and performance characteristics. Modern technologies and emissions
regulations have created demands beyond the capabilities of classical combustion [111]. Of growing social and governmental importance is the control of harmful emissions with regards to both
climate warming and air quality. With our present infrastructure and the energy density provided by
hydrocarbons, synthetic and fossil fuel combustion will remain a dominant energy source well into
the 21st century [121], though present methods result in the production and release of greenhouse

46

Figure 3.1 Experimental explosion limits of stoichiometric H2 -O2 combustion in a spherical KClcoated vessel from [30] with original data credited to [103]. Dotted line corresponds to where
2k2 = k5 nM .

47

Figure 3.2 Schematic of plasma assisted combustion enhancement pathways taken from [80].
gases (e.g CO2 ) and pollutants (e.g. NOx ). The prevalence of internal combustion car engines and
gas turbines for both aircraft propulsion and electricity generation paired with increasing fuel costs
have also driven a demand for better operating efficiency. The quest for efficiency and emission
control has resulted in modern designs that push the boundaries of classical combustion lean and
low-temperature flammability limits.
Beyond efficiency and emission control, interest in hypersonic propulsion calls for ignition
control in extremely difficult conditions. At or near Mach 5 flight, the flow residence time in a
combustor is significantly shorter than ignition time offered by spark sources [35] resulting in an
incomplete burn. The use of specialized cavity geometries initially showed promise with regards
to ignition and flame holding; however, thermoacoustic instabilities result in unstable operation
risking blowout and the altered geometries induced stagnation pressure hampering engine efficiency [13]. As with fuel efficiency and by-product control, this further motivates a search for new
methods of reliable, fast, and robust ignition and flame holding methods.

48

3.3.1

PAC Fundamentals

In the context of using a discharge to initiate combustion or maintain a flame, multiple mechanisms
exist through which the generation of a plasma can affect the gas. Two thermal effects come from
an (nearly-) equilibrium plasma in the form of: 1) local gas heating as a result of energy release
increases chemical reaction rates, and 2) flow perturbations generated from inhomogeneous gas
heating causing turbulence and mixing. In turn, the formation of a non-equilibrium plasma results
in three new mechanisms for combustion modification: 1) excitation, dissociation, and ionization through electron impact leading to non-equilibrium radical species production, 2) momentum
transfer from the electric field to the gas through interaction with space charge (ionic wind), and 3)
ion and electron drift can lead to modified radical fluxes along gas flow gradients [153]. No obvious ways of separating these effects exist due to their inherent coupling, fig. 3.2, in turn requiring
special considerations for systems attempting to analyze a single mechanism.
Of primary interest to this work is the generation of radical species from non-equilibrium
plasma chemistry. Considering the set of reactions presented in the previous section, the introduction of electron impact driven or initiated chain initiation, branching, and propagation pathways
has the potential to completely alter asymptotic behavior. An extremum case would be the bypassing of the extinction limit resulting in a completely monotonic explosion limit curve. Three types
of species excitations are involved in the non-equilibrium generation of radicals: vibrationally and
electronically excited species as well as ions. Unlike the short lived rotational states, vibrationaltranslational relaxation of vibrationally excited states of N2 and O2 occurs over 10-100 µs, though
inversely proportional to hydrocarbon density [129]. This allows for an accumulation of vibrational
states which undergo reactions with other vibrationally excited molecules such as
k (v)

H2 (v) + O −−−→ H + OH(w)

(R.8)

Multiple reaction mechanisms exist using a vibrational temperature which assumes a Boltzmann
distribution over the vibrational levels, however, their validity comes into question in non-equilibrium
conditions [26]. Instead, a state-to-state model employs tracking of individual densities of vibra49

Figure 3.3 Comparison of chain propagation reactions with electron impact dissociation rates
from [149] using an EEDF found using a two term Boltzmann solver assuming different field
strengths, E/N, and a fixed gas composition.
tionally excited states [151]. As shown in [104], rate constants can vary by several orders of
magnitude for different vibrational levels of both reactants and products, requiring a description of
individual vibrationally excited levels.
In addition to the contribution from vibrational levels, four different methods exist through
which to produce radicals from electronically excited molecular states as a result of electron impact [152]:
1. Excitation to a pre-dissociative state resulting in dissociation and the formation of two radicals

3
1
O2 + e −−→ O2 (B3 Σ−
u ) −−→ O( P) + O( D) + e

(R.9)

2. Excitation to an electronic state which produces radicals through chain reactions
O2 + e −−→ O2 (a1 ∆g ) + H −−→ O( 3P) + OH

(R.10)

3. Excitation of the molecule with dissociation via quenching with another molecule
1
O2 + e −−→ O2 (A3 Σ+
g ) + CH4 −−→ O2 + CH3 + H( S)

50

(R.11)

4. Excitation to an electronic state with radiative emission resulting in molecular dissociation
3
H2 + e −−→ H2 (a3 Σ+
g ) −−→ H2 (b Σg ) + hν

(R.12)

O2 + hν −−→ O + O

(R.13)

While the second mechanism only requires a weak electric field due to low energy thresholds of
excited states (≈ 1eV for oxygen), the other mechanisms require a relatively large field strength
in order to have high energy electrons reach threshold energies (

5eV for oxygen). Depending

on both the electron energy distribution function (EEDF) and electron density, these reactions can
occur at a significantly faster rate than thermally driven radical generation. For example, with
E/N of 100-300 Td, the electron impact dissociation occurs two orders of magnitude faster than
classical radical sources at the equivalent temperature, as shown in fig. 3.3. With cross sections
known for electron impact reactions, a Boltzmann equation solver can yield an EEDF for rate
constants; however, the both the distribution and resulting rate depend on both the electric field
strength as well as the ratio of gases involved. While comprehensive kinetic models exist for
different gas mixtures (e.g N2 -O2 [91], H2 -N2 -O2 [150],O2 -Cx Hy [3]), it is agreed that rates are
not necessarily well known in all regimes and cross section data may not be available.
Finally, ions can also participate in dissociation and recombination reactions resulting in radical
generation. While ionization requires a significant amount of discharge energy and ions typically
have a quick recombination time (1-10 µs with Tg of 300-3000 K respectively), it bears responsibility for maintaining the electron population in the system along with introducing new pathways
by which to produce radical species:
1. Electron-ion dissociative recombination
O2 + + e −−→ O( 3P) + O( 1D)

(R.14)

O4 + + e −−→ O2 + O + O

(R.15)

O2 + + O2 − + M −−→ O2 + O + O + M

(R.16)

2. Recombination of cluster ions

3. Ion-ion recombination

51

Compared to the four pathways discussed previously, ion reactions result in higher heat release
and less radical generation [2]. Though not efficient pathways to alter radical production, these
reactions are not so insignificant as to be excluded from considered chemistry mechanisms.

3.3.2

Physical system for study of kinetics

The lack of spatial variation in a global modeling study makes it impossible to study the effects of
flow perturbations and inhomogeneous radical flux, however, it bypasses the need for an extremely
resource-heavy coupling of full reaction chemistry with a fluid or PIC code. As such, the experimental platform for numerical studies should demonstrate a weak dependence on geometry and
flow variation. Studies of plasma discharge modifying ignition characteristics started back in the
early 20th century with the work by Haselfoot and Kirkby [66], demonstrating ignition with electrical discharge in low pressure regimes which would otherwise not ignite. Gorchakov and Lavrov’s
study on the interaction between electrical discharge and autoignition times [56] demonstrated that
this effect was mapped across the entire range of the classical explosion limits, fig. 3.4, and not just
the first ignition limit. While applications of equilibrium plasmas have been extensively studied,
thermal mechanisms from plasma interaction do not present an efficient method for combustion
enhancement due to an inability to selectively choose excitation pathways. In searching for efficient methods of combustion augmentation, initial studies attempted to understand the interaction
between the flame and electric field [102]. While showing increased electron temperature [21]
and improvement of flame extinction limits [75], the strong electric field resulted in an ionic wind
whose effect was not separable from electron and ion chemistry effects. Attempting to isolate the
effects of electric field enhancement from the ionic wind, the use of microwave discharges showed
both an increase in electron and burned gas temperatures [57]. Recent numerical modeling of the
microwave interaction showed accelerated flame propagation with a flame temperature increase of
≈200K [157]. However, in limiting microwave power to avoid breakdown, the electron density
was low, resulting in the burned gas absorbing a large portion of microwave energy and an overall

52

Figure 3.4 Experimental extension of explosion limits of H2 -O2 combustion under electrical discharge from [149] with original reference to [56].
loss in discharge modification efficiency.
Further attempts to isolate flow effects from the plasma combustion chemistry with regards
to ignition and extinction limits, while increasing the strength of the electric field and number
of electrons relative to microwave discharge experiments, resulted in idealized experiments using
gliding arcs in a counter flow diffusion flame [49]. Driven predominantly by thermal effects,
the flame extinction limit increased by a factor of three, while kinetic effects were attributed to
lowering the ignition temperature in hydrogen and methane air systems through reaction pathways
involving NOx species [122]. Followed by experiments using nanosecond discharges showing
kinetic effects to both ignition (predominant) and extinction limits [158], questions arose as to
what extent non-thermal effects can alter extinction limits and whether or not kinetic pathways
could drive ignition speed so low as to bypass the classical flame extinction limit altogether.

53

Figure 3.5 Types of discharges with regards to electron density and electric field values, taken
from [152].
3.3.3

Nanosecond Discharge System

Nanosecond discharge (NSD) presents the best experimental platform for a focused study of
plasma interaction with the combustion process by allowing for selective excitation gas. Relative
to other discharge methods, fig. 3.5, NSD provides a large amount of input energy being consumed
for the production of high energy level states, otherwise not attainable from lower E/N discharges,
without resulting in large thermal effects. As a result of peak field intensities existing on the order of thousands of Townsend for 2-3 ns before decreasing to a few hundred, a uniform discharge
develops in space allowing for modeling in a volume averaged context. The fast ionization wave
created by the peak field intensities develops on a significantly shorter time scale than ignition,
allowing for separation of the discharge and combustion kinetics.
Using a counter flow flame configuration, Sun et. al demonstrated the removal of the standard
flammability limits as a function of input flow fuel mass fraction [159], further motivating the use
of a NSD. The experimental system consisted of two stainless steel nozzles, with inner diameters
of 25.4 mm and flow diffusers, facing each other at a separation of 16 mm. A counter flow system
was chosen due to simplified 1-D flame geometry and a flow velocity gradient (strain rate) along

54

(a) Oxidizer mass fraction XO = 0.34

(b) Oxidizer mass fraction XO = 0.62

Figure 3.6 Experimental demonstration of the bypassing of extinction limits. Plotted against fuel
mass fraction, XF , with different cases of oxidizer mass fraction, XO , left demonstrates the classical
S-curve behavior, while the case with higher oxidizer levels, right, shows monotonic behavior;
taken from [159].
the centerline as the two jets reach a stagnation plane [122]. The strain rate, dependent only on
nozzle separation, gas densities, and flow velocities, is the inverse of the flow residence time and
is used to scale ignition time measurements. Ignition occurring at the same temperature but with a
higher strain rate corresponds to a lower ignition time. The nozzles attached to a high voltage pulse
generator operating with a 12 ns (FWHM) pulse at a constant 7.6 kV and repetition frequency of 24
kHz. Time integration of voltage and current profiles estimated an energy input of 0.73 mJ/pulse
corresponding to a total input power of 17.5 W. With the combustion chamber kept at a constant
pressure of 72 Torr, CH4 and O2 flow through the anode and cathode nozzles respectively with He
as a buffer gas. Oxidizer temperature was found to be 650 ± 20 K and fuel temperature at 600 ±
20 K. Ignition times were observed while changing the fuel mass fraction and fixing a strain rate
of 400 1/s (equivalently a residence time of 2.5 ms), the oxidizer mass fraction, and the driving
frequency.
While cycling the fuel mass fraction both upwards and downwards between XF = [0.1, 0.35],
different configurations of fixed oxidizer mass fractions were used, with XOa = 0.34 and XOb =
0.62, as shown in fig. 3.6. In the case of XOa = 0.34, fig. 3.6a, raising the fuel concentration
(solid markers) did not have a measurable OH emission until XF = 0.265, at which point emission

55

changed abruptly along with a visible flame, demonstrating ignition. Decreasing the fuel (hollow
markers) showed the opposite behavior at XF = 0.2, with visible emission disappearing and OH
emission dropping quickly to noise levels. The second case, with XOb = 0.62, both directions of
changing fuel concentration resulted in remarkably similar curves demonstrating an ignition and
extinction limit of XF = 0.09, shown in fig. 3.6b. A middle case of XO = 0.55 shows shifted limits
with an ignition limit of XF = 0.14 and extinction limit of XF = 0.1. With fixed fuel and oxidizer
mass fractions, similar curves can be created as a function of repetition frequency with a cutoff of
1k Hz, or else radical species begin to quench between pulses.

3.4

KGMf overview with PAC

With NSD systems representing the most appropriate platform for global model studies, this section covers the specifics of the system of equations, considered mechanisms and experimental
conditions replicated.

3.4.1

KGMf analytic additions

In order to adequately describe the kinetics of NSD experiments in the literature, terms for gas
heating during discharge, heat flux lost to the combustor walls, and flow induced species flux
need definitions for inclusion in the KGMf global model. Modeling results from Popov [130],
corroborated by experimental data, demonstrated that gas heating during air plasma discharge
predominantly occurs during the preliminary dissociation of electronically excited states of oxygen
molecules. These states come from either electron impact reactions or via quenching of higher
excited gas species. Further, they showed that over a wide range of reduced electric field values
(E/N), the amount of discharge power involved in heating the gas corresponds to the energy lost
during excitation of electronic degrees of freedom. This agrees with adiabatic theory whereby the
gas heating rate follows from an energy balance. With inspiration from the work by Flitti et. al [47],
the discharge power absorbed by electrons, Pabs from eq. 2.4, distributes between internal potential

56

energy of major species, Pchem , and translational degrees of freedom for the electrons, Pe , and
major species,Pg ,
Pabs = Pe + Pg + Pchem

(3.1)

In a similar short notation, the two temperature equations from the previous chapter can be rewritten as
Pe = Pabs − Sei

and Pg = −Pchem + P˙proc

(3.2)

where Pchem corresponds to enthalpy changes from chemical reactions, P˙proc is the heat flux from
the adiabatic discharge process, and Sei is the second term in eq. 2.4 for energy lost from electron
impact reactions. Solving for P˙proc , the gas heating term from electrical discharge is equal to
the energy lost by electrons in impact processes. The implementation of this term in the KGMf
conservation equations takes advantage of the common-sub-expression algorithm and negligibly
affects evaluation time.
Heat flux to the walls of the combustor is defined through classical conduction heat transfer
theory,
λ (Tg )(Tg − Tw )
Q˙ =
Λ2

(3.3)

where Tw is the ambient wall temperature, λ (Tg ) is the thermal conductivity of the wall material,
and Λ is the characteristic length of the chamber cross section. In a similar approximation, species
flux from gas flow is treated with linearized boundary conditions of inflow and outflow velocities,
Vin and Vout , such that the flow flux for species k in eg. 2.2 becomes
out
in
(G − L)k = nin
k Vin A − nk Vout A = nk (Vin −Vout )A − 2nkVout A

(3.4)

where Vin and Vout are inflow and outflow velocities , A is the cross sectional area of the combustor
chamber, and nin
k is the density of species k within the inflow gas mixture. This assumes a laminar
flow discharge region and relies on spatial uniformity of the flame, commonly demonstrated in
NSD experiments.

57

3.4.2

Mechanisms

In modeling with a volume-averaged framework, a major emphasis is placed on the evolution of
chemical processes as opposed to other physical interactions. As such, the validity of the chemistry
mechanism used is of the utmost importance. Mentioned previously, several attempts have been
made to compile detailed kinetic models for different gas mixtures [91, 3, 150], however, their
validity comes into question under different pressures, temperatures, and electric field strengths.
In the case of gas phase reactions with excited species, many rate coefficients are approximated as
no experimental measurements exist. Electron impact reactions are further plagued by the need to
convolve reaction cross sections with an EEDF. Cross sections are notoriously difficult to measure
experimentally and calculation methods typically do not extend to molecular targets. Reactions
with unknown cross sections are typically provided as experimental fits of modified Arrhenius
rates, however, the EEDF of the measured system is cooked into the result and not extractable.
These rates may be the best option for a given reaction, as not including a pathway would be
worse than misrepresenting it, but the context of their measurement needs consideration when applied to another system. System dependent sensitivity analysis can be used to determine which
rates predominantly determine evolution; specifying which reactions must be represented accurately. Delving deeper, uncertainty quantification methods are used to estimate possible error in
simulation results coming from reaction rate uncertainty bounds.
Further complicating matters, having cross section data for a given reaction does not guarantee
validity, since the rate is also dependent on the EEDF. Outside of particle tracking software, the
EEDF can be found in a variety of methods; such as solving the two-term Boltzmann equation
approximation as a function of electron impact processes and the reduced electric field [60] or
using an isotropic-collision Monte Carlo code to represent electrons [107]. Ignoring the fact that
assumptions taken within most Boltzmann solvers (two-term approximation) do not allow for large
exchanges of mass in a system, their accuracy is also limited by a lack of data. By definition, the
Boltzmann solver should have access to every electron impact reaction cross section, which is not

58

always available. Secondly, due to the computational cost of evaluating the Boltzmann equation,
it is common to evaluate the EEDF as a pre-processing step instead of a complete coupling with
the energy and continuity equations. This reduces EEDF and/or Te dependent rates to constant
values, simplifying the mechanism evaluation at the expense of physical behavior. The majority of
Boltzmann solvers available were not written with combustion in mind, resulting in few warnings
about imposed assumptions. This is not a problem in systems evolving electrically without mass
transfer (e.g pure argon discharge at low pressures such that Ar+
2 can be neglected), but in the
context of non-equilibrium discharges in combustion, the mass fractions change considerably while
electric field densities remain relatively constant. A complete coupling of the Boltzmann equation
solver would provide the most physical results, with a severe increase in the overall evaluation
time.
Of primary interest to this work is the role of plasma-induced excited species in the ignition
process and their susceptibility to variations in the EEDF. Validation studies were performed with
both H2 and CH4 in inert (O2 -Ar) and synthetic (N2 -O2 -Ar) air, however, this analysis is composed
of H2 -O2 -Ar cases. This is done for three reasons: 1) the hydrogen mechanism contains a higher
percentage of electron impact rates defined with cross section data, 2) a large amount of the plasmainduced combustion kinetics are attributed to electronic states of oxygen and nitrogen species
(atomic hydrogen is still a major chain carrier), and 3) as discussed in appendix A.3, the inclusion
of carbonic species greatly increases the computational intensity of global model evaluation.

3.4.2.1

Argon buffer gas

Argon is chosen to be the neutral buffer gas for two reasons: 1) as one of the cheapest noble gases,
argon is commonly used in combustion experiments and 2) argon has a relatively low ionization
energy (∼ 15.7 eV) compared to lighter inert gases and is used as an electron source. Remaining
∗
relatively inert even in excited states, the inclusion of argon brings five species (Ar, Ar+ , Ar+
2 , Ar

Ar∗∗ , the latter two are excited 4s states with Ar∗∗ representing metastable states and Ar∗ the other
two states) and eleven reactions with explicit dependence on an argon species. Six of the eleven
59

Figure 3.7 Electron impact cross sections for Argon from the Morgan LXcat database [119].
reactions have cross sections from the LXcat Morgan database [119], shown in fig. 3.7. Though
not included in this mechanism, systems with ≥ 80% Ar and high field strengths can experience
dissociation of O2 through quenching with excited argon,
Ar( 3P0 , 3P2 ) + O2 −−→ O( 3P) + O( 1D, 1S)

(R.17)

Even without explicit dependence, argon behaves as a collision partner in any reaction with arbitrary impact species.

3.4.2.2

H2 - O2

The H2 -O2 combustion mechanism embodies the core data component of presented analyses and
was created from a variety of sources depending on cross section data availability. A collection
of rate data provided by Igor Adamovich [1], predominantly sourced from Kossyi et. al [91] and
Itikawa et. al [72], which was combined with the GRI-Mech 3.0 [144], giving the former priority if
both databases contained the reaction. This mechanism was then supplemented with cross section
data from the Morgan database [119], replacing other reactions if not defined by cross section data.
60

Figure 3.8 Electron impact cross sections for Argon from the Morgan LXcat database [119].
This combination results in a mechanism with 55 distinct species and 231 reactions. Out of these,
106 reactions are defined with cross section data.
Experimental and modeling work has demonstrated that a primary avenue of PAC/PAI is through
pathways involving excited molecular states of oxygen [142, 148]. Predominately produced through
electron impact,
O2 + e −−→ O2 (a1 ∆g ) + e
O2 + e −−→ O2 (b1 Σ+
g )+e

(R.18)
(R.19)

singlet oxygen species, O2 (a1 ∆g ) and O2 (b1 Σ+
g ) are involved in a good number of dissociation
and branching reactions. With extremely small threshold energies, these singlet states have a huge
overlap with the convolved EEDF, resulting in fast reaction rates. Using a 0.7 m reactor at ∼ 800
K and 10 Torr, with constant discharge power and gas flow of ∼ 25 m/s, Smirnov et. al showed
that a 5% increase of O2 (a1 ∆g ) in the input flow resulted in an induction zone reduction by a factor
of 5 (effectively related to the ignition time) [142]. This is attributed to the creation of ground and
excited states (e.g. O(1 D)) of atomic oxygen
O2 (a1 ∆g ) + e −−→ O + O( 1D) + e
61

(R.20)

Figure 3.9 O generation from gas phase reactions. Reactions involving M collision species are
scaled assuming 100 Torr for matching1units.
O2 (b Σ+
(R.21)
g ) + e −−→ O + O + e
which can drastically alter ignition time with high enough electron and singlet densities. Cross
sections for reactions R.18, R.19, and R.20 are plotted in fig. 3.8 for comparison with argon.

3.4.2.3

Adding N2

In fuel-air mixtures, the dominant avenues of atomic particle production are dissociation of molecules
from electron impact or quenching of electronically excited molecular nitrogen states [129]. This
is partially due to these reactions having extremely short timescales,
N2 (C3 Πu ) + O2 −−→ N2 (ν) + O( 3P, 1D) + O( 3P, 1D)
N2 (a

1 Σ− ) + O
2
u

−−→ N2

(ν) + O( 3P) + O( 3P)

(R.22)
(R.23)

of 3-4 ns and tens of ns respectively. Unlike argon however, nitrogen comes with the formation of
NOx products and significantly more reactions. In the used model, adding nitrogen results in 56
new species and 453 new reactions (68 of which are defined with cross section data).
In low initial temperatures, NO and NO2 have a catalytic effect on hydrogen and hydrocarbon
fuels as demonstrated by Takita et. al [160]. In a plasma system, NO can be generated efficiently
62

through the following pathway
N2 + e −−→ N( 2D) + N + e

(R.24)

N( 2D) + O2 −−→ NO + O( 3P, 1D)

(R.25)

2
N2 (A3 Σ+
u ) + O −−→ NO + N( D)

(R.26)

where N2 (A3 Σ+
u ) is likely also created from an electron impact reaction. The catalytic effect was
found to decrease ignition time as the pressure increased, however, the pathway breaks down at
higher temperatures relative to the effect from O radicals.

3.5

Modeling

The overarching goal of this work, relative to PAI, is to both analyze kinetic characteristics of
deviations from classical combustion chemistry and to validate global models in a multi-physics
scenario. The relative computational efficiency of volume-averaged modeling allows for computationally cheap parametric studies. This efficiency comes at the cost of a high-data dependency
and needing a large amount of physics handled through approximations. Alluded to previously, for
most gas mixtures, there does not exist a kinetic mechanism which the PAC/PAI community collectively agrees will function in all regimes. Even with regards to classical combustion, updates to
mechanisms are regularly published, either for handling new regimes or concentration-dependent
reaction pathways [173, 138]. This is due in part to the regime-dependent experimental nature of
acquiring reaction rates of specific processes.
The validity of global models relies on a combination of: 1) a kinetic mechanism describing
all relevant reactions, 2) a dominance of volumetric processes in conditions of reasonable homogeneity, and 3) models for spatial characteristics such as gas flow, excited species diffusion, and
thermal conduction to better match experimental conditions. The selection of focusing on experimental NSD systems was made to meet the second criterion, as a system reliant on turbulence
will not be well described by a global model. This work predominantly focuses on the first criterion, with the latter used as tuning parameters for matching experimental conditions. The intention
63

Figure 3.10 Plot of EEDF, eq. 2.9, with different values of X and values of Te.
being to match result variation that occurs under parameter changes, not necessarily to yield 1:1
result matches with experimental conditions.
The intention with this work is to study the plasma effects on ignition times. Two H2 -O2 -Ar
cases are considered with total system pressures of 10 and 100 Torr respectively. First, changes
to reaction rates and integrated results, relative to EEDF selection, are quantified with regards to
ignition times. Secondly, following the work of Smirnov et. al [142], an investigation of the effects
of electronically excited oxygen species on the ignition process is performed. This experiment used
a DC discharge for generation of oxygen singlet states, which although not the most power efficient
manner of PAI, has been demonstrated to contribute to both flame-holding and ignition [101].
Finally, the dynamic characteristics of NSD are considered relative to constant power deposition
in the ignition process.

64

(b) Dissociation of excited O2 (a1 ∆g )

(a) Excitation of O2

Figure 3.11 Effect of EEDF selection on different electron impact rates: excitation of molecular
oxygen on the left, O2 → O2 (a1 ∆g ), and dissociation through quenching on the right, O2 (a1 ∆g ) →
O + O(1 D).
3.5.1

Rate Variation From EEDF dependence

Beyond the relatively static nature of gas-phase reactions, the dependence of electron impact processes on EEDF selection can result in non-intuitive transitions of dominant pathway chains (e.g.
the pair of reactions R.18 and R.20 producing atomic oxygen species). The evolution of an EEDF
from a Maxwellian form results in a lower number of higher energy electrons, as shown in fig. 3.10,
involved in collisional processes. This in turn modifies convolved reaction rates in a non-linear
manner due to artifacts in individual reaction cross sections. The overall shifting of reaction rates
is shown in fig. 3.11; however, changes relative to equilibrium distributions highlight the non-linear
transitions occurring, fig. 3.12.
The collective effects of all electron driven reactions being altered is only realizable through
integration of the system of equations as the strength of each reaction is also determined by the
product of reactant species densities, which in turn are based on the evolution of the system. However, with the transitions that occur in the relative difference between rates, fig. 3.12, it is in principle possible to configure a system with a specified EEDF such that certain reaction rates have a
larger effect. The severity of the transitions that occur imply that selection of the wrong EEDF can
result in different pathways being treated preferentially.

65

(a) X2 - X1

(b) X5 - X1

Figure 3.12 Effect of EEDF selection on different electron impact rates, shown through relative
error from rates calculated with a Maxwell-Boltzmann distribution (X=1).
3.5.2

O2 (a1 ∆g ) Injection and Ignition Time

As a validation case, the KGMf is applied to the experimental and modeling work by Smirnov et. al [142]
on the effect of excited singlet oxygen on the combustion process. Molecular oxygen has two lowlying electronically excited singlet states, O2 (a1 ∆g ) and O2 (b1 Σ+
g ), which are predicted to reduce
the ignition time of H2 -O2 combustion [147]. The experiment is composed of a discharge cell
with oxygen throughput piping into a heated flow reactor with hydrogen injection. The system is
maintained at 10 Torr with reactor walls heated to 780 K and a gas mixture of H2 :O2 ::5:2. A dc
discharge, operated with I= 10 mA and V= 4.2kV, is used to excite a portion of injected oxygen
into the O2 (a1 ∆g ) and O2 (b1 Σ+
g ) states prior to mixing with hydrogen. At this temperature and
pressure, the system is already in the auto-ignition region and will eventually ignite without introducing excited states, however reduction of the ignition time relative to O2 (a1 ∆g ) concentration
provide insight into the dominant reactions involved. Absorbed molar input power was measured
to be Ws = 107Jmmol−1 resulting in the singlet state reaching number densities of (1.13±0.3)1019
m−3 . At a flow velocity of 25 ms−1 , it was found that O2 (b1 Σ+
g ) depleted rapidly while traversing
the drift tube due to quenching down to the O2 (a1 ∆g ) state. Ignition times were found as a function
of the singlet state mole fraction. Without excitation, the system achieved ignition in 12.9 ms at a

66

distance of 51 cm from the injection point of the hydrogen gas. As the concentration of O2 (a1 ∆g )
was increased to 1.2%, 4%, and 6% of the total gas, the ignition time reduced to 6.8 ms, 4.7ms,
and 4.2 ms respectively.
This system was replicated using the KGMf in the original pressure of 10 Torr with a few
considerations; namely handling a characteristic of the reported flow and supplementing a small
amount of power to maintain the excited oxygen species. Without adding any external power,
the effective electron temperature equation remained extremely unstable and regularly trended to
just above 0 eV. While the averaged electron energy is expected to be small in this regime, the
zeroing is unphysical behavior brought out by unbalanced momentum transfer and attachment
processes. Though momentum transfer reactions should be the dominant driving factor for lowering electron temperature, electron attachment to atomic oxygen, forming O− , was found to have
an overestimated threshold energy. Using the experimental geometry of a cylinder with a diameter
of D = 5.2mm and length of L = 0.3m, a background input power of 0.01W is required to maintain
the electron temperature and density to reasonable values. As reported in the original modeling of
the experiment, it is necessary to account for inhomogeneous flow from low pressure. At pressures
of 10 Torr, the flow in the system has a parabolic velocity profile relative to the cross section of the
flow regime [142]. Using the linear flow model in the KGMf, symmetric inflow and outflow from
the system is tuned in order to match the baseline case of pure ground state oxygen being injected.
Tuning resulted in an optimized flow of ∼ 32m−1 to match the ignition time of the unsupplemented
oxygen case.
For proof of concept of the previous section and calibration of the O2 (a1 ∆g ) injection effect,
base line temperature dependent ignition times were calculated for the system in the case of both
a Maxwellian EEDF and a Druyvesteyn EEDF, fig. 3.13. In the relatively low pressure of 10 Torr,
the combustion chemistry is dominated by other chain reactions and the effect of electrons is negligible. With an introduction of energy into the system, either through pumping in excited species or
stronger discharge powers, electron impact processes become more dominant. Figure 3.14 presents
the full reaction term (rate and product of reactant species) evaluated in time while integrating
67

(a) Maxwellian EEDF

(b) Druyvesteyn EEDF

Figure 3.13 Ignition time in 10 Torr of H2 :O2 ::2:1 with equilibrium and non-equilibrium EEDFs.

(a) 0% and 0.0W

(b) 6% and 0.0W

Figure 3.14 Rates of radical producing reactions with (a) and without (b) singlet delta oxygen
supplementation. Without external power, the rates appear quite similar but differ in scale. The
total system number density is used for the M collision species.
the global model. Reactions which appear unimolecular are dependent on the electron density,
which in the case of quasi-neutrality is the charge weighted sum of all ionic species in the system.
Though no power is being explicitly passed to the global model, initializing O2 (a1 ∆g ) at such a
high concentration is effectively passed potential energy. The drop in reaction terms in fig. 3.14a
is explained by a drop in electron density from an over production of O− . This is counteracted in
fig. 3.14b by the higher initialized O2 (a1 ∆g ) concentration and its reaction with O− releasing the
captured electron.

68

Figure 3.15 Changes to ignition time due to injection of 6% O2 (a1 ∆g ) at 10 Torr.
With the global model calibrated with background power and flow velocity, ignition times
were found for O2 (a1 ∆g ) initialized, and flowing in, at specific molar fractions of the total feed
gas, fig 3.15. The ignition times found through integration with the KGMf were slightly slower
than reported values. Though remaining within a factor of two of the experimental results under
6% singlet delta oxygen flow, the initial decrease in time is not as extreme. This may be explained
by a lower simulated electron density than necessary to drive the O radical production chain. As the
concentration of O2 (a1 ∆g ) is increased, potential energy can trickle through to electrons after the
excited species undergo quenching. Ignition times were also found for 1.0 W and 10.0 W power
depositions for comparison purposes. The further reduced ignition time in the high power cases
would physically be followed by breakdown resulting in a thermal ignition. The nearly unchanged
ignition time in the 0 W case of fig. 3.15 is explained by the unbalanced drop in Te from O− over
production.
The effect of introducing O2 (a1 ∆g ) in the initial flow is better seen when comparing the evaluated rate terms when using 0.1 W and 0.01 W of external power, fig. 3.16, to the un-powered cases,
fig. 3.14. The rate at which the system consumes the fuel and oxidizer has also increased, resulting
69

(a) 0% O2 (a1 ∆g ) and 0.1W DC power

(b) 6% O2 (a1 ∆g ) and 0.01W DC power

Figure 3.16 Effect of O2 (a1 ∆g ) injection on the rates of reactions involved.
in a fast change of reaction rate dominance to matching the depletion of reactant species.

3.5.3

Nanosecond discharge

While the dc discharge provides insight into the underlying processes that are occurring during PAI,
nanosecond discharge operates on a time scale that directs almost all of its energy into electronically excited species. These short kicks to the system are absorbed by the electrons and transmitted
to the gas through collisions. For validation of the KGMf NSD dynamics, a global model is made
to match experimental conditions by Yin et. [178]; hydrogen-air flows initially heated to 400-500
K, at pressures of 40-150 Torr and equivalence ratios of φ = 0.5 − 1.2 are ignited with 25 ns pulses
at a repetition frequency of 40 kHz. Of interest to this work are the experimental measurements of
power coupling, ignition delay time, and the effect of the stoichiometric ratio on the PAI process.
Analogous to the dc discharge study, baselines are made with a Maxwellian EEDF, X=1, and
normal ignition time behavior is recognized, fig. 3.17. In the unpowered case, the ignition time
vertically asymptotes along an isotherm of the used mechanism. In the over-powered case, 10W,
ignition time plateaus until an unphysical overgrowth of charged species occurs. The electron number density remains within reasonable ranges, however they are only defined through the quasineutrality relationship with charged species. The cause is attributed to an imbalance of energy

70

(a) Maxwellian EEDF

(b) Druyvesteyn EEDF

Figure 3.17 Baseline results at 100 Torr with different EEDFs.

(a) 1000K

(b) 1200K

Figure 3.18 Pulsed power 60 kW 100 Torr. Electron impact reaction terms jump in magnitude due
to temporary increases of Te during pulses.
being deposited into the system and over production of electronic species which are not well represented in the gas phase mechanism. In the Druyvesteyn case, ignition times remain relatively
unchanged, due to the higher pressure driving gas phase reaction dominance, however run stability
limited lower temperature results.
Monitoring NSD effects is significantly more time and computationally intensive than modeling systems with constant or gradual external forces. In order to resolve smooth 25 ns impulses,
time steps need to be on the order of 1 ns; resulting in 107 calls to the integrator evolving the global
model system of ODEs. With a ‘complete’ set of air-fuel chemistry, aside from the function call

71

time, the resulting set of ODEs is extremely stiff. Paired with the intensity on the NSD only operating on a subset of the ODE species, fig. 3.18, requires using stiff integrators which operate with
iterative sub-stepping. An alternating time step scheme is used such that time steps just before and
during a pulse are 10−10 s, increasing in size as the system equilibrates before the next pulse. This
can reduce the number of time steps by over an order of magnitude in lower temperature regimes.
Following from the work by Yi et. al, a 25 ns pulse operating at a repetition frequency of
ν = 40kHz and having a Gaussian envelope, is used to ignite hydrogen with synthetic air. In
this regime, the contributions from nitrogen species are negligibly different from any other large
collision partner, justifying the simplification. Experimentally measured pulse energy coupling to
the plasma at 100 Torr is ∼ 1.6 mJ, corresponding to a 60 kW discharge. Using this discharge,
KGMf ignition times were found to agree with the asymptotic behavior of reported values from
cooler temperatures and lower pressures. In lower temperatures, the ignition time is found to be
five times faster than in the baseline case. However the upper temperature bound of 1600 K is well
beyond the auto-ignition limit, and is ambivalent towards both NSD and introduction of O2 (a1 ∆g )
states, fig. 3.19. Overall agreement is achieved, however the fine tuning to initial conditions and
integrator parameters required in lower temperature simulations hampers the ability to parameter
scan as effectively.

3.6

Conclusion

In this chapter, motivation and fundamentals of PAC are discussed followed by validation examples
with two experimental works. Relatively good agreement is found between KGMf and experimental results, however system calibration/tuning is typically required for simulation times longer than
a few milliseconds. This is primarily associated with unbalanced system initialization and integration time step instabilities. In the first case, the effect of generating O2 (a1 ∆g ) prior to combustion
has a dramatic effect in reducing ignition times. In the second validation case, extrema of the
experimental ranges are matched with KGMf results continuing on to higher temperature regimes.

72

Figure 3.19 Pulsed Ignition time at 100 Torr with a 60 kW 25 ns pulse. Contains modeling results
for stoichiometric combustion, φ = 1, with varied EEDF, and examples of lean and rich behavior.

Figure 3.20 Species densities in 1000K NSD ignition.

73

Figure 3.21 Pulsed 1000 K Temperatures.
The method proposed for handling integration difficulties is a double-mechanism configuration. A reduced chemistry model is integrated between NSD pulses, while a full chemistry mode
is used during the pulse and cool down period thereafter. This is done as the electron temperature,
shown in fig. 3.21, drops to thermal equilibrium between pulses. Complications arise at the interface between integrated models as some electronically excited species have significant lifetimes
in these pressures. In order to significantly drop the evaluation time, the reduced model consists
of only non-excited gas species. The long lifetime of O2 (a1 ∆g ) otherwise requires that the full
chemistry model is used until all excited species have dropped to negligible densities. Transitioning from the reduced model to the full model requires estimation of excited species density ratios
relative to ground states. Since the cool down period can be longer than the pulse repetition rate,
shown in fig. 3.20, excited species are re-initialized at non-negligible values resulting in significant
differences as a function of chosen excited species ratios and switching time between models.

74

CHAPTER 4
RARE GAS METASTABLE LASERS WITH KGMF

4.1

Introduction

Since the conceptual advent of triggered stimulated emission in the middle of the 20th century,
laser technology has expanded to cover a huge range of incorporated physics. Symmetries found
in nature and similarities between systems have led to the development of multiple types of lasers
operating on widely differing physical principles. While the principles of operation differ, the underlying fundamentals remain relatively unchanged. Requirements of different laser applications
have fueled regular development of new technologies, opening the doors to new process optimization and physical manipulation of matter.
Optically-pumped rare gas lasers, first presented in 2012, represents a new hybrid of laser and
plasma technology with novel avenues of laser characteristic optimization [64]. Relying on a combination of plasma discharge and optical pumping creates a system requiring a deep understanding
of underlying electron kinetics. Despite experimental success, numerical study of these systems
remains an immature field with several outstanding questions. This work aims to study the nonequilibrium interaction of plasma kinetics and laser operation. Starting from an overview of laser
development and fundamentals, this work then derives an intracavity intensity model for coupling
with previously discussed KGMf continuity equations. After validation of the chemistry and laser
models, the effect of non-equilibrium plasma on laser properties is analyzed for optimum operating
conditions.

75

4.2
4.2.1

General Laser History and Concepts
Laser Background

Lasers have become common instruments with consumer, industrial, and research applications.
The process of light amplification by stimulated emission of radiation, first given the acronym
“laser” in a notebook by Gordon Gould on the 13th of November in 1957, relies on theory which
took 40 years prove. Albert Einstein’s ground breaking work, On the Quantum Theory of Radiation [44], laid the theoretical foundations for stimulated emission by combining the concept of
probability coefficients for radiative reactions (now called Einstein coefficients) with Max Planck’s
law of radiation, earning him the Nobel Prize in Physics in 1921. Rudolf Ladenburg validated
the concept experimentally in 1928 [90], but the idea of stimulated emission as an amplifier was
not proposed until Valentin Fabrikant’s 1939 doctoral thesis [108]. Bringing the laser theoretical
piece required, Alfred Kastler proposed optical pumping as a stimulation method in 1950 [84] and
demonstrated it experimentally in 1952 [22], earning him the Nobel Prize in Physics in 1966.
Overlapping with the final theoretical framework of lasers, Joseph Weber proposed the concept of a microwave amplifier by stimulated emission of radiation (maser) in 1952 at the Institute
of Radio Engineers Vacuum Tube Research Conference in Ottowa, followed by the first device
demonstration by Charles Townes and his graduate students in 1953 . While novel, this original
implementation could not demonstrate continuous-wave (CW) operation as a result of difficulties
from the two-level excitation structure used. Concurrently in the Soviet Union, Nikolay Basov and
Aleksandr Prokhorov built a CW maser by incorporating a three-level excitation through pumping to a higher energy level of ammonia gas [7]. These combined and independent works earned
Townes, Basov, and Prokhorov the 1964 Nobel Prize in Physics.
The origins of the laser process, specific from masers, are shrouded in controversy as a result
of patent claims and the Iron Curtain of the USSR. While accepted that the earliest reference of
the laser acronym comes from a notebook of Gordon Gould, no single person can be credited for
the overall concept. Gould’s coining of the term came after discussions with Townes, who at the
76

time worked at Bell Labs with Arthur Schawlow. Shortly thereafter in 1958, Bell Labs applied
for a patent on an optical maser, or laser, and submitted theoretical calculations to Physical Review [136]. The patent office denied Gould’s application in 1959, instead awarding it to Bell Labs
in 1960, resulting in several decades of legal battles between the two entities. While Gould won the
patent case in 1987, people familiar with the situation believe that instead of malicious idea theft,
Townes simply had no interest in keeping the idea a secret [106]. Back in the UUSR, Prokhorov
and Basov suggested in 1955 that a viable method for obtaining positive gain could come from optical pumping of a multi-level system. Further complicating matters, Fabrikant submitted a request
for a patent in 1951, equivalently the laser process without the name, which was granted in 1959.
Due to a lack of records on changes to the application, no one knows if Fabrikant took inspiration
from works published in the mid 1950s [15]. Interestingly, Fabrikant developed this work without
progressing through the maser process, but both maser and laser systems had been realized by the
time the patent was published.
Regardless of initial conception, the first laser device was operated by Theodore Mainman
in 1960 [113], using a solid-state ruby crystal as the gain medium. Ali Javan William Bennett
and Donald Herriott demonstrated the first gas laser (helium and neon) in 1961 [77], followed by
Robert Hall demonstrating the first semiconductor laser (gallium arsenide) in 1962. These three
types of lasers, with the addition of liquid/dye lasers, make up the four major classes as defined by
their gain medium. Before describing benefits and disadvantages of these different class, the next
section will cover basic laser fundamentals and the metrics by which to judge different systems.

4.2.2

Laser Fundamentals

The premise of both maser and laser operation relies on having two energy states, as shown in
fig. 4.1a, between which radiative transfer is allowed, and a method to over-populate the higher
energy state to produce a population inversion. Considering a system composed of two-particle
states, |1 and |2 , with population densities N1 and N2 and energies E1 and E2 , a higher energy
particle interacting with an field of frequency ν21 = (E2 − E1 )/h produces stimulated emission at
77

(a) 2-level scheme

(b) 3-level scheme

Figure 4.1 Two and three level laser schemes.
a rate of
dN2
= −B21 ρ(ν)N2
dt

(4.1)

releasing a photon of the same frequency. Similarly, the absorption rate is expressed as
dN2
= B12 ρ(ν)N1
dt

(4.2)

with radiation density ρ(v) and Einstein B coefficients. Einstein showed that the proportionality
constants are equal, B21 = B12 , and therefore the difference between state populations directly
correlates to net photon flux from stimulated emission,
dN2
= B21 ρ(v) [N2 − N1 ] = B21 ρ(v)∆N
dt

(4.3)

where N2 − N1 = ∆N > 0 defines a population inversion. If contained within a chamber having
walls at temperature T and allowed to equilibrate without a source of external power, the relative
number of particles in each energy level comes from Boltzmann statistics,
p21 N1
=
= e(E2 −E1 )/kT
p12 N2

(4.4)

with Boltzmann constant k and transition probabilities, p21 and p12 , for 2 → 1 and 1 → 2 respectively. Since E2 > E1 , the ratio is never less than unity and the higher probability of a downward
transition, p21 , accounts for the difference in populations. With the addition of a radiation source
that triggers the transitions between the two states with probability S and again allowing the system to reach a (non-thermal) equilibrium, fluxes between the two states must balance, resulting in

78

a new relation
N1 S + p21
=
N2 S + p12

(4.5)

which again is greater than unity as p21 > p12 . In the special case where S

p21 , a close enough

equivalence of N1 and N2 results in a saturated state. Though attainable with microwave transitions,
as hν is on the order of 10−5 eV while kT at room temperature is ∼ 0.026 eV, the high energy of
photons in the visible spectrum, hν ∼ 1−3 eV, requires an excessively large input power in order to
surpass thermal equilibration, destroying efficiency of the amplifier. Without a third pumping level
with quenching to the laser level, the energy pumping mechanism contributes to the S term, populating both states and working against itself. In Towne’s original maser experiment, the two energy
states of ammonia were spatially separated with an electric field resulting in a beam formation of
predominantly higher energy particles.
In order to demonstrate an efficient laser with optical photons, the lasing scheme requires at
least one more energy level, as shown in fig. 4.1b. For a three level system undergoing pumping
for the 1 ↔ 3 transition and allowed to equilibrate, the equivalent set of balanced flux equations
are
N1 (p12 + p13 + P + S) = N2 (p21 + S) + N3 (p31 + P)
N2 (p21 + p23 + S) = N1 (p12 + S) + N3 p32
N3 (p31 + p32 + P) = N1 (p13 + P) + N2 p23
with pumping transition probability P. Defining the total particle density as N = N1 + N2 + N3 and
assuming saturated conditions such that P
N1 = N3 =

N
2+λ

(p31 , p13 ), the steady state solutions become

and N2 =

Nλ
2+λ

where

p + p32 + S
λ = 12
p21 + p23 + S

(4.6)

and substituting the solutions into the population inversion condition produces
E2 −E1

p12 e kT

E3 −E2

− 1 > p23 e kT

79

−1 .

(4.7)

Finally, the laser amplification is directly proportional to the population inversion
N2 − N1 = N 1 −

3
2+λ

(4.8)

which becomes larger as λ increases. Maximization of the population inversion and other laser
characteristics greatly depends on the gain medium aside from effects not considered in this simplified case: transition probabilities p∗ are temperature dependent, S may alter other transitions in
the system (e.g. 2 ↔ 3), and all of these considerations assume CW behavior which is not required
for pulsed operation. However, except in the special case of low transition energies, three or more
levels are required for efficient operation [155].
Choice of the gain medium largely defines laser properties and the selection thereof is extremely situationally dependent. Though originally described as a “solution seeking a problem”
by Maiman [114], the ability to generate a narrow, tight bandwidth, and intense beams of light has
found applications in industrial, military, science, and medical fields. In being an amplifier system, a major metric of importance is the amount of power returned for a given pumping intensity,
however, specific applications may have emphasis on different laser properties. For example, spectroscopy requires a tunable monochromatic light source whereas a medical application may only
need a single wavelength and a military application demands only high-power scalability without
a requirement as to output bandwidth.
At the most general level, four categories classify different gain media based on their phase:
solid-state, dye, semiconductor, and gas lasers. Solid-state lasers rely on a crystal structure where a
“dopant” maintains the population inversion and lasing is induced by optical pumping in the form
of a arc lamp, or more recently, diode-pumping. As a result of the nature of the gain medium,
solid-state lasers tend to have long energy storage times and high thresholds for damage making
them appropriate for high power applications. Dye lasers use an optically-pumped liquid gain
medium composed of a carbon-based, soluble, stain which yields a broad output spectrum with
a high Stokes efficiency (difference between pump and output wavelengths) and high gain. The
latter implies both a low pumping threshold and high incident noise. Semiconductor lasers, or

80

laser diodes, are electrically pumped and formed by a p-n junction in the gain material. Thanks
to decreasing manufacturing costs, laser diodes are becoming commonplace in consumer products
(eg. barcode reader) and as pumping mechanisms for other gain materials. Finally, gas lasers use
an electrical or optical discharge in a gas mixture and resonant cavity in order to produce coherent light. Applications range from kilowatt-class industrial CO2 lasers for cutting and welding
applications to He-Ne frequency-tunable lasers.
A variety of methods exist for populating the excited states of gas lasers and all rely on a sealed
cavity containing mirrors forming a Fabry-Pérot resonator. The CO2 laser relies on a electrical
discharge to populate vibrationally excited states of nitrogen which excites CO2 through collisional
processes. Chemical lasers rely on a chemical process (e.g. C2 H4 -NF3 yielding fluorine radicals)
to introduce energy to the system. Excimer lasers rely on electrical discharge and high pressure
to create noble gas molecules, excimers, which undergo spontaneous and stimulated emission to
a repulsive ground state followed by dissociation into unbound atoms. Ion lasers similarly use
an electrical discharge to generate a population density; however, the energy cost is hampered
by a large amount of produced waste heat. Metal-vapor lasers, such as the copper vapor laser,
have extremely narrow linewidths making them excellent candidates for Raman spectroscopy. The
generation of copper vapor requires temperatures on the order of 1200 K, limiting containment
options and energizing pumping for both vapor generation and metal ionization.

4.3
4.3.1

Motivation of this Work
DPAL

Another gas laser of primary interest to this work is the concept of diode-pumped alkali vapor lasers
(DPAL), coined by Krupke with a patent in the early 2000s [93]. Alkali metals present an attractive
gain media due to having a single valence electron in the outer-most occupied shell however their
reactivity can pose problems with regards to experimental configurations and longevity. The idea,
as shown in fig. 4.2, comes from creating a three level laser system with the ground state, n2 S1/2 ,
81

Figure 4.2 Electronically excited levels of a 3-level DPAL system with ground state n2 S1/2 and
further excited levels n2 P1/2 and n2 P3/2 [94]
and two low-lying electronically excited states of an alkali metal, n2 P1/2 and n2 P3/2 , with pumping
transition line D2 and lasing transition D1 . This model has an extremely high quantum efficiency
due to the similarity in energy of the two states. As one of the first tested gain media, CW operation
of an optically pumped Cs vapor laser was demonstrated by Jacobs, Rabinowitz, and Gould in the
early 1960s [73, 134], but did not see continued attention due to a lack of reliable pumping sources.
In the following four decades, multiple experimental works demonstrated the effects of different
buffer gases on the collisional relaxation from n2 P3/2 → n2 P1/2 with different alkali metals. At the
turn of the century, the concept of producing a continuous population inversion through pumping
along the D2 transition was demonstrated, leaving questions as to which pumping source would
spectrally match the centerline of the D2 transition while efficiently coupling to the narrow atomic
transition bandwidth.
The development of broad-area laser diodes (AlGaAs or InGaAs) coupled with techniques
for reducing their spectral bandwidth from a few nm to a fraction of a nm answered the first
question. The second question is handled by the buffer gas accompanying the excited levels [94].
Collisions with the buffer gas transform the D2 line into a Lorentzian lineshape. Assuming equal
transition line-strengths, the Lorentzian transition will have a slightly smaller peak cross section
and a wider cross section than with the Gaussian signal, allowing for off-center pumping. Though

82

tight bandwidth overlap would be the most efficient pumping method, the off-center pumping
allows for flexibility of the pumping hardware. The same effect presents itself on the D1 transition,
however this mitigates spectral hole burning and facilitates single frequency DPAL operation.
Complications still remain with regards to energy scaled DPAL systems, independent on the
reactivity of the alkali vapor. A considerable amount of pumped power is absorbed as heat by the
gain medium resulting in thermal lensing. In a non-flowing system, the buffer gas is responsible
for the thermal conduction to a cooling surface which is limited by the low thermal conductivity
of He. While unique designs have mitigated thermal effects [182], typical DPAL systems employ
a heavier buffer gas species or convective cooling. Secondly, while the He buffer gas is used
to broaden the pumping absorption lines, it has a low cross section for alkali metal spin-orbit
relaxation [65]. Increasing the gas pressure is not a viable means of accelerating the relaxation
due to degradation of the beam quality from a lower system opacity. Instead, the addition of small
chain hydrocarbons (e.g. CH4 ) results in an increase to the rate of spin-orbit relaxation without over
quenching the 2 P levels [69]. Hydrocarbon species also induce the potential for formation of metal
hydrides, or “laser snow,” which result in carbonaceous deposits damaging window transmittance
and mirror reflectivity [82].

4.3.2

OPRGL

Motivated by the difficulties in handling the reactive alkali species, Han and Heaven demonstrated
in 2012 that a similar three level system exists in electronic excitations of noble gases, resulting
in the concept of an optically pumped rare gas laser [64]. Completely circumventing the problem
of chemically reactive species, the proposed scheme mimics the electronic structure of the DPAL
scheme by shifting up in energy. Under electric discharge, the np5 (n + 1)s configuration of noble
gases easily populates excited metastable levels (n + 1)s[3/2]2 and (n + 1)s [1/2]0 [154]. Akin to
alkali metals, strongly allowed optical transitions exist connecting these metastable states to the
np5 (n + 1)p excitation manifold [92]. Using the lower energy metastable state, (n + 1)s[3/2]2 ,
as the bottom level, the laser scheme then pumps to populate the (n + 1)p[5/2]3 state, which
83

Figure 4.3 Excited level diagram for an argon oRGL system.
undergoes spin-orbit relaxation to the (n + 1)p[1/2]1 state followed by lasing emission back to
the metastable state, fig. 4.3. Due to the similar electron excitation orbitals, the excited noble gas
transitions between s and p manifolds are analogous to the DPAL transitions, exhibiting similar
transition probabilities, upper state radiative lifetimes, lasing wavelengths, pressure broadening
coefficients and ionization potentials. Along with the chemically stable nature of noble gases,
heating is no longer required to populate the system (replaced by electrical discharge) and other
heavy noble gases are available for the spin-orbit relaxation [82] without strong interactions to the
rest of the system.
The use of electric discharge as opposed to thermal heating results in both a more efficient
method of populating the bottom level of the laser system as well as longer hardware lifetime.
Experimental DPAL systems require heating to roughly 120-145 ◦ C and deposited thermal energy
is absorbed by every species. By comparison, electric discharge deposits the majority of its energy
into the electron population and the excited metastable level, which forms the bottom laser level,
is the most readily excited of the argon states. Non-equilibrium discharges further offer the opportunity to focus electrical energy predominantly into species of interest. This can be done at room
temperature without requiring the need for heat dissipation.
These first experimental studies demonstrated pulsed lasing with argon, krypton, and xenon
metastable states, using helium as a buffer gas due to the high ionization threshold, operating at
pressures of 0.2 - 2 bar [64]. With a repetition frequency of 10 Hz, a 10 ns discharge pulse is
84

Figure 4.4 Experimental time-resolved traces of 5p[1/2]1 ↔ 5s[3/2]2 gain in Kr, with colors denoting total pressures of 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.4, 1.6, 1.8, and 2 bar [64].
followed by a 10 ns optical pulse with a tuned gap time depending on the partial pressures. The
discharge pulse operated at a minimum voltage of 18 kV with an electrode separation of 2.5 cm
in the form of a commercial excimer laser cavity which was not optimized for RGL transitions.
The tuned optical pulses had a linewidth of 30 GHz, beam diameter of 0.4cm, and approximately
0.25 mJ of energy. Initially tested with a Kr/He system, the population inversion heavily depended
on fast collisional transfers between the two excited p manifold states. In order to explore this
dependence, a fixed partial pressure of 27 mbar of Kr was combined with He to form total pressures
ranging from 0.2 - 2 bar, with results shown in fig. 4.4. The optimized delay between discharge
and optical pumping varied from 5-25 µs.
Pressures above 0.6 bar showed nonzero gain; increasing until a saturation pressure of 1.4 bar
as a result of competition between radiative and collisional relaxation of the top laser level. Other
rare gas mixtures exhibited similar behavior with the additional case of argon lasing without a
helium buffer above pressures of 1 bar. Efficiency characteristics are reported as being similar for
each mixture, having a photon conversion efficiency of 13% and a metastable number density of at
85

least 4 × 1019 m−3 .
Following this work, Demyanov et. al performed initial kinetic studies on CW operation of
the RGL scheme with argon [39]. Assuming the metastable population is maintained by a glow
discharge, laser efficiency, gain, and power output were found as functions of argon number fraction, total pressure, pump intensity, and system size by solving for the steady-state operation of
a reduced reaction mechanism. Considering the four 1s excitations and first three 2p excitations,
where the laser system cycles between 1s5 → 2p9 → 2p10, electron impact reactions are only
considered for excitation of the 1s states, spontaneous emission is only considered along the pump
and lase transitions, collisional relaxation reactions are only included for the species involved in
the laser process, and rates are considered independent of temperature. While the simplification
of the mechanism allows for analytical solutions, the excluded state-to-state collisional excitations
were neglected with an argument of no available rates in the literature and that experimental results
at similar pressures show a larger 1s5 density relative to other states in s manifold, implying fast
and complete energy transfer to the 1s5 state.
From steady-state analytic expressions, Demyanov et. al demonstrated that the discharge power,
required for maintaining the 1s5 state population, and laser output power both increase with pumping intensity. Beyond saturation of the population inversion, the total efficiency of RGL operation
also increases with pumping intensity. The simplified reaction mechanism causes estimates of the
required discharge power to be considered as an upper bound of the physical system. The rates of
both metastable and excimer production are proportional to total argon concentration in the system.
Finally, the threshold pump intensity was found to be on the order of tens of watts with a saturation
intensity an order of magnitude higher and quadratic in pressure. As such, lasing in this scenario
should be achieved relatively easily.
Later in 2013, Han et. al demonstrated another RGL system using CW optical pumping and 1
µs pulsed electrical discharges [63] in an argon helium system. Stable operation was achieved with
a pulse repetition frequency of 1 kHz with pressures at or above 0.5 bar. Whereas the previous work
relied on high instantaneous pump intensities (∼ 0.2 MW cm−2 ) to achieve pulsed laser output,
86

Figure 4.5 Experimental time-resolved traces of 5p[1/2]1 ↔ 5s[3/2]2 intensity in Ar, with colors
denoting total discharge voltages of covering [-1 kV,-2 kV] in steps of 100V [63].
Han et. al showed that the power range of a CW diode pumping source need only be on the order of
5-10 W. The DC discharge, operated between -1-2kV, used electrodes separated by 0.5 cm within
15 cm long glass cell of diameter of 5 cm. A mirror and output coupler were placed on either end
of the glass cell forming a 35 cm long plane parallel resonance cavity. The cavity is end-pumped
with a diode laser centered at 811.754 nm (wavelength of the 4s[3/2]2 → 4p[5/2]3 transition)
with a maximum power of 8 W and approximate beam diameter of 200 µm. Laser power output
was found to be dependent on cavity alignment, proving that lasing was occurring as opposed to
emitting a laser-induced fluorescence. Unlike with the previous system, lasing with pure argon
was limited to pressures of ∼ 90 mbar due to arcing occurring at higher pressures.
Dependence of the lasing system on the discharge voltage was carried out using a mixture of 40
mbar and 420 mbar of Ar and He respectively. While glow discharges were created for lower discharge voltages, lasing was not observed until the voltage was below -1.2 kV. With lower voltages,
both laser peak intensity and laser pulse duration increased; the latter to times of ∼ 0.15µs. With-

87

Figure 4.6 Optical gain along microplasma vs pumping intensity; 2% Ar - 98% He at 769 Torr
with flow of 3.7 mmoles/s and total microwave power of 9W, data adapted from [135].
out optical simulation, Ar metastables have a lifetime of several µs. Lasing threshold was found to
be significantly lower than required intensity from the original RGL experiment at ∼ 11 kW/cm2 .
A concluding remark implied a key requirement for further development is the identification of a
discharge capable of continuous production of the (n + 1)s[3/2]2 near atmospheric pressures.
In attempting to answer the question of continuous discharge without breakdown, Rawlins et. al
demonstrated CW optical gain at atmospheric pressure using a linear micro-discharge array driven
by a 900 MHz microwave [135]. Unlike traditional atmospheric discharge limited by short-term
breakdown, the linear arrays, credited to Hopwood et. al [181, 176], provide a stable and spatially
uniform CW discharge at driving frequencies of ∼ 900 MHz. Though operated at low power (under
30 W), the array resonators have high field strengths and create electron and argon metastable
densities on the order of 1019 m−3 in a 2% Ar-98%He mixture. Arrays composed of 15 micro-strip
resonators, each producing a microplasma, resulted in a 19 mm long microplasma in a width of
25µm from ∼ 9W of discharge power. Optical excitation was end-pumped along the microplasma,
without a resonance cavity, using a CW laser with up to 1 W of power.

88

Noting localized gas temperature increases up to 600 K, Rawlins et. al measured gains G0 up
to ∼ 0.7 cm−1 with pump saturation occurring ∼ 2 kW/cm2 . The temperature increase marks
a larger effect of Doppler broadening, reducing the efficiency of single-frequency pumping but
allowing for a wide pump bandwidth. Unfortunately this also implies that the pumping laser will
also populate other 2p levels. Gain was found to vary significantly in space, with the maximum
gain of ∼ 1 cm−1 located downstream from the discharge region. Following the discussion by
both Han et. al and Rawlins et. al, kinetic modeling needs to include radiative loss from 2p10 →
1s4 as this makes up over 20% of the radiative decay from the 2p10 state. Disagreeing with
their own steady-state analysis (simplified mechanism) and counter to experimental observations
in pulsed operation by Han et. al, Rawlins et. al report relatively small populations of the 1s4
state in CW discharge conditions. This was explained by a temperature dependence on the overall
transition from 1s4 → 1s5 with up to a 10x variation between Han’s 300 K experiment and the 600
K measured by Rawlins’ team. This suggests that scaling metastable lasers to higher gains in CW
operation might require gas heating beyond room temperature. Limits of this analysis and results
were attributed to unknowns in the rate coefficients at higher temperatures.

4.4

Derivation of Lasing Model

In order to describe these laser systems, the KGMf requires additional terms for both the lasing and
pumping processes. Acknowledging the similarity between DPAL and RGL systems, an analytic
model for CW DPAL operation is implemented alongside KGMf continuity and conservation equations. Building from initial modeling attempts of DPAL systems assuming longitudinally averaged
number densities [8, 10], Hager and Perram developed an analytic two-way averaged plane-wave
model assuming single frequency optical pumping of a three-level laser system [62, 61]. Since
extremely narrow band pumping is not required for DPAL systems due to broadening of the D2
absorption line, Zameroski et. al extended this model to account for broadband pumping and spectral overlap of the pump and transition lines in 2011 [180]. While the analytic model presented

89

Figure 4.7 Assumed two-way passing schemes for induced lasing emission, adapted from [62].
only accounts for the three species involved in the laser process, the continuity equation for intracavity intensity and terms corresponding to pumping and stimulated combine with collisional
dynamics to form a complete reaction mechanism for the system.
With spontaneous emission and collisional rates previously defined with regards to KGMf implementation, fine-structure mixing, pumping, and stimulated emission terms are discussed herein.
Consider the continuity equations of an example model with the three laser states: n1 corresponds
to the number density of the lower laser level, n2 the upper laser level, n3 the relaxing level above
n2 :
dn3 Ω(t) ˜ mix
=
− K32 ± ∑ O3
dt
hν31
dn2
mix − K˜ stim ±
= K˜ 32
∑ O2
21
dt
dn1
stim − Ω(t) ±
= K˜ 21
O
dt
hν31 ∑ 1

(4.9)

where Oi represents particle flux from collisional and spontaneous emission losses described in the
first chapter, Ω(t) is the pump power absorbed per unit volume, and K˜ mix and K˜ stim are respectively
the spin-orbit mixing and stimulated emission flux terms. The fine-structure mixing term is written
as


ij 
K˜ imix
j = kmix ni − n j

90



∆Ei j
k
e B Tg 

(4.10)

ij

with i > j, Boltzmann constant, kB , gas temperature, Tg , mixing rate between states i and j, kmix ,
and ∆Ei j = hνi j is the energy difference between states. The latter two have constant values
specific to the gas composition and states involved.

4.4.1

Lasing Model

As opposed to tracking the radiation density within the gain region, an intracavity circulating intensity model is used to account for multiple passes within a resonance cavity in a volume averaged
context. Extending from sec. 4.2.2, the population inversion in the case of different degeneracies,
gi for each level, is written as
g
∆n21 = n2 − 2 n1
g1

(4.11)

which is paired with the stimulated emission cross section for transition i → j
σi j (ν) =

c2 Ai j
8νi2j π

Gi j (ν)

(4.12)

where c is the speed of light, νi j is the frequency of the transition line, Ai j is the Einstein coefficient, and Gi j is the spectral lineshape function. Assuming a Lorentzian lineshape function and
inhomogeneous broadening in a gas [34],
Gi j (ν) =

∆νi j
1
2π (ν − νi j )2 + (∆νi j /2)2

⇒

Gi j =

1
∆νi j

(4.13)

the dependence on ν is removed, relying only on the collisionally broadened line width ∆νi j . With
species dependent Doppler broadening rates as γnb , for each broadening species, nb , the linewidth
is written as
∆νi j =

Tg
γ n
298[K] ∑ b b

(4.14)

Assuming small signal gain resulting in stimulated emission from n2 → n1 , the governing equation
for stimulated emission and corresponding solution are:
dI
= σ21 (ν)∆n21 I(z)
dz

⇒

91

I(z) = Iin eσ21 (ν)∆n21 z

(4.15)

with gain defined as the product of σ21 and ∆n21 . Originally developed by Raymond Beach [9]
with extensions by Hager et. al, the averaged two-way intracavity intensity, Ψ, is defined as the integrated intensities of left and right traveling plane waves, IL− and IL+ , in the gain medium, averaged
over the gain length, lg ,
Ψ=

1 lg + −
IL + IL dz
lg 0

(4.16)

Substituting the solution in eq. 4.15 for each directional intensity in eq. 4.16 and using single pass
intensities from fig. 4.7 for Iin , the expression for Ψ comes to
Ψ=

lg
1 + lg σ ∆n z
I2
e 21 21 dz + I2−
eσ21 ∆n21 (lg −z) dz
lg
0
0

= (I2+ + I2− )
Ψ = (I2+ + I2− )

eσ21 ∆n21 lg − 1
σ21 ∆n21 lg
eα21 − 1
α21

(4.17)

where α21 = σ21 ∆n21 lg is introduced for convenience. The gain medium is assumed to be in a
container with two windows facing each other and having transmission factor t. A perfect mirror
and an output coupler, of reflection factor r, form the resonance cavity by encasing the ends of the
container. With a reference point of I1+ inside the cavity behind the coupler, the intensity is tracked
through a round trip of the cavity, picking up a factor of eα21 through the gain medium and t
through the glass. With relationships between intensities at different points, fig. 4.7, and boundary
conditions across the coupler, four identities are used to relate Ψ to the output lasing intensity Ilase :
I2+

= trI4−

I2−

I4−
= α
te 21

Ilase = (1 − r)I4−

I1+ =

r
I
1 − r lase

(4.18)

Combining these relationships with eq. 4.17 allows for the laser output intensity to be written in
terms of the average intracavity intensity
α21 (1 − r)teα21
Ilase =
Ψ
1 + t 2 reα21 [eα21 − 1]

(4.19)

Depending on the mirrors and glass used, a considerable amount of power may be lost from scattering and transmission loss from imperfect windows. Summing up the difference between intensities
92

at each window transition,
Iscat = (I1+ − I2+ ) + (I3+ − I4+ ) + (I1− − I2− ) + (I3− − I4− )
= Ilase

r(1 − t)
1 + eα21 t + t 2 + t 3 eα21
1−r

(4.20)

arrives at an expression for the scattered laser intensity. The conservation equation for Ψ incorporates gains induced from a round trip in the cavity along with a spontaneous emission noise term
to start the laser oscillation,
dΨ
Ψ
n c2 σ21 hν21
= (rt 4 e2σ21 ∆n21 lg − 1)
+ 2
dt
τRT
lg

(4.21)

where τRT is the residence time in the cavity [62]. Finally, dΨ/dt couples with the system of
equations 4.9 through the stimulated emission rate
stim = σ ∆n Ψ(t)
K˜ 21
21 21
hν21

(4.22)

which remains dependent on species densities n2 and n1 as well as intracavity intensity.

4.4.2

Pumping Model

Hager et. al assumed single-frequency (narrowband) pumping which is not quite representative of
optical pumping at atmospheric pressures due to collisional broadening of both D1 and D2 absorption lines. Zameroski et. al extended the model to account for broadband pumping and spectral
overlap of the pump and D2 transition [180]. Analogously to eq. 4.17, the two-way averaged pump
intensity inside the gain medium is written as
+ + I− )
ϑP = (IP0
P0

eα31 − 1
α31

(4.23)

with the new shorthand term α31 for the D2 transition. Assuming the rear mirror is not perfect with
regards to the pumping wavelength (therefore having reflectivity rP ) and the pump has a spectral
distribution of f (ν) and an amplitude time envelope of P(t), the tracked round trip pumping intensities form the identities
+ = t P(t) f (ν)
IP0
P

− = t 3 r P(t) f (ν)eα31
IP0
P P

and
93

(4.24)

analogous to eq. 4.18, with glass windows having a different transmission coefficient for the pumping wavelength, tP . Plugging the relationships back into eq. 4.23 yields the pumping intensity
ϑP (t, ν) =

f (ν)P(t)tP α
(e 31 − 1)(1 + tP2 rP eα31 )
α31

(4.25)

however this needs to be integrated across the frequency range of the pump for the total energy
absorbed by the upward stimulated emission,
Ω(t) = −

∞

g
σ31 (ν) n3 (t) − 3 n1 (t) ϑP (t, ν) dν
g1
−∞

(4.26)

It can reasonably be assumed that the pump has a normalized Gaussian spectral distribution centered around the transition frequency ν31 , such that
4 ln (2)
exp
π∆νP

fP (ν) =

−4 ln (2)

ν − ν31 2
∆νP

(4.27)

where ∆νP is the spectral pump width (FWHM). Using the pump spectral distribution and the
broadened stimulated emission cross section, eq. 4.26 can be integrated to yield
Ω(t) =

P(t)
tp
lg

∆νP (eα31 − 1)(1 + tP2 rP eα31 )

(4.28)

Finally, imposing a pump power envelope in time
E
PG (t) = P
AP

4 ln (2)
exp
π∆tP

−4 ln (2)

t − t0 2
∆tP

(4.29)

closes the system of equations defined in eq. 4.9 with terms for pump energy, EP , beam area, AP ,
and temporal width, ∆tP .

4.5

Modeling Results

This work was motivated by two main goals: 1) validation of the KGMf implemented laser model
and 2) determining if it is possible to pump the laser scheme electronically as opposed to using
another laser. While diode laser technology has both increased in performance and decreased in
cost, optically-pumped RGL systems still require two different powering schemes for operation.
94

Figure 4.8 Excited levels of argon considered in the model. Levels involved in the laser process
are in red. Data comes from the LXcat Biagi database [16].
Of interest is whether a ‘tuned’ EEDF is capable of selectively pumping the gas to maintain a
population inversion without over-pumping the reversed laser transition. The experimental works
by both Han and Heaven [64] and Rawlins et. al [135] are used for validation of the intracavity
laser model. Both systems operate at atmospheric pressure using transitions between excited states
of argon gas, fig 4.8. Argon was chosen as it was the only reported rare gas to demonstrate lasing
without having a collisional partner. Cases with helium as a buffer gas are also included. Lacking
data for interactions between excited species of argon and helium, only ground and ion species of
helium are considered. The systems differ in the specifics of the optically pumping source and the
method of maintaining the metastable state population. The first case uses a nanosecond discharge
pulse while the latter uses an RF discharge.

4.5.1

Pulsed Discharge

Matching the original RGL example by Han and Heaven [64], the first case assumes a gain length
of 60 cm inside of a resonance cavity. The discharge pulse is operated at 18kV across a distance

95

(a) Species response to pumping

(b) Laser output vs. delay time

Figure 4.9 Optimized delay time between discharge and optical pulses depends on species equilibration.
of 2.5 cm for 10 ns. Optical pumping follows a pressure dependent optimized delay time of 5-25
µs. Centered along the D2 transition of 811.75 nm, the pumping laser has a bandwidth of ∆ν = 30
GHz and deposits 0.25 mJ of energy in a 10 ns timespan. With a beam width of 0.4 cm, peak input
intensity is assumed to be 3.14 kW/cm2 . A repetition frequency of 10 Hz was used experimentally
to separate the effect between pulses, allowing for single pulse simulation. An optimized delay
time of 2.3 µs between discharge and optical pumping was found, however the relative increase in
output intensity is nominal, fig. 4.9.
Introducing a buffer gas to the system provides a faster quenching between the pumped and
lasing states of argon. A ratio of He:Ar::1:1 was found to provide the greatest increase in output
intensity. Keeping the gas temperature fixed removes any pressure dependent (explicit) terms,
therefore comparison simulations varied total pressure to maintain the same argon number density
in each case. With laser gain proportional to the weighted difference of species concentrations,
varying the pressure provides an equal opportunity for each system reach a peak intensity
Overall good agreement was found between reported species densities in the original work [64]
and temporal profiles of pulses match well with measurements in the following publication by
Han et. al [63]. Differences from reported total output intensity are attributed to efficiency factors
in the conversion from intracavity intensity and output intensity

96

(a) Output Intensity

(b) Species during optical pulse

Figure 4.10 Effect of adding helium to the system over a range of ratios. The line type in the right
plot transitions from solid to dotted with the same labels as the colors in the left plot.
4.5.2

Continuous Discharge

The second validation case was in modeling the micro-discharge array presented by Rawlins
et. al [135]. The experimental system is composed of 15 micro-strip resonators generating plasmas
in a linear array. Modeling assumptions used here however, considers a single resonator and lasing
region. This was done primarily to avoid the spatial variation of gain reported in the experiment.
A secondary consideration was to avoid the amplifying effect lasing points releasing photons into
neighboring lasing regions. The geometry of the strip was approximated as a cylinder with radius of R = 0.6 mm to match discharge volume. The gas mixture is composed of 2% Ar in a He
buffer. Pumping intensities were reported to be as high as ∼ 10 kW/cm2 with a spectral width of
∆ν = 2GHz.
Unlike the previous example, this system is not within a resonant cavity, bringing into question
the validity of using the two-way averaged plane wave model. Instead, the averaged model was
reduced to a single-pass configuration [45] with a gain length of 1 cm as reported in the measurements of laser extraction in the original work.
Gain is observed in modeling results in the saturated limit with matching microwave power and
species densities. Pulse lengths of optical stimulation are chosen such that population inversion,
and therefore output intensity, maintains steady state operation. While in the experimental work
97

(a) Relevant species

(b) Laser intensities

Figure 4.11 Species densities and involved intensities as a result of repetitive optical pumping.

(a) Intensity

(b) Laser output vs. delay time

Figure 4.12 Equilibration of species and intensities with 50% helium buffer gas.
the argon is at steady-state densities prior to optical pumping, fig. 4.12, demonstrates pumping
during the equilibration phase of the integrator. The relative strength of the pumping process is
significantly stronger by comparison to the electron impact excitation rates.

4.5.3

EEDF driven

With the start of this work, there was a question as to whether or not non-equilibrium kinetics could
be used to optimize the laser process. This idea hinges on the ability to aid pumping to the upper
level or quenching to the laser level by using a non-Maxwellian EEDF. This was shown to not be

98

Figure 4.13 Electron impact cross sections involved in the laser process from LXcat Biagi [16].
energetically favorable as a result of the sharpness of EEDF required for a non-trivial change to
laser kinetics. In the ideal scenario, an EEDF and electron temperature would be chosen such that
excitation to the (n+1)p excited levels occurred more rapidly than to the lower metastable level.
As can be seen in the electron impact cross sections of argon, fig. 4.13, there is no point of overlap
between the cross sections for the bottom laser level (1s5) and the two 4p states of the laser system
(2p9 and 2p10). This implies that an electronic discharge could not pump the rare gas laser to
the point of achieving a population inversion. Instead, the modified EEDF and system electron
temperature required would need a narrow bandwidth such as to maximize production of the 2p9
and 2p10 excited levels while minimizing excitation to higher 4p and 4s states. This condition
would allow for a more efficient optical pumping scheme. However, the gains from including an
electron beam as well as a optical-pumping source would be significantly less efficient than using
a more powerful diode-laser.

4.6

Conclusion

The implementation of the laser model in the KGMf has been validated, using acceptable tuning
parameters, with regards to rare gas metastable laser operation. Validation cases demonstrating

99

Figure 4.14 Same cross sections fig. 4.13, zoomed in to look at relevant lasing states.
the relatively new process use both constant and pulsed discharge power to generate a plasma and
optical pumping to maintain a lasing population inversion. Though the model takes into account
spectral overlap in pumping intensities due to broadening, this work motivates the use of spectral
tracking within the global model by tracking all spontaneous and stimulated emission processes in
order to have a picture of radiation density in the system.
Studies were performed on pulsed and continuous electric discharge effects on populating the
lower laser level in argon. Pulsed operation of the electric discharge was experimentally performed
in order to avoid localized breakdown effects from the power requirements of exciting a large
volume. Continuous discharge operation was stably achieved at atmospheric pressure using a
micro-array microwave discharge. When comparing the two discharge methods, peak intensities
were found with optical pumping during the electric discharge as opposed to when optical and
electric pulses were offset. A population inversion is readily achieved under the combination of
both discharge and optical pumping. Overlapping the optical pulse with the electric discharge onset
resulted in the largest gain, explained by faster equilibration of the 1s5 state and correspondingly
higher pump efficiency, though near steady state population inversion was observed with optical
100

pumping and steady state species densities.
The optimal EEDF for the laser process in argon would require a bump-on-tail EEDF with
an extremely narrow bandwidth in order to avoid overpopulation of higher energy 4p levels. The
difficulty arises from the proximity of 4p state threshold energies and similar cross sections. While
this EEDF could be achieved with electron beam sources, it is not possible to directly drive a
population inversion with an electric discharge. The equivalent effects of an optimized EEDF can
be achieved with higher powered optical pumping in a more energy efficient manner.

101

CHAPTER 5
CONCLUDING REMARKS AND FUTURE WORK

This dissertation work revolved around the development of a general-purpose global modeling
framework and applications to non-equilibrium plasma chemistry. A good deal of effort was spent
on guaranteeing generality through reaction data flexibility and the ability to add new physics effect
modules (e.g. intracavity intensity). Verified relative to published global model implementations
of both plasma and combustion chemistry, studies were then performed on the non-equilibrium
chemistry found within plasma assisted combustion and the dynamics of optically pumped lasing
between plasma generated excited states of argon after validating each respective global model
with experimental results.
Mechanisms were compiled for plasma assisted combustion of hydrogen and methane in a
variety of oxygen, argon and nitrogen compositions; incorporating a mixture of fitted and EEDF
dependent rates. These mechanisms were used to study the effect of the electron energy distribution on alternative pathways for combustion reaction chains. In order to agree with experimental
measurements, a linear transport model was implemented within the KGMf to incorporate spatial
effects from gas flow. Following analysis of supplemental excited oxygen, O2 (a1 ∆g ), on the H2 -O2
combustion process, the effect of continuous RF discharge on electron-impact chain-branching reactions was studied. Though decreasing the ignition time, a significant portion of discharge power
is wasted as thermal energy within the gas. Nanosecond discharges were shown to reduce ignition
time through radical species production without significantly contributing to the gas temperature.
Integrating different physics into the KGMf, a stimulated emission and pumping model was
developed to study the reaction kinetics of rare gas metastable lasers. Validation of the model was
carried out with both continuous and pulsed electric discharge for populating the bottom laser level.
Driving the population inversion through electric discharge was found to be inefficient relative to
diode laser pumping with regards to direct excitation pathways.

102

5.1

Future KGMf

With any code development or research effort, there is always more that can be implemented or
studied. Plans are in place for future development of the KGMf to: incorporate more self-consistent
behavior, become a platform for sensitivity analysis and uncertainty quantification of reaction data,
and act as a chemistry engine within a spatial framework including transport equations. While
some of this work is being performed by new students within our research group (effectively inhouse), we hope that the open-source nature of KGMf will inspire collaboration for improvements
beyond present design, as has been the case with other PTSG open-source codes.

5.1.1

Self-consistent EEDF

The major limitation of the KGMf is the lack of a self-consistent EEDF. In order to remedy this, a
Boltzmann equation solver or Monte Carlo method could be coupled with KGMf time integration.
This development is multi-faceted and non-trivial in the context of maintaining global model computational efficiency. As opposed to the present version, where EEDF dependent reaction rates are
pre-calculated assuming distribution dependence on a system variable, such as Te , solutions for the
EEDF result in scalar values of rate coefficients for considered densities and system parameters at
a fixed instance in time. With regards to efficiency, a brute force implementation of solving for the
EEDF on each time step reduces the speed advantages of global models. Instead, this implementation requires recognition of the multi-physics time scales involved and a prescription for updating
the EEDF selectively. In the same manner that reaction data can be cycled for comparison, interaction between the KGMf and Boltzmann solvers allows for testing of the downstream effects of
EEDF assumptions. This allows for a comparison of solution methods; for example, checking the
validity of two-term Boltzmann equation approximations in the context of large collisional mass
fluxes.
It should however be noted that solutions to the Boltzmann equation are just as data dependent
as the global models relative to electron impact reactions. Mentioned in the PAC chapter, EEDF

103

solutions are typically found for an initial gas mixture when fields are known to be relatively
fixed. This assumption is questionable when concentrations of gas species no longer match initial
conditions, however electron impact cross section data is often not available for radical and excited
molecular species. Analysis of EEDF variation from considered gas species would lend credence
to the ability to pre-compute the EEDF for a given configuration.

5.1.2

Uncertainty Quantification

With functional assumptions based in theory, the reaction mechanism required for plasma systems
relies on a mix of gas and plasma phase reaction data. These measured and calculated data describe
extremely non-linear systems benefiting from uncertainty quantification (UQ), beyond standard
sensitivity analysis, for uncertainty driven result fluctuation estimates and reaction mechanism
optimization.
In order to limit the UQ parameter space, dominant reaction pathways and parameters are determined using the PumpKin [116] package along with preliminary sensitivity analyses. Analogously
to classical combustion UQ in the literature, uncertainty in the rate parameters, including cross
section data and parameterized EEDF shape, is quantified for forward propagation. Focusing on
ignition delay times, result confidence intervals are calculated using polynomial chaos expansion
and non-intrusive pseudospectral projection tools from the Chaospy [46] package with a wrapper
over KGMf batch methods.

5.1.3

Spatial - Fluid Code

While global models are inherently light with respect to computational intensity, the KGMf was
developed with evaluation time as a dominant focus. Though partly related to developer impatience, the run time optimizations motivate the inclusion of KGMf models as exportable chemistry
engines in spatial resolving codes. The intention here is to create hybrid methods with active chemistry fluids. One major application is related to studying the interaction of plasma-chemistry with

104

turbulence effects in plasma assisted combustion systems. It has already been demonstrated that
both radical and electron species densities can play a large role in the combustion process; spatial
variation of these species can therefore lead to kinetic hot spots analogous to thermal hot spots
from spark ignition. Localized heating, for example at sheath edge regions, and consequent active
radical generation and transport may be key to improving PAC performance.

5.2

Future PAC/RGL

Each of the numerical research studies can be readily expanded for more thorough analyses of
the underlying chemistry and non-equilibrium behavior. In both cases, global models can be created to match different experimental regimes in the literature. Development of the self-consistent
EEDF also provides the opportunity for numerical PAC and RGL studies which have yet to be
demonstrated.
With validated models for both PAC and RGL chemistry, large scale parametric studies can be
performed with regards to system parameters and resulting ignition time or laser gain respectively.
In both cases, sensitivity analysis methods should be used to account for result variation due to
poorly known input parameters. Once the operating regimes are understood well enough for stable
integration, result-search methods can be used to isolate pathways which can be optimized.

5.2.1

PAC

Plasma assisted combustion and ignition represents a breakthrough in the understanding of combustion kinetics. If developed beyond its infancy and incorporated into consumer products, PAC
technology has the potential to dramatically change global resource requirements and power structure. In order for this to come to fruition, a firm understanding of the underlying fundamentals are
necessary for future application development.
The study presented on PAC can be significantly built upon with new regimes, data, and selfconsistent physics. Of primary interest are the non-equilibrium effects on ignition below the au-

105

toignition threshold and contributions of nitrogen species when considering air as the oxidizer
source. While our studies focused on the additional contribution that plasma can provide to the
combustion process, efficiently achieving ignition below the autoignition threshold fundamentally
distorts classical combustion concepts. Vibrational levels of nitrogen have also been shown to
dramatically alter reaction pathways through collisional quenching [1] motivating the inclusion of
more species. It is still not understood whether excited nitrogen or oxygen molecules are the main
drivers behind PAC behavior.
First and foremost, the combustion mechanism including nitrogen species should be validated
for low temperature operation. Once dominant pathways have been found through sensitivity analysis, we propose implementing nonlinear control in order to determine domineering parameters of
the ignition time. If the KGMf has been coupled to the Boltzmann equation solver, a supplementary study should be performed on variation of the EEDF as a function of radical species densities
to check the validity of pre-computing EEDFs. If not yet finished, a more exhaustive set of testing
assumed EEDFs should be performed.

5.2.2

RGL

The metastable rare gas laser was originally demonstrated in a number of noble gases while our
analyses were performed assuming only argon and helium (as a buffer species). With the intracavity laser module validated, recommended studies revolve around testing behavior of other
rare gases and determining if other excitation pathways exist so as to remove the need for optical
pumping. Paired with the use of other lasing species, the effects of the helium buffer and possible
diatomic gas contamination should be explored. Parasitic electron energy absorption by rotational
and vibrational states could alter the discharge driven laser excitation pathways.
Few quenching rates for the upper laser level are available with temperature dependence or
measured at the high micro-array operating temperature. Sensitivity analysis to these rates should
be performed in order to determine the relative importance of having an accurate model for quenching at higher temperatures. Expanding upon the implemented laser model, which only takes pho106

tons of interest into account, an extension would be to include approximate spectral radiation density within the system to determine overlaps between stimulated emission pathways. This would
allow for laser schemes with more than three levels along with an accurate description of broadening effects.

107

APPENDICES

108

APPENDIX A
APPENDIX TO KGMF

A.1

Derivation of Gas Temperature

The amount of heat released during a chemical process is governed by the laws of thermodynamics
and directly results in a change in gas temperature. From the 1st law of thermodynamics, the
balance between different forms of energy in the system is written as
du + pdv = dh − vd p = δ q + δ w

(A.1)

where u is the internal energy, v is the specific volume or inverse mass density, p is the pressure,
δ q is the heat transfer from external sources, and δ w is the frictional work. Introducing the species
mass fraction, Yi , the total specific enthalpy of a multi-species system is related to species’ partial specific enthalpies by h = ∑ Yi hi . Using the temperature dependence of the partial specific
enthalpies
T

hi = hi,ref +

Tref

c pi dT

(A.2)

where c pi is the specific heat capacity at constant pressure, hi,ref is the reference enthalpy, and Tref
is the temperature of the reference, eq. A.1 is differentiated to yield
N

dh = c p dT + ∑ hi dYi

(A.3)

i=1

with i summing over all species and c p = ∑ Yi c pi as the total specific heat capacity. In order to
arrive at an equation for the temperature evolution of the system, mass and energy conservation
equations need to be defined. Assuming a system with both mass and heat flux, the species continuity equation is defined as
ρ

R
R
DYi
= −∇ · ji +Wi ∑ νil ωl = −∇ · ji +Wi ∑ νil Ki ∏ n j
Dt
j
l=1
l=1

(A.4)

where ji is the diffusion flux, Wi is the mass of species i, l sums over all R reactions involving
species i, νil is the net stoichiometric coefficient of species i in reaction l, Ki is the reaction rate,
109

and j creates a product of all reactant species number densities. The last term of eq. A.4 is the
chemical mass source term for species i. Ignoring the effects of thermodiffusion, the diffusion flux
is defined as ji = −ρD∇Yi where D is the binary diffusion coefficient.
Defining energy from external sources as qr and the total energy flux as je , the total energy
equation
∂ (ρe)
= −∇ · je + q˙r
∂t

(A.5)

where e is the total energy per unit mass, is turned into a total derivative of the internal energy
ρ

Du
= −∇ · jq − p∇ · v + q˙r
Dt

(A.6)

after dropping dependence on the viscous stress tensor and substituting
je = ρev + pv + jq

(A.7)

for the total energy flux, where v is the mass average velocity, p is the pressure, and jq is the total
heat flux. Using Fourier’s law of thermal conductivity, the total heat flux is written as
N

jq = −λ ∇T + ∑ hi ji

(A.8)

i

where λ is the thermal conductivity and the Dafour heat flux has been neglected. Using a substitution of ∇ · v = ρD(ρ −1 )/Dt, eq. A.6 is re-written to be in the same form as eq. A.1 after having
neglected the δ w term
D
Du
+p
Dt
Dt

1
ρ

=

1
−∇ · jq + q˙r
ρ

(A.9)

Noting the relationship between enthalpy, internal energy, and pressure, h = u + p/ρ, eq. A.9 can
be re-expressed as the enthalpy balance equation
ρ

Dh Dp
−
= −∇ · jq + q˙r
Dt
Dt

(A.10)

In order to solve for the temperature equation, another expression for the total derivative of enthalpy
is found by multiplying all terms in eq. A.3 by ρ/Dt and substituting eq. A.4 for the DYi /Dt term,
resulting in
ρ

N
Dh
DT
= ρc p
+ ∑ hi (m˙ i − ∇ · ji )
Dt
Dt i=1

110

(A.11)

where m˙ i is the last term in eq. A.4 representing the species mass change due to chemical reactions.
Finally, setting enthalpy derivatives in eqs. A.10 and A.11 equal to each other and replacing the
heat flux with eq. A.8 yields the temperature equation
ρc p

N
Dp
DT
=
+ ∇ · (λ ∇T ) − ∑ (hi m˙ i + ji ∇ · hi ) + q˙r
Dt
Dt
i=1

(A.12)

Prior to solving for a volume-averaged form of the temperature equation, it is worth considering
the chemical source term in terms of heat contribution as opposed to mass changes. Writing out
the third term from eq . A.12 and including the definition of m˙ i , the heat release of reaction l, Ql ,
is found through an exchange of summations
N

R

∑ him˙ i =
i

∑

N

R

∑ Wihiνil ωl =

l=1 i=1

∑ Ql ωl

(A.13)

l=1

Removing terms corresponding to spatial variation and including eq. A.13, the temperature equation comes to the general form
R
DT
1 Dp
=
− ∑ Ql ωl + q˙r
Dt
ρc p Dt l=1

(A.14)

It should be noted at this point that eq. A.14 is extremely coupled to the species continuity equations. While the total specific heat, c p , is dependent on the mass fractions, or effectively densities, of all involved species, the summation over heats of reactions includes scaled forms of every
species continuity equation. Though easier to derive the temperature equation in terms of mass
fractions, KGMf has equations implemented in units of species number density. In order to convert eq. A.14 to eq. 2.3, a conversion from mass fraction to molar concentration must be performed,
followed by simplification of enthalpy and specific heat terms taking the NASA polynomial format
into account.

A.2

Control Files for KGMf

Operation of the KGMf is largely defined by two control files. The simulation input file, fig. A.1
defines system parameters, such as external power, pressure, and EEDF, and contains the flags for
111

# G e n e r a l ( d e f a u l t ) s i m u l a t i o n i n p u t f i l e f o r BASIC_PLAS__Ar
#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
# Model S p e c i e s Name
,
Constant Flag (T/ F)
,
I n i t i a l Value
#−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
PRESSURE
,
T
,
760 [ T o r r ]
T_g
,
T
,
1500 [K]
T_e
,
F
,
1 . 2 [ eV ]
EDF
,
T
,
[ EEDF&{TParam ; x v a l = 1 . 0 [ u n i t ] ; w e i g h t = 1 . 0 [ u n i t ] } ]
#POWER
,
F
,
[ c o l f r e q &{Pmethod0 ; E e f f _ e d i t = 2 . e3 [V/m ] ; f r e q = 2 . 5 4 e9 [ Hz ] } ]
POWER
,
F
,
[ p u l s e &{ s q u a r e ; a m p _ e d i t = 50 0[W] ; w i d t h = 1 . e −8[ s ] ; t s t a r t = 5 . e −5[ s ] ; r e p = 8 ;
→ r e p g a p = 1 . e −4[ s ]}&{ c o n s t ; amp = 0 . 0 1 [W] } ]
#POWER
,
F
,
[ p u l s e &{ c o n s t ; amp =1 0 [W] } ]
============================
GEO
,
cyl
,
R 0 . 0 2 5 [m]
,
L 0 . 6 [m]
#GEO
,
cyl
,
R 0 . 0 1 [m]
,
L 0 . 0 0 4 2 [m]
SHEATH
,
None
RATS
,
Ar_g 1 . 0
GP2PLAS
,
BASIC
,
maj : T
,
chrg :T
INIT_VALS
,
u n i t =m^−3 , i n i t 0 = 1 . 4 e12 , i n i t 1 = 1 . 1 e14 , i n i t 2 = 1 . 1 e9 , i n i t 3 = 1 . e12
CUSTOM_VAR
,
t e s t v a l =0.0
CUSTOM_SET
,
c r o s s a n a _ r e s o l u t i o n = [ geo ; 1 . 4 ; 3 . 0 ; 2 0 0 . 0 ]
============================
Ar_g
,
F
,
[ init3 ]
Ar_exc_1
,
F
,
[ init0 ]
Ar_exc_2
,
F
,
[ init0 ]
Ar_exc_3
,
F
,
[ init0 ]
Ar_exc_4
,
F
,
[ init0 ]
Ar_exc_5
,
F
,
[ init0 ]
Ar_exc_6
,
F
,
[ init0 ]
Ar_exc_7
,
F
,
[ init0 ]
Ar_exc_8
,
F
,
[ init0 ]
Ar_exc_9
,
F
,
[ init0 ]
Ar_exc_10
,
F
,
[ init0 ]
Ar_exc_11
,
F
,
[ init0 ]
Ar_exc_12
,
F
,
[ init0 ]
Ar_exc_13
,
F
,
[ init0 ]
Ar_exc_14
,
F
,
[ init0 ]
Ar_pion_0
,
F
,
[ init1 ]
Ar_pion_1
,
F
,
[ init2 ]
,
t h r e s h = 2 0 . 0 [ eV ]
,
degen =1.0
Ar_pion_2
,
F
,
[ init2 ]
,
t h r e s h = 2 7 . 0 [ eV ]
,
degen =1.0
Ar2_pion_0
,
F
,
[ init3 ]
,
t h r e s h = 0 . 0 [ eV ]
,
degen =1.0

Figure A.1 Example simulation input file for an argon plasma.
determining if a species is allowed to evolve or held fixed in time. Database values for threshold
energy and degeneracy can be overwritten through declaration in this file.
The reaction file, fig. A.2, contains the reaction network and linking information for data corresponding to each reaction. The structure of these files is reduced to hash numbers in order to
determine if a model can be recovered from saved data as opposed to fully recompiled. Reactions
under the ‘REACTIONS’ header are included directly from the local database while reactions under ‘EXTRA_REACTIONS’ are added by the user. These additions can be dependent on numerical
data or passed algebraic expressions. The ‘SPEC_MAPS’ section creates species placeholders for
grouped species.

112

SPECIES
Ar_g , Ar_exc_1 , Ar_exc_2 , Ar_exc_3 , Ar_exc_4 , Ar_exc_5 , Ar_exc_6 , Ar_exc_7 , Ar_exc_8 , Ar_exc_9 , Ar_exc_10 ,
→ Ar_exc_11 , Ar_exc_12 , Ar_exc_13 , Ar_exc_14 , Ar_pion_0 , Ar_pion_1 , Ar_pion_2 , Ar2_pion_0 , A r _ e x c _ 4 s
SPEC_MAPS
A r _ e x c _ 4 s = Ar_exc_1 + Ar_exc_2 + Ar_exc_3 + Ar_exc_4
REACTIONS
Ar_g + E l e c t r o n −> Ar_g + E l e c t r o n
===
DATA( EI_CSmomen__list , 0 )DATA
Ar_g + E l e c t r o n <−> Ar_exc_1 + E l e c t r o n
===
DATA( E I _ C S _ _ l i s t , 1 )DATA
Ar_g + E l e c t r o n <−> Ar_exc_2 + E l e c t r o n
===
DATA( E I _ C S _ _ l i s t , 2 )DATA
#...
Ar_g + E l e c t r o n <−> Ar_exc_13 + E l e c t r o n
===
DATA( E I _ C S _ _ l i s t , 1 3 )DATA
Ar_g + E l e c t r o n <−> Ar_exc_14 + E l e c t r o n
===
DATA( E I _ C S _ _ l i s t , 1 4 )DATA
Ar_g + E l e c t r o n −> A r _ p i o n _ 0 + 2∗ E l e c t r o n
===
DATA( E I _ C S _ _ l i s t , 1 5 )DATA
Ar_exc_4 ~> Ar_g
===
DATA( R_A__const , 1 6 )DATA
Ar_exc_2 ~> Ar_g
===
DATA( R_A__const , 1 7 )DATA
#...
Ar_exc_5 ~> Ar_exc_3
===
DATA( R_A__const , 4 5 )DATA
Ar_exc_5 ~> Ar_exc_4
===
DATA( R_A__const , 4 6 )DATA
EXTRA_REACTIONS
A r _ e x c _ 4 s + E l e c t r o n −> A r _ p i o n _ 0 + 2∗ E l e c t r o n
===
DATA( customEI , ( 1 . 3 7 e −13) ∗( Te_eV ∗ ∗ 0 . 5 )∗exp ( − 4 . 1 1 / Te_eV ) ,
→ u n i t s =m^ 3 / s , t h r e s h = 4 . 1 1 [ eV ] ) DATA
A r _ p i o n _ 0 + 2∗ E l e c t r o n −> Ar_g + E l e c t r o n
===
DATA( customEI , ( 8 . 7 5 e −39) ∗( Te_eV ∗∗( −4.5) ) , u n i t s =m^ 6 / s ,
→ t h r e s h =1.5∗ Te_eV )DATA
A r 2 _ p i o n _ 0 + E l e c t r o n −> Ar_g + A r _ e x c _ 4 s
===
DATA( customEI ,
→ ( 1 . 0 4 e −12) ∗ ( ( Te_K / 3 0 0 . 0 ) ∗∗( −0.67) ) ∗(1− exp ( − 4 1 8 . 0 / Tg_K ) ) / ( 1 − 0 . 3 1 ∗ exp ( − 4 1 8 . 0 / Tg_K ) ) , u n i t s =m^ 3 / s ,
→ t h r e s h =1.5∗ Te_eV )DATA
#...
A r _ e x c _ 4 s + Ar_g −> 2∗Ar_g
===
DATA( custom , ( 3 . 0 e −21) , u n i t s =m^ 3 / s )DATA
# Analytic cross sections
A r _ p i o n _ 1 + E l e c t r o n −> A r _ p i o n _ 2 + 2∗ E l e c t r o n
===
DATA( c u s t o m E I c r o s s a n a ,
→ t e s t v a l ∗ c r _ s c ∗ ( 1 . 0 / ( ( ep−t h r e s h ) +emax ) ) ∗ ( ( ( ep−t h r e s h ) / ( ( ep−t h r e s h ) +emax ) ) ∗ ∗ 1 . 1 2 7 ) , c r _ s c = 2 . 4 e −20[m^ 2 ] ,
→ emax = 9 0 . 0 [ eV ] , u n i t s =m^ 2 , t h r e s h = 5 . 5 [ eV ] , b o u n d s = [ 5 . 6 ; 1 0 0 . 0 ] )DATA
A r _ p i o n _ 0 + E l e c t r o n −> A r _ p i o n _ 1 + 2∗ E l e c t r o n
===
DATA( c u s t o m E I c r o s s a n a ,
→ c r _ s c ∗( Te_eV / ( ( ep−t h r e s h ) +emax ) ) ∗ ( ( ( ep−t h r e s h ) / ( ( ep−t h r e s h ) +emax ) ) ∗ ∗ 1 . 1 2 7 ) , c r _ s c = 2 . 8 3 5 8 e −20[m^ 2 ] ,
→ emax = 7 8 . 9 5 4 [ eV ] , u n i t s =m^ 2 , t h r e s h = 4 . 0 [ eV ] , b o u n d s = [ 5 . 6 ; 1 0 0 . 0 ] )DATA
# Passed numeric c r o s s s e c t i o n
Ar_exc_3 + E l e c t r o n <−> Ar_exc_13 + E l e c t r o n
===
DATA( c u s t o m E I c r o s s , d a t a f r o m p l o t s . t x t , u n i t s =m^ 3 / s , ID = 0 )DATA
Ar_exc_3 + E l e c t r o n <−> Ar_exc_11 + E l e c t r o n
===
DATA( c u s t o m E I c r o s s , d a t a f r o m p l o t s . t x t , u n i t s =m^ 3 / s , ID = 1 )DATA
#...
Ar_exc_1 + E l e c t r o n <−> Ar_exc_7 + E l e c t r o n
===
DATA( c u s t o m E I c r o s s , d a t a f r o m p l o t s . t x t , u n i t s =m^ 3 / s , ID = 7 )DATA
Ar_exc_1 + E l e c t r o n <−> Ar_exc_6 + E l e c t r o n
===
DATA( c u s t o m E I c r o s s , d a t a f r o m p l o t s . t x t , u n i t s =m^ 3 / s , ID = 8 )DATA

Figure A.2 Example reaction input file for an argon plasma.

A.3

Common-Sub-Expression Evaluation and Benefits

As can be seen above, the set of species continuity and conservation equations is extremely coupled. However, not only are they coupled, but there are repeated expressions when considering all
equations together. Representing large repetitions, species continuity equations, eq. 2.2, appear in
the energy equations in full and the electron energy equation, eq. 2.4, can appear in the gas temperature equation, eq. 2.3. A common smaller repetition is the product of reaction rate and reactant
species densities, which appear in the continuity equation of each involved reactant and product
species. While the repeated operations are barely registered in time for smaller systems, the number of species and considered reactions grows exponentially with the gas mixtures and physics
involved, table. A.1, which in turn burdens evaluation and compilation times of the system.

113

Model and Phase

ID

# Var. Spec.

# Rates

Energy Eqns.

Ar Plasma

0

17

61

Te

Ar-H2 -O2 Gas

1

10

64

Tg

Ar-CH4 -O2 Gas

2

31

348

Tg

Ar-H2 -O2 Mix

3

56

169

Tg & Te

Ar-CH4 -O2 Mix

4

142

552

Tg & Te

Ar-H2 -O2-N2 Mix

5

115

691

Tg & Te

Ar-CH4 -O2-N2 Mix

6

213

1305

Tg & Te

Table A.1 Species and reaction counts for a range of chemistry models demonstrating a huge
increase in complexity with included species.

Discussed previously, the algebra of the global model (and Jacobian) is separated from the
evaluation of the spline reaction rates and compiled. This is done to both speed up evaluation
time relative to a pure Python implementation and to support recovery of a previously run global
model instead of recomputing it. Reaction rates defined as analytic functions are included in these
expressions, assuming all used functions exist within the cmath library. This process is achieved
by taking advantage of the symbolic description of equations and transcribing them into C-code.
There is always motivation to have a faster evaluation time, but as the problems get larger, compiler
optimization begins to have trouble with the extremely large and coupled expressions. In particular, the recursive nature of gcc results in astronomical memory usage; it fails to compile the two
example methane cases due to a lack of RAM on a system with 64GB available. Interestingly, the
non-recursive clang compiler finishes compilation in the same cases, however, the objects created
segmentation fault during linking with the python instance.

114

ID

CSE

Setup (s)

F Op. Cnt.

F Call (s)

J Op. Cnt.

J Call (s)

0

False

4.29

518

7.88e-6

1251

1.08e-5

0

True

4.92

354

7.76e-6

807

1.02e-5

1

False

23.86

6139

8.07e-5

81264

8.91e-4

1

True

10.40

725

1.35e-5

3731

2.05e-5

2

False

7059.95

77182

5.17e-4

4704460

2.87e-2

2

True

156.94

4708

6.50e-5s

32295

1.08e-4

3

False

308.07

17630

8.37e-5

776015

3.07e-3

3

True

123.48

4771

2.60e-5

304263

3.01e-4

4

False

X

165612

X

23116971

X

4

True

1554.13

19434

8.96e-5

2596179

3.74e-3

5

False

12238.46

117938

3.47e-4

8328475

2.50e-2

5

True

1209.77

19020

7.37e-5

2132907

3.46e-3

6

False

X

335801

X

45667621

X

6

True

6599.70

42369

1.52e-4

8112852

1.18e-2

Table A.2 GCC-5.4 Compile and evaluation timing. Runs marked with X failed to compile due to
a lack of memory on a machine with 64GB of RAM available.

Instead of simply throwing the system of ODEs at a compiler, it is possible to take advantage
of the known structure of ODE expressions to reduce both the total operation count (linked to
evaluation time) and the necessary compilation time. Inspired by SymPy’s common subexpression
(CSE) detection and collection algorithm, the C-code is written with local variables pre-computing
expressions which are repeated throughout the ODE system. First, any species continuity or conservation equations which exists in another ODE equation are given a symbolic placeholder and
evaluated before the ODE system. Secondly, rate terms and similarly created shared expressions
115

are given placeholders and moved earlier in the evaluation time line. Finally, SymPy’s CSE algorithm is employed to recursively scan over local variable expressions and ODE equations (with
placeholders) such that no operation chain is repeated in a given evaluation of the system of ODEs
or Jacobian. This optimization not only reduces evaluation calls by up to an order of magnitude,
but also reduces compilation time and memory consumption significantly, as shown in tables A.2
and A.3.

ID

CSE

Setup (s)

F Op. Cnt.

F Call (s)

J Op. Cnt.

J Call (s)

0

False

4.64

518

9.35e-6

1251

1.21e-5

0

True

6.07

354

9.32e-6

807

1.20e-5

1

False

20.37

6139

6.65e-5

81264

8.32e-4

1

True

10.71

725

1.40e-5

3731

2.07e-5

2

False

1100.36

77182

5.18e-4

4704460

2.89e-2

2

True

154.65

4708

5.93e-5

32295

1.13e-4

3

False

169.155

17630

9.04e-5

776015

3.21e-3

3

True

113.36

4771

3.19e-5

304263

3.54e-4

4

False

X

165612

X

23116971

X

4

True

1133.84

19434

9.83e-5

2596179

4.04e-3

5

False

2502.094

117938

3.49e-4

8328475

3.00e-2

5

True

819.84

19020

7.19e-5

2132907

3.54e-3

6

False

X

335801

X

45667621

X

6

True

3150.56

42369

1.64e-4

8112852

1.26e-2

Table A.3 Clang-3.8 Compile and evaluation timing. Runs marked with X failed due to a segmentation fault during function loading.

116

APPENDIX B
APPENDIX TO PAC

B.1

Gas Phase Hydrogen-Oxygen Mechanism

The complete set of 16 reactions required to adequately describe H2 -O2 combustion in a classical
context. Double headed arrows imply the reaction is reversible

k

H2 −−0→ 2 H
k

O2 −−0→ 2 O
k

O2 + H ←−1→ OH + O
k

H2 + O ←−2→ OH + H
k

H2 + OH ←−3→ H2 O + H
k

H + wall −−4→ 0.5 H2
k

O2 + H + M ←−5→ HO2 + M
k

2 HO2 ←−6→ H2 O2 + O2
k

HO2 + H2 ←−7→ H2 O2 + H
k

HO2 + H ←−8→ 2 OH
k

HO2 + H ←−9→ H2 + O2
k

10
HO2 + H ←−
→ H2 O + O
k

11
H2 O2 + M ←−
→ 2 OH + M
k

12
H2 O2 + H ←−
→ H2 O + OH
k

13
H2 O2 + H ←−
→ H2 + O2
k

14
H2 O2 + OH ←−
→ H2 O + HO2
k

15
OH + H + M ←−
→ H2 O + M
k

16
H + H + M ←−
→ H2 + M

117

BIBLIOGRAPHY

118

BIBLIOGRAPHY
[1]

I V Adamovich, I Choi, N Jiang, J-H Kim, S Keshav, W R Lempert, E Mintusov, M Nishihara, M Samimy, and M Uddi. Plasma assisted ignition and high-speed flow control: nonthermal and thermal effects. Plasma Sources Science and Technology, 18(3):034018, 2009.

[2]

NL Aleksandrov, SV Kindysheva, MM Nudnova, and A Yu Starikovskiy. Mechanism of
ultra-fast heating in a non-equilibrium weakly ionized air discharge plasma in high electric
fields. Journal of Physics D: Applied Physics, 43(25):255201, 2010.

[3]

NB Anikin, SM Starikovskaia, and A Yu Starikovskii. Oxidation of saturated hydrocarbons
under the effect of nanosecond pulsed space discharge. Journal of Physics D: Applied
Physics, 39(15):3244, 2006.

[4]

Sumio Ashida, C. Lee, and M. A. Lieberman. Spatially averaged (global) model of time
modulated high density argon plasmas. Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films, 13(5):2498–2507, 1995.

[5]

Orbital ATK. Magic tool suite. https://www.orbitalatk.com/magic.

[6]

Peter Banks. Collision frequencies and energy transfer ions. Planetary and Space Science,
14(11):1105–1122, 1966.

[7]

NG Basov and AM Prokhorov. Applications of molecular beams to radio spectroscopic studies of rotation spectra of molecules. Journal of Experimental and Theoretical
Physics,(USSR), 27:431, 1954.

[8]

Raymond J Beach. Cw theory of quasi-three level end-pumped laser oscillators. Optics
Communications, 123(1-3):385–393, 1996.

[9]

Raymond J Beach. Cw theory of quasi-three level end-pumped laser oscillators. Optics
Communications, 123(1):385 – 393, 1996.

[10] Raymond J Beach, William F Krupke, V Keith Kanz, Stephen A Payne, Mark A Dubinskii,
and Larry D Merkle. End-pumped continuous-wave alkali vapor lasers: experiment, model,
and power scaling. JOSA B, 21(12):2151–2163, 2004.
[11] Houshang F Behbahani, ARTHUR FONTIJN, KLAUS MULLER-DETHLEFS, and FELIX J WEINBERG. The destruction of nitric oxide by nitrogen atoms from plasma jets.
Combustion Science and technology, 27(3-4):123–132, 1982.
[12] Stefan Behnel, Robert Bradshaw, Craig Citro, Lisandro Dalcin, Dag Sverre Seljebotn, and
Kurt Smith. Cython: The best of both worlds. Computing in Science & Engineering,
13(2):31–39, 2011.
[13] Adela Ben-Yakar and Ronald K Hanson. Cavity flame-holders for ignition and flame stabilization in scramjets: an overview. Journal of Propulsion and Power, 17(4):869–877,
2001.
119

[14] Carlo Benedetti, Andrea Sgattoni, Giorgio Turchetti, and Pasquale Londrillo. Aladyn: A
high-accuracy pic code for the maxwell–vlasov equations. IEEE Transactions on plasma
science, 36(4):1790–1798, 2008.
[15] Mario Bertolotti. The history of the laser. CRC press, 2004.
[16] Biagi v8.97 database. retrieved on July 4, 2012.
[17] Charles K Birdsall and William B Bridges. Space-charge instabilities in electron diodes and
plasma converters. Journal of Applied Physics, 32(12):2611–2618, 1961.
[18] Charles K Birdsall and A Bruce Langdon. Plasma physics via computer simulation. CRC
Press, 2004.
[19] José A Bittencourt. Fundamentals of plasma physics. Springer Science & Business Media,
2013.
[20] C de Boor. Splines as linear combinations of b–splines, a survey. Approximation Theory, II,
pages 1–47, 1976.
[21] Derek Bradley and Said MA Ibrahim. Electron temperatures in flame gases: experiment and
theory. Combustion and Flame, 24:169–171, 1975.
[22] Jean Brossel, Alfred Kastler, and Jacques Winter. Gréation optique d’une inégalité de population entre les sous-niveaux zeeman de l’état fondamental des atomes. Journal de Physique
Radium, 13:668, 1952.
[23] Oscar Buneman. Dissipation of currents in ionized media. Physical Review, 115(3):503,
1959.
[24] Alexander Burcat. Thermochemical Data for Combustion Calculations, pages 455–473.
Springer US, New York, NY, 1984.
[25] Rodney L Burton and PJ Turchi. Pulsed plasma thruster. Journal of Propulsion and Power,
14(5):716–735, 1998.
[26] Mario Capitelli, Carlos M Ferreira, Boris F Gordiets, and Alexey I Osipov. Rate coefficients
of chemical reactions. In Plasma Kinetics in Atmospheric Gases, pages 167–191. Springer,
2000.
[27] P Chabert, J Arancibia Monreal, J Bredin, L Popelier, and A Aanesland. Global model
of a gridded-ion thruster powered by a radiofrequency inductive coil. Physics of Plasmas,
19(7):073512, 2012.
[28] Liu Chen, A Bruce Langdon, and CK Birdsall. Reduction of the grid effects in simulation
plasmas. Journal of Computational Physics, 14(2):200–222, 1974.
[29] R K Cheng and A K Oppenheim. Autoignition in methane-hydrogen mixtures. Combustion
and Flame, 58(2):125 – 139, 1984.

120

[30] Norman Cohen. Flammability and explosion limits of H2 and H2/CO: a literature review.
Technical report, DTIC Document, 1992.
[31] Andrew Collette. Hdf5 for python, 2008–.
[32] Tech-X Corporation. Usim. https://www.txcorp.com/usim.
[33] Tech-X Corporation. Vsim. https://www.txcorp.com/vsim.
[34] Mark Steven Csele. Laser Modeling: A Numerical Approach with Algebra and Calculus.
CRC Press, 2014.
[35] Edward T Curran. Scramjet engines: the first forty years. Journal of Propulsion and Power,
17(6):1138–1148, 2001.
[36] UC Davis, US Santa Barbara, and UC San Diego. Kepler. https://kepler-project.
org/.
[37] John Dawson. One-dimensional plasma model. The Physics of Fluids, 5(4):445–459, 1962.
[38] Pierre Degond, Lorenzo Pareschi, and Giovanni Russo. Modeling and computational methods for kinetic equations. Springer Science & Business Media, 2004.
[39] AV Demyanov, IV Kochetov, and PA Mikheyev. Kinetic study of a cw optically pumped
laser with metastable rare gas atoms produced in an electric discharge. Journal of Physics
D: Applied Physics, 46(37):375202, 2013.
[40] J Dietz and The Iter Joint Central Team. The iter fusion experiment. Vacuum, 47(6-8):911–
918, 1996.
[41] Rajesh Dorai, Khaled Hassouni, and Mark J. Kushner. Interaction between soot particles
and nox during dielectric barrier discharge plasma remediation of simulated diesel exhaust.
Journal of Applied Physics, 88(10):6060–6071, 2000.
[42] BD Dudson, MV Umansky, XQ Xu, PB Snyder, and HR Wilson. Bout++: A framework for
parallel plasma fluid simulations. Computer Physics Communications, 180(9):1467–1480,
2009.
[43] JW Eastwood and RW Hockney. Computer simulation using particles. New York: Mc
GrawHill, 1981.
[44] Albert Einstein. Zur quantentheorie der strahlung. Physikalische Gesellschaft Zürich, 18,
1916.
[45] Tso Fan and ROBERTL Byer. Modeling and cw operation of a quasi-three-level 946 nm nd:
Yag laser. IEEE Journal of Quantum Electronics, 23(5):605–612, 1987.
[46] Jonathan Feinberg and Hans Petter Langtangen. Chaospy: An open source tool for designing
methods of uncertainty quantification. Journal of Computational Science, 11:46 – 57, 2015.

121

[47] A Flitti and S Pancheshnyi. Gas heating in fast pulsed discharges in n2–o2 mixtures. The
European Physical Journal Applied Physics, 45(2):21001, 2009.
[48] Jeffrey P Freidberg. Ideal magnetohydrodynamic theory of magnetic fusion systems. Reviews of Modern Physics, 54(3):801, 1982.
[49] Alexander Fridman, Alexander Gutsol, Shailesh Gangoli, Yiguang Ju, and Timothy Ombrello. Characteristics of gliding arc and its application in combustion enhancement. Journal of Propulsion and Power, 24(6):1216–1228, 2008.
[50] Alexander Fridman and Lawrence A Kennedy. Plasma physics and engineering. CRC press,
2004.
[51] A Friedman, JJ Barnard, RJ Briggs, RC Davidson, M Dorf, DP Grote, E Henestroza, EP Lee,
MA Leitner, BG Logan, et al. Toward a physics design for ndcx-ii, an ion accelerator for
warm dense matter and hif target physics studies. Nuclear Instruments and Methods in
Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 606(1):6–10, 2009.
[52] Josiah Willard Gibbs. A method of geometrical representation of the thermodynamic properties of substances by means of surfaces. Connecticut Academy, 1873.
[53] V A Godyak. Soviet Radio Frequency Discharge, chapter 5. Delphic Associates, Inc., Falls
Church, VA, 1986.
[54] Valery A Godyak. Nonequilibrium eedf in gas discharge plasmas. IEEE transactions on
plasma science, 34(3):755–766, 2006.
[55] Viktor Evgenevich Golant, Aleksei Petrovich Zhilinskii, and Igor Evgenevich Sakharov.
Fundamentals of plasma physics. John Wiley & Sons, 1980.
[56] G Gorchakov and F Lavrov. Influence of electric discharge on the region of spontaneous
ignition in the mixture 2H2–O2. Acta Physicochim. URSS, 1:139–44, 1934.
[57] Edward G Groff and Mark K Krage. Microwave effects on premixed flames. Combustion
and flame, 56(3):293–306, 1984.
[58] David P Grote, Alex Friedman, Irving Haber, and Simon Yu. Three-dimensional simulations of high current beams in induction accelerators with warp3d. Fusion engineering and
design, 32:193–200, 1996.
[59] ESI Group.
Cfd-ace+.
https://www.esi-group.com/software-solutions/
virtual-environment/cfd-multiphysics/ace-suite/cfd-ace.
[60] G J M Hagelaar and L C Pitchford. Solving the boltzmann equation to obtain electron
transport coefficients and rate coefficients for fluid models. Plasma Sources Science and
Technology, 14(4):722, 2005.

122

[61] Gordon D Hager and Glen P Perram. Extended saturation analysis and analytical model of
diode-pumped alkali lasers. In SPIE LASE, pages 75810J–75810J. International Society for
Optics and Photonics, 2010.
[62] Gordon D Hager and GP Perram. A three-level analytic model for alkali metal vapor lasers:
part i. narrowband optical pumping. Applied Physics B: Lasers and Optics, 101(1):45–56,
2010.
[63] Jiande Han, Leonid Glebov, George Venus, and Michael C. Heaven. Demonstration of a
diode-pumped metastable ar laser. Opt. Lett., 38(24):5458–5461, Dec 2013.
[64] Jiande Han and Michael C Heaven. Gain and lasing of optically pumped metastable rare
gas atoms. Optics letters, 37(11):2157–2159, 2012.
[65] William Happer. Optical pumping. Reviews of Modern Physics, 44(2):169, 1972.
[66] CE Haselfoot and PJ Kirkby. Xlv. the electrical effects produced by the explosion of hydrogen and oxygen. The London, Edinburgh, and Dublin Philosophical Magazine and Journal
of Science, 8(46):471–481, 1904.
[67] E Havlıcková. Fluid model of plasma and computational methods for solution, 2006.
[68] Ady Hershcovitch. High-pressure arcs as vacuum-atmosphere interface and plasma lens for
nonvacuum electron beam welding machines, electron beam melting, and nonvacuum ion
material modification. Journal of Applied Physics, 78(9):5283–5288, 1995.
[69] David A Hostutler, Gordon D Hager, and Michael C Heaven. Collisional quenching and radiation trapping kinetics for rb (5p) in the presence of ethane. In High-Power Laser Ablation
2008, pages 700522–700522. International Society for Optics and Photonics, 2008.
[70] Comsol Inc. Comsol multiphysics modeling software. https://www.comsol.com.
[71] Ist-lisbon database. retrieved on October 15, 2013.
[72] Y. Itikawa, A. Ichimura, K. Onda, K. Sakimoto, K. Takayanagi, Y. Hatano, M. Hayashi,
H. Nishimura, and S. Tsurubuchi. Cross sections for collisions of electrons and photons
with oxygen molecules. Journal of Physical and Chemical Reference Data, 18(1):23–42,
1989.
[73] S Jacobs, G Gould, and P Rabinowitz. Coherent light amplification in optically pumped cs
vapor. Physical Review Letters, 7(11):415, 1961.
[74] Lance Jacobsen, Campbell Carter, Robert Baurle, and Thomas Jackson. Plasma-assisted
ignition in scramjets. In 41st Aerospace Sciences Meeting and Exhibit, page 871, 2003.
[75] HC Jaggers and A Von Engel. The effect of electric fields on the burning velocity of various
flames. Combustion and Flame, 16(3):275–285, 1971.

123

[76] G M Janssen, J van Dijk, D A Benoy, M A Tas, K T A L Burm, W J Goedheer, J A M van der
Mullen, and D C Schram. PLASIMO, a general model: I. applied to an argon cascaded arc
plasma. Plasma Sources Science and Technology, 8(1):1, 1999.
[77] Ali Javan, William R Bennett Jr, and Donald R Herriott. Population inversion and continuous optical maser oscillation in a gas discharge containing a he-ne mixture. Physical Review
Letters, 6(3):106, 1961.
[78] T. Johnson, L. Palumbo, and A. Hunter. Kinetics simulation of high-power gas lasers. IEEE
Journal of Quantum Electronics, 15(5):289–301, May 1979.
[79] Eric Jones, Travis Oliphant, Pearu Peterson, et al. SciPy: Open source scientific tools for
Python, 2001–. [Online; accessed <today>].
[80] Yiguang Ju and Wenting Sun. Plasma assisted combustion: Dynamics and chemistry.
Progress in Energy and Combustion Science, 48:21 – 83, 2015.
[81] Yiguang Ju and Wenting Sun. Plasma assisted combustion: Dynamics and chemistry.
Progress in Energy and Combustion Science, 48:21 – 83, 2015.
[82] Md Humayun Kabir and Michael C Heaven. Collisional relaxation of the kr (4p 5 5p) states
in he, ne, and kr. In Proc. SPIE, volume 8238, page 823807, 2012.
[83] Fumihiko Kannari, Minoru Obara, and Tomoo Fujioka. An advanced kinetic model of
electron-beam-excited krf lasers including the vibrational relaxation in krf*(b) and collisional mixing of krf*(b,c). Journal of Applied Physics, 57(9):4309–4322, 1985.
[84] Alfred Kastler. Quelques suggestions concernant la production optique et la détection optique d’une inégalité de population des niveaux de quantifigation spatiale des atomes. application à l’expérience de stern et gerlach et à la résonance magnétique. Journal de Physique
Radium, 11:255–265, 1950.
[85] Robert J Kee, Fran M Rupley, Ellen Meeks, and James A Miller. Chemkin-iii: A fortran
chemical kinetics package for the analysis of gas-phase chemical and plasma kinetics. Sandia national laboratories report SAND96-8216, 1996.
[86] HC Kim, Felipe Iza, SS Yang, M Radmilovi´c-Radjenovi´c, and JK Lee. Particle and fluid
simulations of low-temperature plasma discharges: benchmarks and kinetic effects. Journal
of Physics D: Applied Physics, 38(19):R283, 2005.
[87] I. Kimura, H. Aoki, and M. Kato. Use of a plasma jet for flame stabilization and promotion
of combustion in supersonic air flows. Combust. Flame; (United States), 42:3, Sep 1981.
[88] Juergen F Kolb, A-A H Mohamed, Robert O Price, R James Swanson, Angela Bowman,
RL Chiavarini, M Stacey, and KH Schoenbach. Cold atmospheric pressure air plasma jet
for medical applications. Applied Physics Letters, 92(24):241501, 2008.

124

[89] C Konz, W Zwingmann, F Osmanlic, B Guillerminet, F Imbeaux, P Huynh, M Plociennik,
M Owsiak, T Zok, and M Dunne. First physics applications of the integrated tokamak
modelling (itm-tf) tools to the mhd stability analysis of experimental data and iter scenarios.
EPS, page O2, 2011.
[90] H Kopfermann and R Ladenburg. Experimental proof ofnegative dispersion.’. Nature,
122:438–439, 1928.
[91] IA Kossyi, A Yu Kostinsky, AA Matveyev, and VP Silakov. Kinetic scheme of the nonequilibrium discharge in nitrogen-oxygen mixtures. Plasma Sources Science and Technology, 1(3):207, 1992.
[92] A. Kramida, Yu. Ralchenko, J. Reader, and and NIST ASD Team. NIST Atomic Spectra
Database (ver. 5.3), [Online]. Available: http://physics.nist.gov/asd [2015, April 16]. National Institute of Standards and Technology, Gaithersburg, MD., 2015.
[93] W. Krupke. Diode pumped alkali-molecular lasers and amplifiers, November 18 2004. US
Patent App. 10/832,490.
[94] William F. Krupke. Diode pumped alkali lasers (dpals)—a review (rev1). Progress in
Quantum Electronics, 36(1):4 – 28, 2012. Special issue in honor of Professor J. Gary Eden
on the occasion of his 60th birthday.
[95] Matthew W Kunz, James M Stone, and Xue-Ning Bai. Pegasus: A new hybrid-kinetic
particle-in-cell code for astrophysical plasma dynamics. Journal of Computational Physics,
259:154–174, 2014.
[96] Mark J Kushner. A model for the discharge kinetics and plasma chemistry during plasma
enhanced chemical vapor deposition of amorphous silicon. Journal of Applied Physics,
63(8):2532–2551, 1988.
[97] Mark J Kushner. Hybrid modelling of low temperature plasmas for fundamental investigations and equipment design. Journal of Physics D: Applied Physics, 42(19):194013, 2009.
[98] Irving Langmuir. Oscillations in ionized gases. Proceedings of the National Academy of
Sciences, 14(8):627–637, 1928.
[99] Antoine Laurent Lavoisier. Traité élémentaire de chimie: présenté dans un ordre nouveau
et d’après les découvertes modernes; avec figures, volume 1. Chez Cuchet, Libraire, 1793.
[100] C. Lee and M. A. Lieberman. Global model of ar, o2, cl2, and ar/o2 high-density plasma
discharges. Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films,
13(2):368–380, 1995.
[101] Sergey B Leonov, Igor V Kochetov, Anatoly P Napartovich, Vladimir A Sabel’Nikov, and
Dmitry A Yarantsev. Plasma-induced ethylene ignition and flameholding in confined supersonic air flow at low temperatures. IEEE Transactions on Plasma Science, 39(2):781–787,
2011.

125

[102] Bernard Lewis. The effect of an electric field on flames and their propagation. Journal of
the American Chemical Society, 53(4):1304–1313, 1931.
[103] Bernard Lewis and Guenther Von Elbe. Combustion, flames and explosions of gases 3rd .
Academic Press, 1987.
[104] G.C. Light and J.H. Matsumoto. The effect of vibrational excitation in the reactions of oh
with h2. Chemical Physics Letters, 58(4):578 – 581, 1978.
[105] F. A. Lindemann, Svante Arrhenius, Irving Langmuir, N. R. Dhar, J. Perrin, and W. C.
McC. Lewis. Discussion on "the radiation theory of chemical action". Trans. Faraday Soc.,
17:598–606, 1922.
[106] David Lindley. Focus: Invention of the maser and laser. Physics, 15:4, 2005.
[107] S Longo and M Capitelli. A simple approach to treat anisotropic elastic collisions in monte
carlo calculations of the electron energy distribution function in cold plasmas. Plasma
Chemistry and Plasma Processing, 14(1):1–13, 1994.
[108] Svetlana G Lukishova. Valentin a. fabrikant: negative absorption, his 1951 patent application for amplification of electromagnetic radiation (ultraviolet, visible, infrared and radio
spectral regions) and his experiments. Journal of the European Optical Society-Rapid publications, 5, 2010.
[109] John L Lumley. Computational modeling of turbulent flows. Advances in applied mechanics, 18:123–176, 1979.
[110] Alejandro Luque. Bolos: An open source solver for the boltzmann equation.
[111] Heather L MacLean and Lester B Lave. Evaluating automobile fuel/propulsion system
technologies. Progress in energy and combustion science, 29(1):1–69, 2003.
[112] Bruce H Mahan. I microscopic reversibility. J. Chem. Educ., 52:299–303, 1975.
[113] T. H. MAIMAN. Stimulated optical radiation in ruby. Nature, 187(4736):493–494, Aug
1960.
[114] Maiman builds first working laser. APS News, 10(5), 2010.
[115] Rudolf Maly. Spark Ignition: Its Physics and Effect on the Internal Combustion Engine.
Springer US, Boston, MA, 1984.
[116] Aram H Markosyan, Alejandro Luque, Francisco J Gordillo-Vázquez, and Ute Ebert. Pumpkin: A tool to find principal pathways in plasma chemical models. Computer Physics Communications, 185(10):2697–2702, 2014.
ˇ
[117] Aaron Meurer, Christopher P Smith, Mateusz Paprocki, Ondˇrej Certík,
Sergey B Kirpichev,
Matthew Rocklin, AMiT Kumar, Sergiu Ivanov, Jason K Moore, Sartaj Singh, et al. Sympy:
symbolic computing in python. PeerJ Computer Science, 3:e103, 2017.

126

[118] W. J. Miloch. Plasma physics and numerical simulations.
studier/emner/matnat/math/MEK4470/h14/plasma.pdf.

https://www.uio.no/

[119] Morgan database. retrieved on October 15, 2013.
[120] S K Nam and J P Verboncoeur. Effect of electron energy distribution function on the global
model for high power microwave breakdown at high pressures. Applied Physics Letters,
92(23):231502, 2008.
[121] Richard G Newell and Stuart Iler. The global energy outlook. Technical report, National
Bureau of Economic Research, 2013.
[122] Timothy Ombrello, Yiguang Ju, and Alexander Fridman. Kinetic ignition enhancement of diffusion flames by nonequilibrium magnetic gliding arc plasma. AIAA journal,
46(10):2424–2433, 2008.
[123] S Pancheshnyi, B Eismann, G J M Hagelaar, and L C Pitchford. Computer code ZDPlasKin.
University of Toulouse, LAPLACE, CNRS-UPS-INP, Toulouse, France, 2008.
[124] Wonchull Park, EV Belova, GY Fu, XZ Tang, HR Strauss, and LE Sugiyama. Plasma
simulation studies using multilevel physics models. Physics of Plasmas, 6(5):1796–1803,
1999.
[125] FI Parra, E Ahedo, JM Fife, and M Martinez-Sanchez. A two-dimensional hybrid model of
the hall thruster discharge. Journal of Applied Physics, 100(2):023304, 2006.
[126] E.L. Petersen, D.M. Kalitan, S. Simmons, G. Bourque, H.J. Curran, and J.M Simmie. Nui
galway. Proc. Combust. Inst., 31:447–454, 2007.
[127] L. R. Peterson and J. E. Allen Jr. Electron impact cross sections for argon. The Journal of
Chemical Physics, 56(12):6068–6076, 1972.
[128] Nuno Rombert Pinhão. PLASMAKIN: A library to evaluate the chemical kinetics in gases,
2001–.
[129] N A Popov. Kinetics of plasma-assisted combustion: effect of non-equilibrium excitation
on the ignition and oxidation of combustible mixtures. Plasma Sources Science and Technology, 25(4):043002, 2016.
[130] NA Popov. Investigation of the mechanism for rapid heating of nitrogen and air in gas
discharges. Plasma physics reports, 27(10):886–896, 2001.
[131] Robert K Porteous and David B Graves. Modeling and simulation of magnetically confined
low-pressure plasmas in two dimensions. IEEE Transactions on Plasma Science, 19(2):204–
213, 1991.
[132] Quantemol-P: Plasma kinetics.

127

[133] M A Raadu, I Axnäs, J T Gudmundsson, C Huo, and N Brenning. An ionization region
model for high-power impulse magnetron sputtering discharges. Plasma Sources Science
and Technology, 20(6):065007, 2011.
[134] Paul Rabinowitz, Steven Jacobs, and Gordon Gould. Continuous optically pumped cs laser.
Applied Optics, 1(4):513–516, 1962.
[135] WT Rawlins, KL Galbally-Kinney, SJ Davis, AR Hoskinson, JA Hopwood, and
MC Heaven. Optically pumped microplasma rare gas laser. Optics express, 23(4):4804–
4813, 2015.
[136] Arthur L Schawlow and Charles H Townes. Infrared and optical masers. Physical Review,
112(6):1940, 1958.
[137] Andrea Sgattoni, Luca Fedeli, Stefano Sinigardi, Alberto Marocchino, Andrea Macchi,
Volker Weinberg, and Anupam Karmakar. Optimising piccante-an open source particlein-cell code for advanced simulations on tier-0 systems. arXiv preprint arXiv:1503.02464,
2015.
[138] Kazuya Shimizu, Atsushi Hibi, Mitsuo Koshi, Youhi Morii, and Nobuyuki Tsuboi. Updated kinetic mechanism for high-pressure hydrogen combustion. Journal of Propulsion
and Power, 27(2):383–395, Mar 2011.
[139] Uri Shumlak. Algorithm development for the multi-fluid plasma model. Technical report,
DTIC Document, 2011.
[140] Uri Shumlak and J Loverich. Approximate riemann solver for the two-fluid plasma model.
Journal of Computational Physics, 187(2):620–638, 2003.
[141] Van Sijde, B Der, JJAM van der Mullen, and DC Schram. Collisional radiative models in
plasmas. Beiträge aus der Plasmaphysik, 24(5):447–473, 1984.
[142] VV Smirnov, OM Stelmakh, VI Fabelinsky, DN Kozlov, AM Starik, and NS Titova. On the
influence of electronically excited oxygen molecules on combustion of hydrogen–oxygen
mixture. Journal of Physics D: Applied Physics, 41(19):192001, 2008.
[143] G P Smith, D M Golden, M Frenklach, N W Moriarty, B Eiteneer, M Goldenberg, C T
Bowman, R K Hanson, S Song, W C Gardiner Jr., V V Lissianski, and Z Qin.
[144] Gregory P. Smith, David M. Golden, Michael Frenklach, Nigel W. Moriarty, Boris Eiteneer,
Mikhail Goldenberg, C. Thomas Bowman, Ronald K. Hanson, Soonho Song, Jr. William
C. Gardiner, Vitali V. Lissianski, and Zhiwei Qin. Gri-mech 3.0.
[145] Josip Z Soln. Finite element, finite difference, and finite volume methods: Examples and
their comparisons. Technical report, DTIC Document, 1996.
[146] Timothy J Sommerer and Mark J Kushner. Numerical investigation of the kinetics and
chemistry of rf glow discharge plasmas sustained in he, n2, o2, he/n2/o2, he/cf4/o2, and
sih4/nh3 using a monte carlo-fluid hybrid model. Journal of applied physics, 71(4):1654–
1673, 1992.
128

[147] AM Starik, BI Lukhovitski˘ı, VV Naumov, and NS Titova. On combustion enhancement mechanisms in the case of electrical-discharge-excited oxygen molecules. Technical
Physics, 52(10):1281–1290, 2007.
[148] AM Starik and NS Titova. Kinetics of detonation initiation in the supersonic flow of the
h2+ o2 (air) mixture in o2 molecule excitation by resonance laser radiation. Kinetics and
Catalysis, 44(1):28–39, 2003.
[149] S M Starikovskaia. Plasma assisted ignition and combustion. Journal of Physics D: Applied
Physics, 39(16):R265, 2006.
[150] SM Starikovskaia, A Yu Starikovskii, and DV Zatsepin. Hydrogen oxidation in a stoichiometric hydrogen-air mixture in the fast ionization wave. Combustion Theory and Modelling,
5(1):97–129, 2001.
[151] Andrey Starikovskiy. Plasma assisted combustion mechanism for small hydrocarbons. In
53rd AIAA Aerospace Sciences Meeting, page 0158, 2015.
[152] Andrey Starikovskiy and Nickolay Aleksandrov. Plasma-assisted ignition and combustion.
In Max Mulder, editor, Aeronautics and Astronautics, chapter 12. InTech, Rijeka, 2011.
[153] Andrey Starikovskiy and Nickolay Aleksandrov. Plasma-assisted ignition and combustion.
Progress in Energy and Combustion Science, 39(1):61 – 110, 2013.
[154] DH Stedman and DW Setser. Chemical applications of metastable rare gas atoms. In Progr.
Reaction Kinetics, volume 6, pages 193–238. Wiley, 1978.
[155] K. W. H. Stevens. Masers. Contemporary Physics, 2(1):1–13, 1960.
[156] P.H. Stewart, C.W. Larson, and D.M. Golden. Pressure and temperature dependence of
reactions proceeding via a bound complex. 2. application to 2ch3 → c2h5 + h. Combustion
and Flame, 75(1):25 – 31, 1989.
[157] Emanuel S Stockman, Sohail H Zaidi, Richard B Miles, Campbell D Carter, and Michael D
Ryan. Measurements of combustion properties in a microwave enhanced flame. Combustion
and Flame, 156(7):1453–1461, 2009.
[158] Wenting Sun, Mruthunjaya Uddi, Sang Hee Won, Timothy Ombrello, Campbell Carter,
and Yiguang Ju. Kinetic effects of non-equilibrium plasma-assisted methane oxidation on
diffusion flame extinction limits. Combustion and Flame, 159(1):221–229, 2012.
[159] Wenting Sun, Sang Hee Won, Timothy Ombrello, Campbell Carter, and Yiguang Ju. Direct
ignition and s-curve transition by in situ nano-second pulsed discharge in methane/oxygen/helium counterflow flame. Proceedings of the Combustion Institute, 34(1):847–855,
2013.
[160] Kenichi Takita, Naoyuki Abe, Goro Masuya, and Yiguang Ju. Ignition enhancement by
addition of no and no 2 from a n 2/o 2 plasma torch in a supersonic flow. Proceedings of the
Combustion Institute, 31(2):2489–2496, 2007.
129

[161] K Tao, D Mao, and J Hopwood. Ionized physical vapor deposition of titanium nitride: A
global plasma model. Journal of applied physics, 91(7):4040–4048, 2002.
[162] Ping King Tien and J Moshman. Monte carlo calculation of noise near the potential minimum of a high-frequency diode. Journal of Applied Physics, 27(9):1067–1078, 1956.
[163] J. Troe. Fall-off curves of unimolecular reactions. Berichte der Bunsengesellschaft für
physikalische Chemie, 78(5):478–488, 1974.
[164] Tamás Turányi. Applications of sensitivity analysis to combustion chemistry. Reliability
Engineering & System Safety, 57(1):41–48, 1997.
[165] Martin A Uman. The peak temperature of lightning. Journal of Atmospheric and Terrestrial
Physics, 26(1):123–128, 1964.
[166] V Vahedi, CK Birdsall, MA Lieberman, G DiPeso, and TD Rognlien. Verification of frequency scaling laws for capacitive radio-frequency discharges using two-dimensional simulations. Physics of Fluids B: Plasma Physics, 5(7):2719–2729, 1993.
[167] Stéfan van der Walt, S. Chris Colbert, and Gaël Varoquaux. The numpy array: A structure
for efficient numerical computation. Computing in Science & Engineering, 13(2):22–30,
2011.
[168] Ayyaswamy Venkattraman and Abhishek Kumar Verma. plasmafoam: An openfoam framework for computational plasma physics and chemistry. In APS Meeting Abstracts, 2016.
[169] John P Verboncoeur. Particle simulation of plasmas: review and advances. Plasma Physics
and Controlled Fusion, 47(5A):A231, 2005.
[170] John P Verboncoeur, Maria Virgínia Alves, V Vahedi, and Ch K Birdsall. Simultaneous
potential and circuit solution for 1d bounded plasma particle simulation codes. Journal of
Computational Physics, 104(2):321–328, 1993.
[171] John P Verboncoeur, A Bruce Langdon, and NT Gladd. An object-oriented electromagnetic
pic code. Computer Physics Communications, 87(1-2):199–211, 1995.
[172] EL Vold, F Najmabadi, and RW Conn. Fluid model equations for the tokamak plasma edge.
Physics of Fluids B: Plasma Physics, 3(11):3132–3152, 1991.
[173] Quan-De Wang. An updated detailed reaction mechanism for syngas combustion. RSC Adv.,
4:4564–4585, 2014.
[174] K-D Weltmann, E Kindel, R Brandenburg, Ch Meyer, R Bussiahn, Ch Wilke, and
Th Von Woedtke. Atmospheric pressure plasma jet for medical therapy: plasma parameters and risk estimation. Contributions to plasma physics, 49(9):631–640, 2009.
[175] John Henry White. The history of the phlogiston theory. Ams Pr Inc, 1932.
[176] Chen Wu, Alan R Hoskinson, and Jeffrey Hopwood. Stable linear plasma arrays at atmospheric pressure. Plasma Sources Science and Technology, 20(4):045022, 2011.
130

[177] Liang Wu, Alexander A Fridman, and Andrey Yu Starikovskiy. Kinetics of plasma assisted
combustion at low reduced electric fields. In 48th AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition, pages 4–7, 2010.
[178] Zhiyao Yin, Keisuke Takashima, and Igor V Adamovich. Ignition time measurements in
repetitive nanosecond pulse hydrogen–air plasmas at elevated initial temperatures. IEEE
Transactions on Plasma Science, 39(12):3269–3282, 2011.
[179] Nathan D. Zameroski, Gordon D. Hager, Wolfgang Rudolph, and David A. Hostutler. Experimental and numerical modeling studies of a pulsed rubidium optically pumped alkali
metal vapor laser. J. Opt. Soc. Am. B, 28(5):1088–1099, May 2011.
[180] Nathan D Zameroski, Gordon D Hager, Wolfgang Rudolph, and David A Hostutler. Experimental and numerical modeling studies of a pulsed rubidium optically pumped alkali metal
vapor laser. JOSA B, 28(5):1088–1099, 2011.
[181] Zhi-Bo Zhang and Jeffrey Hopwood. Linear arrays of stable atmospheric pressure microplasmas. Applied Physics Letters, 95(16):161502, 2009.
[182] BV Zhdanov, MK Shaffer, and RJ Knize. Scaling of diode-pumped cs laser: transverse
pump, unstable cavity, mopa. In SPIE LASE, pages 75810F–75810F. International Society
for Optics and Photonics, 2010.

131