METHODS TO ENSURE SATISFACTION OF CHARGE CONSERVATION FOR FINITE

ELEMENT PARTICLE-IN-CELL

By

Scott Timothy O’Connor

A DISSERTATION

Michigan State University

in partial fulfillment of the requirements

Submitted to

for the degree of

Electrical Engineering – Doctor of Philosophy

Computational Mathematics, Science and Engineering – Dual Major

2021

ABSTRACT

METHODS TO ENSURE SATISFACTION OF CHARGE CONSERVATION FOR FINITE

ELEMENT PARTICLE-IN-CELL

By

Scott Timothy O’Connor

Simulation of space charge and plasma is an integral part of many scientific and
engineering processes, with specific applications of pulsed power, particle accelerators,
integrated chip manufacture, satellites, and medicine, to name a few examples. A number
of methods have been used to simulate the underlying physics. One such method, Particle-
In-Cell (PIC), does this by mapping all particle-field interactions through a mesh, and using
macro-particles as statistically significant markers on the distribution function phase space.
While PIC has been around since the 1960 much of the work has been focused on Finite
Difference Time Domain (FDTD) for solving electromagnetic fields. Recently there has
been a push for using Finite Element Methods due to ease at which it can model complex
geometry. The challenge with using FEM or any method with an unstructured grid is
developing a consistent methodology for evolving fields and particle populations such
that all conservation laws are held. Developing this scheme is the principal contribution of
this thesis.

Specifically, to set the stage for the development, we start with developing a set of
three dimensional validation tests for PIC methods. These tests have either an analytic or
quasi-analytic solutions. Next, we delve into structure preserving geometric algorithms
that inform how different quantities that appear both Maxwell’s equations and Newton’s
laws should be represented. Using these, we prescribe a consistent particle mapping
methodology for an implicit field solver. We prove and demonstrate that conservation laws
are satisfied. Finally, to extend the regime of applicability of this approach, we develop
quasi-Helmholtz projectors that ensures analytical satisfaction of the equation of continuity
and Gauss’ law for magnetic flux density, and discrete satisfaction of Gauss’ law for electric

flux density. The efficacy of these methods and their benefits are presented on benchmark
examples developed earlier.

Copyright by
SCOTT TIMOTHY O’CONNOR
2021

ACKNOWLEDGMENTS

The path of grad school is not taken alone and I would not be here without the support
of many individuals. I would like to thank my advisor Dr. Shanker Balasubramaniam for
his support, encouragement and guidance. He has gone above and beyond what advisors
need to do for their students. I would also like to thanks my committee, Dr. Chahal, Dr.
Verboncoeur, Dr. Christlieb and Dr. Luginsland, for their support over the years. Your
guidance and advice for research, life and a professional career has been valuable. I would
like to thank my fellow graduate students in the EMRG (ElectroMagnetic Research Group)
and PTSG (Plasma Theory and Simulation Group) for there support encouragement and
ideas. I would also like to thank the ECE and CMSE department for there support over the
years. There are countless individuals who have helped me navigate the graduate school
path. I would specifically like to thanks all of Dr. Shanker’s students, this has been a team
effort and would not be here without your support. I would also like to thank my parents
Nancy and Steve, and my brothers Peter and Andrew for their support.

v

TABLE OF CONTENTS

LIST OF TABLES .

. .

LIST OF FIGURES . .

.

.

.

.

.

.

LIST OF ALGORITHMS .

.

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii

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

ix

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii

CHAPTER 1

INTRODUCTION AND BACKGROUND . . . . . . . . . . . . . . .
1.1 Background and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.

1.1.1 Outline of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
8

CHAPTER 2 A SET OF BENCHMARK TESTS FOR VALIDATION OF 3D

.

.

.
.

.
.

.
.
. .

.
Introduction .

.
2.1 Abstract .
2.2
.
2.3 Problem Description .

2.3.1 The Discrete Domain . . . . . . . . . . . . . . . . . . . . . . . . . . .

PARTICLE IN CELL METHODS . . . . . . . . . . . . . . . . . . . . 10
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
. 14
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4.1 Basic Particle Motion . . . . . . . . . . . . . . . . . . . . . . . . . . .
. 16
2.4.2 Expanding Particle Beam . . . . . . . . . . . . . . . . . . . . . . . . . . 21
. 27
2.4.3 Adiabatic Expanding Plasmas . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

2.4 Benchmark Tests .

2.5 Summary . .

. .

.

.

.

.

.

.

.

.

CHAPTER 3 TIME INTEGRATOR AGNOSTIC CHARGE CONSERVING FI-

3.4 Modified TDFEM-PIC .

.

3.1 Abstract .
.
3.2
Introduction .
3.3 Preliminaries

.

.
.
.

.
.
.

.
.
.

.
.
. .
. .

NITE ELEMENT PIC . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3.1 Overview of Method . . . . . . . . . . . . . . . . . . . . . . . . . . .
. 34
3.3.2 Unconditionally Stable Time Marching . . . . . . . . . . . . . . . . . . 37
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
Integrator Agnostic Charge Conserving Scheme . . . . . . . . . . . . 39
3.4.1
3.4.2 Particle Pusher .
. 40
3.4.3 Particle Path and Current Mapping . . . . . . . . . . . . . . . . . . . . 41
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.5.1 Higher Order Particle Motion . . . . . . . . . . . . . . . . . . . . . . . 44
3.5.2 Expanding Particle Beam . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.5.3 Adiabatic Expanding Plasma . . . . . . . . . . . . . . . . . . . . . . . 50
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

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

. .

. .

. .

.

.

.

.

.

.

.

.

.

.

.

.

3.5 Results .

3.6 Summary

CHAPTER 4 QUASI-HELMHOLTZ DECOMPOSITION FOR FINITE ELE-

4.1 Abstract .

MENT PARTICLE-IN-CELL . . . . . . . . . . . . . . . . . . . . . . . 53
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

.

.

.

.

.

.

.

.

.

vi

.
.

.
.
.

.

.

.

.

. .
. .

4.2
Introduction .
4.3 Formulation .

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

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
.
4.3.1 Preliminaries
.
. 57
4.3.2 Helmholtz Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . 58
. 59
4.3.3 Discrete Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Quasi-Helmholtz Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.5 Decomposition .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.5.1 Projector Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
. 64
4.5.2
ns . . . . . . . . . . 65
. . . . . . . . . . .
. 66
4.6 Null spaces, gauging and other remarks . . . . . . . . . . . . . . . . . . . . . 67
4.7 Results .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
Field Solver Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.7.1.1
Eigenvalues . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.7.1.2 Neumann Boundary Condition . . . . . . . . . . . . . . . . . 71
4.7.2 Particle Beam .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.7.3 Adiabatic Expansion of Plasma . . . . . . . . . . . . . . . . . . . . . . 74
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

Solution for non solenoidal component En
Solution to the solenoidal components

System of Equations
4.5.2.1
4.5.2.2

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

4.8 Conclusion .

4.7.1

. .

.

.

.

.

.

.

.

.

.

.

.

.

.

CHAPTER 5 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

APPENDICES .
APPENDIX A

.

.

.

APPENDIX B

.

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

.
A NOVEL PORT/NETWORK PARAMETER EXTRACTION
TECHNIQUE FOR COUPLING CIRCUITS WITH FULL-
WAVE TIME DOMAIN INTEGRAL EQUATION SOLVERS . . 82
SIMPLE CHARGE CONSERVATION EXAMPLE . . . . . . . . 113

BIBLIOGRAPHY .

. .

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

vii

LIST OF TABLES

Table 2.1: Cyclotron Motion .

.

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

. 18

Table 2.2: Expanding Particle Beam Parameters

. . . . . . . . . . . . . . . . . . . . . 22

Table 2.3: Beam Profile

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

Table 2.4: Adiabatic Expanding Plasmas

. . . . . . . . . . . . . . . . . . . . . . . . . 28

Table 2.5: Adiabatic Expanding Plasma . . . . . . . . . . . . . . . . . . . . . . . . . . 29

Table 3.1: Cyclotron Motion .

.

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

. 45

Table 3.2: Expanding Particle Beam Parameters

. . . . . . . . . . . . . . . . . . . . . 47

Table 3.3: Adiabatic Expanding Plasmas

. . . . . . . . . . . . . . . . . . . . . . . . . 50

Table 4.1: Expanding Particle Beam Parameters

. . . . . . . . . . . . . . . . . . . . . 74

Table A.1: Simulation parameters used in the diode mixer example . . . . . . . . . . 101

Table A.2: Parameters used in diode attenuator example . . . . . . . . . . . . . . . . 105

Table A.3: Table of voltages across the 86 Ω output resistor in both TD and FD

codes for various values of the control resistance Rc . . . . . . . . . . . . . 106

Table A.4: Parameters used in Balanced Mixer example . . . . . . . . . . . . . . . . . 108

Table B.1: Values of charge conservation for first order time step.

. . . . . . . . . .

. 117

Table B.2: Current using 1 quadrature point and analytic integral for linear path.

. 117

Table B.3: Non-linear path analytic current, and for 1 and 2 point quadrature rule.

. 118

Table B.4: Current Values for particle crossing cells.

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

. 119

Table B.5: Single quadrature point integration current.

. . . . . . . . . . . . . . . . . 119

viii

LIST OF FIGURES

Figure 1.1: Phase diagram of matter. . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

Figure 1.2: Basic particle in cell (PIC) system. Other details such as collision,

emissions, absorption etc. are left out for simplicity.

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

2

7

Figure 2.1: Particle acceleration due to electric field. . . . . . . . . . . . . . . . . . . . 19

Figure 2.2: Cyclotron particle motion of an electron in a magnetic field.

. . . . . . . 20

Figure 2.3: Cyclotron motion error from just Boris push with several different
time step sizes, as well as influence from the mesh. The time step
∆t (cid:3) 46.1766ps or 11340 steps per cycle.

. . . . . . . . . . . . . . . . . . . 20

Figure 2.4: Error as time step size decreases for Boris push.

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

. 21

Figure 2.5: Drift tube for particle beam test case . . . . . . . . . . . . . . . . . . . . . 22

Figure 2.6: XOOPIC vs FEMPIC particle position at steady state. Approximate

analytic solution profile is shown by the dashed line.

. . . . . . . . . . . 24

Figure 2.7: Particle Beam Energy with a 2ns turn on time . . . . . . . . . . . . . . .

. 25

Figure 2.8: Noise measured from three beam test cases using FEMPIC method

with different random number generator seeds.

. . . . . . . . . . . . . . 25

Figure 2.9: Fourier transform of the electric field from FEMPIC method for three

beam tests that had different random number generator seeds. . . . . . . 26

Figure 2.10: Charge Conservation for particle beam test case measured by |(ρn+1 −

ρn)/∆t − [∇·]Jn+1/2|/|(ρn+1 − ρn)/∆t|/Npart

. . . . . . . . . . . . . . . . . 26

Figure 2.11: Distribution of electrons radially at two instants in time with analytic

solutions. .

. .

. .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

Figure 3.1: Particle path for a single particles with start and location and intersec-

tion point.

. .

. .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

Figure 3.2: Mean relative error in position for Adam-Bashforth Orders 1-5 com-

pared with the Boris push, shown in black.

. . . . . . . . . . . . . . . . . 45

ix

Figure 3.3: Mean relative error in velocity for Adam-Bashforth Orders 1-5 com-

pared with the Boris push, shown in black.

. . . . . . . . . . . . . . . . . 46

Figure 3.4: Expanding particle beam macro particles in the z vs r plan. Particle
locations from both mixed finite element methods and wave equation
versions are compared with XOOPIC beam profile.

. . . . . . . . . . . . 48

Figure 3.5: Electric field values are the radial component half way down the tube
16 mm from the center of the tube. Multiple simulation with different
time steps are performed.

. . . . . . . . . . . . . . . . . . . . . . . . . . . 48

Figure 3.6: Discrete Gauss’s Law error per particle for Newmark-Beta mixed
finite element (NM-MFEM), Newmark-Beta wave equation (WE),
backwards difference mixed finite elements (BD-MFEM) using the
charge conservation technique provided here.

. . . . . . . . . . . . . . . 49

Figure 3.7: Discrete Gauss’s law error per particle various time steps using the

mixed finite element methods using Newmark time stepping.

. . . . . . 49

Figure 3.8: Mixed Finite Element particle beam expansion . . . . . . . . . . . . . .

. 51

Figure 3.9: Wave Equation Adiabatic Expanding Plasma . . . . . . . . . . . . . . .

. 51

Figure 4.1: Tree (solid) and cotree (dashed) edges for a 2D geometry.

. . . . . . . . . 67

Figure 4.2: Tree (solid) and cotree (dashed) faces for a simplified 2D tree structure. . 67

Figure 4.3: Eigenvalues Newmark-beta time stepping of MFEM and solenoidal-

MFEM system for PEC box.

. . . . . . . . . . . . . . . . . . . . . . . . . . 70

Figure 4.4: Total electric field for plane wave through a box with Neumann

boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

Figure 4.5: Non-solenoidal electric field for plane wave through a box with Neu-

mann boundary conditions.

. . . . . . . . . . . . . . . . . . . . . . . . . . 72

Figure 4.6: Solenoidal electric field for plane wave through a box with Neumann

boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

Figure 4.7: Comparison of radial electric field ¯E between MFEM, qausi-MFEM,

wave equation, and quasi wave equation.

. . . . . . . . . . . . . . . . . . 75

Figure 4.8: Comparison of radial rotational electric field ¯Es between MFEM, qausi-

MFEM, wave equation, and quasi wave equation.

. . . . . . . . . . . . . 75

x

Figure 4.9: Comparison of radial irrotational electric field ¯Eir between MFEM,

qausi-MFEM, wave equation, and quasi wave equation. . . . . . . . . . . 76

Figure 4.10: Comparison of radial cotree magnetic field between cotree solve and

newmark solve.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

Figure 4.11: Comparison between regular mixed finite element, wave equation,
and quasi-Helmholtz wave equation error in the discrete continuity
equation. The growing in time DC null space in the wave equation is
remedied by the quasi-Helmholtz decomposition.

. . . . . . . . . . . . . 77

Figure 4.12: Quasi-Helmholtz decomposition adiabatic expansion using quasi-MFEM. 78

Figure A.1: (a) Perfect electric conducting surface with several circuits attached;

(b) An equivalent model–circuit system with a port admittance.

. . . . . 86

Figure A.2: EM-CKT port with current and voltage from circuit interacting with

current density and the electric fields.

. . . . . . . . . . . . . . . . . . . . 92

Figure A.3: (a) Dipole antenna geometry and coupling. (b) Chebyshev low-pass

filter schematic.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

Figure A.4: Comparison of port voltages obtained via the port extraction method
vs. solving either (a) the full coupled TDIE system directly (left) and
(b) a harmonic-balance solution using frequency-domain port extraction. 100

Figure A.5: Port voltages at the bowtie antenna input due to the diode attenuator
circuit obtained via the TD port extraction method. The solution is
compared with a Harmonic Balance (HB) solution using 45 harmonics. . 102

Figure A.6: Port extracted signal from bowtie using third order circuit time basis

function .

.

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

Figure A.7: (a) Diode mixer circuit schematic. (b) Layout of the bowtie antenna

and coupling with circuit.

. . . . . . . . . . . . . . . . . . . . . . . . . . . 103

Figure A.8: Comparison of the number of non-linear iterations and time for the
port extraction and full coupled method simulating a dipole with a
non-linear mixer.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

.

.

Figure A.9: Comparison between port extraction and full coupled methods for
dipole with diode mixer. The L2 − norm of the difference between the
two transient solutions is 1.171 · 10−16

. . . . . . . . . . . . . . . . . . . . 104

xi

Figure A.10:(a) Circuit diagram of the two-port microstrip interconnect coupled to
the diode attenuator circuit; (b) port voltage waveforms computed in
both TD and FD codes for Rc (cid:3) 10 kΩ. . . . . . . . . . . . . . . . . . . . . 105

Figure A.11:(a) Circuit diagram of the four-port hybrid-90◦ structure coupled with
the diode mixer circuit; (b) voltage waveform across the 86 Ω output
resistor calculated in both the TD and FD codes.

. . . . . . . . . . . . . . 106

Figure A.12:Balanced Mixer geometry . . . . . . . . . . . . . . . . . . . . . . . . . .
Figure B.1: Basis functions N1,2(r), N1,3(r), N2,3(r) for cell 1 with cubic path of the

. 107

particle. .

.

. .

. .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

Figure B.2: Two points are chosen to show the standard linear charge conservation

(left). Zoomed in showing curved analytic path and linear path (right).
Figure B.3: Basis functions N1,2(r), N1,3(r), N2,3(r) for cell 2 with cubic path of the

particle .

.

. .

. .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

117

xii

LIST OF ALGORITHMS

Algorithm 3.1: Particle Path Finding Algorithm . . . . . . . . . . . . . . . . . . . 42

Algorithm 3.2: Non-Linear Bi-Section Method . . . . . . . . . . . . . . . . . . . . 43

xiii

CHAPTER 1

INTRODUCTION AND BACKGROUND

1.1 Background and Motivation

The role of computational models or any other model is to be able to make predictions
about the world. Often these predictions are needed when experimental observations are
not feasible, too expensive, or take too long to perform. Formulating robust models of
physical phenomena can inform experiments, explain results and build insight into how
the world works. Plasma physics has utilized computational modeling for decades, arising
in the 1950’s [1]. The need to study plasma stems from it’s abundance in the universe and
role in various application in both physics and engineering. The computational challenges
that the plasma community face are non-linear, multi-physics and multi-scale. Various
methods have been developed to model plasma behavior, each aimed at different regimes
of plasma physics.

Here, we lay the groundwork that will put the chapters of this thesis in perspective. We
provide background on what a plasma is, where it can be found, what the applications are
in science and technology, and what kind of modeling has already been done. We will
explore a certain type of model called Particle-In-Cell and discuss what has been done in
this area and what specifically on the electromagnetic side can be improved.

A plasma, the fourth state of matter, is a gaseous type state where electrons disassociate
from ions and can move around. They can be formed in a variety of different manners
and can exhibit a variety of behaviors as well as various different chemistry. To provide
a general perspective, Fig. 1.1, shows that solids have rather strong chemical bonds that
result in a fixed shape and volume. Liquids have weaker bonds that conserve its volume
but will form to its container. Gases have very little, to no molecular bonds, have variable
volumes and will conform to the container. Plasma is similar to gases, in that the volume

1

is variable and will fit to the container that is being held. The exception is that there is
enough kinetic energy such that the electrons that disassociate from the ions move around
freely, not bound to any ion or molecule.

Deposition

Melting

Vaporization

Ionization

Solid

Liquid

Gas

Plasma

Freezing

Condensation

Recombination

Sublimation

Figure 1.1: Phase diagram of matter.

As an object changes from solid to plasma, the amount of kinetic energy needed
to change each phase is increased. One of the main differences between plasma and
the other three states of matter is that the interaction between particles in the system is
electromagnetically driven. In [2], this is described as collective effects. A stationary charged
particle produces an electric radial field. Other charged particles in the vicinity will feel a
force from this field.

Now that we know what plasma is, why do we care about it? One reason is that it is the
most abundant form of matter in the universe. One example of plasma that occurs in nature
is solar wind. The solar wind consists charged particles released from the atmosphere
of the sun, mostly in the form of electrons, protons, and alpha particles. These particles
are ejected from the sun’s atmosphere and thrown across the solar system, and can be
especially strong when solar flares occur. As the earth moves through the solar system, it
is constantly bombarded by the solar wind and solar flares. The Earth’s magnetic field acts
as a defense from this activity, protecting the lifeforms from the destructive potential of the
solar wind and flares. These charged particles interact with earth’s magnetic field and can

2

cause disturbances in it. The most visible results of these disturbances, is the formation
of the aurora borealis. A more subtle but equally important consequence of solar wind
is the affect on satellites. The solar wind can affect satellite orbits, shorten their lifetime,
and interfere with on-board electronics. Having a thorough understanding of how the
solar wind works and how it affects satellites is important for satellite protection. As the
world becomes more and more dependent on digital technology, satellites become essential
infrastructure. The stock market, airline travel and shipping are a few of the industries that
depend on satellites for daily functions.

Plasma not only occurs in nature, but can be harnessed to solve complicated engineering
challenges. One area that plasma research is constantly found, is in the background of
the communication system infrastructure. When communication system need to generate
large amounts of electromagnetic power, devices such as magnetrons, klystrons and other
vacuum tube devices, are used. The klystron, for example, is used to amplify a radio
frequency signal in communication relays, television transmitters, satellite communication,
radar and has been used in particle accelerators. Klystron utilizes an electron beam in a
vacuum tube that is perturbed by an the incoming electromagnetic signal. Understanding
how charged particles interact with electromagnetic waves is essential to understanding
how these devices work.

Ion thrusters is another interesting avenue of research for the past several decades. Ions
ejected from a spacecraft can be used to provide a small amount of thrust for communication
satellites. This plasma is controlled via electromagnetic fields and allows the space craft
to overcome small amounts of atmospheric drag as well as reorient themselves. These
systems originally were conceptualized in the 1960’s and first put into use in the 1970’s.
Since then they have been used on hundreds of satellites.

Fusion research, which hopes to one day become a sustainable source of energy, relies
on confining plasma using magnetic fields. Various forms of reactor designs (tokamak
and stellerators) rely on computational tools to model the physics involved before building

3

these expensive devices.

Another application is sputtering devices that are used in the chip and film industry.
They allow for doping various metals onto a substrate material. The manufacturing of
computer chips and applying film coatings on to many optical devices relies on these
devices. Designing a sputtering device relies on computational tools. Modeling chip
etching is important in the design and manufacturing of computer chips. How wide traces
can be made and how well they are etched relies on fundamentally understanding plasmas.
These are only a few examples of how plasma affects our daily lives. Research and
development of plasma has had a far reaching impact on our knowledge of the physical
world and development of various technologies.

Plasma, being a state of matter, can come in various forms and behave differently
depending on the density, temperature, chemical make, etc. Plasma is often characterized
by a few properties which dictate which techniques are used to model the physics. Consider
for a moment a bulk neutral plasma that consists of electrons and a heavy ion. Suppose
we perturb the electrons in such a way that there is a region of electrons gathered in one
location forming a slightly negative area. The ions where the electrons came from form a
slightly positive charged area. Upon release, the electrons would repel each other while
trying to restore equilibrium. Due to momentum, they overshoot and oscillations occur.
This plasma frequency, fp is a function of the species number density ηe, which for the
example presented here are electrons. The other parameters are charge q and mass me of
the charges species (electrons), and vacuum permitivity 
0.

(cid:19)2

(cid:18) ηeq2

me 
0

fp (cid:3)

1
2π

.

(1.1)

Another parameter that is import to characterizing plasma is Debye length. Consider
again a bulk neutral plasma consisting of electrons and heavy ions. Suppose an extra
electron is placed in the center of the neutral region. The electrons would repel the other
electrons around it, creating a slightly positive region around the electron, effectively

4

shielding the electron. This shielding limits the length of interactions of that particle. This
length is referred to as the Debye length and is a function of the temperature T and density,
with k being the boltzman constant.

(cid:19)2

(cid:18) 
0kT

ηeq2

λD (cid:3)

(1.2)

There are many regimes in which plasma can be modeled. Most of the decisions on how
each regime is modeled depends on the underlying assumptions. Fluid models for example
treat plasma as a fluid and solve Naiver-Stokes equation as well as Maxwell’s equations. On
the other end of the spectrum, molecular dynamics problems model individual particles
and the evolution of the system is dependant on pair-wise interactions. Particle-in-Cell
(PIC) treat plasma as macro-particles that are essentially a discrete approximation of the
phase space distribution function (PSDF) and interact through a mesh. Many hybrid
models will incorporate a combination of multiple methods to simulate different parts of
the physics. The focus of this thesis will be on collisionless plasma, modeled by PIC.

Collisionless plasma’s behavior can be described by the Vlasov equations,

∂t fα(t, r, v) + v · ∇ fα(t, r, v)+
[E(t, r) + v × B(t, r)] · ∇v fα(t, r, v) (cid:3) 0.,

q
m

(1.3)

where fα(t, r, v) is the PSDF for a given species α. Source such as current density and
charge density are described using the moments of the PSDF,

ρα(t, r) (cid:3) qα

fα(t, r, v)dv

(1.4a)

Jα(t, r) (cid:3) qα

v fα(t, r, v)dv.

(1.4b)
The fields, E(t, r) and B(t, r), in Vlasov equation are solutions to the Maxwell’s equations
Eq. 4.2.

Þ
Þ

− ∂B(t, r)
∂t

(cid:3) ∇ × E(t, r)

5

(1.5a)

∂D(t, r)

∂t

(cid:3) ∇ × H(t, r) − J(t, r)
∇ · B(t, r) (cid:3) 0

∇ · D(t, r) (cid:3) ρ(t, r)

(1.5b)

(1.5c)

(1.5d)

Solving the Vlasov equation directly as shown in Eq. 4.1 is extremely expensive mainly due
to the fα(t, r, v) being a 6 dimensional quantity. Instead to solve this problem, phase space
is sampled using a discrete set of "macro particles". These particles are then pushed around
using Newton’s law and the Lorentz force Fα (cid:3) qα(E + v × B) due to the electromagnetic
fields acting on the particles. The moving particles are mapped as sources for field
quantities and the cycle continues. These are the basis set of steps are the basic step for PIC.
PIC has been used extensively since the 1960’s [1]. It has been used in a variety of
applications and has many variations. The method is a time domain simulation of moving
charge (discrete phase space) which interacts electromagnetically, through the use of a grid
or mesh. The system consists of four parts that repeat each time step, particle push, mapping
current to the grid, field solver, and interpolation. The first step, particle push, is where
particles are moving using Newton’s equations of motion due to the force applied by
the Lorentz force. It important to note that the particles (discrete charge) represent a
collection of particles (phase space) and are often represented in the discrete world as delta
function, but in some cases shape function. After the particle has moved, it is mapped
backed to the grid as either charge or current depending on whether one use electrostatic
or electromagnetic solvers. Once the sources (charge or current) have been mapped to the
grid, the electromagnetic system is solved to find field values. The methods used in the
field solver have varied widely from electrostatic to solving Maxwell’s equations, the wave
equation and even potential based methods. Traditionally, finite difference time domain
(FDTD) has been the solver of choice, although as of late extension to finite element and
discontinuous Galerkin have been made. Once field values are known, the electromagnetic

6

fields can be interpolated at the location of the macro particle for use in the particle push
and the system repeats.

Particle Push

Map Fields to

Particles

Map particle in current

and charge

Update Maxwell’s

Equations

Figure 1.2: Basic particle in cell (PIC) system. Other details such as collision, emissions,
absorption etc. are left out for simplicity.

Certain features such as collisions and secondary emission are commonly incorporated
into PIC models but their discussion is beyond the scope of this thesis. For a robust look at
all the methods and development see [3]. For a general state of the art as of 2004, see [4].
While traditional FDTD has been used for PIC there has been a recent push to develop
finite element based PIC methods. The reasons are many-fold. FDTD while efficient has
trouble modeling complicated object and suffers from the stair stepping approximations
of curved surfaces. Often small mesh sizes needed describe the geometry results in
prohibitively small time step sizes due to the Courant–Friedrichs–Lewy (CFL) condition.
The (CFL) condition ties the smallest element in a mesh to the time step size that can be
taken [5]. To over come this methods such as cut cells have been developed that help model
boundaries.

Finite element representation of the fields allows for better geometry representation.
Tetrahedral elements can easily allows for various cell sizes to resolve finer features without
the need for the entire mesh to be high resolution. Finite element also can be extended
to higher order basis functions in time, space and geometry. Implicit solvers such as
Newmark-Beta alleviate the mesh dependant nature of explicit solver.

Finite elements method (FEM), while a robust candidate for adaption to PIC, has a few
challenges that must be addressed. When mapping sources in PIC one must take special

7

care to properly ensure charge is conserved. Charge conservation is ensuring that Gauss’s
law holds and that spurious charge doesn’t exist on the mesh. Special care must be taken
when developing a PIC methods such that this constraint is held to machine precision.
Non-charge conserving schemes can quickly build up fictitious charge on the mesh that can
impact particle trajectory. The enforcement of the continuity equation is the key to being
able to conserve charge. Divergence cleaning is a common method to ensure that fictitious
charge does not build up on the mesh by adding correcting term by solving Poisson’s
equation and removing the charge. Another way to properly ensure charge conservation
is by properly mapping current onto the mesh consistent with the underlying basis sets.
Charge conserving schemes for FEM-PIC have been implemented [6], but the mostly focus
on explicit schemes. This brings us to the underlying problems addressed in this thesis.

1.1.1 Outline of Thesis

The over arching theme of this thesis is to address the challenges in enforcing Gauss’s

law with 3D finite element particle in cell methods.

• Chapter 2: focuses on developing a test bank that allows for methods developed to
be validated and verified. One of the first challenges we ran into in developing a
3D particle-in-cell codes is the lack of robust simple simulation that have analytic or
quasi-analytic solutions that can test us of the code is correct. This paper add to the
literature of validation and verification methods but focuses on 3D. The chapters that
follow refer to the test presented here for comparison between methods.

• Chapter 3: Building upon existing charge conserving scheme as we develops a
methodology for enforcing charge conservation that is agonist the time stepping
methods. We apply the methodology to both Maxwell and the wave equation and
demonstrate the effectiveness on a non-neutral beam with space charge and an
adiabatic expansion of a neutral plasma.

8

• Chapter 4: Expanding on chapter 3, implicit FEM transient solvers support an
underlying null space that corresponds to a gradient of a scalar potential ∇Φ(r)
(t∇Φ(r) in the case of wave equation solvers). While explicit schemes do not, they are
only conditionally stable, time step sizes are mesh dependent, and very small. A way
to overcome this bottleneck, and indeed, better satisfy all four Maxwell’s equation
is to use a quasi-Helmholtz formulation on a tessellation. In the re-formulation
presented, we strictly satisfy both the equation of continuity and Gauss’laws for both
the electric and magnetic flux densities. Results demonstrating the efficacy of this
scheme will be presented.

9

CHAPTER 2

A SET OF BENCHMARK TESTS FOR VALIDATION OF 3D PARTICLE IN CELL

METHODS

2.1 Abstract

While the particle-in-cell (PIC) method is quite mature, verification and validation of
both newly developed methods and individual codes has largely focused on an idiosyncratic
choice of a few test cases. Many of these test cases involve either one- or two-dimensional
simulations. This is either due to availability of (quasi) analytic solutions or historical
reasons. Additionally, tests often focus on investigation of particular physics problems,
such as particle emission or collisions, and do not necessarily study the combined impact
of the suite of algorithms necessary for a full featured PIC code. As three dimensional
(3D) codes become the norm, there is a lack of benchmarks test that can establish the
validity of these codes; existing papers either do not delve into the details of the numerical
experiment or provide other measurable numeric metrics (such as noise) that are outcomes
of the simulation. This paper seeks to provide several test cases that can be used for
validation and bench-marking of particle in cell codes in 3D. We focus on examples that
are collisionless, and can be run with a reasonable amount of computational power. Three
test cases are presented in significant detail; these include, basic particle motion, beam
expansion, adiabatic expansion of plasma. All presented cases are compared either against
existing analytical data or other codes.

We anticipate that these cases should help fill the void of bench-marking and validation

problems and help the development of new particle in cell codes.

10

2.2 Introduction

Simulation of space charge and plasma is an integral part of many scientific and
engineering processes, with specific applications of pulsed power, particle accelerators,
integrated chip manufacture, satellites, and medicine, to name a few examples. A number
of methods have been used to simulate the underlying physics; these range from molecular
dynamics, many-body, global-model, 2-fluid, magneto-hydrodynamics (MHD), etc. [7].
The key with all of these methods is the self-consistent evolution of the fields and plasma.
In the area of kinetic computational plasma physics, this means evolving the particle
distribution function in a manner consistent with dynamic electric and magnetic fields.
One such method, the Particle-in-Cell (PIC), does this by (a) mapping all particle-field
interactions through a mesh, and (b) using macro-particles as statistically significant
markers of the distribution function’s phasespace [1, 8]. PIC has been around for several
decades, originating in the 1960’s. Much of the early work focused on using finite-difference
time-domain (FDTD) for electromagnetic fields solutions on a regular grid, more specifically,
Yee Cells [9]. This method has become popular due to its simplicity and scalabilty [10],[11].
Additionally, these methods are local, which makes parallel implementation relatively
straight forward. While efficient, curved surfaces, or any geometry that is not aligned with
a regular structured grid, is poorly represented. Significant work has been done to adapt
Cartesian grids with better model boundaries. The most popular is the Dey-Mittra adaption
of FDTD to PIC [12]. It should be noted that there have been a number of other advances
in conformal FDTD[3], but they have not made their way into PIC codes. An alternative is
to use an axisymmetric methods when possible by treating a three dimensional problem in
a two-dimensional setting. An excellent example is XOOPIC, which is open source and has
been used extensively [13]. For an extensive review, see [3] and references therein.

Since the introduction of FDTD, the electromagnetic fields community has invested
time and resources to developing high fidelity field solvers. These include finite element
method (FEM) and discontinuous Galerkin methods (DG). Both can provide higher order

11

representation of fields and geometry. Additionally, one has better understanding of
efficient time marching schemes [14]. Here, our focus is a finite element based Maxwell
solver, that has been a) adapted to PIC and b) shown to be devoid of several problems
associated with charge conservation.

The challenge with using unstructured grids, using FEM or DG, is developing a
consistent scheme to solve the coupled problem of evolving fields and particle population.
Given the relative lack of maturity of FEM-PIC, it is important to understand how the
nuances of these algorithms affect the accuracy of the physics that is predicted. Small
variations in method can have significantly different prediction; for instance, the time
stepping method, solving the first order Maxwell’s equations vs. solving the second order
wave equation, order of spatial basis functions, order of time basis function, etc., can yield
different predictions. Perhaps, the most critical of these nuances is charge conservation in
that it explicitly ties together the evolution of the fields and the particles. Additionally,
while there has been extensive efforts for FEM in the context of source-free electromagnetics,
FEM-PIC is much less developed. The most recent of these is the work in [6] where a
charge conserving scheme for Whitney basis was presented. This has been developed
further by Teixeira [15].

This begs the question as to how PIC codes are validated. Correctness of a simulation
has typically been done through a number of methods such as simple property tests such
as energy conservation, to more rigorous tests such as rates of convergence. These test
are described in more rigor in [16]. The most rigorous tests require a comparison with
analytical solutions. But, this is much harder to come by in PIC settings, as analytic solution
for most problems do not exist. Most tests rely either demonstrating proper behavior, or
comparison to another code. As a result, there is a need for a set of test examples that can be
used to validate/benchmark and see how these variations affect the solution of the solver.
Validation and verification methods play import roles in quantifying error, numerical noise,
and comparison of method performance. Of late, in the PIC community, there has been

12

a push for more robust bench marking and validation of (PIC) codes [17, 18, 16, 19, 20].
Much of this current literature focuses on discharge and particle emission and collisions.
Less attention is paid to test examples where the grid-particle interaction dominate the
numerical effects. This is the void we seek to fill.

In collision-less PIC, there are a number of metrics one may want to measure such
as field noise, charge conservation and energy fluctuation. A set of benchmark tests are
needed to measure each of these metrics. In this paper, we have devised a set of test cases
to measure various metrics of particle in cell codes to both verify the code and allow for
comparisons between codes. The four cases presented here have (quasi)-analytic solutions,
and can be verified against other methods such as an electrostatic PIC codes. We provide
details of each case such that comparison and validation is easy. Each of these cases are
a 3-D examples that can be set up with simple boundary conditions – perfect electric
conductor (PEC) or first order absorbing boundary condition (ABC). The domain size to
capture the physics of the problem is limited such that large amounts of computational
resources are not needed, but can be scaled in a parallel context, to test algorithm scaling.

2.3 Problem Description

Next, we present a succinct outline of the problem description and implementation
details for the test cases presented. Consider a domain Ω enclosed by a boundary ∂Ω. The
domain consists of a free space permittivity and permeability defined by 
0 and µ0, and
can contain multiple charge species. Both the fields and charge distribution evolve over
time self consistently with accordance to both the Vlasov Eq. (4.1) and Maxwell’s Eq. (4.2),

∂t f(t, r, v) + v · ∇ f(t, r, v)+
[E(t, r) + v × B(t, r)] · ∇v f(t, r, v) (cid:3) 0.

(2.1)

q
m

Here, f(t, r, v) is the phase space distribution function. Typically, one does not solve Eq.
(4.1) directly. The usual approach, that we follow as well, is to define charge and current
Ω f(t, r, v)dv

via the first and second moment of this distribution function; i.e., ρ(t, r) (cid:3) qÞ

13

and J(t, r) (cid:3) qÞ

Ω v(t) f(t, r, v)dv. These are then solved self-consistently with Eq. (4.2) and

the particle location are updated using the Lorentz force F (cid:3) q(E + v × B).

∂D(t, r)

(cid:3) ∇ × E(t, r)

− ∂B(t, r)
∂t
(cid:3) ∇ × H(t, r) − J(t, r)
∇ · B(t, r) (cid:3) 0

∂t

∇ · D(t, r) (cid:3) ρ(t, r)

(2.2a)

(2.2b)

(2.2c)

(2.2d)

We assume that the boundary can be subdivided into either Dirichlet ΓD, Nuemann ΓN or
impedance ΓI, boundary conditions as defined by,

(2.3a)

(2.3b)

ˆn × E(r, t) (cid:3) ΨD(r, t) on ΓD,
ˆn × B(r, t)
(cid:3) ΨN(r, t) on ΓN ,
− Y ˆn × ˆn × E(r, t) (cid:3) ΨI(r, t) on ΓI.

µ

µ

ˆn × B(r, t)

where ˆn is an outward pointing normal to ΓD, ΓN or ΓI. We define Y (cid:3) (cid:112)
0/µ0 as

the free space surface admittance, and ΨD(r, t), ΨN(r, t), ΨI(r, t) are boundary condition
functions. Initial conditions must satisfy Gauss’s electric and magnetic law, such that
∇ · D(0, r) (cid:3) ρ(0, r) and ∇ · B(0, r) (cid:3) 0. The algorithm’s preservation of these conditions
can then be tested by our problems.

(2.3c)

2.3.1 The Discrete Domain

To solve the system, we use a 3D finite element Maxwell solver. Consider a domain Ω
broken into an irregular non-overlapping tetrahedral elements with boundaries ∂Ω. In
each tetrahedral element the field quantities E and B can be represented using Whitney
spaces. The moments of the distribution function can be represented using using delta

14

functions, shown in Eq. (2.4a),(2.4b),

ρα(t, r) (cid:3) qα

δ(r − rp(t))

Np

Np


p

p

Jα(t, r) (cid:3) qα

vp(t)δ(r − rp(t))

(2.4a)

(2.4b)

(1)
i

for each charge species α, where qα is the charge of the macro-particle and Np is the number
of macro-particles. The system proceeds as follows. The Whitney basis functions are used
to convert Maxwell’s continuous equation to discrete ones. Current and charge are mapped

the solution to Faraday’s and Ampere’s law uses a leap frog scheme. Details to pertinent
(r) and

(cid:1) vector basis set and Whitney-0(cid:0)W

to the grid via Whitney-1(cid:0)W
to each stage of the scheme can be found in [6, 15]. Using E(t, r) (cid:3)
Ne
B(t, r) (cid:3)
N f
En+1 (cid:3) En + ∆t[(cid:63)
]−1 ·(cid:18)[Dcurl]T · [(cid:63)µ−1] · Bn+1/2 − Jn+1/2

(r) and Galerkin testing results in,

Bn+1/2 (cid:3) Bn−1/2 − ∆t[Dcurl] · En

i(cid:3)1 bi(t)W

(cid:1) scalar basis set, and

i(cid:3)1 ei(t)W
(cid:19)

(2.5)

(0)
j

(2)
i

(1)
i

where Ne and N f are the number of edges and faces, ∆t is the time step size,

1

En (cid:3) [e n

, e n
2
Bn+1/2 (cid:3) [bn+1/2
Jn+1/2 (cid:3) [jn+1/2

],
, ..., e n
Ne
, e n+1/2
, ..., bn+1/2
],
, ..., jn+1/2
, jn+1/2
]

N f

1

2

1

2

Ne

with En, Bn+1/2, and Jn+1/2 being the field coefficients at time step n and n + 1/2. The
Hodge and discrete curl matrices in Eq. (3.2a) and Eq. (3.2b) are,

(1)
(1)
[(cid:63)
]i,j (cid:3) (cid:104)W
(r), ε · W
(r)(cid:105)
i
j
(2)
(2)
[(cid:63)µ−1]i,j (cid:3) (cid:104)W
(r), µ−1 · W
(r)(cid:105)
j
i
(1)
(r)(cid:105)
[Dcurl]i,j (cid:3) (cid:104)ni, ∇ × W
j

15

(2.6)

(2.7)

(2.8)

where, [Dcurl] is the discrete curl matrix, [(cid:63)µ−1] is the Hodge star matrix, [(cid:63)
] is the Hodge
epsilon matrix. The current is mapped using a charge conserving method presented in
[6, 15],

Þ rn+1

p

jn+1/2

i

(cid:3)

qp
∆t
qp
∆t

(1)
i

W
(rn
p)W

(r) · dl
(0)
(rn+1
j

p

rn
p

(0)
i

(2.9)

(2.10)

(cid:3)

[W

(rn
p)],
(r) is the barycentric coordinate for node i within a cell.

) − W

(rn+1

)W

(0)
i

(0)
j

p

where W

(0)
i

2.4 Benchmark Tests

In this section, we present four cases that could be used for validation. In each of these
cases, we will present a number of metrics and discuss what would happen if some of these
metrics were not to be satisfied. Each test is designed to highlight an aspect of the method.
The first test case, basic particle motion, highlights the particle push’s performance. The
expanding particle beam is setup to verify current mapping, discrete Gauss’s magnetic
and electric law, and the discrete continuity equation. The adiabatic expanding plasma
test case allows for multi-species evolving distributions and can be scaled up to test the
scalability of methods.

2.4.1 Basic Particle Motion

The first set of tests of any PIC method should start with validating the particle motion.
This is done by subjecting a single particle to external electric and magnetic fields to
determine if the motion calculated is correct. To this end, we set up two tests, 1) a simple
acceleration test by a uniform electric field and 2) a rotation test in a uniform magnetic field.
We use Newton’s law and Lorentz Force to obtain position and velocity at every time point.
Such a test verifies that the trajectory of the particle — as mapped on a unstructured mesh
— is computed correctly. Both cases presented here have analytic solutions that allows for
error in particle trajectory to be related to time step size. It is advisable for both particle

16

motion tests to be run twice. Once with just the external background fields, and a second
time assuming radiation, that is solving Newton and Maxwell’s equation self consistently,
which we refer to as being fully coupled. The second run should use a macro-particle with
a small amount of charge such as an electron such that the radiated fields are significantly
smaller than the constant background fields. This choice causes the main source of error in
the particle trajectories to be dominated by the particle push and not the field solver. The
solutions between the first and second run should be approximately identical and allow
for verification that the machinery of the method is correct. This setup also presents an
opportunity to verify current continuity and the discrete Gauss’s law in a setup where
the particle path is for all intents and purposes is known. Although the field values are
small, they are non-zero. This case is simple enough that it can be performed on a course
unstructured mesh, where large amount of computational resources are not necessary.
Both the external E-field and B-field cases presented below are performed on a course
unstructured meshed 1m cube with a minimum edge length of 0.11 m and max edge length
0.33 m

External E-Field The first term of the Lorentz force can be checked by placing a particle
at rest in an uniform background electric field of 1 V/m and letting the particle accelerate
as shown in Fig. 2.1. The macro-particle has the mass m and charge q of a single electron.
We calculate the analytic motion in x direction and compare it with the computed push,

x (cid:3)

q
2m

Ext2.

For the sake of brevity we only show error analysis and convergence for the magnetic field
motion below.

External B-Field To check the v × B term, we apply a uniform external magnetic flux
density B (cid:3) 6.82272 · 10−5ˆz Wb/m2 with no electric field. A single particle with the mass

17

and charge of an single electron is given an initial velocity v (cid:3) 3 · 106 ˆy m/s and position
r (cid:3) [0.75, 0.5, 0] which results in a gyro-radius of 0.25m, as shown in Fig. 2.2.

Table 2.1: Cyclotron Motion

Parameter

B
Q
m
v0
r0

Value

6.822756 · 10−5ˆz T
−1.60217646 · 10−19 C
9.10938370 · 10−31 kg
3 · 106 ˆy m/s
[0.75, 0.5, 0.0] m

We analyze the error for this test case under four scenarios shown in Fig. 2.3. The four
cases are particle cyclotron motion due to a fully coupled system labeled mesh and the
motion due to only external background fields at three different time step sizes ∆t, 0.1∆t
and 10∆t. These four cases were run to test that our full coupled systems mechanics are
working. The error is calculated by finding the distance between the analytic position
ra(i∆t) and the simulated position ri at time i∆t, and distance di (cid:3) ri − ra(i∆t),

error (cid:3)(cid:112)di · di.

The analytic position for our set up is [21],
cos(ωct) +
sin(ωct) +

ra,x(t) (cid:3)
ra,y(t) (cid:3)

vy(0)
ωc
vy(0)
ωc

vx(0)
ωc
vx(0)
ωc

sin(ωct) − vy(0)
cos(ωct) − vx(0)

ωc

ωc

(2.11)

(2.12a)

(2.12b)

+ ra,x(0)
+ ra,y(0).

Results from this test are shown in Fig. 2.3. The two dips in the error are from the particle
starting at the same location — zero error — and from when the particle rotates back to
the starting point. In this specific case the volume preserving nature of the Boris push is
what causes this second dip. We can also look at the convergence rate of relative error.
This becomes important when higher order methods are used. In Fig. 2.4, we compare the
relative error over one cycle of cyclotron motion vs. the number of step per cycle. This
was setup in the same manor as Fig. 2.3, but we vary the time step size such that between

18

Figure 2.1: Particle acceleration due to electric field.

100 and 20000 steps are taken to complete a single rotation. The error is calculated by
taking the relative error for all the position. This is calculated by taking the L2-norm of all
√
the distance errors ¯d (cid:3) [√
dn · dn], divided by the L2-norm of the
d1 · d1, ··· ,
analytic positions ¯ra (cid:3) [ra(0), ra(∆t), ··· , ra(n∆t)]

d0 · d0,

√

relative error (cid:3)

|| ¯d||2
||¯ra||2

.

(2.13)

The slope of the graph is first order, which is also reflected in Fig. 2.3 as a time step 10 time
smaller produced an error that is 10 times smaller.

19

Figure 2.2: Cyclotron particle motion of an electron in a magnetic field.

Figure 2.3: Cyclotron motion error from just Boris push with several different time step
sizes, as well as influence from the mesh. The time step ∆t (cid:3) 46.1766ps or 11340 steps per
cycle.

20

Figure 2.4: Error as time step size decreases for Boris push.

2.4.2 Expanding Particle Beam

In the second test, a beam of electrons are injected into a cylindrical drift tube and is
evolved over time. The reasons for this test are, first it has an approximate solution found
in [22]. This is important because if the method being developed does not come close to
the approximate solution then we know something is wrong. The second reason is charge
conservation. It is easy to measure charge conservation in a beam case since any error
tends to accumulate due to the charge motion being one directional. Any issue also shows
up in the solution as higher noise levels or spatial striations in the beam [23].

We create a cylindrical drift tube geometry as shown in Fig. 2.6 and 2.5, with a radius
of 0.02 m and a length of 0.1 m. We set all exterior surfaces to perfect conductors (PEC)
with the interior volume set to vacuum. The drift tube is discretized such that the average
edge length is 4.11 mm. The beam has a voltage and current of 7.107 kV and 0.25 A. This
setup amounts to a beam radius of 8 mm centered in the cylinder with a starting velocity of
v0 (cid:3) [0.0, 0.0, 5 · 107] m/s. The beam is ramped up by changing the macroparticle weight
from 0 to 52012 electrons per macroparticle linearly over the first 2 ns of the simulation.
The parameters are listed in Table 4.1. The current in this example was chosen such that a
converged solutions of the expansion could be represented with a coarse mesh to reduced

21

Figure 2.5: Drift tube for particle beam test case

Table 2.2: Expanding Particle Beam Parameters

vp
vp/c

beam radius rb

Number particles per time step

Parameter

Cavity Radius
Cavity Length

Boundary Conditions

species

Turn on time
beam current

macro-particle size
min edge length
max edge length

Mesh Type

Value
20 mm
100 mm
5 · 107 m/s
0.16678
8.00 mm

PEC

10

electrons

2 ns
0.25 A
52012.58
1.529 mm
6.872 mm

Unstructured

the simulation time needed to solve the problem. We also compare the results with the
2-D FDTD code XOOPIC in r-θ mode. As is evident, the agreement between the two codes
— despite obvious advantages afforded by asymmetry in XOOPIC — agree well with each
other. The trajectories of the particle is shown in Fig. 2.6, and there is excellent agreement

22

between FEMPIC and XOOPIC. In Table 2.3, we show the beam profile for both FEMPIC as
well as XOOPIC. The profile is found by binning the particle positions such that each 10
mm section of the beam in the z direction is a different bin. In each bin the particle with the
largest radius is found. This particle is going to be very close to the downstream end of the
bin range. We then compare these radius values between both methods to find the relative
error in beam radius. While FEMPIC and XOOPIC agree with each other, with a relative
error in the radius of 0.0063 as the beam exits, they differ slightly from the approximate
analytic data. The approximate analytic solution derived in [22] is also shown in Fig. 2.6, it
assumes an idealized beam, with uniform density radially and ignores the perfect electric
conductor (PEC) boundary conditions on the entrance and exit surfaces of the beam. We
postulate that this leads to the slight disagreement. Note, it is difficult, if not impossible
to create the exactly identical conditions for either the numerical simulation or analytical
analysis. Nerveless, we include this data to show that the same trends are captured.

Next, the kinetic and potential energy in both the electric and magnetic field obtained
by FEMPIC and XOOPIC agree as shown in Fig. 2.7. We also examined field noise. Note,
the boundaries are perfect electrically conducting, and there are no losses in this problem
with the exception of particles leaving the domain. This results in cavity modes that are
excited. The radial component of the electric field, Er, at halfway through the tube at a
radial distance of 18 mm from the center of the tube is shown for three cases in Fig. 2.8.
Each case uses a different random number generator seed for the location of where the
particle in the beam is injected. This allows us to see if the noise is random or does it fall
into certain modes. Various modes of the cylindrical cavity are shown in Fig. 2.9. These
show up upon taking the Fourier transform of the noise shown in Fig. 2.8, after the system
reaches steady state at 6.66ns and the DC component had been removed. Non-random
injection does artificially reduces the noise of the system.

Discrete charge continuity is an import metric to measure in a PIC method. Current
that doesn’t satisfy the continuity equation allows for build up of fictitious charge on the

23

Figure 2.6: XOOPIC vs FEMPIC particle position at steady state. Approximate analytic
solution profile is shown by the dashed line.

Table 2.3: Beam Profile

z (mm) FEMPIC r (mm) XOOPIC r (mm) Relative Error

10
20
30
40
50
60
70
80
90
100

8.0119
8.0673
8.2028
8.3484
8.6291
8.9395
9.3018
9.7719
10.1852
10.6962

8.0012
8.0432
8.1622
8.3284
8.5577
8.8722
9.2343
9.6360
10.0949
10.6294

0.0013
0.0029
0.0050
0.0024
0.0083
0.0076
0.0073
0.0141
0.0089
0.0063

mesh. For details on how charge conservation works and how it measured see [15]. The
discrete continuity equation for this problems was measured as a function of time in Fig.
2.10, it is evident charge is conserved to machine precision.

24

Figure 2.7: Particle Beam Energy with a 2ns turn on time

Figure 2.8: Noise measured from three beam test cases using FEMPIC method with
different random number generator seeds.

25

Figure 2.9: Fourier transform of the electric field from FEMPIC method for three beam
tests that had different random number generator seeds.

Figure 2.10: Charge Conservation for particle beam test case measured by |(ρn+1− ρn)/∆t −
[∇·]Jn+1/2|/|(ρn+1 − ρn)/∆t|/Npart

26

2.4.3 Adiabatic Expanding Plasmas

In the third test, we simulate an adiabatically expanding plasma ball with a Gaussian
distribution in the radial direction. This experiment has an approximate analytic solutions
[24] that can used for validation. These analytic solutions have been compared with
experimental results in [25]. In earlier examples, the analytical examples are obtained in
1-D and 2-D, making it difficult to show a proper expansion with 3-D numerical code [15].
While many of the goals of this test are similar to the beam case, having this example in
mind is important. This test case will also be important in testing multi-scale codes as well.
As the sphere expands in time beyond what is shown here more computational power
will be needed to resolve the solution. As analytical data is available in 3-D, we design
the test to reflect the experiment done in [25] but change the density such that the Debye
length can be fully resolved with modest computational resources. The expanding plasma
has two species, Sr+ ions with an initial temperature of 1 K and electrons at 100 K. The
initial number density of both species is 5e8 particle per m−3, and the initial σ (cid:3) 0.01m
such that the system is initial charge neutral everywhere. Other details of the simulation

Figure 2.11: Distribution of electrons radially at two instants in time with analytic solutions.

27

Table 2.4: Adiabatic Expanding Plasmas

Boundary Conditions

First order ABC

Parameter

Mesh Radius

Tion

Telectron

Species

Number Particles

Macro-Particle Size
Min Edge Length
Max Edge Length

Mesh Type

Value
6mm

1K
100K
8000

Electrons and Sr +

52012.58
1.529mm
6.872mm

Unstructured

parameters can be found in Table 3.3. When the simulation starts, the electrons expand
outward faster than ions, due to the higher temperature, creating a radial electric field.
After the initial expansion outward, the growing electric field slows the expansion of
the electrons. At the same time the ions expand outward pulled by the electrons. This
example can be scaled – by changing the density of charges – such that various levels of
computational power are needed thus making it a good test of the scalability of a code.
We use an first order absorbing boundary condition (ABC) to truncate the domain and
place the boundary conditions sufficiently far from the initial distribution to capture the
initial expansion of the electrons but not exit the mesh. Early electron expansion can be
validated in small geometries, which we do here. The analytic distribution of electrons at
initial conditions and at time 184 ns, are show in Fig. 2.11. This was derived in [25], and is
succinctly presented here. The density function fα of species α can be expanded as,

(cid:20)

(cid:21)3/2

(cid:26)− r2

(cid:27)

fα(t, r, v) (cid:3) n0

mα

2πkbTα(0)

exp

2σα(t)2

− mαv2
2kbTα(t)

.

(2.14)

The radial distribution is a symmetric Gaussian distribution function with a deviation
σα(t) is a function of time,

α(t) (cid:3) σ2
σ2

α(0)(1 + t2/τ2

exp).

(2.15)

28

Table 2.5: Adiabatic Expanding Plasma
r (mm) FEMPIC ηe m−3 Analytic ηe m−3
2.6559 · 108
0.8692
2.5551 · 108
3.1874
2.3445 · 108
5.5055
2.0517 · 108
7.8236
1.7125 · 108
10.1418
1.3632 · 108
12.4599
1.0350 · 108
14.7780
7.4947 · 107
17.0962
5.1761 · 107
19.4143
3.4094 · 107
21.7324
2.1419 · 107
24.0506

3.0924 · 108
2.5357 · 108
2.0866 · 108
1.9706 · 108
1.6154 · 108
1.1742 · 108
8.5786 · 107
6.1211 · 107
4.3343 · 107
2.6577 · 107
1.6639 · 107

The velocity distribution depends on the temperature Tα(t) of species α,

Both σα(t) and Tα(t) depend on the characteristic expansion time τexp

,

(2.16)

Tα(t) (cid:3)

Tα(0)
1 + t2
τ2
exp

(cid:115)

mασ2α(0)

(cid:21)3/2

(cid:20) Tα(t)

Tα(0)

τexp (cid:3)

(2.17)
Integrating velocity space, we can write the number density η(t, r) as a function of time
and radius from center of the sphere,

kb[Te(0) + Ti(0)] .

η(t, r) (cid:3) η0

−r2
2σ2α .

e

(2.18)

The results Fig. 2.11. are show in Table 3.3 for easy reference.

2.5 Summary

In this paper, we have provided three tests that can be used to validate 3D PIC codes.
These tests have been developed as (a) there is significant recent interest in developing
PIC codes to fully resolve complex boundaries, as properties of these boundaries, and (b)
as there is a void in examples with details that can be used to validate these codes. In
all examples that we have developed, there exist quasi-analytic solutions, as well as open

29

source codes that can be used for benchmarking results. We have also provided sufficient
detail so as to facilitate the user to examine all possible metrics that may subtly affect the
physics.

30

TIME INTEGRATOR AGNOSTIC CHARGE CONSERVING FINITE ELEMENT PIC

CHAPTER 3

3.1 Abstract

Developing particle-in-cell (PIC) methods using finite element basis sets, and without
auxiliary divergence cleaning methods, was a long standing problem until recently. It
was shown that if consistent spatial basis functions are used, one can indeed create a
methodology that was charge conserving, albeit using a leap-frog time stepping method.
While this is a significant advance, leap frog schemes are only conditionally stable and time
step sizes are closely tied to the underlying mesh. Ideally, to take full advantage of advances
in finite element methods (FEMs), one needs a charge conserving PIC methodology that
is agnostic to the time stepping method. This is the principal contribution of this paper.
In what follows, we shall develop this methodology, prove that both charge and Gauss’
laws are discretely satisfied at every time step, provide the necessary details to implement
this methodology for both the wave equation FEM and Maxwell Solver FEM, and finally
demonstrate its efficacy on a suite of test problems. The method will be demonstrated by
single particle evolution, non-neutral beams with space-charge, and adiabatic expansion of
a neutral plasma, where the debye length has been resolved, and real mass ratios are used.

3.2 Introduction

Simulation of space charge and plasmas is critical to a number of areas in science and
engineering. These range from, applications of pulsed power to particle accelerators to
satellites and medicine [26, 27, 28]. The means to do so has largely relied on Particle-
in-cell (PIC) methods. PIC has been around since the 1950’s and is a popular methods
of modeling plasma and space charge due to its simplicity and ease of use [29]. PIC
enables a self consistent solution to Maxwell’s equation and equations of motion for

31

charged species. Traditionally, PIC is based on finite difference time domain to evolve
fields [4]. The use of regular cubical grids presents challenges, especially in modeling
complex geometry. Modeling curved features requires small cell sizes, and this results
in a stair-stepped approximation of the desired geometry as well as small time steps in
keeping with the Courant–Friedrichs–Lewy condition. Using cut-cells has improved the
geometry representation by allowing boundaries to cut across cells [30]. Complex and
fine features, as well as multi-scale objects, require the use of a prohibitively expensive
number of small cells for high fidelity simulations. As a result of these challenges, there has
been persistent investigation into the use of more sophisticated field evolution techniques
[31, 32, 33, 3]. A natural choice is using time domain finite-element method (TDFEM)
due to (a) unconditionally stable time stepping methods, (b) ability to model complex
geometries, and (c) well developed extensions to higher order (both in representation of
fields and geometry) [14].

While TDFEM can be thought of as a panacea for modeling complex geometries, it is
not so for crucial quantities that must be conserved. These include Gauss’ law and charge
conservation. Indeed, developing a numerical scheme that implicitly conserved charge
was an unsolved problem until [6, 34]. Prior to this development, one used divergence
cleaning methods to remove spurious charge accumulation [35]. The key to realizing
charge conservation relied on (a) following the de-Rham sequence to represent physical
quantities on a mesh and (b) use explicit time stepping methods. A more recent paper
prescribes three conditions must be satisfied by self-consistent charge conserving schemes
[36]; this assertion is proved and illustrated for different PIC schemes. The TDFEM-PIC
method relies on Maxwell solvers, in that one solves Maxwell’s first order equations as
opposed to the wave equation, and leap-frog time stepping. The structure of the solver is
such that one avoids a time growing null space corresponding to DC modes. Unfortunately,
leap frog is only conditionally stable. As a result, there is a limit on the time-step sizes that
one can take, and this closely tied to the underlying discretization. In classical TDFEM,

32

this has been overcome using Newmark-beta time stepping, which is second order and
unconditionally stable. Unfortunately, implicit time stepping poses a number of challenges
to satisfaction of conservation laws that must be satisfied and is an open problem[37]. This
paper provides the theoretical framework for resolving this bottleneck.

Implicit time stepping permits taking significantly larger time steps, un-constrained
by the mesh; and unconditional stability is an added bonus. Unfortunately, as will be
evident in the paper, applying these directly to TDFEM-PIC violates both Gauss’ law and
the equation of continuity. In addition, in solving the field equations, one needs to evolve
the locations of particles over time via Newton’s laws. A larger time step size, implies that
additional infrastructure needs to be in place to accurately compute all aspects of particle
trajectory (including information necessary to map it back on the mesh). Resolution to
these challenges associated implicit time stepping with a TDFEM framework will be the
main contribution of this paper. We will

1. Develop the methods to ensure that both Gauss’ law and equation of continuity is

satisfied for implicit methods. The methods rely on insight provided in Ref. [36].

2. We will show that the proposed method is agnostic to time stepping schemes.

3. We will develop methods to evolve particle parameters (path, velocity along the path,

and mapping path to the mesh).

4. Finally, we will present results validating these methods for both the Maxwell and

wave equation TDFEM solvers.

Our hope is to present the technique with sufficient lucidity such that they can be retrofitted
with existing codes.

The rest of this paper is organized as follows: In the next Section, we present an overall
rubric of implicit TDFEM solvers (both Maxwell and wave), and why direct application
of implicit time stepping fails to conserve quantities. Next, in Section 3.4, we present

33

details on how these may be modified so as to conserve charge, satisfy Gauss’ law, and be
independent of time stepping approach. In addition, we present details of the method used
to evolve particle parameters. In Section 4.7, we present a number of results that validate
our claims. Finally, we conclude this paper in Section 3.6 outlining future directions of
research.

3.3 Preliminaries

Consider a domain Ω whose boundaries are denoted by ∂Ω. It is assumed that the
domain comprise charged species that exist in a background medium defined by ε0 and
µ0, the permittivity and permeability of free space, and the speed of light denoted using
c (cid:3) 1/√
µ0ε0; for simplicity of the exposition, we consider only one species. It is also
assumed that there exists an electromagnetic field, both impressed and arising from motion
of the charged species. Both the fields and the charged species evolve in time. The
distribution of charge can be represented by a phase space distribution function (PSDF)
f(t, r, v) that satisfies the Vlasov equation

∂t f(t, r, v) + v · ∇ f(t, r, v)+
[E(t, r) + v × B(t, r)] · ∇v f(t, r, v) (cid:3) 0.

q
m

(3.1)

While we do not solve this equation directly, our approach is conventional in that we make
a particle approximation for the PSDF in (4.1).

3.3.1 Overview of Method

defined as ρ(t, r) (cid:3) qÞ

Ω f(t, r, v)dv and J(t, r) (cid:3) qÞ

Using this PSDF, we follow the conventional definition of the charge and current density
Ω v(t) f(t, r, v)dv as moments of the
PSDF. The fields, E(t, r) and B(t, r), in the Vaslov equation are solutions to Maxwell’s curl
equations with the sources (charge and currents) defined earlier

− ∂B(t, r)
∂t

(cid:3) ∇ × E(t, r)

34

(3.2a)

∂D(t, r)

∂t

(cid:3) ∇ × H(t, r) − J(t, r)

(3.2b)

and boundary conditions. These can be either Dirichlet or impedance boundary conditions
on ∂ΩD or ∂ΩI, to bound the domain,

ˆn × E(r, t) (cid:3) ΨD(r, t) on ∂ΩD,
− Y ˆn × ˆn × E(r, t) (cid:3) ΨI(r, t) on ∂ΩI.

µ

ˆn × B(r, t)
∇ ×(cid:18) 1

(cid:19)

Instead of using (3.2), the wave equation
∇ × E

µr

+


r

∂2E
∂t2

1
c2
0

(cid:3) −µ0

∂J
∂t

(3.3a)

(3.3b)

(3.4)

can be used instead. The magnetic field can be obtained from (3.2a) and the impedance
boundary condition is defined using a time derivative on (3.3b) and using (3.2a). The fields
should also satisfy Gauss’ laws

∇ · D(t, r) (cid:3) ρ(t, r)

∇ · B(t, r) (cid:3) 0

(3.5)

(3.6)

though they are not explicitly solved.

As alluded to earlier, we use the moments of PSDF to find the fields generated and
then evolve their position using Newton’s equations and Lorentz force, viz., F(t, r) (cid:3)
q(t, r)(E(t, r) + v(t, r) × B(t, r)), and so on, for the duration of the simulation. Thus far, our
description has been in continuous world. To perform an actual simulation, we would need
to represent all the quantities involved in terms of functions defined on a discretization of
space and time. This is typically referred to as a particle in cell (PIC) approach and is the
subject of our next discussion.

Our starting point is the representation of both Ω and ∂Ω in terms of a finite set of
tetrahedra or a mesh that contains Ns nodes, Ne edges and N f faces. On these tetrahedra,
we define basis functions that follow the de-Rham sequence, enabling us to represent fields,
fluxes and sources []. But before proceeding too far ahead, note that we are going to

35

follow the usual PIC cycle; (a) map charges and currents on the mesh, (b) solve for electric
and magnetic fields on the mesh, (c) move particles due to Lorentz force and find the
current due to this motion, and (d) find the fields due the updated sources. The cycle then
continues.

With no loss of generality, we follow the usual procedure such that ρ(t, r) (cid:3) qα
Np
and J(t, r) (cid:3) qα
Np

The starting point of the simulation is to define the charge and currents due to PSDF.
δ(r−rp)
p(cid:3)1 vp(t)δ(r − rp). This implies that PSDF is sampled with Np shape
functions, each being a delta function. Generalization to other shape functions is possible
[36] and is agnostic to the crux of this paper.

p(cid:3)1

Specifically, the electric fields using Whitney edge basis functions, E(t, r) (cid:3)
Ne
The electric and magnetic fields are represented using Whitney basis functions[32, 14, 6].
(r).

N f
The magnetic flux density is represented using Whitney face basis function, B(t, r) (cid:3)
(2)
i(cid:3)1 bi(t)W
(r). Here, Ne are the number of edges and N f are the number of faces in
i
the mesh. Two different approaches can be used to solve Maxwell’s equations; (a) either
solve them in the coupled form or (b) solve the wave equation for the electric field and
then obtain the magnetic field. To set the stage for both these solvers, we introduce the
following Hodge matrix operators

(1)
i(cid:3)1 ei(t)W
i

(1)
[(cid:63)
]i,j (cid:3) (cid:104)W
i
(2)
[(cid:63)µ−1]i,j (cid:3) (cid:104)W
i

(1)
(r)(cid:105)
(r), ε · W
j
(2)
(r), µ−1 · W
(r)(cid:105),
j

the surface impedance matrix

(1)
[(cid:63)I]i,j (cid:3) (cid:104)ˆni × W
i

(1)
(r), µ−1 · ˆn j × W
j

(r)(cid:105)

and discrete curl operator

(1)
[∇×]i,j (cid:3) (cid:104)ni, ∇ × W
j

(r)(cid:105).

(3.7)

(3.8)

(3.9)

(3.10)

36

These matrices are used to build the semidiscrete Maxwell system

[I]

(cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125)

0
[(cid:63)
0
¯¯CM

0

]

∂t ¯B

∂t ¯E

 +



(cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125)

[∇×]T[(cid:63)µ−1]

[∇×]
0

0

 ¯B

¯E

 (cid:3)

¯¯KM

 0
(cid:124)(cid:123)(cid:122)(cid:125)

− ¯J

¯¯FM

(3.11)

where the degree of freedom vectors ¯E (cid:3) [e1(t), e2(t), . . . , eNe
and ¯J (cid:3) [j1(t), j2(t), ...jNe
becomes

(1)
(t)] with ji(t) (cid:3) (cid:104)W
i

(t)],
(t)], ¯B (cid:3) [b1(t), b2(t), . . . , bN f
(r), J(t, r)(cid:105). For the wave equation, the system

¯E (cid:3) −∂t ¯J

(3.12)

[(cid:63)
0
¯¯MW
where [(cid:63)S] (cid:3) [∇×]T[(cid:63)µ−1][∇×].

](cid:124)(cid:123)(cid:122)(cid:125)

t ¯E + c[(cid:63)I](cid:124)(cid:123)(cid:122)(cid:125)

∂2

¯¯CW

(cid:124)(cid:123)(cid:122)(cid:125)

∂t ¯E + c2[(cid:63)S]
¯¯KW

3.3.2 Unconditionally Stable Time Marching

The mixed finite element system in (3.2) is discretized in time using Newmark-Beta,
an unconditionally stable time stepping method. This method has been extensively used
in for the wave equation [14] and examined for the mixed finite element method in [5],
allowing for much larger time step sizes than the traditional leapfrog method. In this
method, the fields in time are represented by three temporal basis functions

corresponding to i ∈ [1, 3] and weighting function

Nn+1−i(t) (cid:3)

t − tn+1−j
tn+1−i − tn+1−j

W(t) (cid:3)

tn−t
∆t
t−tn
∆t

0

t ∈ [tn−1, tn]
t ∈ [tn, tn+1]
otherwise

.

2


j(cid:3)0
j(cid:44)i

37



(3.13)

(3.14)

This combination of basis function and weighting function creates a non-disappative,
unconditionally stable time marching scheme, which can be written as recurrence formula

provided in [38], corresponding to parameters γ (cid:3) 0.5 and β (cid:3) 0.25. When applied to (3.2),
this becomes

(0.5 ¯¯CM + 0.25∆t
+(0.5 ¯¯CM + .025∆t

¯¯KM) ¯Xn+1 − 0.5∆t
¯¯KM ¯Xn
¯¯KM) ¯Xn−1 + 0.25∆t ¯Fn+1
M + 0.25∆t ¯Fn−1
+0.5∆t ¯Fn

M
M (cid:3) 0

where ¯Xi (cid:3) [ ¯Bi,T , ¯Ei,T] and ¯Fi

M

(cid:3) [0,−
−1 ¯Ji,T]. Likewise, (4.11) becomes
( ¯¯MW + 0.5∆t

t

¯¯CW + 0.25∆2
+(−2 ¯¯MW − 0.5∆2
¯¯CW + .025∆2
t
t ¯Fn
W + 0.25∆2

¯¯KW) ¯En+1
¯¯KW) ¯En
¯¯KW) ¯En−1
t ¯Fn−1
W (cid:3) 0

t

( ¯¯MW + 0.5∆t
t ¯Fn+1
W + 0.5∆2

+0.25∆2

(3.15)

(3.16)

However, treating the current as written in (3.2) will not preserve the necessary link
between Ampere’s law and Gauss’ law needed to create a charge conserving scheme. This
is immediately apparent after applying a discrete divergence operator to semidiscrete
Ampere’s law. After using the identity that [∇·][∇×]T (cid:3) 0, this operation yields

] ¯En+1 − 0.5[∇·][(cid:63)
0

0.5[∇·][(cid:63)
0
−0.25∆t[∇·] ¯Jn+1 − 0.5∆t[∇·] ¯Jn − 0.25∆t[∇·] ¯Jn−1.

] ¯En−1 (cid:3)

When the same operator is applied to the semidiscrete wave equation, one gets

[∇·][(cid:63)
0

] ¯En+1 − 2[∇·][(cid:63)
0

] ¯En−1
−0.25∆t[∇·] ¯Jn+1 − 0.5∆t[∇·] ¯Jn − 0.25∆t[∇·] ¯Jn−1.

] ¯En + [∇·][(cid:63)
0

(3.17)

(3.18)

] ¯Ei, then it becomes clear that neither (3.17) nor
Making the substitution of ¯ρi (cid:3) [∇·][(cid:63)
0
(3.18) satisfy Gauss’ law or the continuity equation. Instead, a different treatment of the
right hand side, the particle current density, is needed in order to create a charge conserving
scheme.

38

3.4 Modified TDFEM-PIC

3.4.1

Integrator Agnostic Charge Conserving Scheme

It is apparent that, as written, time conservation fails for both Maxwell solver and the
wave equation. The reasons are two fold: (a) the order of time derivatives on the current
(on the right hand side) and those on the electric field are off by one; (b) this requires the
discrete time integrator to remember initial conditions. The latter holds the key to solving
the puzzle. Newmark time stepping schemes are, in effect, stable time integrators. The
crux of our approach is to correctly evaluate the time integral of the current. As elucidated
in [36], the time integral of the current is readily obtainable, and indeed a part of the PIC
scheme. Specifically, starting with the definition of the PDSF,
¯Jn(r(τ))dτ.

(3.19)

where tn (cid:3) n∆t. As shown in [36] this equation can be rewritten as,

0

− ¯ρn(r(tn)) (cid:3) [∇·]Þ tn
¯ρn(r(tn)) (cid:3) −[∇·]Þ r(tn)
Þ rp(tn)

dτvp(τ)δ(r − rp(τ)) (cid:3)

r(0)

¯Ji(˜r)d˜r.

Þ tn

0

Following the details in Ref. [36], it is immediately apparent that for any particle p

(3.20)

(3.21)

d˜rδ (r − ˜r)

rp(0)

Note, for each charge, its trajectory is determined by the solution to Newton’s equations.
The integration along a particle path can be computed to a high degree of accuracy. To
develop a charge conserving methodology, we define
¯Jn(τ)dτ

Þ tn

¯Gn (cid:3)

(3.22)

0

This equation is readily evaluated using (3.20). It follows that instead of using ¯J in (3.2)
(and therefore, in (3.17)) one can instead use ∂t ¯G. Discrete implementation with a Maxwell
equation solver results in in the divergence of Ampere’s law to be

0.5[∇·][(cid:63)
0

] ¯En+1 − 0.5[∇·][(cid:63)
0
−.5[∇·] ¯Gn+1 + 0.5[∇·] ¯Gn−1.

] ¯En−1 (cid:3)

(3.23)

39

Examining (3.24) term by term reveals that both sides of the equation are identical given
that [∇·] ¯Gn (cid:3) − ¯ρn. In a similar manner, one can use ∂2
¯G instead of ∂t ¯J in (4.11) to yield

t

¯ρn+1 − 2 ¯ρn + ¯ρn−1 (cid:3) −[∇·] ¯Gn+1 − 2[∇·] ¯Gn−1 + [∇·] ¯Gn−1.

(3.24)
Here, we have taken the liberty of substituting, ¯ρn (cid:3) [∇·][(cid:63)
0
] ¯En. At this point, we note
that the proposed approach is agnostic to the time stepping scheme (or integrator) used to
solve Maxwell’s equations; both the equation of continuity and Gauss’ laws are satisfied by
design.

A word of caution is in order before we proceed. While, the method developed is exact,
it should be noted that to obtain ¯En+1, one needs to solve either (3.17) (or (3.18)) with the
appropriate substitutions for ¯Gn+1 instead of ¯Jn+1. Obviously, the solution to these sets
of equations is subject to errors that arise due to vagaries of a linear algebraic solution
(tolerances, excitation of null-spaces, etc). As a result, as will be seen in the results section,
our errors are small but not identically zero. Next, we discuss a higher order particle
pusher to solve the equations of motion consistently.

3.4.2 Particle Pusher

Using an implicit time stepping scheme has advantages as well as challenges. The
principal advantage is an unconditionally stable time step size independent scheme as
opposed to a conditionally stable scheme like leap frog whose stability depends on the
time step size. The downside of using large time steps is that one must capture the nuances
of both the path and velocity of the particle. Thus, solving the equations of motion using
a Boris push [39], with its linear representation of the particle position and velocity can
introduce large errors into the system. Our goal is to develop a higher order scheme.

As is well known, the particle positions and velocities are updated by solving Newton’s
equations via the Lorentz force, giving us the following set of coupled first order ODEs for

40

each particle,

∂tvp(t, rp) (cid:3) ap(t, rp) (cid:3)
∂trp(t, rp) (cid:3) vp(t, rp).

qα
mα

(cid:0)E(t, rp) + vp(t, rp) × B(t, rp)(cid:1)

(3.25)

(3.26)

These form a pair of first order ODEs in time, and there are a number of methods that
can be applied. Our choice is to use a higher order order Adams-Bashforth scheme. An
exemplar recursion relation for v and r for a 4th order Adam’s-Bashforth method is as
follows:

p +

vn+1
p (cid:3) vn
rn+1
p (cid:3) rn

p +

(55an
∆t
24
(55vn
∆t
24

p − 59an−1
p − 59vn−1

p + 37an−2
p + 37vn−2

p − 9an−3
)
p − 9vn−3
).

p

p

(3.27)

(3.28)

where ∆t is the time step size. Given that the Newmark scheme is second order, we choose
the Adams-Bashforth scheme to be at least two orders higher so as to accommodate a
second time derivative in on ¯Gn. The path used for interpolating the position is a fourth
order Lagrange polynomial k (cid:3) 4 + 1 is defined as,
(cid:96) j(t)
t − tn+1−m
tn+1−j − tn+1−m

rp(t) (cid:3)
(cid:96)(t) (cid:3)

k




(3.29b)

(3.29a)

rn−j

j(cid:3)0

p

,

0≤m≤k
m(cid:44)j
where rp(t) is the position at time t and rn

p is the location of particle p at the t (cid:3) n∆t.

3.4.3 Particle Path and Current Mapping

The final piece of the puzzle is mapping the path to the underlying tesselation. In order
to do so, we note that that the integrator used to solve the equation of motion implicitly
assumes a Lagrange polynomial interpolant. As a result, the order of the method used
maps to order of the interpolant. This information needs to be used to find out where the
particle enters and leaves the cell.

41

Once the particle locations at each time step are known from the particle push, the path
through the unstructured mesh needs to be found. This includes finding the locations of
where a particle enters a cell and where it leaves and is detailed in Algorithm 3.1. Since
we are using a higher order representation of a particle path, finding these entry and
exit points of the cell with the tetrahedron becomes a non-linear problem and is detail
in Algorithm 3.2. Assume that we are given the normal to surface ˆn, and vertices of the
triangle rv,1, rv,2, rv,3. The intersection between the trajectory rp(t) and the plane can be
obtained by solving

ˆn ·(cid:2)(cid:0)rp(t) − rv,1

(cid:1) × (rv,2 − rv,1)(cid:3) (cid:3) 0

(3.30)

if All quadrature points between are in same cell then

Integrate using a quadrature rule along path.
Return and go onto next particle

if All quadrature points are in same cell then

Integrate using a quadrature rule along path.
Path is complete
Return and go onto next particle.

end if

Algorithm 3.1: Particle Path Finding Algorithm
1: Push particle finding rp, f
2: if rp, f is in same cell as rp,s then
3:
4:
5:
6:
7: end if
8: Find exit point rp,i of path in cell
9: Integrate from rp,s to rp,i
10: while Path not complete do
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
end if
21:
22: end while

Find next cell that path travels through
if r f is in same cell as rp,i then

Find exit point of path in cell.
Integrate from rp,is to rp,i f .

end if

else

42

t f (cid:3) tq

Algorithm 3.2: Non-Linear Bi-Section Method
1: ts (cid:3) 0,t f (cid:3) 1,th (cid:3) 0.5
2: if Any quadrature points are outside of the cell. then
3:
4: end if
5: while |ts − t f | < tol do
6:
7:
8:
9:
10:
11:
12: end while

if th is in same cell as ts then

end if
th (cid:3) 0.5(ts + t f)

else

ts (cid:3) th

t f (cid:3) th

Note, the path rp(t) can be parameterized using (3.29). Using this parameterization, one
can use a non-linear iteration (such as Newton-Raphson) to solve (3.30). For convinience,
we take a simpler approach by implementing a bi-section method that moves along the
path checking whether candidate points are inside or outside of the cell. For test cases
presented in this paper, this method converges rather robustly. Every step takes around
47 steps to converge below a tolerance of 1 · 10−15 (0.547 (cid:3) 7.1 · 10−15). Once the method
converges, we then compute the integral along each path segment in each cell using a set
of quadrature points. To illustrate this, consider Fig. 3.1 containing an example particle
starting position rp,s, and finishing position rp, f with the intersection point being rp,i. The
quadrature points would lie along the path between the rp,s and rp,i, then another set of
quadrature points between rp,i and rp, f .

Before we discuss results obtained using the above approach, a few points are in
order; to evaluate ¯Gn, (a) the integral over the path be evaluated using quadrature rules
to very high precision as the order of the path is know; (b) when the path passes though
multiple cells, the integration is broken up into pieces over each cell; (c) one can save on
computational cost of by updating the integral.

43

Figure 3.1: Particle path for a single particles with start and location and intersection point.

3.5 Results

In this Section, we present a number of results demonstrating the efficacy of the
proposed scheme with respect to conservation laws, as well accuracy of key steps that are
integral to the process.

3.5.1 Higher Order Particle Motion

One of the key advantages in using implicit time stepping is the possibility of using
much larger time step sizes. Unfortunately, this also implies that one needs higher order
methods to capture both the path as well as velocity. In this section, we demonstrate
convergence of our algorithm for particle motion using various orders of Adams-Bashforth
integrator and compare these to standard non-relativistic Boris push.

To do so, we set up a classic cyclotron [40] motion test where a single particle was given
an initial velocity in a constant magnetic field resulting in circular motion as shown in
Fig. 3.3. The parameters are shown in Table 3.1 with a particle’s initial velocity v0 with a
background magnetic fields B with a given mass m and charge q. The particle will move in
a circle due to the Lorentz force as shown in Fig. 3.2. The relative error in both position
and velocity for various time step sizes with multiple order of Adams-Bashforth and Boris
are is shown in Fig. 3.3. The average error is calculated by taking the norm of the distance
errors of each point r and dividing by the normal of the analytic positions ra (see Ref. [40]

44

Figure 3.2: Mean relative error in position for Adam-Bashforth Orders 1-5 compared with
the Boris push, shown in black.

for details),

error (cid:3)

||r − ra||2
||ra||2

(3.31)

The slopes for each of the Adams-Bashforth methods match its order. Boris on the other
hand has a second order velocity update with a first order positional update. This test
essentially validates out pusher as well as helps correlate error (or approximately so) in
particle motion with time step size.

Table 3.1: Cyclotron Motion

Parameter

B
Q
m
v0
r0

Value

6.822756 · 10−5ˆz T
−1.60217646 · 10−19 C
9.10938370 · 10−31 kg
3 · 106 ˆy m/s
[0.75, 0.5, 0.0] m

45

Figure 3.3: Mean relative error in velocity for Adam-Bashforth Orders 1-5 compared with
the Boris push, shown in black.

3.5.2 Expanding Particle Beam

Next, we consider an expanding beam test [40]. An expanding particle beam is injected
into a cylindrical cavity with an initial velocity of magnitude v0. As the beam travels down
the tube, the electrons repel each other causing the beam beam to expand. This expansion
rate can be compared with other codes to validate the solution. The detail of the mesh and
beam parameters used are shown in Table 4.1.

Both the wave equation and mixed finite element trajectories are compared in Fig.
3.4 and show good agreement with XOOPIC [4] (an extensively used and well validated
quasi-2D FDTD code). We sample the electric field half way down the tube 16 mm from
the center of the tube. The radial field values are plotted over time shown in Fig. (3.5) for
simulations with different time steps. We compare four runs with time steps of α∆t where
α is scale factor and ∆t (cid:3) 0.333 ps is the largest stable step size in a leap frog time marching
method for the given mesh. Note, 2 ns corresponds to 1 transit of the tube. It is evident
from this figure that the proposed method provides stable results; indeed, as is evident
from this figure, the data at 7.5∆t, 15∆t and 30∆t are almost identical to each other, where

46

Table 3.2: Expanding Particle Beam Parameters

Parameter

Cavity Radius
Cavity Length

Boundary Conditions

Value
20 mm
100 mm
5 · 107 m/s
0.16678
8.00 mm

PEC

electrons

2 ns
0.25 A
52012.58
1.529 mm
6.872 mm

ns

Number particles per time step

10

v0
v0/c

beam radius rb

species

Turn on time
beam current

macro-particle size
min edge length
max edge length

∆t

as the one at 1485∆t is slightly different. This points to significant gains that can be made
with Newmark time stepping (provided the method is charge conserving).

This leads to the next argument. Shown in Fig. 3.6 is data from two different methods
for the same set up run using MFEM with backward difference at ∆t, MFEM with Newmark
at 7∆t and the wave equation (WE) at 7∆t. As evident, all three methods conserve charge
to almost machine precision. It should be noted that both MFEM and WE have a null
space. In the case of the former, it is fields that behave like ∇φ(r), and the latter, as t∇φ(r).
However, as is evident from these results, our mapping on to these null spaces is small and
behaves as expected.

To further illustrate the robustness of the method to time step sizes, in Fig. 3.7 we
compare the satisfaction of Gauss’ law for all four time steps used in Fig. 3.5. As is evident
from here, charge is again conserved almost to machine precision (around 10−18 for all
with slight difference evolution of trajectory).

47

Figure 3.4: Expanding particle beam macro particles in the z vs r plan. Particle locations
from both mixed finite element methods and wave equation versions are compared with
XOOPIC beam profile.

Figure 3.5: Electric field values are the radial component half way down the tube 16 mm
from the center of the tube. Multiple simulation with different time steps are performed.

48

Figure 3.6: Discrete Gauss’s Law error per particle for Newmark-Beta mixed finite element
(NM-MFEM), Newmark-Beta wave equation (WE), backwards difference mixed finite
elements (BD-MFEM) using the charge conservation technique provided here.

Figure 3.7: Discrete Gauss’s law error per particle various time steps using the mixed finite
element methods using Newmark time stepping.

49

Table 3.3: Adiabatic Expanding Plasmas

Parameter

Mesh Radius

Boundary Conditions

First order ABC

Value
6mm

1K
100K
8000

Electrons and Sr +

52012.58
1.529mm
6.872mm

Tion

Telectron

Species

Number Particles

Macro-Particle Size
Min Edge Length
Max Edge Length

3.5.3 Adiabatic Expanding Plasma

Finally, for a third validation case we simulate an adiabatic expansion of a plasma ball
with radial Gaussian distribution in the radial direction. This case has an analytic solutions
[24] and allows for good comparison and validation. We change some of the parameters
from the original numerical experiments [40] such that the Debye length can be fully
resolved. This example is described in more detail in[40]. We simulate this example both
MFEM and WE. For both examples we get excellent agreement in the expansion rate with
both the wave equation, Fig. (3.9), and the mixed formulation, Fig. 3.8, when compared
with analytic densities.

3.6 Summary

In this paper, we have presented a solution to a problem that has been long-standing–
charge conserving FEM-PIC methods for implicit time stepping systems without the need
to adopt divergence cleaning. In other words, rubrics have been developed such that
conservation laws are implicitly obeyed. Indeed, the method presented is agnostic to any
time stepping scheme. We have demonstrated the efficacy of this approach for a set of test
problems, using different time step sizes and different time stepping schemes, as well as
both MFEM and WE solvers. The results reliably attest our claims. The above approach
opens multiple doors that will further the state of art of FEM-PIC; these include higher

50

Figure 3.8: Mixed Finite Element particle beam expansion

Figure 3.9: Wave Equation Adiabatic Expanding Plasma

51

order schemes in both space and time, quasi-Helmholtz decomposition to get a better
handle on null-spaces, and domain decomposition to effect rapid solution by parallelizing
the scheme. Papers on these will be presented soon in other forums

52

CHAPTER 4

QUASI-HELMHOLTZ DECOMPOSITION FOR FINITE ELEMENT

PARTICLE-IN-CELL

4.1 Abstract

Development of particle in cell methods using finite element based methods (FEMs)
have been a topic of renewed interest; this has largely been driven by (a) the ability of finite
element methods to better model geometry, (b) better understanding of function spaces that
are necessary to represent all Maxwell quantities, and (c) more recently, the fundamental
rubrics that should be satisfies in space and time so as to satisfy all conservation laws. In
that vein, methods have been developed that satisfy conservation laws and are agnostic to
time stepping methods. While these methods are indeed a significant advance, it should be
noted that implicit FEM transient solvers support an underlying null space that corresponds
to a gradient of a scalar potential ∇Φ(r) (or t∇Φ(r) in the case of wave equation solvers).
While explicit schemes do not, they are only conditionally stable, time step sizes are mesh
dependent, and very small. A way to overcome this bottleneck, and indeed, better satisfy
all four Maxwell’s equation is to use a quasi-Helmholtz formulation on a tesselation. In the
re-formulation presented, we strictly satisfy both the equation of continuity and Gauss’
laws for both the electric and magnetic flux densities. Results demonstrating the efficacy of
this scheme will be presented.

4.2 Introduction

Simulation of plasma and space charge has many applications in science and technology
ranging from particle accelerators to satellites and medicine[26, 27, 28]. A popular method
for simulation of plasma is particle-in-cell (PIC) which self-consistently evolves particle
motion for charge species with Maxwell’s Equations. Traditionally PIC has been based

53

on finite difference time domain with the domain represented by Yee cells[9]. The stair
stepping nature of a structured grid presents challenges on geometry modeling especially
fine features and multi-scale objects. Methods such as cut-cells[41] have been developed
to alleviate this challenge. Another fundamental bottleneck is that these methods are
conditionally stable, with time step size of the field solver dictated by the mesh, and is
relatively small so as to satisfy the Courant-Friedrich-Levy condition. As a result, there has
been renewed interest in developing methods that not only better represent the geometry
but better capture the underlying physics.

The route taken is to use finite element methods (FEMS); they have an extensive
history of use in field solvers due to their ability to better represent the geometry and
physics. But more importantly, there is a well developed rigorous body of literature
[32, 42, 43] on the mathematical foundation of FEM. Indeed, as is to be expected, FEM
has made inroads into PIC modeling more of late with some excellent work rubrics of
this approach[44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60]. Of particular
note in seminal work in [6], wherein rubrics for charge conservation for FEM-PIC where
developed. These have been refined and developed further in a series [15, 61, 62]. By
and large, these papers follow a similar rubric, use a leap-frog scheme to step through a
solution for the field, update position and velocity of the particle and solve for the fields
again. In this setting, it was shown that the methodology is charge conserving. While this
the classical PIC cycle, a more recent trend has been to solve both the electromagnetic and
Newton’s equations self-consistently at the same time step[53]. This approach necessary
when the motion is close to being relativistic.

In the problems that we address in this paper, our motion is non-relativistic. As a
result, we are in the regime where the classic PIC is valid. Even so, being restricted to a
time step size that is tightly related to the mesh is undesirable, especially, when there are
a number of well established methods that make the field solver unconditionally stable
[38]. Unfortunately, developing a method where Gauss’ laws and equation of continuity

54

is satisfied is not a trivial task. As an aside, in the rest of the paper, we will use the term
conservation laws to denote these equations and admit to an abuse of terminology! This
problem was solved relatively recently [63, 36] with results demonstrating satisfaction of
conservation laws to almost machine precision. But this is not really the end of the story.
While leap frog FEM is known to not have a null space, implicit time stepping methods do.
There are two flavors of time domin FEM-those that discretize Maxwell’s equations directly
(MFEM) and those that solve the wave equation (WE-FEM). These equation suffer from a
null space that behave as a gradient of a potential and a product time with a gradient of a
potential, respectively. As a result, it is possible to add spurious charge into the system
corresponding to the divergence of these null-spaces. This is especially true when one
uses iterative methods to solve the system of equations where level of excitation of the
null space depends on the error threshold. Given that satisfaction conservation laws is
paramount in a PIC solution, we prescribe an alternate solution.

Traditional solution relies on solving for fields via the two curl equations. Gauss’ laws
will be satisfied by these fields provided equations of continuity hold and there are no null
spaces excited. But this is easier said than done. It follows that perhaps an alternative is to
use a Helmholtz decomposition of fields in the construction of field solvers. This would
permit solution to all four equations. Means to effect this is the overall thrust of this paper.
It is well known that an exact Helmholtz decomposition is not possible in a discrete
setting. Indeed, using an analytical decomposition as a starting point prior to discretization
will lead to erroneous results as one needs to ensure that proper mapping on de-Rham spaces
are satisfied. What one can achieve is a quasi-Helmholtz decomposition where one of the
two conditions are strongly satisfied. The development of quasi-Helmholtz decomposition
dates back to the early days of development vector basis for modeling electromagnetic
fields [64, 43, 65, 66]. Seminal literature was devoted to understanding the necessary
mathematics in the language of differential forms and was driven by the need to solve
magnetostatic problems. Development of quasi-Helmholtz decomposition can be traced

55

back to tree-cotree decomposition for magnetostatics [67, 68, 69, 70, 71]. Development
of of this body of work has progressed from simply connected to multiply-connected
objects [72, 73]. As an aside, this period was rich with development of state of the art FEM
techniques, including MFEM leap frog methods [74].

This body of work is the starting point of the work presented in this paper. We will
restrict ourselves to simply connected objects and zero-th order basis function; extension
to multiply connected objects and higher order basis is underway and will be presented
elsewhere. In this paper, our principal contributions of this paper are as follows:

1. We will develop quasi-Helmholtz projectors to partition fields into solenoidal and

non-solenoidal components for both the electric field and magnetic flux density.

2. In specializing these to simply connected systems, we will use a co-tree identification
of degrees of freedom such that the total number of degrees of freedom is identical
to the original system.

3. We will show that the resulting system of equations map exactly to Gauss’ law and

the equation of continuity is exactly satisfied.

4. We apply these to both MFEM and WE-FEM, discuss ramification of null spaces and

show connections/variations from the well known A − Φ formulation.

5. Finally, we will present numerous results that validate the arguments presented in

the paper.

The rest of this paper is organized as follows: Next, in Section 4.3 we present the
desired formulation of the PIC system in the continuous domain. Section 4.4 presents
the qausi-Helmholtz decomposition in a discrete setting for both the electric field and
magnetic flux density, and derivation of the revised sets of equation. We also show how
these equations satisfy conservation laws. Next, in Section 4.6, we discuss null spaces that

56

arise out of these equations as well as gauge considerations. Finally, results are presented
in in 4.7 that demonstrates the advantages of the proposed method.

4.3 Formulation

4.3.1 Preliminaries

Consider a domain Ω whose boundaries are denoted by ∂Ω and enclose the domain.
The domain consists of free space and permittivity ε0 and permeability µ0 with the speed
of light denoted using c (cid:3) 1/√
µ0
0. The enclosing boundaries are either Nuemann or
Dirichlet. The excitation of the system comes from electromagnetic field due to moving
charges whose position (and velocity) evolve over. The moving charges are assumed to be
collisionless and follows the Vlasov equation,

∂t f(t, r, v) + v · ∇ f(t, r, v)+
[E(t, r) + v × B(t, r)] · ∇v f(t, r, v) (cid:3) 0,

q
m

definition of the charge and current density defined as ρ(t, r) (cid:3) qÞ
J(t, r) (cid:3) qÞ

where f(t, r, v) is a phase space distribution function (PSDF). We follow the conventional
Ω f(t, r, v)dv and
Ω v(t) f(t, r, v)dv as moments of the PSDF [2]. The fields, E(t, r) and B(t, r), in

Vaslov equation are solutions to the Maxwell’s equations Eq. 4.2.

(4.1)

(4.2a)

(4.2b)

(4.2c)

∇ · D(t, r) (cid:3) ρ(t, r)
∇ · B(t, r) (cid:3) 0
∇ × E(t, r) (cid:3) −∂tB(t, r)
∇ × H(t, r) (cid:3) ∂tD(t, r) + J(t, r)

(4.2d)
where D(t, r) (cid:3) εE(t, r) and H(t, r) (cid:3) µ−1B(t, r). The boundary conditions enclosing the
domain are either Neumann ∂ΩN or Dirichlet ∂ΩD,

ˆn × E(t, r) (cid:3) ΨD(t, r) on ∂ΩD,
ˆn × µ−1B(t, r) (cid:3) ΨN(t, r) on ∂ΩN ,

(4.3a)

(4.3b)

57

where ˆn is an outward pointing normal to ∂ΩN and ∂ΩD and ΨN(t, r) and ΨD(t, r) are
Neumann and Dirichlet boundary condition functions. In addition, it can be shown from
Maxwell’s equations that the equation of continuity

∇ · J(t, r) + ∂t ρ(t, r) (cid:3) 0.

particles. We follow the usual representation where ρ(t, r) (cid:3) q
Np
J(t, r) (cid:3) q
Np

(4.4)
holds. While f(t, r, v) is not solved directly, as it is computationally prohibitive, we
use an approach where one approximates the PSDF as a collection of discrete point
δ(r − rp(t)) and
p(cid:3)1 vp(t)δ(r − rp(t)). Generalization to other shape function (instead of delta
functions) is certainly possible[36], but is beyond the scope of this paper and not pertinent
to the central thesis of this paper. As in PIC schemes, the particles are moved using the
Lorentz force F(t, r) (cid:3) q(E(t, r) + v(t, r)× B(t, r)) and Newton’s equations. A self-consistent
solution to both evolution of charges and fields constitute a PIC methodology.

p(cid:3)1

4.3.2 Helmholtz Decomposition

It is well established via Helmholtz theorem that any sufficiently smooth vector field
Λ(t, r) (cid:3) Λir(t, r) + Λr(t, r) can be represented into a rotational (divergence-free) Λr(t, r)
and irrotational (curl-free) Λir(t, r) components. This implies that ∇ × Λir(t, r) (cid:3) 0 and
∇ · Λr(t, r) (cid:3) 0. Along this vein, it follows that both the electric and magnetic fields can be
decomposed in a similar manner. Parenthetically, we note here that in the PIC context, our
point charges are measured on the mesh, thereby smoothing the sources and the resulting
vector field. Given this decomposition we can rewrite Maxwell’s equations as,

∇ · Dir(t, r) (cid:3) ρ(t, r)
∇ · Br(t, r) (cid:3) 0
∇ × Er(t, r) (cid:3) −∂tBr(t, r)
∇ × Hr(t, r) (cid:3) ∂tDr(t, r) + ∂tDir(t, r) + J(t, r)

(4.5a)

(4.5b)

(4.5c)

(4.5d)

58

Additionally the continuity equation is updated as, ∇ · Jir(t, r) + ∂t ρ(t, r) (cid:3) 0 The main
takeaways are the following; (a) the three components of fields that one needs to determine
are Eir(t, r), Er(t, r) and Br(t, r), (b) Maxwell’s equations are sufficient to determining
these components, and (c) solving these would ensure that null-spaces do not corrupt
conservation laws. Of course, it is well known that it is not possible to construct a complete
Helmholz decomposition in a discrete setting. As result, in what follows, we will develop
a quasi-Helmholtz decomposition on discretization of the domain. For simplicity of
exposition, we assume that the domain is simply connected [43] .

4.3.3 Discrete Framework

E(t, r) (cid:3) span(cid:8)W1(r)(cid:9) and B(t, r) (cid:3) span(cid:8)W2(r)(cid:9). Specifically, E(t, r) (cid:3)
Ne
and B(t, r) (cid:3)
N f

The analysis framework starts with the de-Rham complex[43] to represent fields and
fluxes. To that end, we assume that the domain is represented using a collection of finite
elementsK (cid:3) {N , E, F , T } defined using Nn nodes, Ne edges, N f faces and Nt tetrahedron.
On these elements one can define Whitney spaces [43, 32] W0(r), W1(r), W2(r) and W3(r) to
represent quantities on the primal grid. In what follows, we represent electric fields using
(r)
(r). It follows that these coefficients can be arranged into vector
of the form ¯E(t) (cid:3) [e1(t), e2(t), . . . , eNe
(t)]; the value at any
time n∆t is denoted using ¯E(n∆t) (cid:3) ¯En (likewise for ¯B(n∆t) (cid:3) ¯Bn) where ∆t is the time
step size. Similar notation is adopted for the coefficients of both the magnetic field and
electric flux density.

(1)
i(cid:3)1 ei(t)W
i

(2)
i(cid:3)1 bi(t)W
i

(t)], ¯B(t) (cid:3) [b1(t), b2(t), . . . , bN f

In keeping with this framework, it follows that the electric flux densities and magnetic
field lies in the dual space, but can be presented using primal space quantities. But before
we do so, we define the following Hodge matrices [(cid:63)
]i,j and [(cid:63)µ−1] that map primal grid

59

quanties to the dual grid:

(1)
[(cid:63)
]i,j (cid:3) (cid:104)W
i
(2)
[(cid:63)µ−1]i,j (cid:3) (cid:104)W
i
(3)
[(cid:63)ρ]i,j (cid:3) (cid:104)W
i

(1)
(r), ε · W
(r)(cid:105); i, j ∈ E
j
(2)
(r), µ−1 · W
(r)(cid:105); i, j ∈ F
j
(r), W

(r)(cid:105); i, j ∈ T

(3)
j

In addition, we define the following
(1)
[Mg]i,j (cid:3) (cid:104)W
i
(2)
[Mc]i,j (cid:3) (cid:104)W
i
(3)
[Md]i,j (cid:3) (cid:104)W
i

(0)
(r)(cid:105); i ∈ E, j ∈ N
(r), ∇W
j
(1)
(r)(cid:105); i ∈ F , j ∈ E
(r), ∇ × W
j
(2)
(r)(cid:105); i ∈ T , j ∈ F
(r), ∇ · W
j
[∇] (cid:3) ε[(cid:63)
]−1[Mg]
[∇×] (cid:3) µ−1[(cid:63)µ−1]−1[Mc]
[∇·] (cid:3) [(cid:63)ρ]−1[Md]

(4.9f)
where [∇×] is the discrete curl operator, [∇] is the discrete gradient operator, and [∇·] is
the discrete divergence operator. We note that using this notation, one can obtain the
coefficients for the electric flux density as ¯Dn (cid:3) [(cid:63)
] ¯En. Using this framework, one can
write the discretized equations in space as

[I]

0



∂t ¯B

 +



0

0
[(cid:63)ε]
where ¯J (cid:3) [j1(t), j2(t), ...jNe
one can arrive at a similar discrete system for the wave equation as

−[∇×]T[(cid:63)µ−1]
∂t ¯E
(1)
(t)] with ji(t) (cid:3) (cid:104)W
i



 ¯B

¯E

 (cid:3)

 0

− ¯J



[∇×]
0

(4.6)

(4.7)

(4.8)

(4.9a)

(4.9b)

(4.9c)

(4.9d)

(4.9e)

(4.10)

(r), J(t, r)(cid:105), and ¯Jn (cid:3) ¯J(n∆t). Likewise,

[(cid:63)ε]∂2

t ¯E + [(cid:63)S] ¯E (cid:3) −∂t ¯J

(4.11)
where [(cid:63)S] (cid:3) [∇×]T[(cid:63)µ−1][∇×]. Using these equations as a backdrop, we next discuss the
steps necessary to effect a discrete quasi-Helmholtz decomposition.

60

4.4 Quasi-Helmholtz Decomposition

As noted earlier, it is well known that one cannot develop a complete Helmholtz
decomposition of fields. It follows, that we may be able to achieve one condition, viz.,
either divergence or curl free fields. To this end, prior to developing a numerical framework,
we ask what is necessary for the fields to satisfy. In plasma physics, it follows that one
must satisfy the condition ∇ · Dr(t, r) (cid:3) 0 and Bns(t, r) (cid:3) 0 strongly; in other words, we
cannot necessarily ensure ∇ × Dir(t, r) (cid:3) 0. As eloquently pointed out by Vecci [75], the
condition that we are strongly imposing is the solenoidal nature of Dr(t, r); i.e., Dir(t, r) is
non-solenoidal but not irrorational, i.e., ∇×Dir(t, r) (cid:44) 0. Henceforth, we will use subscripts
“s” and “ns” to denote solenoidal and non-solenoidal components, respectively.

To discuss this decomposition, we start with the electric field. The usual starting
∇−→ W1. It is
point of the analysis is de-Rham diagrams. We focus on the sequence W0
apparent that in the framework being pursued we seek to exploit the representation of

the potential φ(t, r) (cid:3) span(cid:8)W0(r)(cid:9). As ∇W0 ∈ W1, and we seek a decomposition of the

form E(t, r) (cid:3) Ens(t, r) + Es(t, r) it follows that these statements can be made in the discrete
setting via

Ne


i(cid:3)1

(r) (cid:3) ∇ Nn


(r) + Es(n∆t , r)

i W1
j,nsW0
e n
e n
i
j
[(cid:63)ε] ¯En (cid:3) ε[Mg] ¯En

j(cid:3)1

ns + Λ ¯En
s

(4.12a)

(4.12b)

¯Dn (cid:3) Σ ¯En

ns + Λ ¯En
s

(4.12c)
We have chosen not to represent Es(n∆t , r) explicitly, leaving it as the remainder. It follows,
that one can choose basis sets (the so called co-tree basis), but that is not quite important
right now. It is sufficient to say that there exists a matrix Λ that is an inner product of these
basis with W1(r).

In a similar manner, we seek a representation of the magnetic flux density such that
B(t, r) (cid:3) Bns(t, r) + Bs(t, r). To do so we take recourse again to the de-Rham picture and note

61

that the magnetic field, H(t, r), lies the dual space(cid:101)W1(r), where tilde denotes quantities
associated with the dual space. As before, we take recourse to the sequence(cid:101)W0
∇−→(cid:101)W1
and note that(cid:101)∇(cid:101)W0(r) ∈(cid:101)W1. In spirit of (4.12), it follows that one can write

N f


i(cid:3)1

(r) (cid:3)(cid:101)∇ Nt


i(cid:101)W1

i

hn

j,ns(cid:101)W0

j

hn

j(cid:3)1

[ ˜(cid:63)µ] ¯Hn (cid:3) µ[∇·]T ¯Hn

(r) + Hs(n∆t , r)

(cid:48)
m ¯Hn
ns + µΛ
s
ns + Λm ¯Bn
s

¯Bn (cid:3) Σm ¯Bn

(4.13a)

(4.13b)

(4.13c)

where [ ˜(cid:63)µ] ¯Hn (cid:3) ¯Bn. In all both (4.12) and (4.13), we require ΣTΛ (cid:3) 0 and ΣT
mΛm (cid:3) 0. This
imposes the necessary solenoidal nature of the fields. Parenthetically, we note that Σ and
Σm can have a one dimensional null space that corresponds to the zero mean constraint on
the potentials; as will be evident in the next Section, this depends on boundary conditions.

4.5 Decomposition

As was evident from the discussion thus far, the decomposition that we espouse relies
on imposing solenoidal nature of a component. To do so, we have specified Σ and Σm.
Mapping to approaches used in the literature is apparent; Σ and Σm are associated with
tree’s formed either out of edges or faces respectively. The complement of these are the
co-trees (of edges and faces). While we shall return to this concept later, the approach we
take is more generalizable along the lines of Ref. [68]. This follows the observation that in a
Helmholtz decomposition in multiply connected domains, both the irrotational component
and the harmonic component are divergence free, and it is not important to distinguish
between the two. To this end, we introduce projector that effect this decomposition.

4.5.1 Projector Matrices

Development of projectors, given (4.12) and (4.13), is relatively straightforward. In
what follows, we will illustrate this for the electric field. As alluded to earlier, ΣTΛ (cid:3) 0

62

implies ΣTEs(n∆t , r) (cid:3) 0. It follows from (4.13c) that

ns (cid:3)(cid:0)ΣTΣ(cid:1)†

¯En

ΣT ¯Dn

(4.14)
where the † denotes a Morse-Penrose psuedo inverse. We note though that abstraction
hides insight; in this case, ΣT ¯Dn is related to the Gauss’ law for electric fields and can
be used to solve for ¯En
ns directly, a tactic we will flesh out in detail in next subsection.
Completing the train of thought, one can now define projectors

[P]Σ
[P]Λ

e (cid:3) Σ(ΣTΣ)†
ΣT
e (cid:3) I − [P]Σ

e

Using these projectors, it follows that (4.13c) can be rewritten as

¯Dn (cid:3) Σ ¯En

ns + [P]Λ

e

¯Dn

(4.15a)

(4.15b)

(4.16)

If, ¯En
ns has already been computed, what remains is to compute the rest of the contributions.
Similar approach is possible for the magnetic flux density. Here, the projectors take the
form

(4.17)

(4.18)

b (cid:3) I − Σm
and the magnetic flux density can be written as

[P]Λ

mΣm

ΣT
m

(cid:0)ΣT

(cid:1)†

¯Bn (cid:3) Σm ¯Bn

ns + [P]Λ

b ¯Bn

It follows that the magnetic Gauss’ law is trivially satisfied provided we choose ¯Bn
ns (cid:3) 0.
As an aside, the motivation to present ideas in terms of projectors is to lay the foundation
for extending these to both higher order and make the decomposition genus free. This
extension will be the thrust of another forthcoming paper. The main goal of this paper,
is to fully develop and demonstrate these ideas for 0th order Whitney basis and simply
connected structures.

63

4.5.2 System of Equations

Identification of both the non-solenoidal and the solenoidal component has been studied
for (non)simply connected structures for a long time []. In this paper, this forms our starting
point. Given any network, one can create a minimum spanning tree as shown in Fig.
4.1. From this minimum spanning tree, one can identify both the tree and co-tree edges;
specifically, we note that Σ maps N −→ E and Λ maps Ec −→ E where Ec ⊂ E. In a similar
manner, one can construct a tree (and co-tree) of faces; this is illustrated in two-dimensions
in Fig. 4.2. Here, Σm is a map from T −→ F , and Λm is a map from Fc −→ F where
Fc ⊂ F . Here, and henceforth, the subscript “c” will denote co-tree quantities. It follows
that one can define matrices

Likewise,

[Ce
c]i,j (cid:3)

[Cb
c]i,j (cid:3)

 1,
 1,

i ∈ E, j ∈ Ec
0, otherwise

i ∈ F , j ∈ Fc
0, otherwise

(4.19)

(4.20)

A word of caution; if Dirichlet boundary conditions are applied the tree must be grounded
to the PEC surface as in Fig. 4.1. We need one more piece of the puzzle prior to setting up
our solution system. As we had noted earlier, the both Σ and Σm have a 1-D null space
that arises as a result of boundary condition that need to be imposed; this can be related
to the appropriate Whitney spaces being closed as well. For ¯En
ns it manifest itself as zero
mean constraint when Neumann boundary conditions are used and for ¯Bn
c when Dirichlet
boundary conditions as used. To effect these, we introduce a zero mean constraint matrices
[Cq
z]

[Cq
z] (cid:3)

(4.21)

 .



1
0
− αq
1
αq
n

···

0
. . .
··· − αq
n−1
αq
n

64

αk

W k
i

i (cid:3)

(4.22)
for q (cid:3) {e, b} and k (cid:3) 0 when q (cid:3) e and k (cid:3) 3 when q (cid:3) b. Note, these matrices become
identity if the constraint is not necessary; likewise these need be modified when the net
charge is not zero. Using this framework, we now examine projectors and solution to
our field equations. We note the following: [∇]TΛ (cid:3) 0 and [∇·]Λm (cid:3) 0 and more to the
point, [∇]T ¯Dn and [∇·] ¯Bn are the discrete analogs of the continuous Gauss’ laws ([∇]T is a

discrete divergence [(cid:101)∇·] in the dual grid). This fact allows use to redefine projectors; viz,

with a slight abuse of notation

z](cid:0)[Ce
z](cid:16)[Cb

[P]Σ
[P]Σ

e (cid:3)Σ[Ce
(cid:3)Σm[Cb

b

z](cid:1)−1 [Ce
z](cid:17)−1 [Cb

z]T[∇]TΣ[Ce
z]T[∇·]Σm[Cb

z]T[∇]T
z]T[∇·]
e and [P]Λ

(4.23a)

(4.23b)

Using these in (4.15b) and (4.17) enables us to derive [P]Λ
b and then using (4.16)
and (4.18) we have a complete description of the electric and magnetic fields/flux densities.
We now have the necessary quantities to prescribe a PIC scheme to consistently evaluate
¯En
ns, ¯En
ns (cid:3) 0. Just to reiterate, there is a slight
abuse of notation. Quantities are assumed to evaluated at n∆t. We have not specified
temporal basis sets used for discretization or how functions are measure in time. We will
leave that to later in this section when it is more appropriate.

s given ¯Jn. We will, a-priori, set ¯Bn

s and ¯Bn

4.5.2.1 Solution for non solenoidal component En
ns

The starting point of the solution is the discrete Ampere’s law at time point n∆t which
z, operating both sides
e (cid:3) 0 (cid:3) [∇]T[∇×]T) results in

reads ∂t ¯Dn − [∇×]T[(cid:63)µ−1] ¯Bn (cid:3) − ¯Jn. Using (4.16). the constraint [C]e
by [∇]T, and using the identities prescribed earlier ([∇]T[P]Λ

Þ

Ω

(r)dr

[∇]TΣ[Ce
[∇]T[(cid:63)ε][∇][Ce

z]∂t ¯En
z] ¯En

Þ n∆t
ns (cid:3) −[∇]T ¯Jn
ns (cid:3) −[∇]T

0

65

¯J(τ)dτ (cid:3) ¯ρn

(4.24)

It is apparent that (4.24) is the discrete Gauss’ law for the electric fields. The left hand side
of this equation is the discrete Laplacian, and coefficients ¯En
ns are the appropriate values
of the potential due to the charge density on the right hand side at that instance of time.
Further note, as shown in [36] this equation can be rewritten as,

Þ n∆t

Þ r(n∆t)

ji(t)dτ (cid:3)

(1)
(cid:104)W
i

(˜r), J(t, ˜r)(cid:105)d˜r

¯Gn
i (cid:3)

0

r(0)

(4.25)
It follows that, ¯ρn (cid:3) −[∇]T ¯Gn. From the equation above, it follows that the integration over
time is subsumed in the path integral used to evaluate the charge density[36, 63]. As a
result, one obtains the solution to the Laplacian directly. Two points that should be noted:
(a) the equation of continuity is analytically satisfied, and (b) the electric flux density is
obtained by explicitly solving Gauss’ law. Having obtained ¯En
ns, obtaining the rest is fairly
straightforward for both MFEM and WE-FEM.

4.5.2.2 Solution to the solenoidal components

Having obtained the solution to the non-solenoidal components, the solution to the rest

is relatively straightforward. We begin with rewriting Maxwell’s equations as

[Z]11∂t ¯Bn
[Z]21∂t ¯En

s + [Z]12 ¯En
s − [Z]22 ¯Bn

s (cid:3) −[Z]13 ¯En
s (cid:3) −∂t ¯Gn − [Z]23∂t ¯En

ns

ns

where

b

[Z]11 (cid:3)[Cb
[Z]12 (cid:3)[Cb
[Z]13 (cid:3)[Cb
[Z]21 (cid:3)[Ce
[Z]22 (cid:3)[Ce
[Z]23 (cid:3)[Ce

c]
[Cb

c]T[P]Λ
c]T[∇×][(cid:63)ε]−1[P]Λ
e [(cid:63)ε][Ce
c]
z]
c]T[∇×][(cid:63)ε]−1Σ[Ce
c]T[P]Λ
c]T[∇×]T[(cid:63)µ−1][P]Λ
c]TΣ[Ce
z]

c]
e [(cid:63)ε][Ce

[Cb
c]

b

66

(4.26a)

(4.26b)

(4.27a)

(4.27b)

(4.27c)

(4.27d)

(4.27e)

(4.27f)

Figure 4.1: Tree (solid) and cotree (dashed) edges for a 2D geometry.

Figure 4.2: Tree (solid) and cotree (dashed) faces for a simplified 2D tree structure.

Since Σ (cid:3) [(cid:63)ε][∇], it is trivial to show that [Z]13 (cid:3) 0. The above is a re-writing of Maxwell’s
equation on the tetrahedral grid. Using the above, the wave equation can be trivially
derived using (4.26a) in the time derivative of (4.26b). Specifically, using the above notation,
these can be written as

[Z]21∂2
s + [Z]22[Z]−1
t ¯En
− [Z]23∂t ¯En

11

ns

¯Gn

(4.28)

[Z]12 ¯En

s (cid:3) −∂2

t

These equations form the backbone of both our Maxwell and wave equation solvers. As
presented we have not presented any specific scheme to evolve both the field equation or
the particles. In the results presented, we have used Newmark-β to evolve the fields and a
fourth order Adams method to map the trajectory and velocity of the particles. These are
developed detail in Ref. [63] and are not repeated here.

4.6 Null spaces, gauging and other remarks

The methodology presented in this paper has several advantages over conventional
FEM-PIC; one of the fundamental challenges with time domain FEM is the presence of

67

a null-space. It has been shown the leap-frog mixed FEM (MFEM) does not suffer from
a null space when its initial conditions are specified. But the downsides of the leapfrog
MFEM is that it is conditionally stable with a mesh-dependent time step size. MFEM with
Newmark overcomes this fundamental bottleneck; time step size can be arbitrary large
(albeit it needs to be chosen to capture the physics). Unfortunately, there is a null space that
behaves as ∇φ(r) where φ(r) is some scalar function. The level of excitation of this null
space depends on methods used to invert the matrices in the time marching system. In
both the quasi-static limit and high tolerance of the iterative solver, this is small. Unlike
MFEM, the null space for the wave equation FEM (WE-FEM) grows like t∇(r). The effect
of this on charge conservation is evident in Ref. [63].

What is abundantly clear from the above, is that null space of both MFEM and WE-FEM
play a significant role on the spurious numerical artifacts that manifest themselves in
conservation laws. From this perspective, we note the following:

1. The solutions for ¯En

ns does not depend on a time integrator on the left hand side. This

is effected on the RHS when integrating over the path.

2. The resulting equation for ¯En

ns is exactly the discrete Gauss’ law with the equation of

continuity used to create the RHS.

3. The null space of these equation are harmonic solutions that are not excited provided

boundary conditions are properly imposed.

4. The solution to both ¯En
5. Despite the presence of the null space, by design, [∇·] ¯Bn (cid:3) 0 and the divergence of the

s have null spaces that are of the form ∇φ(r).

s and ¯Bn

solenoidal component to the electric flux density [∇]T[P]Λ

e [(cid:63)ε][C]e

c ¯En

s (cid:3) 0.

6. As a result, conservation laws are satisfied and not corrupted by null spaces.

68

s and ¯Bn

7. Parenthetically, we note that the number of degrees of freedom of ¯En

s is
identical for simply connected domains (via Euler’s relations). These may be useful
for some symplectic PIC schemes.
Next, we note that the above scheme bears a strong resemblance to the well known A− Φ
formulation. But we note some not-so-subtle differences. Note standard decomposition
relies on the electric field. We note that the decomposition in this paper relies on the
electric flux density, and is not a complete Helmholtz decomposition (which is not possible
in a discrete setting). Decomposition of the flux density induces a mesh based metric in
mapping between the primal and dual grids (can thought of simple mesh based anisotropy).
Using this is imperative in solving Maxwell’s equations. Gauge condition (Coulomb)
are, however, identical; likewise, final system of equations bears a strong similarity to
those obtained in an A − Φ formulation. By virtue of our accuracy in maintaining the
continuity equation, the use of the Coulomb gauge does not adversely impact our numerical
description of the space-charge and space-current sources[76].

4.7 Results

In what follows, we present a sequence of results to verify the methodology presented
earlier. We will show that (a) the solution to (4.24) does not have a null space provided one
imposes a zero mean constraint and the solution to (4.26) has a null space corresponding
to a gradient of a time invariant scalar potential. In the course of this test, we will also
show that the solution to (4.26) is unconditionally stable. Next, we will test the proposed
methodology against two test problems and validate the behaviour of the solution.

69

Figure 4.3: Eigenvalues Newmark-beta time stepping of MFEM and solenoidal-MFEM
system for PEC box.

4.7.1 Field Solver Validation

4.7.1.1 Eigenvalues

The temporal solution of (4.26) utilizes Newmark-beta as a time stepping method. To
verify that the resulting system is unconditionally stable and examine its null space, we
find the eigenvalues of the discrete system as described in Ref.[5]. Our discrete system is
obtained from discretizing a cube of side length 1m using tetrahedra with average edge
length of 0.214m. Dirichlet boundary condition (ˆn × E(t, r) (cid:3) 0) is imposed on the walls.
The eigenvalues for reduced system in (4.26) is compared against those obtained for a full
MFEM solve. As is evident from Fig. 4.3 eigenvalues for both systems lie on the unit
circle making the system energy conserving and unconditionally stable. Furthermore, the
presence of eigenvalues at zero indicate a null space that behaves as a gradient of potential,
a feature that both methods share.

70

4.7.1.2 Neumann Boundary Condition

Next, we examine the behavior of the system under Neumann boundary conditions.
Recall that one has impose a zero mean constraint on En
ns. This test validates proper
imposition of this condition as well as the solver as a whole. As before, our domain is a
cube with side length 1m. The mesh is discretized with an average edge length of 0.429m
yielding 133 tetrahedral elements. Neumann boundary conditions are equivalent to a plane
wave polarized along −ˆz and propagating along ˆx direction, with temporal variation being
the derivative of modulated Gaussian with center frequency of 10 MHz and bandwidth of
5 Mhz. The explicit formulae is fairly standard and can be obtained from say[5].

ns and ¯En

s . Alternatively, we can solve (4.24) for ¯En

We note the following. Given results from a standard MFEM solve, we can use (4.12) to
partition into solenoidal and non-solenoidal components. That is, given coefficients ¯En we
can obtain ¯En
ns and either (4.26) or (4.28)
for ¯En
s . So we compare can both the component wise solution from the reduced system
as well as the full solution. First at a random unknown in computational domain, we
compared the total electric field from both the quasi-Helmholtz decomposition and MFEM
system in Fig. 4.4. Next, the solutions from the standard MFEM system is decomposed
solenoidal and non-solenoidal components and compared with the solendoidal and non-
solenodial components of the quasi-Helmholtz system. In Fig 4.5 and 4.6 the non-solenoidal
components and solendodial components are compared.

A couple of points to note: (a) from Fig. 4.5 the solution of classical MFEM has a
null space (as expected) whereas that of (4.24) does not. (b) The solution of solenoidal
components of all three have a null space, in that they do not smoothly go to zero once the
excitation vanishes. (c) The agreement between all three is excellent (approximately single
precision) in the regime of interest despite all three being different systems.

71

Figure 4.4: Total electric field for plane wave through a box with Neumann boundary
conditions.

Figure 4.5: Non-solenoidal electric field for plane wave through a box with Neumann
boundary conditions.

72

Figure 4.6: Solenoidal electric field for plane wave through a box with Neumann boundary
conditions.

4.7.2 Particle Beam

Having validated the field solution, we turn our attention to inclusion of particles in
the system. We leverage the test bank developed in Ref. [63] for examples. In the first
test case electron macro particles are injected into a cylindrical cavity. The electrons repel
away from each other and expand. The profile of the trajectories can be compared with
other verified methods as well as quasi-analytic solutions[40]. Note, in this case ¯En
s will
correspond to cavity modes and ¯En
ns that will be the dominant component. Details of the
experiment are provided in 4.1.

In Fig. (4.7), the electric field in the radial direction half way down the tube 16 mm in
the radial direction is compared with the same run using a Newmark-beta formulation
with and without the quasi-Helmholtz decomposition (data for without from[63]). The
agreement between all three sets of data is excellent.

Next, we can examine each of the components. To wit, as before, we decompose ¯En
s and compare this data against that obtained from

from a full MFEM solve into ¯En
the reduced systems presented here (both (4.26) and (4.28)).

ns and ¯En

73

First, as is evident from Fig. 4.8 compares the ¯En
s components, they contain no DC
components and contain all the cavity modes. The ¯En
ns component can be seen in 4.9 is the
DC component shown in the total field Fig. 4.7. Next, the magnetic field shown in Fig.
4.10 shows excellent agreement between the magnetic field obtained from solution with
and without quasi-Helmholtz decomposition.

The removal of the growing in time t∇Φ(r) null space for the wave equation is of
considerable interest. In Fig. 4.11 we compare the MFEM and wave equation error in
discrete Gauss’s law without quasi-Helmholtz projectors. As evident, the error grows with
time for WE-FEM and thresholds at some value for MFEM. With the quasi-Helmholtz
decomposition, the story is very different. As is to be expected, the error for MFEM and
WE-FEM is identical and about a order below the results obtained earlier.

Table 4.1: Expanding Particle Beam Parameters

vp
vp/c

beam radius rb

Number particles per time step

Parameter

Cavity Radius
Cavity Length

Boundary Conditions

species

Turn on time
beam current

macro-particle size
min edge length
max edge length

Value
20 mm
100 mm
5 · 107 m/s
0.16678
8.00 mm

PEC

10

electrons

2 ns
0.25 A
52012.58
1.529 mm
6.872 mm

4.7.3 Adiabatic Expansion of Plasma

The final validation case is the adiabatic expansion of a plasma ball. This test case has
analytic solutions and allows for good comparison and validation[24]. In this example,
electrons and ions are placed in the center of a mesh with a Gaussian radial distribution.
Ions are given a low temperature 1k and the electrons 100k. The system is initially charge

74

Figure 4.7: Comparison of radial electric field ¯E between MFEM, qausi-MFEM, wave
equation, and quasi wave equation.

Figure 4.8: Comparison of radial rotational electric field ¯Es between MFEM, qausi-MFEM,
wave equation, and quasi wave equation.

75

Figure 4.9: Comparison of radial irrotational electric field ¯Eir between MFEM, qausi-MFEM,
wave equation, and quasi wave equation.

Figure 4.10: Comparison of radial cotree magnetic field between cotree solve and newmark
solve.

76

Figure 4.11: Comparison between regular mixed finite element, wave equation, and quasi-
Helmholtz wave equation error in the discrete continuity equation. The growing in time
DC null space in the wave equation is remedied by the quasi-Helmholtz decomposition.

neutral. The electrons initially expand outward creating an electric field that pulls the
ions outward. The density over time has analytic solutions and is presented along with
the measured density from the simulation in Fig. 4.12. Details of this test example can be
found in [40] and not repeated here. As is evident, the proposed solution agrees well with
analytic data.

4.8 Conclusion

In this paper, we have taken a step forward in building a FEM-PIC scheme that uses
implicit field solvers and is robust to corruption due to null spaces. Our approach has
been to define quasi-Helmholtz projectors, and in doing so, we show (a) satisfaction of
conservation laws as well as (b) correctness against data obtained earlier. An important
item to note is that the dimension of our system, in terms of number of degrees freedom,
is slightly smaller than earlier. Of course, this is a foundational paper in that we have
not addressed other pressing issues such as complexity, extension to higher order basis,
multiply-connected domain, sympletic methods, and so on. Work on some of these topics

77

Figure 4.12: Quasi-Helmholtz decomposition adiabatic expansion using quasi-MFEM.

is underway and will be presented elsewhere.

78

CHAPTER 5

CONCLUSION

This work presents issues arising from enforcing Gauss’s law in 3D finite element
particle-in-cell methods. To start things off, we provided three tests that can be used to
validate 3D PIC codes. These tests have been developed as (a) there is significant recent
interest in developing PIC codes to fully resolve complex boundaries, as properties of these
boundaries, and (b) as there is a void in examples with details that can be used to validate
these codes. In all examples that we have developed, there exist quasi-analytic solutions,
as well as open source codes that can be used for benchmarking results. We have also
provided sufficient detail so as to facilitate the user to examine all possible metrics that
may subtly affect the physics.

Next, we presented a solution to a problem that has been long-standing; charge
conserving FEM-PIC methods for implicit time stepping systems without the need to adopt
divergence cleaning. The method presented is agnostic to time stepping scheme and we
demonstrated the efficacy of this approach for a set of test problems for both MFEM and
WE.

Finally, null spaces that exist in the implicit system still need to be handled. To do, this
we developed a quasi-Hehlmholtz decomposition to remove the DC null space from both
the MFEM and wave equation. In doing so, we developed a quasi-Helmholtz projectors
that partitioned fields into non-solenoidal components for both the electric fields and
the magnetic flux density. We mapped the solenoidal component to a cotree system that
allowed the number of unknowns to remain the same. The method was demonstrated
for both the MFEM and WE-FEM systems showing Gauss’s Law is satisfied as well as
removing the DC null space.

Much work still needs to be done. Now that Gauss’s Law can be satisfied to higher order
in time, extension to higher order space are necessary. Utilization of domain decomposition

79

for the electromagnetic system is necessary for electrically large systems. How does
this affect the particle system? Other methods that can be easily parallelized such as
discontinuous Galerkin testing can Extension of the Quasi-Helmholtz decomposition to
multiply-connected domains. Current state of the art presented in this work only works for
simply connected domains. Extension to multiply connected domains require the use of
tearing and other techniques in the cotree. While they have been performed before adaption
to the current methods would be needed. Adding more test examples but in relativistic
scenarios and tightly coupled highly non-linear scenarios would provide a robust test
case for particle stepping and field stepping methods. This would have application in
relativistic tube devices such as magnetrons.

In terms of the physics, most of this thesis presented basic particle mesh interactions.
Extension of more complicated physics such as multi-pactor, collisions, emission and
chemistry would have impacts in areas such as pulsed power, material doping and micro-
plasmas. Developing modules that account for this physics and can be bolted on top of the
underlying PIC algorithm is highly desirable.

80

APPENDICES

81

APPENDIX A

A NOVEL PORT/NETWORK PARAMETER EXTRACTION TECHNIQUE FOR

COUPLING CIRCUITS WITH FULL-WAVE TIME DOMAIN INTEGRAL

EQUATION SOLVERS

Fully coupled transient analysis of electromagnetic-circuit systems is highly desirable
for a number of applications. Over the last decade or so, there have been several papers
that present effective solutions to this problem within the context of full-wave differential
or integral equation solvers. The method typically espoused is as follows; discretize
the underlying electromagnetic (EM) system, couple with the circuit subsystem self-
consistently, and then solve the coupled system. Within this framework, non-linearities are
easily accounted for in the solution process. This is in direct contrast to frequency domain
analysis, wherein one defines ports that couple the electromagnetic and circuit subsystems,
i.e., representation of the EM system using an effective impedance/admittance as seen at
the circuit terminals. This approach permits (a) solution only at the circuit level and (b)
one can readily incorporate different circuits without having to solve the coupled system
repeatedly. Achieving the same functionality in the time domain is a challenge. In this
paper, we present a solution to this problem. We will prescribe the means to effectively
represent the EM subsystem in terms of a transient admittance, thus facilitating a solution
to the coupled system via a circuit-level solve. Several results are presented that serve to
validate the proposed approach and demonstrate its effectiveness.

A.1

Introduction

The co-simulation of circuits together with full-wave electromagnetics is of paramount
importance to radio frequency (RF) device design, understanding and mitigation of
electromagnetic interference, signal integrity, and so on. This is a consequence of rapid
decrease in feature sizes of microwave circuit components, increased complexity, higher
frequencies, and tighter integration of passives and actives that permits close physical

82

proximity of these components. As a result, there is significant mutual dependence
between what can be considered as lumped elements and those that need a full wave
analysis. Not taking this coupling into account in a self-consistent manner results in designs
that can be wildly inaccurate. This need for accurate, rigorous self-consistent modeling
tools is positioned to take a higher importance given current trends in technology. The
approach that one takes is to develop modeling tools that permit both a lumped element
representation for circuits and distributed element representation for the fields, and solve
for the mutual coupling self consistently.

To a large extent, this problem has been visited several times in the recent past. An
extensive body of work exists in both frequency and time domain analysis; see [77, 78] and
references therein. In the frequency domain, the analysis of coupled circuit-electromagnetic
(CKT-EM) analysis is relatively straightforward. The incorporation of weak non-linearities
is effected via combining with harmonic balance (HB) [77, 79]. To incorporate more strongly
non-linear elements, typically one employs transient analysis.

The most extensive work in this vein relies on coupling a typical circuit solver with a finite
difference time domain (FDTD) EM field simulator [80, 81]. Several applications ranging
from analysis of packaging [82, 83] to high-speed interconnects [84, 85] to microwave
devices [86], small signal analysis [87], and active antennas [88] have been explored.
Similarly, there is extensive work on coupled full wave-circuit simulations using finite
elements [89, 90, 91, 92, 93] and time domain integral equations [78, 94, 95].

Given the depth of the existing literature, it is natural to question the nature of open
problems that this paper seeks to address. Common to all techniques is the following
process: (a) partition the CKT system into linear and non-linear contributions, (b) create a
time marching system for the coupled linear CKT and EM system, (c) impose the non-linear
contribution as a right hand-side, (d) solve for the local linearized system for each time
step, (e) use a non-linear solver to update the non-linear component, (f) check convergence
of the non-linear component and update steps (d) and (e) as necessary and (g) march

83

forward in time to obtain the entire time history. This process is efficient in that it does not
use a non-linear update of the combined CKT and EM system [92]. However, the challenge
lies in the computational cost incurred by steps (d) and (e) for each time step. In the
frequency domain analysis of the coupled CKT-EM system, this bottleneck is overcome by
extracting an equivalent port response of the linear EM system using a Schur complement
approach [77] together with harmonic balance to analyze non-linear CKT systems; this
is equivalent to obtaining the impulse response of the linear EM system as seen by the
circuit system. The advantages of this approach are (i) the non-linear iterations are done
only at the circuit level, (ii) it is efficient when the number of iterations required by the
non-linear solver is larger than the number of ports, (iii) it makes the analysis CKT agnostic,
and (iv) the equivalent currents in the EM system can be trivially reconstructed while post
processing at low cost. However, obtaining a transient impulse response is challenging for
a number of reasons. Using an impulse to excite the transient system is untenable as one
needs the excitation to be bandlimited. Alternatively, one could use a Gaussian instead
and then deconvolve to obtain the impulse response. Unfortunately it is well known that
deconvolution is unstable due to inaccuracies at the edges of the band. Another approach is
to use Hermite polynomials as a basis set (as their transforms [96] are analytically available),
but this procedure is limited to very narrow-band signals. A final, not so obvious, hurdle
is that the methods used for solving the transient EM system must be late time stable
for broadband signals. In this paper, we will prescribe an approach to overcome these
challenges.

Specifically, the contributions of this paper are as follows; we will (a) develop a
methodology for extracting the transient port admittance matrices of an EM system
with multiple ports; (b) integrate this with a CKT solver; (c) show that the responses
of a fully coupled system to those obtained using the aforementioned procedure are
identical; and (d) outline a frequency domain CKT-EM approach for extracting port
responses together with harmonic balance for analysis non-linear circuits. The crux to the

84

method relies on obtaining the electromagnetic response due a single basis function. This
implies that if certain properties of the system are satisfied, then the results obtained are
identical to that obtained using classical fully coupled systems. Furthermore, provided the
stated constraints are satisfied, the proposed approach can readily integrated into more
sophisticated fully coupled solvers that have been developed by others [94, 77, 95].

The rest of paper is organized as follows; Section II, presents a thorough formulation of
the problem. Section III presents the proposed method for extracting port admittance. A
number of numerical results are presented in Section IV that demonstrate the efficacy of
the proposed approach. Finally, Appendix I presents the corresponding frequency-domain
analogue of this approach. In what follows, the notation Z is used to denote matrices,
denotes Z vectors, and Z is used to denote a matrix of matrices.

A.2 Formulation of the Coupled EM-CKT System

To set the stage for the analysis, assume that a perfectly conducting object occupies
a domain D− residing in R3. Henceforth, this object is called the EM system. Attached
to the EM system are Np ports, each of which contain circuit sub-systems. With no
loss of generality, assume that that the CKT-EM system is excited by a field EEM,i(r, t) (cid:3)
Ei(r, t) + Ec(r, t) where Ei(r, t) denotes the incident electromagnetic field and Ec(r, t) is
due to coupling with source currents at the ports. Furthermore, it is assumed that the
excitation is temporally band-limited to fmax and is identically zero for t < 0. To analyze
this problem, we will assume the excitations induce voltages in the circuit system and an
electric current density J(r, t) in the EM system. These currents in turn radiate electric
and magnetic fields Es(r, t), Hs(r, t)}. The solution to the overall system is obtained by
self-consistent coupling between field and circuit solvers. Our field solver is based on time
domain integral equations (TDIEs), and has been shown to be late time stable. In what
follows,we describe each component briefly, culminating in the details of self consistent
coupling. Next, we discuss TDIEs used in the analysis; details can be obtained [97, 98].

85

A.2.1 TDIE solver

(a)

(b)

Figure A.1: (a) Perfect electric conducting surface with several circuits attached; (b) An
equivalent model–circuit system with a port admittance.

To formulate the electric field equation, we start with boundary conditions on the total

electric field, viz.,

ˆn׈n×Etotal(r, t)(cid:12)(cid:12)r∈D− (cid:3) 0,
(cid:3)⇒ ˆn׈n×(cid:16)

EEM,i(r, t) + Es(r, t)(cid:17)(cid:12)(cid:12)(cid:12)r∈D−

(cid:3) 0,

(A.1)

Casting Es(r, t) in terms of J(r, t), as shown Fig. A.1a, yields the time-domain electric field
integral equation (TD-EFIE)

where

ˆn׈n×EEM,i(r, t) (cid:3) LT{J(r, t)},

LTJ(r, t)} (cid:3) ˆn׈n×(cid:26) µ0∂t

Þ
dr(cid:48) δ(t − R/c)

Þ

− ∇
4π
0

4π

D−

R

(A.2)

(A.3)

(cid:63)t J(r(cid:48), t)

(cid:48))(cid:27)

dr(cid:48) δ(t − R/c)
Þ t
R
(cid:48)∇(cid:48)J(r(cid:48), t

dt

D−

(cid:63)t

0

Here, ∗t denotes the convolution operation in time t, R (cid:3) |r − r(cid:48)|, ∇(cid:48) denotes that the
derivatives are taken with respect to r(cid:48), c is the speed of light in free space, and ε0 and µ0
are the permittivity and permeability of free space, respectively. The solution to the above

86

equation is obtained by the following discretization scheme. First, the unknown surface
current is represented using a space-time basis as

Ns


Sm(r) Nt


P


m(cid:3)1

j(cid:3)0

p(cid:3)0

J(r, t) (cid:3)

j,pTEM,P
Im

j,p

(t),

(A.4)

where Im
j,p are the unknown coefficients to be determined at time step j, Ns is the number
of spatial unknowns, and Nt is the number of time steps. Further as in [98], the spatial
variation is represented using the standard Rao-Wilton-Glisson (RWG) function Sm(r), and
the temporal variation is represented using Lagrange polynomials of order P within a time
step j as

(t) (cid:3) TEM

(t − j∆t),

TEM,P
j,p

t ∈ [−∆t, 0]
otherwise

,

p

(cid:96)p(t)
P


0

i(cid:3)0
i(cid:44)l

t − ti
tp − ti

,

(t) (cid:3)

TEM
p

(cid:96)p(t) (cid:3)

tp (cid:3) (p/P − 1)∆t, p (cid:3) 0, 1, ..., P,

(A.5)

with ∆t is the time step size. Discretizing (A.2) using a Galerkin testing procedure results
in a system that can be represented as

with



ZEM
0
...

ZEM
k
...

0

(A.6)

(A.7)



IEM
0
IEM
1
...
...

IEM
Nt

(cid:3)



0

VEM
VEM
1
...
...

VEM
Nt



,





ZEMIEM (cid:3) VEM,

. . .

. . .

. . .
··· ZEM

k

. . .
··· ZEM

0

87

with k < Nt and

(cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171)

ZEM,11
j,00
ZEM,11
j,10
...

ZEM,Ns1
j,P0

ZEM

j

(cid:3)

ZEM,11
j,01
. . .

··· ZEM,1Ns

j,0P

...

··· ZEM,Ns Ns

j,PP

(cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)

,

ZEM,mm(cid:48)
j,pp(cid:48)

(cid:3) (cid:104)Sm(r)TEM,P

j,p

(t), LT{Sm(cid:48)(r)TEM,P
j,p(cid:48)

(t)(cid:105),

(cid:3) [IEM,1

j,0

IEM,1
j,1

··· IEM,Ns

j,P

]T ,

IEM
j

VEM

i

(cid:3) [vEM,1

i,0

vEM,1
i,1

··· vEM,Ns

i,P

]T ,

vEM,m
i,p

(cid:3)

1
η0

(cid:104)Sm(r)TEM,P

i,p

, n × n × EEM,i(r, t)(cid:105).

(A.8)

(A.9)

(A.10)

(A.11)

(A.12)

The EM system ZEM is represented in (A.6) by a block-banded lower triangular system of
equations of dimension NtNs(P +1)× NtNs(P +1), with IEM and VEM being the total degrees
s(P + 1)2
of freedom and excitation for all time. The entries of Eq. (A.8) are sparse with NtN2
non-zero entries, a consequence of the retarded potential. An apt illustration of such a
system can be found in [99, 94, 98, 100]. Here, IEM
contains the degrees of freedom for all
Ns(P + 1) spatial degrees of freedom associated with the time step at time t (cid:3) j∆t, and
VEM
represents the discretization of the incident field in (A.2) at time step i. The solution
to this lower triangular system is typically effected via a marching-on-in-time approach

j

i

i


j(cid:3)0

ZEM
i−j IEM

j

(cid:3) VEM

i

,

(A.13)

88

where i is the current time step. Solving for the currents at IEM
previous time steps IEM
are the lowest order basis in [98].

requires knowledge of the
, n (cid:3) 0, 1, ··· , i − 1. Note, the space-time basis used in this study

n

i

A.2.2 Circuit System

Next, we briefly describe a circuit system solver used in the analysis. As is typically
done, the analysis of each of the sub-systems is performed using a SPICE-like framework
that can incorporate an arbitrary network of linear and non-linear elements [101]. In
general, the system of equations representing the circuit physics can be written as

YCKT(t) ∗t eCKT(t) + iCKT,nl(cid:16)

eCKT(t)(cid:17)

(cid:3) iCKT,i(t),

(A.14)

where YCKT denotes the linear modified nodal admittance matrix (MNA), eCKT represents
the circuit unknowns, iCKT,nl denotes branch currents of non-linear elements, and iCKT,i(t) (cid:3)
ic,i(t) + iEM,i(t) denotes the system excitation due to current sources ic,i(t) and the external
electromagnetic fields iEM(t).

The MNA method poses the circuit problem in terms of node voltages and branch
currents (where applicable) as degrees of freedom by enforcing Kirchhoff’s current law
at each node. The circuit system is represented by an admittance matrix with some
modifications to allow for branch elements to be modeled. This approach allows for
both inductors and capacitors to be modeled using derivatives. To discretize the degrees
of freedom in the MNA system, we employ shifted piece-wise Lagrange interpolation
polynomials, representing e(t) as

eCKT(t) (cid:3)

eCKT
j

TCKT,Q
j

(t),

(A.15)

Nt


j(cid:3)0

89

where eCKT

j

represents the value of eCKT(j∆t), and

TCKT,Q
j

(t) (cid:3)



fq(t)gQ−q(t)

(q + j − 1)∆t ≤ t ≤ (q + j)∆t
q (cid:3) 0, . . . , Q

0,

elsewhere

q

Q−q


i(cid:3)1

i(cid:3)1

(t − i∆t)
−i∆t
(t + i∆t)

i∆t

,

.

fq(t) (cid:3)

gQ−q(t) (cid:3)

(A.16)

(A.17)

(A.18)

The circuit basis function consists of Q + 1 piece-wise elements, where Q is the polynomial
order to describe the circuit system. These basis functions have long been used for
simulating MOT systems [102]. While interpolatory basis functions are used together
with point testing to create circuit solvers, much remains to be explored regarding using
transient basis (and testing) like the one used in the EM system. Our approach has been
to restrict ourselves to classical methods, as using basis functions that are used for the
EM system would be tantamount to developing a Discontinuous Galerkin formulation for
circuits. The analysis of such an approach, study of its robustness and accuracy constitutes
a different body of work that lies outside the scope of this paper.

A circuit system of order Q is described in (A.19), similar to the EM system. YCKT, a
(Nt × NCKT) × (Nt × NCKT) sized matrix, is the global block modified nodal admittance
system where NCKT is defined as the number of spatial unknowns, eCKT is the vector of
unknown circuit quantities, and iCKT are the circuit excitations, for all time

YCKTeCKT + iCKT,nl (cid:3) iCKT,

(A.19)

90

where



YCKT (cid:3)

eCKT (cid:3)

iCKT (cid:3)

iCKT,nl (cid:3)

(cid:104)
(cid:104)
(cid:104)

. . .

YCKT

0
...

YCKT

Q
...

0

··· YCKT

Q

. . .
··· YCKT

0

(cid:105) T ,
(cid:105) T ,
(cid:105) T.

Nt

, ··· , eCKT
, ··· , iCKT,i
, ··· , iCKT,nl

Nt

Nt

eCKT
0
iCKT,i
0
iCKT,nl
0

, eCKT

1
, iCKT,i

1
, iCKT,nl

1



,

(A.20)

(A.21a)

(A.21b)

(A.21c)

0

The submatrix YCKT
contains all the purely resistive admittances plus partial expressions
for the derivatives at the current time step. The submatrices YCKT
contain
the terms of the derivative expressions pertaining to time history. As with the EM system,
the circuit system is solved using an marching-on-in-time (MOT) framework

through YCKT

Q

1

YCKT

0

eCKT
i

(cid:3) iCKT,i

i

− iCKT,nl

i

YCKT

q

eCKT
i−q

.

(A.22)

− Q


q(cid:3)1

A.2.3 Coupling the EM-CKT System

While the coupling between the two systems was implicit in the above two sections, our
scheme contains minor variations from what is commonly known in the literature [103].
The key difference is due to EM system that uses temporal basis functions defined within a
time step with Galerkin testing whereas the circuit system uses backward looking Lagrange
polynomials with delta testing (which is more commonly used in the literature). All other
existing methods use temporal basis functions that are interpolatory, i.e., coincident with
a time index. The consequence of this testing procedure is shown in (A.29), expanded

91

Figure A.2: EM-CKT port with current and voltage from circuit interacting with current
density and the electric fields.

out to include all the relevant matrix elements. The off diagonal blocks CCKT and CEM
are connectivity matrices between the circuit node voltages and RWG basis function edge.
Their values are made up of 1’s,-1’s and 0’s multiplied by different coefficients to take into
account the difference in time basis functions and unknown quantities. The matrix CCKT
couples the circuit node voltages to the voltage jump across the gap of the RWG basis
function as show in Fig. A.2. The coefficients apq handle the inner product between the
two time basis sets

q

TCKT,Q
j+q

(t) · TEM,P

(t)dt.

j,p

(A.23)

Þ j∆t

apq (cid:3)

(j−1)∆t

q

q

contains spatial coupling matrix CCKT

. The coupling matrices CCKT
, where q is the circuit
index in the coupling coefficient apq. The EM index p will vary with in CCKT
. The indexes
of CCKT
due to some port k require information on the connectivity of the positive and
negative nodes,+N(k),−N(k) in the CKT system, with the the edge number mk in the EM
in connectivity, being the transpose, but each
system. The matrix CEM is similar to CCKT
entry is multiplied by the edge length of the of basis function such that the surface current

q

q

q

92

density is coupled to the lumped circuit current.

. . .
··· CCKT

0

··· CCKT

Q



apq,

−apq,

0,

i,(P + 1)mk + p
j, N +(k)
i,(P + 1)mk + p
j, N−(k)
elsewhere

CEM
0
...

0

0
CEM

···

. . .

0

CEM





CCKT

0
...

CCKT

Q
...

0

CCKT (cid:3)

[CCKT

q

]i,j (cid:3)



CEM (cid:3)

[CEM]i,j (cid:3)

(cid:96)m,

−(cid:96)m,



0,

ZEM CCKT

CEM YCKT



,



,

(A.24)

(A.25)

,

(A.26)

i, N +(Port)
j,(P + 1)mPort + p
i, N−(Port)
j,(P + 1)mPort + p
elsewhere

 IEM

eCKT

 (cid:3)

VEM

iCKT

 ,

93

,

(A.27)

(A.28)

ZEM
−
i−1
−
Q

0

CCKT
CEM YCKT

0

0

 (cid:3)


 IEM
−
Q

i
eCKT
i

j(cid:3)0 ZEM
q(cid:3)1 YCKT

i−j IEM
j
eCKT
i−q

q

q

eCKT
q(cid:3)1 CCKT
i−q
− iCKT,nl(eCKT
)

i

VEM

i
iCKT
i

(A.29)

 ,

Since the circuit system is delta tested, the coupling from the EM system to the CKT system
only depends on the present time step.

Non-linear elements are incorporated by adding their I-V relationship into the MNA
system of equations and solved using the Newton-Raphson method. A direct approach to
solving this is to use a global minimization approach [94], that has a benefit of the Jacobian
available analytically. Alternatively, as the total number of linear unknowns is significantly
larger than the number of non-linear unknowns, we apply the Newton method in a manner
similar to [78],[92]. This approach formulates the Jacobian for each non-linear element as
opposed to the whole system.

A.3 Transient Port Extraction

The computational bottleneck in solving (A.29) is (a) the evaluation of the convolution
on the right hand side, and (b) inversion of (A.29) for each Newton iteration. The cost
associated with convolution has been solved with the introduction of the plane wave
time domain and the adaptive integral methods. To amortize costs associated with each
Newton step, we develop a port model. Extracting a port response is tantamount to
extracting the impulse response of the linear system (in this case, the EM system). It is
shown in the Appendix that performing a similar analysis in the frequency domain is
tantamount to using a Schur complement of the EM system. The advantages of such a
formulation are many-fold; (a) the non-linear solve is restricted to the circuit portion of the
simulation domain, which is typically small; (b) the non-linear solver is not subject to an
entire domain global minimization (if it is so desired); and (c) since the circuit and EM
system are effectively decoupled, it is possible to simulate different circuit topologies that

94

could be paired with the same EM system. However, translating these ideas to the time
domain presents a challenge.

To set the stage for the solution to this problem, consider the frequency domain
equivalent method shown in the Appendix. To obtain the equivalent port admittance, we
excite the system with unitary voltage while shorting the rest to obtain a column of YEM;
this procedure is then repeated for all ports. The time domain equivalent of this procedure
would to obtain an impulse response; it is obvious from a theoretical perspective that
excitation via a unit impulse is untenable. A potential indirect approach is (a) obtain the
response of the EM system due to a Gaussian excitation and (b) using Fourier transform to
obtain the impulse response. However, as has been noted in a number of papers [104], this
“deconvolution” is inherently unstable. Our approach in solving this problem is to develop
a framework that depends only on time domain.

The key ideas pertinent to the development of the port extraction are as follows; (a)
extraction response due to one temporal circuit basis and (b) use the response due to the
temporal basis to construct the effective transient admittance due the EM system. The
rationale for this approach is as follows; as any excitation can be represented in terms of
these temporal basis sets, the response to one basis can be used to recover the EM response
to this excitation (or construct the effective admittance due to the EM system). However,
obtaining the response due to one temporal basis function places additional constraints; (i)
the testing scheme should be such that it would distinguish the system from a unit-impulse,
and (ii) a stable TDIE solver. If these conditions are satisfied, then the results of the port
extracted system would be indistinguishable from the full-solve, to machine precision. Both
conditions are satisfied by the TDIE (EFIE) solver in [98, 97]. Note, the matrix CCKT
contains
the requisite inner product of temporal basis sets that couple the EM and circuit system.
Exciting one of the coefficients of eCKT
that correspond to the nodes associated with the

q

j

95

port, n and time step j results in a right hand side

ZEM
i−j IEM

j

CCKT

q

eCKT,n
i−q

(cid:3) 0,

(A.30)

− Q

1

q(cid:3)0

0

(cid:3)

eCKT,n
j

j (cid:3) 1, k (cid:3) n

j (cid:44) 1, k (cid:44) n

.

The solution to the above system can be effected using a MOT scheme. After finding IEM
,
for port n for all time, we store only the coefficients corresponding to currents at observer
ports n(cid:48) to formulate a column of the reduced order model for the EM system ˜YEM
j which
is a matrix whose dimensions are N ports × N ports × Nt. An entry of ˜YEM
n,n(cid:48),j corresponds to
the time history at time step j of port n(cid:48) excited by port n is given by ˜YEM
. As a
nn(cid:48),j
result, the system ˜YEM

can be written as

(cid:3) IEM,n(cid:48)

j

j

j

˜YEM
j

(cid:3)

˜YEM
11,j
˜YEM
21,j
...
˜YEM
Nports1,j

˜YEM
12,j
. . .

···

···

˜YEM
1Nports ,j
...

˜YEM
Nports Nports ,j

(cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)

.

(A.31)

The entries of ˜YEM constitute the impulse response of the EM system, which is agnostic to
the circuit system. To use this information within the circuit solver we need to ensure that
it is mapped to the correct nodes of the circuit system. In effect, the result of this mapping
is an equivalent transient admittance that can be incorporated into the MNA system. The

i


j(cid:3)0

(cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171)

96

mapping from ˜YEM to YEM can be described by



Q


q(cid:3)0

j

i


j(cid:3)0

Nports


Nports


l(cid:3)1

k(cid:3)1

YEM
m,n,j

(cid:3)

− ˜YEM

l,k,j

˜YEM
l,k,j

,

˜YEM
l,k,j

,

− ˜YEM

l,k,j

0,

, m (cid:3) N +(l)
n (cid:3) N +(k)
m (cid:3) N +(l)
n (cid:3) N−(k)
m (cid:3) N−(l)
n (cid:3) N +(k)
, m (cid:3) N−(l)
n (cid:3) N−(k)
elsewhere

,

(A.32)

where l, k are port indices and N +(k) and N−(k), as defined in section A.2.3, are the positive
and negative node indices associated with port k. Once YEM has been obtained, the circuit
unknowns at each time step eCKT

can be obtained using

YEM
i−j eCKT

j

+

YCKT

q

eCKT
i−q

(cid:3) iCKT

i

,

(A.33)

The proposed methodology bears a stark similarity to its frequency-domain counterpart,
presented in the Appendix. Algorithmically, we summarize the time domain procedure as
follows:

1.

a) For each port we use eCKT,n

j

where in we set all the entries to zeros (shorted),

except only the entry associated with the port m.

b) Find ˜YEM the time series matrix (A.31), that contains admittances seen by the
ports. This is mathematically equivalent to finding IEM (cid:3) ZEM−1CCKTeCKT for
each port n and storing only the component of IEM seen by all ports.

97

2. Find YEM by reorganizing ˜YEM with the circuit information needed to solve the

system.

3. Solve the reduced system of equations using the Newton procedure specified in

Section A.2.3.

A.3.1 Computational Cost

The benefits of this method come to play when solving a non-linear circuit system
with a very large linear electromagnetic system. The cost of solving a time dependent
non-linear system of equation comes from three components; at any time step (a) the
number of non-linear iterations to converge a solution, (b) the cost for the computation the
right hand side for each iteration, and (c) the cost of solving the system for each iteration.
We note that the matrix system to be solved is sparse, and therefore, the cost of an iterative
solver is directly proportional the number of degrees of freedom; in the case fully coupled
system it would be O(Ns + NCKT), that for port extraction scales as O(NCKT). Further,
given the advances in fast methods [105, 94], we will assume that cost of computation of
the right hand side scales as O(NtNs log2 Ns). Next, we denote the number of non-linear
interactions/per time step as N NL
iter. Using these, it follows that the total cost of for a fully
coupled system scales as

CFC (cid:3) O(NtN NL

iter

System for each NL step

(Ns + NCKT))

(cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125)
(cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125)

(NCKT))

+ O(NtNs log2 Ns)

(cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125)
(cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125)

+ O(NtNPNs log2 Ns)

RHS cost

iter

whereas that of port-extraction scales as
CPE (cid:3) O(NtN NL

(A.34)

(A.35)

System for each NL step
In a typical simulation, the N NL
iter is approximately identical for both port extraction and
fully coupled solves. As a result, it is evident from the above cost complexity that when
N NL
iter is larger than NP (it is in most cases), the port extraction scheme is significantly faster.

PE cost

98

However, we note that there is significant potential for gain as port extraction is circuit
agnostic. As a result, the effective admittance of the EM system can be used with multiple
circuit configurations without the need to simulate the EM system.

A.4 Numerical Results

In this section, we present a number of numerical results that demonstrate the validity
of the proposed method for both linear and non-linear circuit systems. We also present
comparisons with the corresponding frequency domain method presented in the Appendix.

A.4.1 Dipole Antenna with Chebyshev Filter

0.5 m

VPort(t)

−

+

0.5 m

50 Ω

16.308 nH 16.308 nH

1.34 cm

Vs(t)

9.05 pF

13.48 pF

9.05 pF

+
Vport(t)
−

Circuit

(a)

(b)

Figure A.3: (a) Dipole antenna geometry and coupling. (b) Chebyshev low-pass filter
schematic.

In this first example we validate the port extraction procedure by comparing it to the
results of the full coupled system (A.29) as well as a frequency domain harmonic balance
solution. Consider the center-fed dipole Fig. A.3a loaded with a fifth order Chebyshev
low-pass filter illustrated in Fig. A.3b. The chosen width and length, 1 m by 1.34 cm, of the
dipole causes the fundamental mode to be at 150 MHz. The dipole is discretized using
40 RWG basis functions with the delta-gap feed placed at the center edge. A modulated
Gaussian excites the Chebyshev low-pass filter with a center frequency of f0 (cid:3) 300 MHz
and a bandwidth of BW (cid:3) 250 MHz (A.36)

v(t) (cid:3) cos(cid:0)2π f0t(cid:1)e

99

−t2/(2σ2),

(A.36)

σ (cid:3)

3

2πBW

.

We set the time step size to ∆t (cid:3) 90.909 ps to capture the spectral content of the modulated
Gaussian. A third order circuit time basis function is chosen to properly capture the circuit
physics. Lower order circuit time basis functions attenuate the higher circuit frequencies.
A zeroth order Q (cid:3) 0 EM circuit time basis function captures the EM physics with adequate
precision.

The L2-norm of the difference between the two transient solutions is 6.7588× 10−15, and

the solution shows good agreement with the frequency-domain solution.

(a)

(b)

Figure A.4: Comparison of port voltages obtained via the port extraction method vs.
solving either (a) the full coupled TDIE system directly (left) and (b) a harmonic-balance
solution using frequency-domain port extraction.

A.4.2 Bowtie Antenna with Diode Mixer

In this next example, we demonstrate the application of port extraction to analyze a
non-linear system. We compare the port voltage in time at steady state, with an inverse
Fourier-transformed solution to the harmonic balance (HB) solution. Consider a bowtie
antenna coupled to a diode mixer, as shown in Fig. A.7b. The spatial discretization of the
bowtie antenna results in 7767 unknowns to capture frequencies up to 1.85 GHz, allowing

100

00.20.40.60.81·10−7−0.4−0.200.20.4Time(s)PortVoltage(V)FullSolvePortExtraction010020030040050060000.20.40.60.81Frequency(MHz)PortVoltage(V)TDIE-CKTHBfor proper representation of several harmonics. The port extraction results were obtained
using first order temporal basis and third order basis were used for the circuit solver. The
circuit design in Fig. A.7a mixes two signals, the RF and local oscillator (LO) input, so that
the output Vport contains integer multiples of RF and LO signal frequencies. We set the RF
source to 900 MHz, the LO to 800 MHz, and the time step to 27.777 ps to properly capture
frequencies supported by the mesh. Setting Vbias to 0.7 V turns on the diode. The RF and
LO amplitudes are set to 0.4 V. In our analysis, we employ the classic large-signal model,
characterized by the following i − v relationship:

I (cid:3) Isat(e

q

ηKb T V − 1).

(A.37)

Table A.1: Simulation parameters used in the diode mixer example

Parameter
Spatial Unknowns
fRF
fLO
∆t
P
Q
Isat
η
kBT
q
R
θ
a
b

Value
7767
900 MHz
800 MHz
27.777 ps
0
3
2.0 nA
2.0
0.02565 V
0.38
90◦
0.01m
0.01m

The diode parameters, along with other simulation parameters, are given in Tab. A.1.
As seen in Fig. A.5, there is agreement with the time harmonic result. Furthermore, it is
evident from Fig. A.6 that the port extracted system is stable.

A.4.3 Dipole Antenna with Diode Mixer

To illustrate the computational cost savings of the port extraction methods as well as the
accuracy with a non-linear circuit, we examine a dipole A.3a with 39 unknowns connected

101

Figure A.5: Port voltages at the bowtie antenna input due to the diode attenuator circuit
obtained via the TD port extraction method. The solution is compared with a Harmonic
Balance (HB) solution using 45 harmonics.

Figure A.6: Port extracted signal from bowtie using third order circuit time basis function

to a non-linear mixer Fig. A.7a with 15 unknowns. We set the RF source to 300 MHz and
LO to 200 MHz and the time step to 100 ps. The Vbias, RF, and LO voltages are set to 0.4 V.
The lowest order EM P (cid:3) 0 and circuit Q (cid:3) 1 basis functions are used. We measured the
number of non-linear iterations and the time to solve the system Fig. A.8. The ratio of the

102

1 pF

4

5

7

LO

.1 pF

8

1

RF

10 Ω

1 pF

100 Ω 2

1 M Ω

100 Ω

3

6

9

6 nH

3 nH

1kΩ

Vport

10 Ω

−+ Vbias

(a)

R

a

θ

b

Circuit

(b)

(b) Layout of the bowtie antenna and

Figure A.7: (a) Diode mixer circuit schematic.
coupling with circuit.
slopes is proportional to (Ns + NCKT)/NCKT following the scaling laws presented in section
A.3.1. The slopes for the full coupled mixer and port extraction mixer are 3.200 · 10−6 and
8.606· 10−7 seconds per iteration respectively resulting in the unknown ratios and the slope
ratios as 3.6 and 3.719.

Figure A.8: Comparison of the number of non-linear iterations and time for the port
extraction and full coupled method simulating a dipole with a non-linear mixer.

103

Figure A.9: Comparison between port extraction and full coupled methods for dipole
with diode mixer. The L2 − norm of the difference between the two transient solutions is
1.171 · 10−16

A.4.4 Diode Attenuator

In this third example, we demonstrate the method’s ability to handle multiple ports.
Consider the λ/4 microstrip line excited with a fr f (cid:3) 4 GHz source as illustrated by Fig.
A.10a. On the opposite end we place a diode attenuator circuit which attenuates the RF
signal depending on the resistance Rc. The microstrip geometry is 2 cm by 1 mm sitting
0.5 mm above a 1 cm by 2 cm ground plane. The ground plane and strip are connected on
their edges by a 0.5mm tall and 1mm wide vertical tabs. The spatial discretization of the
microstrip geometry results in 683 unknowns. We set the time step to 0.4545 ps with the
EM time basis order set to two, P (cid:3) 2, and the circuit basis order set to three, Q (cid:3) 3.

We measured the amplitude of the sinusoid for the voltage across the 86 Ω resistor
for varying values of Rc. The results displayed in Fig. A.10b show excellent agreement
between the HB and port extraction methods.

104

Table A.2: Parameters used in diode attenuator example

Simulation Parameter Value
Spatial Unknowns
fRF
∆t
P
Q
Isat
η
kBT
q

683
4 GHz
0.4545 ps
2
3
2.0 nA
1.8
0.02565 V

λ/4

Rc
−+ Bias
1 nH
1 pF

1 pF

86 Ω

1 86 Ω

RF

(a)

(b)

Figure A.10: (a) Circuit diagram of the two-port microstrip interconnect coupled to the
diode attenuator circuit; (b) port voltage waveforms computed in both TD and FD codes
for Rc (cid:3) 10 kΩ.

105

-1.5-1-0.5 0 0.5 1 1.5 2 0 2e-10 4e-10 6e-10 8e-10 1e-09Output voltage (mV)Time (s)HB+FDIECKT+TDIETable A.3: Table of voltages across the 86 Ω output resistor in both TD and FD codes for
various values of the control resistance Rc

Rc(Ω) HB (V)
86
3.59e-4
3.661e-4
100
7.045e-4
1000
10000
1.201e-3
1.119e-3
100000

TDIE (V)
3.238e-4
3.314e-4
6.697e-4
1.209e-3
1.116e-3

86 Ω

1 nH

−+ Vbias

1 pF

86 Ω

1 pF

1pF

1 nH

86 Ω

RF

86 Ω

LO

(a)

(b)

Figure A.11: (a) Circuit diagram of the four-port hybrid-90◦ structure coupled with the
diode mixer circuit; (b) voltage waveform across the 86 Ω output resistor calculated in both
the TD and FD codes.

106

-2 0 2 4 6 0 5e-10 1e-09 1.5e-09 2e-09 2.5e-09 3e-09 3.5e-09 4e-09Output voltage (mV)Time (s)HB+FDIECKT+TDIE60 mm

b

c

a

c

c

40 mm

(a)

Figure A.12: Balanced Mixer geometry

A.4.5 Balanced Mixer

Finally, we consider the case of a four-port microstrip circuit coupled with a two-diode
mixer circuit. Voltage sources of fRF (cid:3) 500 MHz and fLO (cid:3) 1 GHz are connected at the
input ports, while the other two ports couple into the mixer, as shown in Fig. A.12a. The
time step is set to 4.545 ps with a EM time basis order of 2 and a circuit time basis order of 3.
The EM structure comprises a number of microstrip traces connected at their ends via tabs
of length 0.5 mm to a 30 mm × 20 mm ground plane. The geometry is discretized using
994 RWG basis functions. The diode saturation current is Isat (cid:3) 2 nA and the efficiency
parameter η is set at 1.5. The thermal voltage, being composed of the Boltzmann constant
kB, temperature T and electron charge q, are set at kBT

q (cid:3) 0.02565 V.

Fig. A.11b shows the output voltage waveforms produced by both the TD and FD
(via inverse FFT) codes. The waveforms show good agreement with each other, and the
nonlinearity of the circuit system is pronounced.

A.5 Conclusion

In this paper, we have presented a method to perform port extraction in the time domain
for efficient analysis of coupled EM-CKT systems. The approach relies on two aspects; (a)

107

Table A.4: Parameters used in Balanced Mixer example

Simulation Parameter Value
Spatial Unknowns
fRF
fLO
∆t
P
Q
Isat
η
kBT
q
a
b
c

994
500 MHz
1 GHz
4.545 ps
2
3
2.0 nA
1.5
0.02565 V
1 mm
1.79 mm
19.5 mm

the use of an highly stable transient integral equation analysis methodology, (b) using a
transient CKT basis function to excite the system, and (c) using temporal Galerkin testing
that couples CKTs to the EM system. As is evident from the result, as is to be expected,
the port extracted system produces equivalent admittances that are stable, results that are
identical to the fully coupled system, and sufficiently captures higher order non-linearities.
Several results are presented that validate the presented approach. Extension of this
approach to more complex TDIE system involving composite objects is currently under
study as is coupling with finite element solvers; these will be published in subsequent
papers.

A.6 Acknowledgement

This work was supported in part by Grant NSF ECCS-1408115 and in part by the
Department of Defense High Performance Computing and Modernization Program. We
would also like to thank HPCC facility at Michigan State University.

Harmonic balance simulation of EM-CKT systems in the frequency domain

In what follows, we briefly present a coupled frequency domain method of moment
solver that is coupled with a circuit solver within a Schur complement framework. The

108

presented framework is used to validate the transient analysis framework presented in this
paper.

A.6.1 Frequency-Domain EM Solver

The frequency-domain analogue of (A.6) is given by

ZEMIEM (cid:3) VEM,

(A.38)
where ZEM is the Ns × Ns impedance matrix, IEM is the Ns × 1 vector of unknown RWG
coefficients Im, and VEM is the Ns × 1 vector representing the discretized incident field.
The elements of the impedance matrix ZEM are given by
(cid:41)
(cid:48) e−jk0|r−r(cid:48)|
|r − r(cid:48)| Sm(r) · Sn(r(cid:48))

Þ
dS
dS
(cid:48) e−jk0|r−r(cid:48)|
|r − r(cid:48)| ∇ · Sm(r)∇(cid:48) · Sn(r(cid:48))

(cid:26)Þ
Þ

jk0η0
4π

∂D

dS

Þ

ZEM

mn (cid:3)
− 1
k2
0

dS

∂D

∂D

and the elements of the incident field vector VEM are given as

∂D

Þ

∂D

Sn(r) · Ei(r) dS.

(A.39)

(A.40)

VEM
m (cid:3)

A.6.2 Circuit formulation

The circuit portion of the overall problem is cast within the framework of the modified
nodal admittance (MNA) method [101] and solved via the harmonic balance (HB) method
[106, 107, 108, 79]. HB is often used to find steady-state solutions to nonlinear differential
equations by assuming the solution has the form of a finite power series of a complex
exponential with fundamental frequency f0. Harmonic generation and intermodulation
effects are accounted for in the dual domain through Fourier transforms. We partition
the circuit into a subcircuit containing all of the linear elements and another subcircuit
containing all of the nonlinear devices, i.e. those with nonlinear i − v relations, connected
via ports consisting of node pairs. The full problem is described by the action of the

109

block-diagonal harmonic MNA matrix YCKT and the nonlinear operator YNL acting on the
circuit unknowns eCKT in response to the excitation iCKT:

(cid:3) iCKT.

(A.41)

YCKTeCKT + YNL(cid:16)

eCKT(cid:17)

Each sub-block corresponds to an independent linear MNA system at each of the N f + 1
simulation frequencies.

The HB algorithm proceeds as follows. Given an initial guess of the voltages in the

system:

1. Calculate port voltages V and compute global vector of linear-side port currents ILIN

2. At each port, obtain time-signature of port voltage via inverse Fourier transform (IFT)

and compute currents through nonlinear devices via i − v relations

3. Fourier transform (FT) each resulting nonlinear port current and assemble global

vector of nonlinear-side port currents INL

4. Evaluate harmonic error 
 (cid:3)(cid:12)(cid:12)(cid:12)(cid:12)ILIN(V) + INL(V)(cid:12)(cid:12)(cid:12)(cid:12)2

5. If 
 is satisfactory, exit; else, update V according to iterative scheme and go to step 2

For solving this system of equations, we employ a Newton-Raphson method equipped
with backtracking line search [109].

A.6.3 EM-CKT formulation

With formulations in hand for both the EM and CKT systems, we now describe our
frequency-domain EM-CKT co-simulation approach. This discussion follows much along
the lines of [77, 78, 94], but it is presented here for the sake of completeness.

As in the time-domain system, we model the coupling between the EM and CKT
systems using a delta-gap style port. In the interest of coupling systems (A.38) and (A.41) at
a single frequency, we re-use our result from (A.27) to define the coupling matrix C (cid:3) CEM.

110

For coupling the CKT degrees of freedom into the EM system, we use a sparse matrix
B (cid:3) −CT of size Ns × (Nnv + Nbc). Using the block structure to represent a system of
independent linear systems which interact via a nonlinear function YNL as in (A.41), the
full coupled EM-CKT system is then posed as

ZEM

B

C YCKT



 IEM

eCKT

 +

(cid:3)

0

YNL(cid:0)eCKT(cid:1)

 ,
VEM

iCKT

(A.42)

where, as before, the boldface quantities represent block-diagonal and block versions of the
previously-introduced matrices and vectors, respectively, with each block corresponding to
a different simulation frequency. Note that the size of the EM system is N f Ns × N f Ns, since
the formulation (A.38) does not facilitate DC analysis. This system of equations may be
solved as (A.41) using the HB approach; however, because the EM system ZEM comprises
a block diagonal set of dense matrices that may be large in dimension, obtaining the initial
guess and finding port currents at each iteration may be quite expensive. To circumvent
this issue, we employ a technique based on Schur complement solution of (A.42) [77]. This
technique allows the EM system to be solved once, and the result to be recycled for use
with arbitrary circuit systems so long as the port edges in the EM structure remain the
same.

By making the back-substitution IEM (cid:3)

row, we obtain

(cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171)YCKT −C
(cid:104)
ZEM(cid:105)−1
(cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125)

(cid:17)YEM

B

(cid:104)

ZEM(cid:105)−1(cid:16)
VEM − BeCKT(cid:17)
(cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)eCKT (cid:3) iCKT
ZEM(cid:105)−1
(cid:104)
(cid:124)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:123)(cid:122)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:125)

−C

VEM

.

(cid:17)iEM

111

in the second block

(A.43)

The nonzero elements of YEM may be stored along with the associated port numbers
and frequencies, and the matrix can then be easily assembled for inclusion with arbitrary
circuit topologies. This procedure reduces the system size from N f (Ns + Nnv + Nbc) to
N f (Nnv + Nbc), drastically decreasing the computational cost of the Newton-Raphson
iteration within the context of repeated EM-CKT co-simulation at the expense of a one-time
EM solve. It is also clear from (A.43) that an EM excitation may also be incorporated in the
form of the vector iEM.

112

APPENDIX B

SIMPLE CHARGE CONSERVATION EXAMPLE

B.1 Appendix: Charge Continuity Example

A contextual example of how charge continuity works can be shown through a simple
analytic test case. Here we provide several examples that illustrate charge conservation for
linear path, higher order path, crossing multiple cells, and numerical integration.

To set up this problem we use two 2-D cells comprised of Whitney-1 basis functions as
show in Fig. B.1 and B.3. We define the Whitney-1 basis function for the edge between
node i and j as,

Ni,j(r) (cid:3) λi∇λ j − λ j∇λi,

(B.1)

where λi is barycentric coordinate for node i. The node locations for the cells are
r1 (cid:3) [0, 0],r2 (cid:3) [1, 0],r3 (cid:3) [0, 1] and r4 (cid:3) [1, 1], with cell 1 containing nodes 1, 2, and 3, and
cell 2 containing nodes 2, 3, and 4. For simplicity we choose a path y (cid:3) x3 for the particle
to travel through. We parameterize the path by t,

r(t) (cid:3) x(t)ˆx + y(t) ˆy

(cid:3) t ˆx + t3 ˆy.

We can also write the velocity of the particle as the following,

v(t) (cid:3) x(t)ˆx + y(t) ˆy

(B.2)

(B.3)

(B.4)

(B.5)
To measure the charge at each node we use Whitney-0 W0(r) (cid:3) λi basis function, which

(cid:3) ˆx + 3t2 ˆy.

is just the barycentric coordinate of the location of the particle,

113

Figure B.1: Basis functions N1,2(r), N1,3(r), N2,3(r) for cell 1 with cubic path of the particle.

λ1(x, y) (cid:3) 1 − x − y
λ2(x, y) (cid:3) x
λ3(x, y) (cid:3) y

(B.6)

(B.7)

(B.8)

Since we have the path of the particle defined in terms of t, we can find the charge

measured at each node at time t.

λ1(t) (cid:3) 1 − t − t3
λ2(t) (cid:3) t
λ3(t) (cid:3) t3

(B.9)

(B.10)

(B.11)

B.1.1

Instantaneous Charge Conservation

We the analytic path described, we can look at charge conservation at any point in time

t0 along the path.

114

To do this we first start by taking the time derivative of the charge at each node,

∂λ1(t)
∂t
∂λ2(t)
∂t
∂λ3(t)
∂t

(cid:3) −1 − 3t2
(cid:3) −3t2

(cid:3) 1.

(B.12)

(B.13)

(B.14)

To help formulate the Whitney-1 basis function we need the gradient of the barycentric

coordinates.

∇λ1 (cid:3) −ˆx − ˆy
∇λ2 (cid:3) ˆy
∇λ3 (cid:3) ˆx

We can then define each basis function, in space, then in time,

W1
12

W1
13

W1
23

(r(t)) (cid:3) λ1∇λ2 − λ2∇λ1
(cid:3) y ˆx + (1 − x) ˆy
(cid:3) t3 ˆx + (1 − t) ˆy
(r(t)) (cid:3) λ1∇λ3 − λ3∇λ1
(cid:3) (1 − y)ˆx + x ˆy
(cid:3) (1 − t3)ˆx + t ˆy
(r(t)) (cid:3) λ2∇λ3 − λ3∇λ2

(cid:3) −y ˆx + x ˆy
(cid:3) −t3 ˆx + t ˆy.

(B.15)

(B.16)

(B.17)

(B.18)

(B.19)

(B.20)

(B.21)

(B.22)

(B.23)

(B.24)

(B.25)

(B.26)

Now that we have the path of the particle measured at each node in time, and the value
of the basis function at the particle over time, we can check to see if the continuity equation
holds for node 1.

115

0 (cid:3)

0 (cid:3)

+ ∇ · J

+ J1,2 + J1,3
+ Q(cid:104)W1

∂ρ
∂t
∂ρ1
∂t
∂λ1
0 (cid:3) Q
∂t
0 (cid:3) −1 − 3t2 + [t3 ˆx + (1 − t) ˆy][ˆx + 3t2 ˆy] + [(1 − t3)ˆx + t ˆy][ˆx + 3t2 ˆy]
0 (cid:3) −1 − 3t2 + [−2t3 + 3t2] + [1 + 2t3]
0 (cid:3) 0

(t), v(t)(cid:105) + Q(cid:104)W1

(t), v(t)(cid:105)

1,3

1,2

(B.27)

(B.28)

(B.29)

(B.30)

(B.31)

(B.32)

This shows that when charge is measured at some node i with W0
i

(B.33)
(r(t)) basis function,
(r(t)), the discrete (in space) continuity equation is held.

and the current is tested with W1
i,j

B.1.2 Linear Path

Standard PIC schemes separate charge and current by a half time step. This means
you have positions at integer time steps and velocity at half time steps. To rectify charge
conservation in this setup you assume the path in linear. To illustrate this two points on
the cubic path are chosen at tn (cid:3) 0.45 and tn+1 (cid:3) 0.55 and shown in Fig. B.2. To find
the current we need the integral of the path from tn (cid:3) 0.45 to tn+1 (cid:3) 0.55. This has been
extensively discussed in [6] and an analytic expression for the linear integral was presented
in [15].

k

(cid:3)

jn+1/2

Nk(r) · qv(t)δ(t − r(t))dt (cid:3)

qNk(r(t)) · v(t)dt
If we assume a linear path then v(t) (cid:3) vn+1/2 and eq. (B.34) should become,
− λ f

Þ r2

Þ t2

qNk(r(t)) · v(t)dt (cid:3)

Nk(r) · dL (cid:3)

jn+1/2

(λs

t1

t1

λ f
j

i

λs
j

i

(B.34)

)

(B.35)

Þ t2

k

(cid:3)

t1

q
∆t

r1

q
∆t

Þ t2

To compute jn+1/2

k

for a linear path we have several options; quadrature points, Teixeria’s

methods and analytic integration.

116

Figure B.2: Two points are chosen to show the standard linear charge conservation (left).
Zoomed in showing curved analytic path and linear path (right).

Table B.1: Values of charge conservation for first order time step.

Qn+1

0.55

Qn

0.45

0.283624

0.458875

0.166375

0.091125

Node 1
Node 2
Node 3

Qn+1−Qn
-1.7525

∆t

1.0

0.7525

∇ · J

1.2475 + 0.5050
-1.2475 + 0.2475
-0.5050 - 0.2475

∇ · J
1.7525
-1.0

-0.7525

∆t

Qn+1−Qn
+ ∇ · J
0.0

-4.440892e-16

0.0

Table B.2: Current using 1 quadrature point and analytic integral for linear path.

Edge 1
Edge 2
Edge 3

Analytic
1.2475
0.5050
0.2475

1 pt Quad

1.2475
0.5050
0.2475

Error

-4.4409e-16

0.0
0.0

B.1.3 Non-linear Path

Þ tn+1
Þ tn+1

tn

∆t

The non-linear path if integrated by hand can be shown to conserve charge.
[−2t3 + 3t2] + [1 + 2t3].
[3t2 + 1].

0 (cid:3) −[1 − tn+1 − (tn+1)3] − [1 − tn − (tn)3]
0 (cid:3) −[1 − tn+1 − (tn+1)3] − [1 − tn − (tn)3]
0 (cid:3) −[1 − tn+1 − (tn+1)3] − [1 − tn − (tn)3]
[1 − tn+1 − (tn+1)3] − [1 − tn − (tn)3]

1
∆t
1
∆t
[t3 + t]|tn+1
1
∆t
[(tn+1)3 + tn+1 − (tn)3 − tn].
1
∆t

0 (cid:3)

∆t

∆t

∆t

tn

tn

+

+

+

+

.

0 (cid:3) 0

117

(B.36)

(B.37)

(B.38)

(B.39)

(B.40)

Table B.3: Non-linear path analytic current, and for 1 and 2 point quadrature rule.

Edge 1
Edge 2
Edge 3

Analytic
1.25254
0.5000
0.2525

1.25
0.5
0.25

1 pt Quad 1 pt error
-0.0025

2 pt Quad

1.2525

0.5

0.2525

0.0

0.0025

2 pt error

-6.661338147750939e-16
2.220446049250313e-16

Analytical curved path integrals for current
[t3 ˆx + (1 − t) ˆy] · [ˆx + 3t2 ˆy] (cid:3)

j1 (cid:3)

1
∆t

j2 (cid:3)

j3 (cid:3)

1
∆t
1
∆t

[(1 − t3)ˆx + t ˆy] · [ˆx + 3t2 ˆy] (cid:3)
[−t3 ˆx + t ˆy] · [ˆx + 3t2 ˆy] (cid:3)

1
∆t

tn

Þ tn+1
Þ tn+1
Þ tn+1

tn

tn

tn

1
∆t

Þ tn+1
Þ tn+1
Þ tn+1

1
∆t

tn

−2t3 + 3t2dt (cid:3)

(cid:20)
(cid:21) tn+1

(cid:20)1

1 + 2t3dt (cid:3)

t +

2t3dt (cid:3)

2 t4

tn

tn

0.0

(cid:20) − 1
(cid:21) tn+1

1
2 t4

tn

2 t4 + t3

(cid:21) tn+1

tn
(B.41)

(B.42)

(B.43)

(B.44)

B.1.4 Non-Linear Path - Cell Crossing

Non-linear cell crossing follows the same rules as the linear cell crossing where the path
needs to be broken into two segements, one in each cell. For this setup an example with a
particle moving on a cubic path (same as example above) only tn (cid:3) 0.6 and tn+1 (cid:3) 0.75

Figure B.3: Basis functions N1,2(r), N1,3(r), N2,3(r) for cell 2 with cubic path of the particle

118

Table B.4: Current Values for particle crossing cells.

Current
Cell 1
Cell 2

Edge 1

0.0

Edge 2

0.0

1.5293264723380642 0.7056414418355151 0.5293264723380637

Edge 3

Edge 4

0.0

Edge 5

0.0

1.2173492419884955 1.7362675048664082 0.7362675048664085

Table B.5: Single quadrature point integration current.

Node 1

0.0
0.0

0.1840

Node 2
0.578125

0.6000

0.68232780382801

Node 3

0.25

0.2160

0.3176721961719803

-2.234967914173579

1.0

1.234967914173579

Node 4
0.171875

0.0
0.0
0.0

-1.2266666666666668
2.234967914173579

2.234967914173579
1.0083012475069124

0.0

0

0.0

-1.5398170847477821
-0.14583333333333318
-1.0000000000000004
1.5398170847477821
0.5398170847477817
0.3939837514144485

-0.23496791417357876
-0.00830124750691194
-4.440892098500626e-16 2.220446049250313e-16

0.22666666666666682
-1.2349679141735788

-1.0

1.0

2.539817084747782
1.1458333333333333

0.0

-2.5398170847477823
-2.5398170847477823
-1.393983751414449

-4.440892098500626e-16

∆t1

∆t2

Qn+1
Qn+α
Qn
Qn+α−Qn
Qn+1−Qn+α
Qn+1−Qn
Cell 1 ∇ · J
Cell 2 ∇ · J
∇ · J
Qn+1−Qn
+ ∇ · J
+ Qn+α−Qn

∆t

∆t

∆t1

Qn+1−Qn+α

∆t2

+ ∇ · J

B.2 Analytic E-Field for Adiabatic Expanding Plasmas

The density of the ions and electrons follow a Gaussian distribution shown as,

−r2
2σ(t)

η(t, r) (cid:3) η0

T0(t)(cid:3)3/2e
(cid:2)Tα(t)
σ2(t) (cid:3) σ2(0)(cid:18)
(cid:19)
(cid:20) Te(t)
(cid:20) Tion(t)
(cid:21)3/2

t2
τ2
exp

1 +

−r2
2σ2(t) − η0

e

Tion(0)

Te(0)

ρ(r) (cid:3)(cid:2)η0

(cid:21)3/2

2σ2(t)(cid:3)

−r2

e

(B.45)

(B.46)

(B.47)

We can find the E-field by taking Gauss’s law with the density of the ions and electrons.

∇ · E (cid:3) ρ/


Þ

ˆndsE (cid:3)

s
4πr2E (cid:3)

Þ
Þ

v

v

ρ

 dv
ρ

 dv

(B.48)

(B.49)

119

Er(r) (cid:3)

(cid:3)

(cid:3)

(cid:3)

(cid:3)

Using the identity,

4πr2

4πr2

1

1

1

0

0

0

0

0

0


 dv

Þ ρ(r(cid:48))
Þ π
Þ 2π
Þ r
Þ 2π
Þ π
Þ r
Þ r
(cid:20) Tion(t)
Þ r
(cid:2)η0
Þ b

Tion(0)

ρ(r

(cid:48))r

−αr2r2dr (cid:3)

e

0

4πr2
2π

4πr2

1
2


0

0

We can finally write Er(r) as

Er(r) (cid:3)(cid:2)η0

α (cid:3)

(cid:20) Tion(t)

Tion(0)

1
2σ2(t)

(cid:21)3/2√

π erf(√
4α3/2

αb)

− be−αb2
2α

ρ(r(cid:48))

 dv
ρ(r(cid:48))



r

(cid:48)2dr

(cid:48)

(cid:21)3/2

(cid:48)2sin(φ)dr

(cid:48)

dθdφ

(cid:21)3/2

(cid:20) Te(t)

Te(0)

2σ2(t)(cid:3) dr

−r(cid:48)2

(cid:48)

e

−r(cid:48)2
2σ2(t) − η0

e

√

π erf(√
4α3/2

αb)

− be−αb2
2α

(B.50)

(B.51)

(B.52)

(B.53)

(B.54)

(B.55)

120

BIBLIOGRAPHY

121

BIBLIOGRAPHY

[1]

[2]

[3]

[4]

[5]

Charles K Birdsall and A Bruce Langdon. Plasma physics via computer simulation. CRC
press, 2018.
José A Bittencourt. Fundamentals of plasma physics. Springer Science & Business
Media, 2013.
Collin S Meierbachtol, Andrew D Greenwood, John P Verboncoeur, and Balasub-
ramaniam Shanker. Conformal electromagnetic particle in cell: A review. IEEE
Transactions on Plasma Science, 43(11):3778–3793, 2015.
John P Verboncoeur. Particle simulation of plasmas: review and advances. Plasma
Physics and Controlled Fusion, 47(5A):A231, 2005.
Zane Crawford, Jie Li, Andrew Christlieb, and Balasubramaniam Shanker. Uncon-
ditionally stable time stepping method for mixed finite element maxwell solvers.
Progress In Electromagnetics Research, 103:17–30, 2020.

[7]

[8]

[9]

[6] Martin Campos Pinto, Sébastien Jund, Stéphanie Salmon, and Eric Sonnendrücker.
Charge-conserving fem–pic schemes on general grids. Comptes Rendus Mecanique,
342(10-11):570–582, 2014.
Francis F Chen. Introduction to plasma physics. Springer Science & Business Media,
2012.
Roger W Hockney and James W Eastwood. Computer simulation using particles. crc
Press, 1988.
Kane Yee. Numerical solution of initial boundary value problems involving maxwell’s
equations in isotropic media. IEEE Transactions on antennas and propagation, 14(3):302–
307, 1966.
Ji Qiang, Robert D Ryne, Salman Habib, and Viktor Decyk. An object-oriented
parallel particle-in-cell code for beam dynamics simulation in linear accelerators.
Journal of Computational Physics, 163(2):434–451, 2000.
Julien Derouillat, Arnaud Beck, F Pérez, Tommaso Vinci, M Chiaramello, A Grassi,
M Flé, G Bouchard, I Plotnikov, Nicolas Aunai, et al. Smilei: A collaborative, open-
source, multi-purpose particle-in-cell code for plasma simulation. Computer Physics
Communications, 222:351–373, 2018.

[10]

[11]

[12] S Dey, R Mittra, and S Chebolu. A technique for implementing the fdtd algorithm on

a nonorthogonal grid. Microwave and optical technology letters, 14(4):213–215, 1997.
John P Verboncoeur, A Bruce Langdon, and NT Gladd. An object-oriented electro-
magnetic pic code. Computer Physics Communications, 87(1-2):199–211, 1995.

[13]

122

[14]
Jian-Ming Jin. The finite element method in electromagnetics. John Wiley & Sons, 2015.
[15] Haksu Moon, Fernando L Teixeira, and Yuri A Omelchenko. Exact charge-conserving
scatter–gather algorithm for particle-in-cell simulations on unstructured grids: A
geometric perspective. Computer Physics Communications, 194:43–53, 2015.

[16] Fabio Riva, Carrie F Beadle, and Paolo Ricci. A methodology for the rigorous

verification of particle-in-cell simulations. Physics of Plasmas, 24(5):055703, 2017.

[17] Miles M Turner, Aranka Derzsi, Zoltan Donko, Denis Eremin, Sean J Kelly, Trevor
Lafleur, and Thomas Mussenbrock. Simulation benchmarks for low-pressure plasmas:
Capacitive discharges. Physics of Plasmas, 20(1):013507, 2013.

[18] Miles M Turner. Verification of particle-in-cell simulations with monte carlo collisions.

Plasma Sources Science and Technology, 25(5):054007, 2016.

[19] Patrick Kilian, Patricio A Muñoz, Cedric Schreiner, and Felix Spanier. Plasma waves

as a benchmark problem. arXiv preprint arXiv:1611.01127, 2016.

[20] AD Greenwood, JF Hammond, P Zhang, and YY Lau. On relativistic space charge
limited current in planar, cylindrical, and spherical diodes. Physics of Plasmas,
23(7):072101, 2016.

[21] Michel Moisan and Jacques Pelletier. Individual motion of a charged particle in
electric and magnetic fields. In Physics of Collisional Plasmas, pages 101–202. Springer,
2012.

[22] Martin Reiser and Patrick O’Shea. Theory and design of charged particle beams, volume

312. Wiley Online Library, 1994.

[23] Robert J Barker and Edl Schamiloglu. High-power microwave sources and technologies.

Wiley-IEEE Press, 2001.

[24] VF Kovalev and V Yu Bychenkov. Analytic solutions to the vlasov equations for

expanding plasmas. Physical review letters, 90(18):185004, 2003.

[25] S Laha, P Gupta, CE Simien, H Gao, J Castro, T Pohl, and TC Killian. Experimental
realization of an exact solution to the vlasov equations for an expanding plasma.
Physical review letters, 99(15):155001, 2007.

[26] Richard Marchand. Ptetra, a tool to simulate low orbit satellite–plasma interaction.

IEEE Transactions on Plasma Science, 40(2):217–229, 2011.

[27] RW Lemke, TC Genoni, and TA Spencer. Three-dimensional particle-in-cell simula-

tion study of a relativistic magnetron. Physics of Plasmas, 6(2):603–613, 1999.

[28] E Fourkal, B Shahine, M Ding, JS Li, T Tajima, and C-M Ma. Particle in cell
simulation of laser-accelerated proton beams for radiation therapy. Medical Physics,
29(12):2788–2798, 2002.

123

[29] Charles K Birdsall and A Bruce Langdon. Plasma physics via computer simulation. CRC

press, 2004.

[30] Chet Nieter, John R Cary, Gregory R Werner, David N Smithe, and Peter H Stoltz.
Application of dey–mittra conformal boundary algorithm to 3d electromagnetic
modeling. Journal of Computational Physics, 228(21):7902–7916, 2009.
Jonathan Squire, Hong Qin, and William M Tang. Geometric integration of the
vlasov-maxwell system with a variational particle-in-cell scheme. Physics of Plasmas,
19(8):084501, 2012.

[31]

[32] Peter Monk. Finite element methods for Maxwell’s equations. Oxford University Press,

2003.

[33] Alexander S Glasser and Hong Qin. The geometric theory of charge conservation in

particle-in-cell simulations. arXiv preprint arXiv:1910.12395, 2019.

[34] Haksu Moon, Fernando L Teixeira, and Yuri A Omelchenko. Exact charge-conserving
scatter–gather algorithm for particle-in-cell simulations on unstructured grids: A
geometric perspective. Computer Physics Communications, 194:43–53, 2015.

[35] C-D Munz, Pascal Omnes, Rudolf Schneider, Eric Sonnendrücker, and Ursula Voss.
Divergence correction techniques for maxwell solvers based on a hyperbolic model.
Journal of Computational Physics, 161(2):484–511, 2000.

[36] Zane D Crawford, Scott O’Connor, John Luginsland, and B Shanker. Rubrics for
charge conserving current mapping in finite element particle in cell methods. arXiv
preprint arXiv:2101.12128, 2021.

[37] Guangye Chen, Luis Chacón, and Daniel C Barnes. An energy-and charge-conserving,
Journal of Computational Physics,

implicit, electrostatic particle-in-cell algorithm.
230(18):7018–7036, 2011.

[38] Olgierd Cecil Zienkiewicz. A new look at the newmark, houbolt and other time
stepping formulas. a weighted residual approach. Earthquake Engineering & Structural
Dynamics, 5(4):413–418, 1977.
Jay P Boris. Relativistic plasma simulation-optimization of a hybrid code. In Proc.
4th Conf. Num. Sim. Plasmas, pages 3–67, 1970.

[39]

[40] S O’Connor, ZD Crawford, J Verboncoeur, J Lugisland, and B Shanker. A set
of benchmark tests for validation of 3d particle in cell methods. arXiv preprint
arXiv:2101.09299, 2021.

[41] David P Grote, Alex Friedman, Jean-Luc Vay, and Irving Haber. The warp code:
modeling high intensity ion beams. In AIP Conference Proceedings, volume 749, pages
55–58. American Institute of Physics, 2005.

[42] A Bossavit. for the student of numerical methods in electromagnetism. 1991.

124

[43] Alain Bossavit. Whitney forms: a class of finite elements for three-dimensional

computations in.

[44] Benedikt Perse, Katharina Kormann, and Eric Sonnenucker. Geometric particle-in-
cell simulations of the vlasov–maxwell system in curvilinear coordinates. SIAM
Journal on Scientific Computing, 43(1):B194–B218, 2021.

[45] Lee F Ricketson and Luis Chacón. An energy-conserving and asymptotic-preserving
charged-particle orbit implicit time integrator for arbitrary electromagnetic fields.
Journal of Computational Physics, 418:109639, 2020.

[46] Yang He, Yajuan Sun, Hong Qin, and Jian Liu. Hamiltonian particle-in-cell methods

for vlasov-maxwell equations. Physics of Plasmas, 23(9):092108, 2016.

[47] Eero Hirvijoki, Katharina Kormann, and Filippo Zonta. Subcycling of particle orbits
in variational, geometric electromagnetic particle-in-cell methods. Physics of Plasmas,
27(9):092506, 2020.

[48] Evstati G Evstatiev and Bradley A Shadwick. Variational formulation of particle
algorithms for kinetic plasma simulations. Journal of Computational Physics, 245:376–
398, 2013.
Jianyuan Xiao, Hong Qin, Philip J Morrison, Jian Liu, Zhi Yu, Ruili Zhang, and
Yang He. Explicit high-order noncanonical symplectic algorithms for ideal two-fluid
systems. Physics of Plasmas, 23(11):112107, 2016.

[49]

[50] Martin Campos Pinto, Katharina Kormann, and Eric Sonnendrücker. Variational
framework for structure-preserving electromagnetic particle-in-cell methods. arXiv
preprint arXiv:2101.09247, 2021.

[51] Philip J Morrison. Structure and structure-preserving algorithms for plasma physics.

Physics of Plasmas, 24(5):055502, 2017.

[52] XIAO Jianyuan, QIN Hong, and LIU Jian. Structure-preserving geometric particle-in-
cell methods for vlasov-maxwell systems. Plasma Science and Technology, 20(11):110501,
2018.

[53] Katharina Kormann and Eric Sonnendrücker. Energy-conserving time propaga-
tion for a structure-preserving particle-in-cell vlasov–maxwell solver. Journal of
Computational Physics, 425:109890, 2021.
Joshua W Burby. Finite-dimensional collisionless kinetic theory. Physics of Plasmas,
24(3):032101, 2017.

[54]

[55] Hong Qin, Jian Liu, Jianyuan Xiao, Ruili Zhang, Yang He, Yulei Wang, Yajuan Sun,
Joshua W Burby, Leland Ellison, and Yao Zhou. Canonical symplectic particle-in-cell
method for long-term large-scale simulations of the vlasov–maxwell equations.
Nuclear Fusion, 56(1):014001, 2015.

125

[56] Yang He, Hong Qin, Yajuan Sun, Jianyuan Xiao, Ruili Zhang, and Jian Liu. Hamilto-
nian time integrators for vlasov-maxwell equations. Physics of Plasmas, 22(12):124503,
2015.

[57] Bradley A Shadwick, Alexander B Stamm, and Evstati G Evstatiev. Variational
formulation of macro-particle plasma simulation algorithms. Physics of Plasmas,
21(5):055708, 2014.

[58] Nicolas Crouseilles, Lukas Einkemmer, and Erwan Faou. Hamiltonian splitting for
the vlasov–maxwell equations. Journal of Computational Physics, 283:224–240, 2015.
[59] Alexander B Stamm, Bradley A Shadwick, and Evstati G Evstatiev. Variational
formulation of macroparticle models for electromagnetic plasma simulations. IEEE
Transactions on Plasma Science, 42(6):1747–1758, 2014.
Jianyuan Xiao, Jian Liu, Hong Qin, and Zhi Yu. A variational multi-symplectic
particle-in-cell algorithm with smoothing functions for the vlasov-maxwell system.
Physics of Plasmas, 20(10):102517, 2013.

[60]

[61] Dong-Yeop Na, Yuri A Omelchenko, Haksu Moon, Ben-Hur V Borges, and Fernando L
Teixeira. Axisymmetric charge-conservative electromagnetic particle simulation
algorithm on unstructured grids: Application to microwave vacuum electronic
devices. Journal of Computational Physics, 346:295–317, 2017.

[62] D-Y Na, H Moon, YA Omelchenko, and FL Teixeira. Relativistic extension of a charge-
conservative finite element solver for time-dependent maxwell-vlasov equations.
Physics of Plasmas, 25(1):013109, 2018.

[63] Scott O’Connor, Zane D. Crawford, O. H. Ramachandran, John Luginsland, and
B. Shanker. Time Integrator Agnostic Charge Conserving Finite Element PIC. arXiv
e-prints, page arXiv:2102.06248, February 2021.

[64] Peter Robert Kotiuga. Hodge decompositions and computational electromagnetics. McGill

University Montreal, Canada, 1984.

[65] L Kettunen, K Forsman, and A Bossavit. Discrete spaces for div and curl-free fields.

IEEE Transactions on Magnetics, 34(5):2551–2554, 1998.

[66] Alain Bossavit. On the geometry of electromagnetism. J. Japan Soc. Appl. Electromagn.

& Mech, 6:17–28, 1998.

[67] Francesco P Andriulli. Loop-star and loop-tree decompositions: Analysis and
efficient algorithms. IEEE Transactions on Antennas and Propagation, 60(5):2347–2356,
2012.

[68] Francesco P Andriulli, Kristof Cools, Ignace Bogaert, and Eric Michielssen. On a
well-conditioned electric field integral operator for multiply connected geometries.
IEEE transactions on antennas and propagation, 61(4):2077–2087, 2012.

126

[69]

[70]

John B Manges and Zoltan J Cendes. A generalized tree-cotree gauge for magnetic
field computation. IEEE Transactions on Magnetics, 31(3):1342–1347, 1995.
JB Manges and ZJ Cendes. Tree–cotree decompositions for first-order complete tan-
gential vector finite elements. International journal for numerical methods in engineering,
40(9):1667–1685, 1997.

[71] Rui Wang, Douglas J Riley, and Jian-Ming Jin. Application of tree-cotree splitting
IEEE

to the time-domain finite-element analysis of electromagnetic problems.
transactions on antennas and propagation, 58(5):1590–1600, 2010.

[72] Saku Suuriniemi, Timo Tarhasaari, and Lauri Kettunen. Generalization of the

spanning-tree technique. IEEE Transactions on Magnetics, 38(2):525–528, 2002.

[73] Lauri Kettunen, Kimmo Forsman, and Alain Bossavit. Formulation of the eddy
current problem in multiply connected regions in terms of h. International journal for
numerical methods in engineering, 41(5):935–954, 1998.

[74] Man-Fai Wong, Odile Picon, and V Fouad Hanna. A finite element method based on
whitney forms to solve maxwell equations in the time domain. IEEE Transactions on
Magnetics, 31(3):1618–1621, 1995.

[75] Giuseppe Vecchi. Loop-star decomposition of basis functions in the discretization of

the efie. IEEE Transactions on Antennas and Propagation, 47(2):339–346, 1999.

[76] YY Lau and D Chernin. A review of the ac space-charge effect in electron–circuit

interactions. Physics of Fluids B: Plasma Physics, 4(11):3473–3497, 1992.

[77] Yong Wang, Dipanjan Gope, Vikram Jandhyala, and C-JR Shi. Generalized kirchoff’s
current and voltage law formulation for coupled circuit-electromagnetic simulation
with surface integral equations. IEEE transactions on microwave theory and techniques,
52(7):1673–1682, 2004.

[78] K. Aygun, B. C. Fischer, Jun Meng, B. Shanker, and E. Michielssen. A fast hybrid
field-circuit simulator for transient analysis of microwave circuits. IEEE Transactions
on Microwave Theory and Techniques, 52(2):573–583, 2004.

[79] Stephen A Maas. Nonlinear microwave and RF circuits. Artech House, 2003.
[80] Kuang-Ping Ma, Min Chen, Bijan Houshmand, Yongxi Qian, and Tatsuo Itoh. Global
time-domain full-wave analysis of microwave circuits involving highly nonlinear
phenomena and emc effects. IEEE transactions on microwave theory and techniques,
47(6):859–866, 1999.

[81] Wenquan Sui, Douglas A Christensen, and Carl H Durney. Extending the two-
dimensional fdtd method to hybrid electromagnetic systems with active and passive
lumped elements. IEEE Transactions on Microwave Theory and Techniques, 40(4):724–730,
1992.

127

[82] B. Houshmand Chien-Nan Kuo and T. Itoh. Full-wave analysis of packaged microwave
circuits with active and nonlinear devices: an fdtd approach. IEEE Transactions on
Microwave Theory and Techniques, 45(5):819–826, May 1997.

[83] Paolo Mezzanotte, Mauro Mongiardo, Luca Roselli, Roberto Sorrentino, and Wolfgang
IEEE

Heinrich. Analysis of packaged microwave integrated circuits by fdtd.
transactions on microwave theory and techniques, 42(9):1796–1801, 1994.

[84] Yui-Sheng Tsuei, AC Cangellaris, and JL Prince. Rigorous electromagnetic modeling
of chip-to-package (first-level) interconnections. IEEE transactions on components,
hybrids, and manufacturing technology, 16(8):876–883, 1993.

[85] Melinda Piket-May, Allen Taflove, and John Baron. Fd-td modeling of digital signal
IEEE Transactions on

propagation in 3-d circuits with passive and active loads.
Microwave Theory and Techniques, 42(8):1514–1523, 1994.

[86] Chien-Nan Kuo, Vincent A Thomas, Siou Teck Chew, Bijan Houshmand, and Tatsuo
Itoh. Small signal analysis of active circuits using fdtd algorithm. IEEE Microwave
and Guided Wave Letters, 5(7):216–218, 1995.

[87] B Toland, B Houshmand, and T Itoh. Modeling of nonlinear active regions with the

fdtd method. IEEE Microwave and Guided Wave Letters, 3(9):333–335, 1993.

[88] B Toland, J Lin, B Houshmand, and T Itoh. Fdtd analysis of an active antenna. IEEE

Microwave and Guided Wave Letters, 3(11):423–425, 1993.

[89] Karine Guillouard, Man-Fai Wong, V Fouad Hanna, and Jacques Citerne. A new
global finite element analysis of microwave circuits including lumped elements. IEEE
transactions on microwave theory and techniques, 44(12):2587–2594, 1996.

[90] Sung-Hsien Chang, Roberto Coccioli, Yongxi Qian, and Tatsuo Itoh. A global
finite-element time-domain analysis of active nonlinear microwave circuits. IEEE
Transactions on Microwave Theory and Techniques, 47(12):2410–2416, 1999.

[91] K Guillouard, MF Wong, V Fouad Hanna, and J Citerne. A new global time-domain
electromagnetic simulator of microwave circuits including lumped elements based
on finite-element method.
IEEE transactions on microwave theory and techniques,
47(10):2045–2049, 1999.

[92] Rui Wang and Jian-Ming Jin. A symmetric electromagnetic-circuit simulator based
on the extended time-domain finite element method. IEEE Transactions on Microwave
Theory and Techniques, 56(12):2875–2884, 2008.

[93] Qing He and Dan Jiao. Fast electromagnetics-based co-simulation of linear network
IEEE

and nonlinear circuits for the analysis of high-speed integrated circuits.
Transactions on Microwave Theory and Techniques, 58(12):3677–3687, 2010.

[94] Ali E Yilmaz, Jian-Ming Jin, and Eric Michielssen. A parallel fft accelerated transient
field-circuit simulator. IEEE transactions on microwave theory and techniques, 53(9):2851–
2865, 2005.

128

[95] Hakan Bagci, Ali E Yilmaz, Jian-Ming Jin, and Eric Michielssen. Fast and rigorous
analysis of emc/emi phenomena on electrically large and complex cable-loaded
structures. IEEE transactions on electromagnetic compatibility, 49(2):361–381, 2007.

[96] Ta-Ming Fang, Sun-Sheng Shei, RJ Nagem, and G v H Sandri. Convolution and
deconvolution with gaussian kernel. Il Nuovo Cimento B (1971-1996), 109(1):83–92,
1994.

[97] Andrew J Pray, Naveen V Nair, and B Shanker. Stability properties of the time domain
electric field integral equation using a separable approximation for the convolution
with the retarded potential. IEEE Transactions on Antennas and Propagation, 60(8):3772–
3781, 2012.

[98] Andrew J Pray, Yves Beghein, Naveen V Nair, Kristof Cools, Hakan Bağcı, and
Balasubramaniam Shanker. A higher order space-time galerkin scheme for time
domain integral equations. IEEE Transactions on Antennas and Propagation, 62(12):6183–
6191, 2014.

[99] Daniel S Weile, Greeshma Pisharody, Nan-Wei Chen, Balasubramaniam Shanker,
and Eric Michielssen. A novel scheme for the solution of the time-domain inte-
gral equations of electromagnetics. IEEE Transactions on Antennas and Propagation,
52(1):283–295, 2004.

[100] Yves Beghein, Kristof Cools, Hakan Bagci, and Daniël De Zutter. A space-time mixed
galerkin marching-on-in-time scheme for the time-domain combined field integral
equation. IEEE Transactions on Antennas and Propagation, 61(3):1228–1238, 2013.

[101] Chung-Wen Ho, A Ruehli, and Pierce Brennan. The modified nodal approach to

network analysis. IEEE Transactions on circuits and systems, 22(6):504–509, 1975.

[102] Balasubramaniam Shanker, Mingyu Lu, Jun Yuan, and Eric Michielssen. Time domain
integral equation analysis of scattering from composite bodies via exact evaluation of
radiation fields. IEEE Transactions on Antennas and Propagation, 57(5):1506–1520, 2009.
[103] Michel Arts Rob Maaskant. Reconsidering the voltage-gap source model used in
moment methods. IEEE Antennas and Propagation Magazine, 52(2):120–125, July 2010.
[104] Paule Charland, Ellen El-Khatib, and John Wolters. The use of deconvolution and
total least squares in recovering a radiation detector line spread function. Medical
physics, 25(2):152–160, 1998.

[105] B. Shanker, A.A. Ergin, Mingyu Lu, and E. Michielssen. Fast analysis of transient
electromagnetic scattering phenomena using the multilevel plane wave time domain
algorithm. Antennas and Propagation, IEEE Transactions on, 51(3):628–641, March 2003.
[106] Rowan J Gilmore and Michael B Steer. Nonlinear circuit analysis using the method
of harmonic balance—a review of the art. part i. introductory concepts. International
Journal of Microwave and Millimeter-Wave Computer-Aided Engineering, 1(1):22–37, 1991.

129

[107] Rowan J Gilmore and Michael B Steer. Nonlinear circuit analysis using the method
of harmonic balance—a review of the art. ii. advanced concepts. International Journal
of Microwave and Millimeter-Wave Computer-Aided Engineering, 1(2):159–180, 1991.

[108] Carlos E Christoffersen, Michael B Steer, and Mark A Summers. Harmonic balance
analysis for systems with circuit-field iterations. In Microwave Symposium Digest, 1998
IEEE MTT-S International, volume 2, pages 1131–1134. IEEE, 1998.

[109] WH Press, SA Teukolsky, WT Vetterling, and BP Flannery. Numerical Recipes in C++.

Cambridge University Press, 2002.

[110] Martin Campos Pinto, Marie Mounier, and Eric Sonnendrücker. Handling the diver-
gence constraints in maxwell and vlasov–maxwell simulations. Applied Mathematics
and Computation, 272:403–419, 2016.

[111] James W Eastwood. The virtual particle electromagnetic particle-mesh method.

Computer Physics Communications, 64(2):252–266, 1991.

[112] A Candel, A Kabel, L Lee, Z Li, C Limborg, C Ng, E Prudencio, G Schussman,
R Uplenchwar, and K Ko. Parallel finite element particle-in-cell code for simulations
of space-charge dominated beam-cavity interactions. In 2007 IEEE Particle Accelerator
Conference (PAC), pages 908–910. IEEE, 2007.

[113] Matthew R Gibbons and Dennis W Hewett. The darwin direct implicit particle-in-cell
(dadipic) method for simulation of low frequency plasma phenomena. Journal of
Computational Physics, 120(2):231–247, 1995.

[114] Johan Carlsson, Alexander Khrabrov, Igor Kaganovich, Timothy Sommerer, and
David Keating. Validation and benchmarking of two particle-in-cell codes for a glow
discharge. Plasma Sources Science and Technology, 26(1):014003, 2016.

[115] M Surendra. Radiofrequency discharge benchmark model comparison. Plasma

Sources Science and Technology, 4(1):56, 1995.

[116] Su Yan and Jian-Ming Jin. A continuity-preserving and divergence-cleaning algorithm
based on purely and damped hyperbolic maxwell equations in inhomogeneous media.
Journal of computational physics, 334:392–418, 2017.

[117] Martin Campos Pinto, Jakob Ameres, Katharina Kormann, and Eric Sonnendrücker.
On geometric fourier particle in cell methods. arXiv preprint arXiv:2102.02106, 2021.
[118] Neelakantam V Venkatarayalu and Jin-Fa Lee. Removal of spurious dc modes in
edge element solutions for modeling three-dimensional resonators. IEEE transactions
on microwave theory and techniques, 54(7):3019–3025, 2006.

[119] Yan-Lin Li, Sheng Sun, Qi I Dai, and Weng Cho Chew. Finite element implementation
of the generalized-lorenz gauged a-φ formulation for low-frequency circuit modeling.
IEEE Transactions on Antennas and Propagation, 64(10):4355–4364, 2016.

130

[120] Nathan M Newmark. A method of computation for structural dynamics. Journal of

the engineering mechanics division, 85(3):67–94, 1959.

[121] Balasubramaniam Shanker, A Arif Ergin, Kemal Aygun, and Eric Michielssen. Anal-
ysis of transient electromagnetic scattering from closed surfaces using a combined
field integral equation. IEEE Transactions on Antennas and Propagation, 48(7):1064–1074,
2000.

[122] Robert Kipp, Chi H Chan, Andrew T Yang, and JTW Yao. Simulation of high-
frequency integrated circuits incorporating full-wave analysis of microstrip dis-
continuities.
IEEE transactions on microwave theory and techniques, 41(5):848–854,
1993.

[123] Albert E Ruehli. Equivalent circuit models for three-dimensional multiconductor

systems. IEEE Transactions on Microwave theory and techniques, 22(3):216–221, 1974.

[124] Rui Wang and Jian-Ming Jin. Incorporation of multiport lumped networks into the
hybrid time-domain finite-element analysis. IEEE Transactions on Microwave Theory
and Techniques, 57(8):2030–2037, 2009.

[125] Rui Wang and Jian-Ming Jin. A flexible time-stepping scheme for hybrid field-
circuit simulation based on the extended time-domain finite element method. IEEE
Transactions on Advanced Packaging, 33(4):769–776, 2010.

[126] Eli Chiprout and Michel S Nakhla. Asymptotic waveform evaluation. In Asymptotic

Waveform Evaluation, pages 15–39. Springer, 1994.

[127] L. J. Jiang P. Li and H. Baˇgcı. Cosimulation of electromagnetics-circuit systems ex-
ploiting dgtd and mna. IEEE Transactions on Components, Packaging and Manufacturing
Technology, 4(6):1052–1061, June 2014.

[128] David M Sheen, Sami M Ali, Mohamed D Abouzahra, and Jin Au Kong. Application
of the three-dimensional finite-difference time-domain method to the analysis of
planar microstrip circuits. IEEE Transactions on Microwave Theory and Techniques,
38(7):849–857, 1990.

131