ROBUST MAXWELL SOLVERS FOR LARGE SCALE PARTICLE-IN-CELL SIMULATIONS By Zane Daniel Crawford A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Electrical Engineering – Doctor of Philosophy Computational Mathematics, Science and Engineering – Dual Major 2022 ABSTRACT ROBUST MAXWELL SOLVERS FOR LARGE SCALE PARTICLE-IN-CELL SIMULATIONS By Zane Daniel Crawford The design of modern devices is impacted heavily by the use and availability of robust, accurate, and efficient computational tools. This includes modeling devices that exploit plasma physics like particle accelerators, klystrons, ion thrusters, and micro-plasma generators among many other applications. While there are a number of current and emerging applications, the common thread between all is the need to accurately and efficiently capture all the relevant physics in geometrically intricate structures. The holy grail is to enable topology optimization to explore the design space. But all this requires rigorous translation from the continuous to the discrete world, while capturing all the underlying physics and not adding spurious artifacts due to discretization. A common computational model to perform this analysis is the particle-in-cell (PIC) method. It provides a straightforward paradigm to self-consistently solve for the distri- bution of the plasma as a collection of particles. The prevailing approach to solve for the fields in PIC is the finite difference time domain method (FDTD), or EM-FDTDPIC. But this effort leaves much to be desired, given the leaps that have been made in the finite element method; indeed, the latter is the method of choice for most commercial tools that that have become the de-facto workhorse in RF design industry. As a result, in the past decade, considerable effort has been expended in developing finite element (FEM) based PIC schemes, EM-FEMPIC. But we are still not there. One major concern of utilizing EM-FEMPIC over EM-FDTDPIC is the computational cost of FEM, which is greater than FDTD, despite the advantages of field and geometry accuracy FEM affords. This dissertation seeks to develop (i) a theoretically rigorous means to translate from the continuous to the discrete world while ensuring that there are no spurious artifacts, (ii) develops a higher order accurate method in both space and time, and (iii) overcomes cost complexity by introducing a linear scaling domain decomposition scheme. In all of these, the methods developed ensure that the necessary conservation properties are satisfied to machine precision. Numerous examples developed demonstrate these claims. Copyright by ZANE DANIEL CRAWFORD 2022 To my parents, Zane and Deborah, and all my family. To God be the glory. v ACKNOWLEDGMENTS The space in this acknowledgment is not sufficient to express the breadth nor the depth of my appreciation of all of those who helped me complete this work. Nevertheless, I must attempt to express my sincerest gratitude. First, I would like to thank my advisor, Dr. Shanker Balasubramaniam, for his encouragement and guidance. I appreciate all of the time and effort spent pushing me to become a better, more complete researcher. I also would like to thank the rest of my thesis committee, Dr. John Luginsland, Dr. Leo Kempel, Dr. John Verboncoeur, Dr. Edward Rothwell, Dr. Andrew Christlieb, and Dr. Daniel Dault. From the time spent at whiteboards going over my research to the helpful conversations that helped give me insight into the larger problems my work impacts, I enjoyed and appreciated every moment. I have been fortunate to have had many excellent colleagues. First, I thank the many members of the EMRG (Electromagnetics Research Group) with whom I had the pleasure of working alongside. I especially want to thank all of the members of the Computational Electromagnetic Laboratory members. A special thank you to those who worked alongside me on the finite element code, namely Dr. Scott O’Connor and Omkar Ramachandran. I am also grateful for the help from other CEML members, especially Luke Baumann, Abdel Alsnayyan, and Jacob Hawkins. I also would like to express my gratitude towards the Department of Electrical and Computer Engineering and the Department of Computational Science, Mathematics, and Engineering, whose students, faculty, and staff provided great support over the years. I would like to acknowledge the support that funded me through graduate school, but also introduced me to outstanding communities of researchers. First, the Department of Energy Computational Graduate Fellowship, which provided outstanding professional development and a peer group of amazing researchers. I also am grateful for the Engineering Workforce Development program, which helped guide me not only in career choices, but vi gave me more opportunities to present my work and receive valuable feedback. I also acknowledge the Sloan program at MSU, which included me in an important community of underrepresented engineers that allowed for excellent, candid conversations that kept me grounded throughout my graduate work. Finally, I would like to extend my thanks to the non-academic support system. I am so grateful for my many friends and well-wishers who encouraged me over the years. I am grateful for my family, the countless aunts, uncles, cousins, grandparents who cheered me on, waiting for this day. I thank my siblings, Aaron, Zaneta, and Angela Grace, for their love and support. And I am eternally thankful for my parents, Zane and Deborah, to whom words are not enough to convey my gratitude. vii In reference to IEEE copyrighted material which is used with permission in this thesis, the IEEE does not endorse any of Michigan State University’s products or ser- vices. Internal or personal use of this material is permitted. If interested in reprint- ing/republishing IEEE copyrighted material for advertising or promotional purposes or for creating new collective works for resale or redistribution, please go to http: //www.ieee.org/publications_standards/publications/rights/rights_link.html to learn how to obtain a License from RightsLink. viii TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii CHAPTER 1 INTRODUCTION AND BACKGROUND . . . . . . . . . . . . . . . 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.3 Particle-in-Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 Finite-Difference Time-Domain . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.5 Finite Element Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.6 Open Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.7 Goals and Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 CHAPTER 2 RUBRICS FOR CHARGE CONSERVING CURRENT MAPPING IN FINITE ELEMENT PARTICLE IN CELL METHODS . . . . . . . 13 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 Problem Statement and Formulation . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 Self consistent modeling framework . . . . . . . . . . . . . . . . . . . . . . . 18 2.3.1 Sampling the distribution function, and equation of continuity . . . . 18 2.3.2 Self-consistent solution of Maxwell’s Equations with active sources . 20 2.3.2.1 de-Rham complex . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.2.2 Discrete and self consistent Gauss’ and Ampere’s Laws . . . 22 2.4 Satisfaction of Conditions with Various EM-PIC Methods . . . . . . . . . . . 24 2.4.1 FDTD-PIC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.4.1.1 Satisfaction of Condition 1 . . . . . . . . . . . . . . . . . . . 25 2.4.1.2 Satisfaction of Condition 2 . . . . . . . . . . . . . . . . . . . 27 2.4.1.3 Satisfaction of Condition 3 . . . . . . . . . . . . . . . . . . . 28 2.4.2 EB-MFEM-PIC Formulation . . . . . . . . . . . . . . . . . . . . . . . . 28 2.4.2.1 Satisfaction of Condition 1 . . . . . . . . . . . . . . . . . . . 29 2.4.2.2 Satisfaction of Condition 2 . . . . . . . . . . . . . . . . . . . 29 2.4.2.3 Satisfaction of Condition 3 . . . . . . . . . . . . . . . . . . . 30 2.4.3 DH-MFEM-PIC Formulation . . . . . . . . . . . . . . . . . . . . . . . 31 2.4.3.1 Satisfaction of Condition 1 . . . . . . . . . . . . . . . . . . . 32 2.4.3.2 Satisfaction of Condition 2 . . . . . . . . . . . . . . . . . . . 32 2.4.3.3 Satisfaction of Condition 3 . . . . . . . . . . . . . . . . . . . 33 2.5 Consequences of Inconsistency . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.5.1 Violation of Condition 1 . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.5.2 Violation of Condition 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.5.3 Violation of Condition 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.5.4 Results for Consistent DH-MFEM-PIC . . . . . . . . . . . . . . . . . . 38 ix 2.6 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 CHAPTER 3 HIGHER ORDER CHARGE CONSERVING ELECTROMAG- NETIC FINITE ELEMENT PARTICLE IN CELL METHOD . . . . 42 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.3 Overview of Discrete Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.3.1 Discretization in Space . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3.2 Evolution of Particle Path and Current Mapping . . . . . . . . . . . . 49 3.3.3 Satisfaction of Gauss’s Magnetic Law . . . . . . . . . . . . . . . . . . . 49 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.4.1 Higher Order Cost . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4.2 Higher Order Particle Motion . . . . . . . . . . . . . . . . . . . . . . . 52 3.4.3 Plasma Ball . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.4.4 Expanding Particle Beam . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.4.5 Pafonsky Quadrupole . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 CHAPTER 4 DOMAIN DECOMPOSITION FRAMEWORK FOR MAXWELL FINITE ELEMENT SOLVERS AND APPLICATION TO PIC . . . . 61 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.2.1 Current Mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.2.2 Field discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.2.2.1 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . 65 4.2.2.2 Temporal Discretization . . . . . . . . . . . . . . . . . . . . . 66 4.3 Domain Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.3.1 MFEM-DP 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.3.2 MFETI-DP2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.4.1 Field Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.4.2 Particle Beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4.3 Plasma Ball . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 CHAPTER 5 CONCLUSION AND FUTURE WORK . . . . . . . . . . . . . . . . . 86 APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 APPENDIX A Supplementary Information for "Rubrics for Charge Con- serving Current Mapping in Finite Element Particle in Cell Methods" . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 APPENDIX B Higher Order Basis Definitions . . . . . . . . . . . . . . . . . . 94 APPENDIX C Unconditionally Stable Time Stepping Method for Mixed Finite Element Maxwell Solvers . . . . . . . . . . . . . . . . . . 97 x BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 xi LIST OF TABLES Table 2.1: First four rectangular cavity modes for consistent and inconsistent FEM solver in MHz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 Table 2.2: Continuity equation error for linear mapping . . . . . . . . . . . . . . . . 37 Table 2.3: Error of Analytic Integration vs Quadrature Rule for Linear Path . . . . . 37 Table 2.4: Continuity Equation Error for Linear Path . . . . . . . . . . . . . . . . . . 37 Table 2.5: Continuity Equation Error for Multiple Cells . . . . . . . . . . . . . . . . . 38 Table 2.6: Error of Analytic Integration vs Quadrature Rule for Quadratic Path . . . 38 Table 3.1: Satisfaction of Gauss’ Magnetic Law for Linear Field . . . . . . . . . . . . 51 Table 3.2: Satisfaction of Gauss’ Magnetic Law for Sinusoidal Field . . . . . . . . . . 51 Table 3.3: Plasma Ball expansion for different basis function orders and geometry radii 56 Table 3.4: Expanding Particle Beam Parameters . . . . . . . . . . . . . . . . . . . . . 56 Table 3.5: Error in impressed Quadrupole Fields . . . . . . . . . . . . . . . . . . . . . 57 Table 4.1: Average number of unknowns per volume and ratio of interior to boundary unknowns for different decompositions of Ω𝑟𝑒 𝑐𝑡 . . . . . . . . . 73 Table 4.2: Ratio of interior to boundary unknowns for different decompositions for particle beam test. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Table C.1: Leapfrog maximum time-step vs mesh edge lengths . . . . . . . . . . . . 104 xii LIST OF FIGURES Figure 1.1: Sample Volume Ω with particles . . . . . . . . . . . . . . . . . . . . . . . 2 Figure 2.1: Flow of typical EM-PIC simulation. © 2021 IEEE . . . . . . . . . . . . . . 18 Figure 2.2: Diagram of electromagnetic quantities as differential forms. © 2021 IEEE 21 Figure 2.3: Two dimensional diagram of where field and source quantities are defined on a simplical mesh. © 2021 IEEE . . . . . . . . . . . . . . . . . . 21 Figure 2.4: FDTD 2D grid. Solid cell is the primal grid. Dotted cells are the dual grid. © 2021 IEEE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Figure 2.5: Discrete curl and divergence on the dual grid. The curl is defined using the edges 𝑙ˆ𝑖 that come out of the page. The divergence is defined using the face normals 𝑛ˆ 𝑖 . © 2021 IEEE . . . . . . . . . . . . . . . . . . . 27 Figure 2.6: FDTD 2D particle current mapping with particles of finite size. © 2021 IEEE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Figure 2.7: FEM 2D curl and divergence mesh. The curl is defined using the edges 𝑙ˆ𝑖 around face center 𝑛ˆ 𝑖 . The divergence is defined with edges that point into node 𝑔𝑖 . © 2021 IEEE . . . . . . . . . . . . . . . . . . . . . . . . 30 Figure 2.8: FEM 2D particle current mapping with point particles. The particle current is mapped to each edge that encloses a cell along the particle path. © 2021 IEEE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Figure 2.9: Discrete curl mesh on the dual grid. The curl is defined with the outward pointing edge 𝑙ˆ𝑖 . © 2021 IEEE . . . . . . . . . . . . . . . . . . . . 32 Figure 2.10: FEM divergence mesh. The divergence is defined using face normals 𝑛ˆ 𝑖 that point out of cell 𝑐 𝑖 . © 2021 IEEE . . . . . . . . . . . . . . . . . . . 33 Figure 2.11: FEM particle mapping with finite shape particles. The particle is mapped to cell 𝑐 4 and 𝑐 3 at first and ends in cell 𝑐1 with current mapped at on the marked faces. The particle shape does not affect the conservation properties of this method. © 2021 IEEE . . . . . . . . . . . 34 Figure 2.12: Particle beam comparison of DH-MFEM-PIC and FDTD-PIC. © 2021 IEEE 39 xiii Figure 2.13: Error in the continuity equation (2.13a) and Gauss’ law (2.13b) for the particle beam case. © 2021 IEEE . . . . . . . . . . . . . . . . . . . . . . . . 39 Figure 2.14: Particle Beam for FDTD, 𝐸𝐵, and 𝐷𝐻 formulation © 2021 IEEE . . . . . . 40 Figure 2.15: Electron density of expanding adiabatic plasma ball using DH-MFEM- PIC compared to analytic expected densities. © 2021 IEEE . . . . . . . . 40 Figure 3.1: Relationship of fields and fluxes with respect to the de-Rham complex. 47 Figure 3.2: Relative Error in E𝑡,r for field propagating through free space. . . . . . . 53 Figure 3.3: Relative Error in E𝜌 for orbital motion . . . . . . . . . . . . . . . . . . . . 54 Figure 3.4: Relative Error in E𝑧 for orbital motion. . . . . . . . . . . . . . . . . . . . . 54 Figure 3.5: Relative Error in r for orbital motion. . . . . . . . . . . . . . . . . . . . . . 55 Figure 3.6: Relative Error of 𝜌(𝑡, ˜ r) vs −[∇]G(𝑡, r) for particle beam. . . . . . . . . . . 57 Figure 3.7: Relative Error of [∇·]D(𝑡, r) vs−[∇]G(𝑡, r) for particle beam. . . . . . . . . 58 Figure 3.8: Relative Error of [∇·]D(𝑡, r) vs−[∇]G(𝑡, r) for particle beam with dif- ferent time marching schemes. . . . . . . . . . . . . . . . . . . . . . . . . 58 Figure 3.9: Schematic for perfect electrical conducting box with Pafonsky quadrapole (shaded). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Figure 3.10: Particle beam traveling in 𝑥 − 𝑧 plane in Pafonsky quadrapole with background B 𝑦 (𝑡, r). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 Figure 3.11: Particle beam traveling in 𝑦 − 𝑧 plane in Pafonsky quadrapole with background B𝑥 (𝑡, r). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 Figure 4.1: Region Ω partitioned into 𝑁𝑣 subdomains. . . . . . . . . . . . . . . . . . 67 Figure 4.2: Condition number of system matrix. . . . . . . . . . . . . . . . . . . . . 75 Figure 4.3: Convergence of history of GMRES for monolithic, MFETI-DP1, and MFETI-DP2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 Figure 4.4: System matrix eigenvalues. . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Figure 4.5: Relative error in electric field for MFETI-DP. . . . . . . . . . . . . . . . . 78 xiv Figure 4.6: Timing data for MFETI-DP. . . . . . . . . . . . . . . . . . . . . . . . . . . 79 Figure 4.7: Error in Gauss’ law for particle beam. . . . . . . . . . . . . . . . . . . . . 81 Figure 4.8: Error in tracked particle positions. . . . . . . . . . . . . . . . . . . . . . . 82 Figure 4.9: Particle beam expansion for enlarged PEC cavity. . . . . . . . . . . . . . . 83 Figure 4.10: Adiabatic plasma ball expansion after 150, 300, and 800 ns. . . . . . . . . 84 Figure C.1: Sample Volume Ω with Different Boundary Conditions. Reproduced courtesy of The Electromagnetics Academy . . . . . . . . . . . . . . . . . 99 Figure C.2: Eigenvalues of A−11 A0 . Reproduced courtesy of The Electromagnetics Academy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 Figure C.3: Space Convergence. Reproduced courtesy of The Electromagnetics Academy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Figure C.4: Time Convergence. Reproduced courtesy of The Electromagnetics Academy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 Figure C.5: Runtime Comparison. Reproduced courtesy of The Electromagnetics Academy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 Figure C.6: Leapfrog Eigenvalues. Reproduced courtesy of The Electromagnetics Academy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 Figure C.7: Newmark-Beta Eigenvalues for Uniform Airbox. Reproduced courtesy of The Electromagnetics Academy . . . . . . . . . . . . . . . . . . . . . . 113 Figure C.8: Newmark-Beta Eigenvalues for Nonuniform Airbox. Reproduced courtesy of The Electromagnetics Academy . . . . . . . . . . . . . . . . . 114 Figure C.9: Newmark-Beta Eigenvalues for Uniform Airbox with First Order Basis functions. Reproduced courtesy of The Electromagnetics Academy . . . 115 Figure C.10: Newmark-Beta Eigenvalues for Uniform Airbox with Second Order Basis Functions. Reproduced courtesy of The Electromagnetics Academy116 Figure C.11: Plane Wave Simulation for 100,000 Time Steps. Reproduced courtesy of The Electromagnetics Academy . . . . . . . . . . . . . . . . . . . . . . 116 xv CHAPTER 1 INTRODUCTION AND BACKGROUND 1.1 Introduction Plasma is the most abundant state of matter in the universe, defined as ionized gas of equal parts positively and negatively charged particles, namely electrons and ions. Plasma is the bulk constituent of stars, solar winds, and natural phenomena such as lightning. It is also the key component in television screens, neon signs, and nuclear fusion. Understanding how plasma interacts with its environments informs the design of novel accelerator technology and beam wave interaction devices. The volume of plasma is approximately charge neutral, though regions may be more positively or negatively charged. This has some interesting repercussions. The charges in plasma can radiate when accelerating, causing self-interactions in the plasma, as well as be affected by other electromagnetic sources. Plasma generally is a good electrical and thermal conductor and even exhibit wave phenomena. The combination of all these features makes plasma very difficult to predict through simple models for the vast number of complicated uses and places it manifests. Thus, computational modeling is often used to study plasma and gain insight into how these behaviors manifest themselves in nature and man-made environments. Many computational tools have been developed in order to accelerate the design of new devices and better understand nature. This is not to say that modeling plasma is a simple task, even after decades of research and development of these tools. Modeling plasma requires one to solve a multiscale, multiphysics problem to simulate the behavior of many, many particles interacting with themselves and potentially other outside sources. This is an inherently nonlinear problem that has been an open area of research, with many paradigms to solve it. This can be simplified by modeling plasma as the interaction of currents and 1 the resulting electromagnetic fields. This would require a robust electromagnetic solver at a minimum to model plasma. It would also be wise to consider methods that can leverage high performance computing architecture, considering that modeling electromagnetic devices alone, without sources, can require vast amounts of computational resources. In order to accomplish this feat of modeling plasma, it is necessary to improve and create computational tools that specifically accommodate the physics one wishes to model. This work specifically seeks to improve the state of the art for numerically finding dynamic electromagnetic fields with charges. 1.2 Problem Statement We begin by simplifying the physics of the problem by following a kinetic description of the plasma. In this model, the representation of the particles in the system move as result of electromagnetic forces, which may come from external forces or interactions between particles. This work does not include the possibility of collisions between particles. Consider a region in space Ω which contains the plasma of interest that is bounded by 𝜕Ω with an outward pointing normal vector 𝑛ˆ that points into a vacuum Ω𝑉 . The plasma Figure 1.1: Sample Volume Ω with particles is represented with particles described with a phase space distribution function 𝑓 (𝑡, r, v) 2 (PSDF) that satisfies the Vlasov equation 𝜕𝑡 𝑓 (𝑡, r, v) + v · ∇ 𝑓 (𝑡, r, v) + a · ∇𝑣 𝑓 (𝑡, r, v) = 0, (1.1) where v is the particle velocity and the particle acceleration is defined by the Lorentz force 𝑞 a= (E(𝑡, r) + v × B(𝑡, r)) . (1.2) 𝑚 In general, the plasma can be composed of many species of particles, but for simplicity, only a single species is considered in this description. The particles are used to define electromagnetic sources, the charge density 𝜌(𝑡, r) and current density J(𝑡, r), which in terms of the PSDF are defined as ∫ 𝜌(𝑡, r) = 𝑞 𝑓 (𝑡, r, v)𝑑v (1.3a) Ω ∫ J(𝑡, r) = 𝑞 v(𝑡) 𝑓 (𝑡, r, v)𝑑v (1.3b) Ω and are subject to Newton’s equations of motion v(𝑡) = 𝜕𝑡 r(𝑡) (1.4a) a(𝑡) = 𝜕𝑡 v(𝑡). (1.4b) The electric field E(𝑡, r) and magnetic flux density B(𝑡, r) satisfy Maxwell’s equations 𝜕𝑡 B(𝑡, r) = ∇ × E(𝑡, r) (1.5a) 𝜕𝑡 D(𝑡, r) = ∇ × H(𝑡, r) − J(𝑡, r) (1.5b) ∇ · B(𝑡, r) = 0 (1.5c) ∇ · D(𝑡, r) = 𝜌(𝑡, r), (1.5d) where the electric flux density D(𝑡, r) and magnetic field H(𝑡, r) are related to the electric field and magnetic flux density through the the constitutive parameters D(𝑡, r) = 𝜀(𝑡, r) ∗𝑡 E(𝑡, r) (1.6a) 3 B(𝑡, r) = 𝜇(𝑡, r) ∗𝑡 H(𝑡, r) (1.6b) where ∗𝑡 is a convolution in time. The permitivity 𝜀(𝑡, r) = 𝜀0 𝜀𝑟 (𝑡, r) and permeability 𝜇(𝑡, r) = 𝜇0 𝜇𝑟 (𝑡, r). This work mostly has Ω composed of free space, implying 𝜀𝑟 (𝑡, r) = √ 𝜇𝑟 (𝑡, r) = 1, with 1/ 𝜇0 𝜀0 = 𝑐, 𝑐 being the speed of light. At the boundary 𝜕Ω, the fields obey the boundary conditions 𝑛ˆ × (EΩ𝑉 (𝑡, r) − EΩ (𝑡, r)) = 0 (1.7a) 𝑛ˆ × (HΩ𝑉 (𝑡, r) − HΩ (𝑡, r)) = J𝑠 (𝑡, r) (1.7b) 𝑛ˆ · (DΩ𝑉 (𝑡, r) − DΩ (𝑡, r)) = 𝜌 𝑠 (𝑡, r) (1.7c) 𝑛ˆ · (BΩ𝑉 (𝑡, r) − BΩ (𝑡, r)) = 0 (1.7d) where the subscript denotes whether it is within Ω or in Ω𝑉 . Additionally, it is necessary to also satisfy the continuity equation 𝜕𝑡 𝜌(𝑡, r) = −∇ · J(𝑡, r) (1.8) which can be derived from (1.5) or from (1.1). The fields which satisfy (1.5) do not have to be obtained directly from those equations. At times, the wave equation is used instead ∇ × 𝜇−1 −2 2 𝑟 (𝑡, r) ∗𝑡 ∇ × E(𝑡, r) + 𝑐 𝜕𝑡 𝜀𝑟 (𝑡, r) ∗𝑡 E(𝑡, r) = −𝜕𝑡 J(𝑡, r) (1.9) or ∇ × 𝜀𝑟−1 (𝑡, r) ∗𝑡 ∇ × H(𝑡, r) + 𝑐 −2 𝜕𝑡2 𝜇𝑟 (𝑡, r) ∗𝑡 H(𝑡, r) = ∇ × 𝜀𝑟−1 (𝑡, r) ∗𝑡 J(𝑡, r). (1.10) Either both equations can be used to obtain both the electric and magnetic fields, or use one and then use one of the curl equations, Faraday’s law (1.5a) or Ampere’s law (1.5b), to obtain the other field. Regardless of which equation or set of equations are used, the main difficulty remains: to choose a scheme which self-consistently yields the fields and particle positions and velocities that satisfies the Vlasov-Maxwell equation, Newton’s equations, Maxwell’s equations, and the appropriate boundary conditions. 4 1.2.1 Preliminaries The first step in solving the proposed the problem is to discretize problem such that it can be described numerically to be solved on computers. First, divide Ω in to a number of non-overlapping cells, also referred to as elements. The shape of the cells are typically hexahedra or tetrahedra in three dimensions, or rectangles or triangles in two dimensions. The collection of these cells are called a mesh or grid. In some methods, the hexahedral (rectangular) mesh will be made regular, such that it will be cubic (square), which is called a structured mesh. Otherwise the mesh is called unstructured. The trade-off is that structured grids allow one to take advantage of the orthogonality of the grid, but lose potential accuracy to staircasing errors as the volume of interest Ω is not usually a square or rectangular object. Regardless of the mesh style that is chosen, the mesh will be composed of nodes, edges, faces, and cells. Associated with some feature of the mesh, there is a basis set used to represent or reconstruct some quantity. The discretization of the equations will require approximation of the vector quantities as 𝑁𝛾 Õ X(𝑡, r) ≈ 𝑥 𝑖 (𝑡)s𝑖 (r) (1.11) 𝑖=1 where 𝑁𝛾 is number of nodes, edges, faces, or cells, and s(r) is some vector basis set. Similarly, scalar quantities are approximated as 𝑁𝛾 Õ 𝑌(𝑡, r) ≈ 𝑦 𝑖 (𝑡)𝑔𝑖 (r) (1.12) 𝑖=1 where 𝑔(r) is a scalar basis set. Often, a shorthand will be used to collect the coefficients of the basis sets, or degrees of freedom, into a single vector, e.g. x = [𝑥1 (𝑡), 𝑥2 (𝑡), · · · , 𝑥 𝑁𝛾 (𝑡)]. The choice of the basis sets depends on the quantity being discretized due to conditions on that quantity. For example, a field quantity like the electric field does not require normal continuity across an interface. In fact, in the presence of a surface charge or change in material, the electric field will be normally discontinuous. However, the tangential component of the electric field is always continuous (as we have discounted the possibility 5 of a fictitious magnetic monopole or current). Furthermore, not every equation defined in the previous section will be turned into an equation that will be solved. However, they still must be satisfied. This will lead to particular choices of basis sets that is explained in detail in Chapter 2. There are many approaches to to obtain the values of the degrees of freedom that best approximate the quantities of interest. The methods that will be presented later in this chapter and in this work, the finite-difference time domain (FDTD), the finite element method (FEM), and the mixed finite element method (MFEM), are all examples of the method of weighted residuals. Consider a spatial linear operator ℒ◦ that operates on a function X(𝑡, r) such that ℒ ◦ X(𝑡, r) = U(𝑡, r), ∀r ∈ Ω (1.13) subject to boundary condition ℒ 𝜕Ω ◦ X(𝑡, r) = W(𝑡, r), ∀r ∈ 𝜕Ω (1.14) with boundary operator ℒ 𝜕Ω ◦. The method of weighted residuals defines a residual function from (1.13) ℛ(𝑡, r) = ℒ ◦ X(𝑡, r) − U(𝑡, r) (1.15) which is has applied inner product, or tested, with a trial basis set s̃ 𝑗 (r) that has compact support, typically over one cell on the mesh defined in Ω, Ω𝑐 . Ideally, the residual would be equal to zero. However, only an approximation is available through the measurement of residual with a compact support basis set through an inner product that is set to zero ∫ 𝑑r̃ s̃(r̃) · ℛ(𝑡, r̃) = 0. (1.16) Ω𝑐 Provided the residual is defined using (1.11), a system of linear equations can be defined as 𝑁𝛾 Õ 𝑥 𝑖 (𝑡)hs̃𝑖 (r), ℒ ◦ s 𝑗 (r)i = hs̃𝑖 (r), U(𝑡, r)i (1.17) 𝑗=1 6 where h·, ·i denotes a standard inner product. This can be more simply written as a matrix vector product Ax = u (1.18) where A ∈ R𝑁𝛾 ×𝑁𝛾 with elements A𝑖𝑗 = hs̃𝑖 (r), ℒ ◦ s 𝑗 (r)i and u ∈ R𝑁𝛾 with elements u 𝑗 = hs̃ 𝑗 (r), U(𝑡, r)i. Often, the choice of s̃(r) is made to be the same as s(r), which would result in a Galerkin method. It is stressed again that s(r) itself should satsify the differential equations and boundary conditions [1]. It is also noted that (1.18) could have been derived through other approaches, such as Rayleigh-Ritz method [2] or through the use of calculus of variations [3]. In later sections and chapters, (1.5) and (1.7) will be used in place of (1.13) and (1.14) to derive FDTD, FEM, and MFEM. Section 2 in particular demonstrates how FDTD, FEM, and MFEM are all related through this framework. 1.3 Particle-in-Cell One of the most used methods to model plasma is particle-in-cell (PIC). It follows the kinetic approach to modeling plasma by tracking the motion of particles, as each particle has a position and velocity with which it is associated. Often, the particles are representative of groups of particles called macroparticles. Because the interest is most often in the bulk behaviour of the particles, reducing the number of particles modeled this way reduces computational burden without greatly changing the overall behaviour of the plasma. This representation of the plasma then needs to be incorporated with the discretization of the electromagnetic system, a process called mapping the plasma to the mesh. Again, a plasma is a collection of charged particles, and is interpreted as charge and current that is being included in the electromagnetic system. This provides a complete description of the discretized electromagnetic system that can be solved for the fields, which are self-consistent and can be used to accelerate the particles. The PIC method is a cycle that repeats this process of mapping the position and velocity of particles under the influence of the fields found through some field solver. The basic 7 steps of the PIC method are as follows 1. Map charge and current to the mesh 2. Solve for the fields including the macroparticles 3. Interpolate the fields at the macroparticle positions 4. Push the macroparticles to new positions. The motion of the macroparticles is determined by (1.4), typically discretized by staggered finite differences in time such that 𝑞 𝑛 v𝑛+ 2 = v𝑛− 2 + (E + v𝑛 × B𝑛 ) 1 1 (1.19a) 𝑚 r𝑛+1 = r𝑛 + Δ𝑡 v𝑛+ 2 . 1 (1.19b) The preceding equation is typically evolved using what is called the Boris push [4], which preserves the phase of the system while being an explicit method. Other explicit formulations tend to have errors in phase space and implicit methods are typically deemed too expensive to calculate given the number of particles that are in the system at any given time. Additionally, modifications have been made to the Boris push to be valid in the relativistic range, where particles move close to the speed of light [5]. The fields are obtained by any number of methods, depending on the relative physics. When the particles are traveling much slower than the wave speed, one can solve Poisson’s equation to yield an electrostatic PIC (ES-PIC) method [6, 7, 8, 9]. However, we focus on the more general case of electromagnetic PIC (EM-PIC) formulations, where the magnetic field (or flux) can be significant. Most commonly, FDTD is used as the field solver, a FDTD-PIC formulation [10]. This formulation not only has been used to codify EM-PIC by the likes of Birdsall, Langdon, Eastwood, and others [11, 12], but has continued to see extensive use [13, 14, 15, 16, 17]. This is due to some of the benefits of FDTD which are highlighted in the next section. However, there has been significant development of FEM-PIC which uses FEM as the field solver in EM-PIC. 8 1.4 Finite-Difference Time-Domain The earliest approach of FDM was the finite-difference time-domain (FDTD) or Yee- method, developed by Yee in 1966 [18]. In this method, finite differences are used to approximate the time derivatives and curl operators of (1.5a) and (1.5b). The degrees of freedom are the electric field and magnetic field defined on the edges of staggered cubical or rectangular mesh, often called a structured grid. The two equations are also staggered in time to make an explicit Leapfrog time marching scheme. The time step size is limited by a Courant–Friedrichs–Lewy condition determined by the size of the mesh cells. Unconditionally stable time marching schemes have been developed for FDTD, but most commonly focus on the Leapfrog method. The simplicity of FDTD makes it a powerful tool for analysis, implemented very simply with the recurrence formulations that define the method using a simple stencil. Its popularity over the decades has created a wealth of implementations, including on high performance architectures. It has been used to analyze and model many problems including electromagnetic scattering, antennas, photonics, and [19, 20, 21, 22]. It has also been used extensively to model kinetic plasmas, which will be discussed at length in Section 2. However, FDTD is not without its limitations. The structured grids used results in geometric errors in representing more complex geometry. This can be remedied somewhat by sub-cell and sub-grid formulations which increase the resolution of regions of the mesh, sacrificing stability and control of the accuracy [23]. Higher order spatial representations are typically limited to expanding the stencil and are not as accessible as other methods. Therefore, often the solution is to increase the resolution of the mesh by creating more cells. Furthermore, these approaches to improve geometric and field accuracy can exacerbate another limitation: the time step size. As stated, the time step size is limited by the smallest features in the mesh, making analyzing certain systems very time consuming from the sheer number of time steps needed to observe the desired phenomena. Some of these problems were remedied by the adaption of FEM for electromagnetics. 9 1.5 Finite Element Method The earliest works in FEM were outside of the electrical engineering discipline. It is credited to the appendix of a paper by Courant [24] in 1943, but it would not be until the late 1960s when Arlett, Bahrani, and Zienkiewicz [25] and Silvester [26] would apply the method to electromagnetics examining waveguides and cavities. Despite studying the time harmonic Helmholtz equation, the process for formulating the FEM system remains much the same. Here, Ω can be subdivided with triangles or tetrahedra in addition to rectangles/rectangular prisms where degrees of freedom and associated basis functions are defined. This generates a sparse, but typically not diagonal, system of equations that can be solved. Importantly, in the earliest works, scalar basis functions associated with the nodes of the mesh were used. However, when extended to the vector wave equation, spurious modes were found which were caused by not properly satisfying Gauss’ law and the associated boundary conditions [27]. In order to overcome this problem, many approaches were developed [28, 29, 30], but the most popular method is using vector basis functions [31, 32]. Importantly, these basis sets, first proposed by [33], only have tangential continuity at the element boundaries, which correctly satisfies the boundary condition for fields (1.7a) without violating the boundary condition for fluxes (1.7d). With this development, accurate results could be reported for not only the wave equation, but Maxwell equations as well. The introduction of vector basis functions for FEM allowed for the exploration of more problems in EM. FEM for the wave equation has been used to model many electromagnetic problems such as circuits, antennas, scattering, and other plasma-related phenomena [34, 35, 36, 37, 38]. In order to do large scale analysis, different approaches have been developed to avoid having to solve all the equations simultaneously. Some approaches like discontinuous Galerkin [39, 40, 41, 42] or element level dual field [43] approaches use flux or field continuity conditions to enable solving (1.9) on an element level. Other domain decomposition methods (DDM) [44, 45] are designed to solve the problem in larger 10 groupings of elements, and using various boundary conditions to ensure satisfaction of boundary conditions. One drawback of solving (1.9) and (1.10) specifically is that it does not implicitly constrain Gauss’ law. Therefore, some began solving (1.5a) and (1.5b), by introducing another basis function with normal continuity across cell faces, to create a mixed finite element method (MFEM) [46, 47, 48, 49]. This method is the primary focus of this thesis. 1.6 Open Problems The creation of more powerful EM-PIC methods requires advancing the underlying field solver technology. EM-PIC is by no means new and has been thoroughly explored using FDTD , as well has implementations in FEM [50], as well as MFEM [51, 52, 53]. However, examining each case reveals deficiencies. The shortcomings of FDTD, namely greater accuracy without greatly increasing spatial discretization, has driven the push toward using FEM and MFEM. At the time of writing, FEM-PIC approaches cannot using exact current mapping, instead relying on pseudocurrents to enforce the continuity equation [54]. MFEM implementations have not utilized unconditionally stable time marching, nor taken advantage of higher order spatial basis sets. Attempts to perform large scale analysis have been restricted to FDTD-PIC and discontinuous Galerkin FEM-PIC methods which still requires divergence cleaning to maintain charge conservation. The approach for large scale analysis using MFEM (and by extension, MFEM-PIC) relies on performing a sparse approximate inverse. However, as will be explored in Chapter 2, there are connections between FDTD, FEM, and MFEM that can be explored to make more robust field solvers, and by extension, better EM-PIC methods. This work uses improvements to MFEM as a means to define and formalize the guiding principles to define robust EM-PIC methods. Furthermore, the improvements and ideas presented here can be extended to modeling other sources beyond kinetic plasmas. 11 1.7 Goals and Outline In this dissertation, I present methods to reduce the cost of MFEM-based EM-PIC. The principle contributions of this work are as follows: 1. Apply various advances in FEM to MFEM, namely • unconditionally stable time marching, • higher order spatial basis sets, • domain decomposition, and 2. Apply these advances to modeling plasmas and space charge with charge conserving PIC schemes 3. Explore the uses and benefits of an underutilized variant of MFEM, DH-MFEM, particularly with respect to DH-MFEM-PIC. The outline of the dissertation is as follows. Chapter 2 presents a rubric to develop charge conserving current mapping schemes in FEM. It examines what properties are necessary to create a charge-conserving current mapping scheme for EM-PIC. Comparisons are made between FDTD and two formulations of MFEM. Chapter 3 presents a EM-MFEMPIC scheme using higher order spatial basis functions. The use of higher order basis functions allows for a reduction of degrees of freedom for a desired level of accuracy compared to previous lowest order descriptions. It also presents a scheme to satisfy Gauss’ magnetic law that is valid for higher order schemes. Chapter 4 presents two domain decomposition schemes for MFEM and applies them to EM-FEMPIC. The cost of solving large problems is reduced by solving a dual problem on fictitious, non-overlapping boundaries. This allows for a speed up in computational time while maintaining the required charge conservation properties. 12 CHAPTER 2 RUBRICS FOR CHARGE CONSERVING CURRENT MAPPING IN FINITE ELEMENT PARTICLE IN CELL METHODS 2.1 Introduction Modeling plasmas and other charged media as sources is a needed capability to better understand and design many real-world devices. Of the many techniques to model this phenomena, the particle-in-cell (PIC) method has been a popular approach. A wide variety of fields depend on PIC as a computational method, ranging from fundamental science of space weather, nuclear fusion, atomic processes on one hand, to high technology applications like particle accelerators, coherent electromagnetic radiation sources, and plasma processing devices on the other. PIC, at a glance works as follows; charged media, represented as a collection of particles, move in a domain under the influence of fields that are due to motion of these particles. This calls for self consistent solution of Maxwell’s equations in conjunction with equation of motion. Given application drivers stated earlier, there has been extensive development of EM-PIC methods [55, 56, 57, 58]. Several of these have been used extensively in the design of devices [59, 60, 61]. While traditional application were electrically large albeit geometrically simple, there is burgeoning interest in devices with significantly geometric complexity but smaller footprint. As a result, there is need for EM-PIC methods with greater capabilities [50]. But before we delve extensively into nuances and needs of modern EM-PIC codes, it helps to understand a bit of history of challenges that have been overcome in developing EM-PIC methods. The principal challenge in EM-PIC modeling is self-consistency of solving Maxwell’s equation together with equations of motion so as to actually capture dynamics of the distribution function. Challenges lies in incorporation of sources (and their variation over time) consistently into Maxwell’s equation. As a result, there has been 13 extensive work in understanding how one may develop methods that conserve charge [62, 63, 64]. Early development of EM-PIC methods largely relied on finite difference time domain (FDTD-PIC), as did addressing many of the aforementioned challenges [65, 66]. More recent PIC development has seen extensive development mixed finite element methods (MFEM-PIC) and discontinuous-Galerkin (DG-PIC) as the EM field solver [53, 67, 68]. These developments are specially pertinent as EM-PIC models that can more accurately capture geometry and physics with greater computational efficiency would soon be neces- sary to modeling physics ranging from novel accelerators to high frequency semiconductor devices. As the complexity of system to be modeled increases one needs EM-PIC suites that leverage advances in both finite element and particle methods. But as is to be expected, integrating the two calls for ensuring that self-consistency is preserved. This implies that discrete solutions of Maxwell’s equations together with Newton’s laws satisfy all of Maxwell’s equations together with the equation of continuity at every time step. This is easier said than done. Errors can manifest themselves as spurious charges and currents. To overcome these, it is possible to use procedures such as divergence cleaning to correct the erroneous fields that arise. Approaches to correct the fields include using Lagrange multipliers [54], solving Poisson’s equation at multiple time steps [4, 66], or changing the field solver [69]. However, there can be a significant computational cost to using these methods. Of course, a more desirable alternative is to prescribe a rubric in which charge is conserved as well as there is no space charge/current. Developing computational methods that are inherently charge conserving has been a topic of interest as long ago as there have been EM-PIC codes [70, 71]. The benefits of having a charge conserving scheme are apparent; they avoid additional work necessary to effect divergence cleaning. Unfortunately, their development has always been closely tied to a particular EM-solution technique, be in FDTD-PIC [65, 72] or MFEM-PIC [51, 52]. The latter is novel in that it develops a method that relies on de-Rham maps to define the 14 proper basis sets for representing the electric field and magnetic flux density. In effect, the primal quantities are defined based on Faraday’s laws. The current mapping schemes demonstrate how they recover the results found in the early literature [65, 72]. As methods developed have been solver specific, significant attention has not been paid to examining the underlying mathematical structure and why that is necessary for charge conservation. In this work, we seek to do just that. While unusual, we start our analysis at the particle approximation of the distribution function and develop a simple albeit rigorous approach to demonstrate charge continuity relationships between charge density and currents. The exposition in this paper straight- forward and uses the sifting property of the delta function. This is unlike that in presented in [72] that relies on an discrete approximation. As an added benefit, this approach lays the groundwork for developing conditions for mapping particles and fields on any discrete grid, and consistently evolving their position and strength. In effect, we create a uniform framework and criterion that is necessary to solve all of Maxwell’s equation consistently, thereby satisfying charge conservation as well. The rubrics prescribed are in agnostic to the method used to solve Maxwell’s equations. Indeed, we demonstrate how these are applicable to FDTD and (M)FEM. We also use a these rubrics to develop a novel MFEM-PIC scheme that uses electric flux density and the magnetic field (DH-MFEM-PIC) and validate its properties including charge conservation, satisfaction of Gauss’ law as well as comparison against quasi-analytic data. The rest of the paper is organized as follows: After defining the EM-PIC problem in Section 2.2 and its discetization using particle approaches, we define the three rules in Section 2.3 in the context of preserving the spatial and temporal connections between Gauss’ law and Ampere’s law. In Section 2.4, we provide examples on how FDTD and two formulations of a variant of FEM, the mixed finite element method (MFEM), satisfy the three conditions. One of these is novel in that it relies on the electric flux density and the magnetic field (DH-MFEM). Next, in Section 2.5, we provide results that show 15 what happens when the conditions are violated as well as affirm the correct behavior of MFEM when the conditions are met. In addition, several results are provided for the DH-MFEM-PIC that demonstrate necessary features for a number of examples (charge conservation, satisfaction of Gauss’ laws, correctness against quasi-nalaytical data). In an Appendix, we present details of the MFEM discretizations and other proofs. 2.2 Problem Statement and Formulation Consider a domain Ω whose boundaries are denoted by 𝜕Ω. It is assumed that domain comprises of charged species that exist in a background medium defined by 𝜀0 and 𝜇0 , the permittivity and permeability of free space; 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 evolution of the charged species is modeled via a fluid approximation and those of fields using Maxwell’s equations. To this end, we assume that the particle distribution function (PDF) 𝑓 (𝑡, r, v) must satisfy the Vlasov-Maxwell equation 𝑞 𝜕𝑡 𝑓 (𝑡, r, v) + v · ∇ 𝑓 (𝑡, r, v) + [E + v × B] · ∇𝑣 𝑓 (𝑡, r, v) = 0. (2.1) 𝑚 The electric field E(𝑡, r) and magnetic flux density B(𝑡, r) satisfy Maxwell’s equations − 𝜕𝑡 B(𝑡, r) = ∇ × E(𝑡, r) (2.2a) 𝜕𝑡 D(𝑡, r) = ∇ × H(𝑡, r) − J(𝑡, r) (2.2b) ∇ · B(𝑡, r) = 0 (2.2c) ∇ · D(𝑡, r) = 𝜌(𝑡, r), (2.2d) where fields and fluxes are related through constitutive parameters B(𝑡, r) = 𝜇0 H(𝑡, r) (2.3a) D(𝑡, r) = 𝜀0 E(𝑡, r). (2.3b) 16 The moments of the particle distribution function yield the the charge density 𝜌(𝑡, r) and current density J(𝑡, r), defined as ∫ 𝜌(𝑡, r) = 𝑞 𝑓 (𝑡, r, v)𝑑v (2.4a) Ω ∫ J(𝑡, r) = 𝑞 v(𝑡) 𝑓 (𝑡, r, v)𝑑v. (2.4b) Ω It follows from the above equations that charge and current densities are related via the equation of continuity − 𝜕𝑡 𝜌(𝑡, r) = ∇ · J(𝑡, r). (2.5) Further, assume that these equations satisfy the initial conditions arising from ∇ · D(0, r) = 𝜌(0, r) (2.6a) and ∇ · B(0, r) = 0 (2.6b) To complete the description of the problem, we will assume that the boundary of the domain 𝜕Ω may be partitioned into several subdomains, Γ𝑖 , such that Γ := ∪𝑖 Γ𝑖 . Each Γ𝑖 has a boundary condition, either Dirichlet (Γ𝐷 ), Neumann (Γ𝑁 ), or impedance (Γ𝐼 ). The boundary conditions are defined as 𝑛ˆ × E(𝑡, r) = 𝚿𝐷 (𝑡, r) on Γ𝐷 (2.7a) B(𝑡, r) 𝑛ˆ × = 𝚿𝑁 (𝑡, r) on Γ𝑁 (2.7b) 𝜇0 B(𝑡, r) 𝑛ˆ × − 𝑌 𝑛ˆ × 𝑛ˆ × E(𝑡, r) = 𝚿𝐼 (𝑡, r) on Γ𝐼 (2.7c) 𝜇0 where 𝑛ˆ the vector normal to the surface 𝜕Ω, 𝑌 is the surface admittance defined as 𝜀0 /𝜇0 , p and Ψ𝐷 (r, 𝑡), Ψ𝑁 (r, 𝑡), andΨ𝐼 (r, 𝑡) represent some function. These equations describe a self consistent set that must be solved numerically. In what follows, we describe the numerics necessary to do so. 17 Map Push Sources Particles to Mesh Interpolate Solve for Fields Fields Figure 2.1: Flow of typical EM-PIC simulation. © 2021 IEEE 2.3 Self consistent modeling framework While we do not solve the Vlasov system directly, we follow the usual practice of representing the particle distribution function using a collection of basis sets and then evolving the trajectory of these in using Newton’s law in concert with Maxwell equations. The basis sets are defined on a cell or simplical complex which discretizes a volume Ω bounded by 𝜕Ω. The mesh consists of 𝑁 𝑔 nodes, 𝑁𝑒 edges, 𝑁 𝑓 faces, and 𝑁𝑣 volumes (in three dimensions). This process of modeling the charged media can be described as the EM-PIC cycle as seen Fig. 2.1. It begins with mapping the charged particles on the grid, forming a discretized representation of 𝜌(𝑡, r) and J(𝑡, r). These discretized currents and charges act as sources in discretized Maxwell solvers, and are used to evolve the fields for one time step. The “new” fields and flux densities are then used update the locations of the charges and thereby obtain currents (and charges). The cycle then repeats. What is apparent from the the above are the the nuances of what ensures self consistent modeling of Maxwell and Vlasov systems discrete grids; this is the focus of this Section. 2.3.1 Sampling the distribution function, and equation of continuity To begin, we start with representing the particle distribution function using 𝑁𝑃 Õ 𝑓 (𝑡, r, v) = 𝑆 r − r𝑝 (𝑡) 𝛿 v − v𝑝 (𝑡)   (2.8) 𝑝=1 18 where shape functions, 𝑆(r), represent the spatial support of particles. It follows from the definitions in (2.4) that 𝜌(𝑡, r) and J(𝑡, r) can be represented as 𝑁𝑝 Õ 𝜌(𝑡, r) = 𝑞 𝑆(r − r𝑝 (𝑡)) (2.9a) 𝑝=1 𝑁𝑝 Õ J(𝑡, r) = 𝑞 v𝑝 (𝑡)𝑆(r − r𝑝 (𝑡)). (2.9b) 𝑝=1 Given the discrete representation of the distribution function, it can be shown that the equation of continuity is still holds in this discrete/particle setting. While this has been proven by Eastwood [72], the approach we take is straightforward, and exploits the sifting ∫𝑡 property of the delta function. Specifically, we use the definition r𝑝 (𝑡) = r0𝑝 (𝑡) + 0 𝑑𝜏v𝑝 (𝜏) and the identity  ∫ 𝑡  𝛿 r − r𝑝 (𝑡) = 𝛿 r − r0𝑝 (𝑡) ★𝑠 𝛿 r − 𝑑𝜏v𝑝 (𝑡)   (2.10) 0 Using this identity, it can be shown that (see Appendix) 𝜕𝑡 𝛿 r − r𝑝 (𝑡) = −∇ · v𝑝 (𝑡)𝛿 r − r𝑝 (𝑡)    (2.11) The equation of continuity follows via 𝑁𝑝 𝑁𝑝 Õ Õ −𝜕𝑡 𝜌(𝑡, r) = −𝑞 𝜕𝑡 𝑆(r − r𝑝 (𝑡)) = −𝑞 𝑆(r) ★ 𝜕𝑡 𝛿 r − r𝑝 (𝑡)  𝑝=1 𝑝=1 (2.12) 𝑁  𝑝  Õ   = ∇ ·  𝑞 v𝑝 (𝑡)𝑆(r − r𝑝 (𝑡)) = ∇ · J(𝑡, r).  𝑝=1    Of course, this is not the only identity that can be derived. As shown in the Appendix, another relationship pertinent to this paper that follows from the spectral representation of (2.10) is 𝑁𝑝 ∫ r𝑝 (𝑡) Õ 𝜌(𝑡, r) = 𝜌(0, r) − 𝑞∇ · 𝑑r̃𝑆 (r − r̃(𝑡)) (2.13) 𝑝=1 r𝑝 (0) As shown in the Appendix, this can be transformed into a line integral along the path taken by each particle. Our analysis, thus far, has dwelt on the representation of the distribution 19 function and its moments, viz., the charge and current. As is evident from (2.12) that the equation of continuity is satisfied when the distribution function is represented by arbitrarily shaped particles. Next, we address mapping/measuring these charges and currents on grid (and how they interact with Maxwell’s equations). 2.3.2 Self-consistent solution of Maxwell’s Equations with active sources For self consistency, one needs to solve Maxwell’s equations with sources defined by (2.9). In this pursuit, we discretize Maxwell’s equations in space with a collection of cells on which we define basis functions to represent fields and sources. Additionally, a temporal discretization is chosen such that a consistent solution for Newton’s laws is obtained to update the position and trajectory of the particle. It is well known that de-Rham diagrams [73, 74] provide appropriate map between fields, fluxes and sources in Maxwell’s equations. In addition, one needs additional framework to ensure that charges and currents are correctly measured and integrated in time. In what follows, we discuss each in turn. 2.3.2.1 de-Rham complex To begin, de-Rham diagrams establish relations between field, fluxes and sources. These diagrams replicate Maxwell’s equations in the discrete world; there exists extensive literature on their use in developing finite elements [75, 76]. In a nutshell, these permit the definition of Whitney spaces that can be used to represent primary quantities (𝜙, E, B) or their duals (𝜌, H, D). It should be noted that the choice of fields on the primal or dual grids are not rigid; in other words, (𝜌, H, D/J) can be defined on the primal grid and (𝜙, E, B) on the dual grid. These relationships are shown in Figs. 2.2 and 2.3. Using Whitney spaces to solve Maxwell’s equations has its antecedents in Bossavit’s seminal work [46] and those of several others [49, 38, 47, 77, 78]; our presentation is along these lines. As has been shown in these papers, these basis sets and their function spaces are defined as follows: the 0-form defined on nodes, the 1-form defined on edges, the 2-form defined on faces, and the 3-form 20 Figure 2.2: Diagram of electromagnetic quantities as differential forms. © 2021 IEEE ρ B E, J Figure 2.3: Two dimensional diagram of where field and source quantities are defined on a simplical mesh. © 2021 IEEE defined on volumes. In this work, we choose to use the Whitney basis sets that 0-form := 𝑊 (0) ∈ 𝐻 1 (2.14a) 1-form := W(1) (r) ∈ 𝐻(𝑐𝑢𝑟𝑙) (2.14b) 2-form := W(2) (r) ∈ 𝐻(𝑑𝑖𝑣) (2.14c) 3-form := 𝑊 (3) (r) ∈ 𝐿2 (2.14d) The Hodge operator relates the corresponding field and flux quantities using the correct constitutive parameters, as well as provides a mapping between the primal and dual grids [79, 80]. 21 2.3.2.2 Discrete and self consistent Gauss’ and Ampere’s Laws In typical solutions to Maxwell equations, there is an implicit assumption that the solution to the two curl equation enforces Gauss’ laws together with the equation of continuity. While this is true in the continuous domain, one needs to handle things with care in the discrete domain. Specifically, the basis functions used to represent or measure the fields, fluxes, currents, and charges need to be consistent with each other. Let ℒ∇· ◦ denote a discrete divergence operator and ℒ∇× ◦ denote a discrete curl operator, whose definitions are dependent on the basis functions. With no loss in generality, let’s denote the basis functions, N𝑖 (r) and 𝜉𝑖 (r), used to represent and measure the fields and fluxes, and charges, respectively. They belong to the appropriate space of functions depending on which quantity is being measured as determined by (2.14) and have local support over a simplex or cell. They also have relationships determined by the de-Rham complex shown in figure 2.2. With no loss of generality, let us define the vectors of quantities that are measured by our candidate basis functions. For instance, on can define 𝑑 𝑖 (𝑡) = hN𝑖 (r), 𝜀E(𝑡, r)i (2.15a) (ℒ∇× ◦ h(𝑡))𝑖 = hN𝑖 (r), ∇ × 𝜇−1 B(𝑡, r)i (2.15b) (ℒ∇· ◦ d(𝑡))𝑖 = h𝜉𝑖 (r), ∇ · 𝜀E(𝑡, r)i (2.15c) 𝑗𝑖 (𝑡) = hN𝑖 (r), J(𝑡, r)i (2.15d) Next, let us define d(𝑡) = [𝑑1 (𝑡), 𝑑2 (𝑡), . . . , 𝑑 𝑁𝛾 (𝑡)], h(𝑡) = [ℎ 1 (𝑡), ℎ2 (𝑡), . . . , ℎ 𝑁𝛼 (𝑡)], and j(𝑡) = [𝑗1 (𝑡), 𝑗2 (𝑡), . . . , 𝑗𝑁𝛽 (𝑡)], which can be used to approximate the continuous functions, Í 𝑁𝛾 such as D(𝑡, r) = 𝑖=1 𝑑 𝑖 (𝑡)N𝑖 (r), where 𝑁𝛾 is the number of degrees of freedom; as we will see, depending on the formulation it can either the number of edges or faces. This definition is the same for 𝑁𝛼 and 𝑁𝛽 . The degree of freedom vectors may require a Hodge operator to represent the function in the correct space. For solutions to be consistent, the systems of equations should be such that the divergence of Ampere’s law yields the 22 time derivative of Gauss’ law, which together with initial conditions will yield Gauss’ law. While this is simply stated, there are a number of subtle nuances. In the discrete world, it follows that the same should be valid. Specifically, the spatial basis sets and the time integrator used has to be such that these relationships are satisfied. Note, these conditions must be satisfied exactly because only the curl equations, Faraday’s law and Ampere’s law, are explicitly solved. Gauss’ law must be satisfied indirectly through the correct choice of spatial basis sets and time integration. As our goal is to ensure the both Gauss’ laws and the equation of continuity are satisfied, it follows that there has to be a relationship between ∇𝜉𝑖 (r) and N𝑖 (r). We begin by considering the divergence of discrete Ampere’s law; specifically, 𝜕𝑡 ℒ∇· ◦ d(𝑡) = ℒ∇· ◦ ℒ∇× ◦ h(𝑡) − ℒ∇· ◦ j(𝑡) (2.16) From the above, one can deduce the following: Condition 1: Consistent spatial basis sets The basis sets should be such that the discrete divergence of the discrete rotation of the magnetic field is zero. In effect, the discrete differential operators should maintain the same properties as their continuous counterpart. To serve as the link between Ampere’s law and the continuity equation, ℒ∇· ◦ ℒ∇× ◦ h(𝑡) = 0. (2.17) Though the discrete differential operators can be defined using properties of the mesh, they must be consistent with the discretized Faraday’s and Gauss’ law. Condition 2: Discrete divergence Let 𝜉(r) denote functions used to measure the charge on the mesh. Then ∇𝜉𝑖 (r) should span the irrorational component of both D(𝑡, r) and J(𝑡, r). In effect, span {∇𝜉𝑖 (r)} = D(𝑡, r) or J(𝑡, r) for ∀𝑡. In other words, this representation is complete. The need for this condition arises from measuring both 23 Gauss’ law and the equation of continuity. Consider measuring Gauss’ law via ∫ ∫ (ℒ∇· ◦ d(𝑡))𝑖 = 𝑑r𝜉𝑖 (r)𝜌(r) = 𝑑r𝜉𝑖 (r)∇ · D(𝑡, r) Ω Ω ∫ ∫ ∫ (2.18) 𝑑r𝜉𝑖 (r)𝜌(r) = − 𝑑r∇𝜉𝑖 (r) · D(𝑡, r) + 𝑑r𝜉𝑖 (r)𝑛ˆ · D(𝑡, r). Ω Ω 𝜕Ω The preceding equation mathematically exemplifies this condition. An almost identical equation can be derived for the equation of continuity. It can be shown boundary terms in vanish interior to the domain and provided the correct boundary condition is imposed on 𝑛ˆ · D(𝑡, r) on 𝜕Ω, Gauss’ law is implicitly satisfied together with condition 3. Condition 3: Discrete Time Integration Finally, discrete time time stepping used to up- date the electric flux density in Ampere’s equation and the path integral to used to evaluate the total current at any instant of time should be consistent with each other. This will ensure that the solution to Ampere’s law agrees with Gauss’ law. To gain more insight, we start by noting that when evolving Maxwell’s equations over-time, we are in effect integrating the current; however, the time integral of the current is obtained via a path integral in (2.13). It follows, that if one were to take a discrete divergence of Ampere’s laws, Condition 1 would yield 𝜕𝑡 ℒ∇· ◦ d(𝑡) = −ℒ∇· ◦ j(𝑡) (2.19) Provided Condition 2 is satisfied, it follows that choosing a consistent rule for integrating the charge along the particle path and integrating Ampere’s law will ensure conservation of charge. 2.4 Satisfaction of Conditions with Various EM-PIC Methods Next, we examine three different formulations of EM-PIC; FDTD-PIC, EB-MFEM-PIC, and DH-MFEM-PIC in light of the conditions discussed thus far. While FDTD-PIC is well known, EM-MFEM-PIC and DH-MFEM-PIC place dominant importance on Faraday’s and 24 Ampere’s laws, respectively. As a result, the dominant variables are represented on the primal grid, the others reside on dual grid. In what follows, we seek to show how each formulation satisfies the three conditions. 2.4.1 FDTD-PIC FDTD-PIC is the combination of the FDTD method [18] with PIC. It is discretized using staggered, structured grids in time and space, as seen in Figure 2.4 for two dimensions. For our discussion, the primal grid contains degrees of freedom for the electric field and magnetic flux density. Likewise, the dual grid has the degrees of freedom defined for the magnetic field and the electric flux density [73]. Particles are mapped to the grid using the 3-forms that measure the charge density. In what follows, we will show that the particle weighting scheme defined by Langdon [66] maps trivially onto the three conditions defined earlier. 2.4.1.1 Satisfaction of Condition 1 The discrete differential operators ℒ∇× ◦ and ℒ∇· ◦ on the dual grid are defined as  ±1 if 𝑛ˆ 𝑗 × 𝑙ˆ𝑖 points into (out of) the face     ℒ∇× ◦ = (2.20)  0  otherwise  and   ±1 if 𝑛ˆ 𝑗 points out of (into) cell 𝑐 𝑖    ℒ∇· ◦ = (2.21)  0  otherwise  The operator ℒ∇× ◦ ∈ R𝑁 𝑓 ×𝑁𝑒 , where 𝑁 𝑓 is the number of faces and 𝑁𝑒 is the number of edges. The operator ℒ∇· ◦ ∈ R𝑁𝑐 ×𝑁 𝑓 where 𝑁𝑐 is the number of volumes. Consider an example mesh in Figure 2.5. Assume the edge on the corners of the cell all point out of the page and the normal to the edges point out of the cell. The corresponding 25 Figure 2.4: FDTD 2D grid. Solid cell is the primal grid. Dotted cells are the dual grid. © 2021 IEEE curl matrix definition would be  1 −1 0 0      −1 0 1 0    [ℒ∇× ◦] =  , (2.22)   0 1 0 −1      0 0 −1 1    The divergence matrix would be h i [ℒ∇· ◦] = 1 1 1 1 (2.23) It is evident that in ℒ∇· ◦ ℒ∇× ◦, each entry in the resulting vector is zero, satisfying (2.17). 26 Figure 2.5: Discrete curl and divergence on the dual grid. The curl is defined using the edges 𝑙ˆ𝑖 that come out of the page. The divergence is defined using the face normals 𝑛ˆ 𝑖 . © 2021 IEEE 2.4.1.2 Satisfaction of Condition 2 The basis function 𝜉(r) in (2.18) is a 3-form, a constant in this formulation. Using the properties of 3-forms, (2.18) reduces to ∫ ∫ 𝑑r𝜉(r)∇ · X(r) = 𝑑r𝑛ˆ · X(r)𝜉(r). (2.24) Ω 𝜕Ω Consider the matrix form of ℒ 𝑑𝑖𝑣 ◦ based on Figure 2.5 1 0 1 1 0 1 0   ℒ∇· ◦ =  . (2.25)  0 1 0 −1 1 0 1   A volume integral would be a summation over the cells of a mesh, 𝑁𝑣 Õ h i [ℒ∇· ◦]𝑖𝑗 = 1 1 1 0 1 1 1 (2.26) 𝑖=1 This would effect the surface integral in (2.24). Therefore, the discrete continuity equation is satisfied exactly and consistency condition is fully satisfied. 27 Figure 2.6: FDTD 2D particle current mapping with particles of finite size. © 2021 IEEE 2.4.1.3 Satisfaction of Condition 3 Consider (2.24) in conjunction with (2.19). Using a leapfrog time marching scheme, the discrete continuity equation would be r𝑛+1 d𝑛+1 − d𝑛 𝑞 ∫ ∫ [ℒ∇· ] = − [ℒ∇· ] 𝑑r 𝑑r̃𝑛ˆ · p(r̃), (2.27) Δ𝑡 Δ𝑡 𝜕Ω r𝑛 where Δ𝑡 is the time step size, p(r̃) is the path of the particle. The integral computes the amount of charge that passes through a face. The divergence operator collects the contributions of each face, giving the net change of charge in each cell, as demonstrated in figure 2.6. This is equivalent to the current mapping in [65] for particles with finite shape. If the particle shape is infinitesimally small, it is equivalent to nearest-grid-point particle weighting. 2.4.2 EB-MFEM-PIC Formulation Next, we analyze EB-MFEM-PIC as defined in [46]. This finite element method can be seen as a generalization of FDTD. Like FDTD, it discretizes Maxwell’s equations using a primal-dual grid system, with the field and flux quantities defined in the same manner. However, the grid itself can be unstructured, using non-cubic cells. The formulation for 28 EB-MFEM can be found in appendix A.0.3.2. We have chosen to demonstrate the validity of the arguments in Section 2.3 using brick elements as it is akin to FDTD-PIC. We defer the discussion on tetrahedral elements to the next section. 2.4.2.1 Satisfaction of Condition 1 The discrete differential operators ℒ∇× ◦ and ℒ∇· ◦ are defined on the primal grid as  ±1 if 𝑛ˆ 𝑗 × 𝑙ˆ𝑖 points into (out of) the face     ℒ∇× ◦ = (2.28)  0  otherwise  and   ±1 if edge 𝑗 points into (out of) node 𝑖    ℒ∇· ◦ = (2.29)  0  otherwise  The operator ℒ∇× ◦ ∈ R𝑁𝑒 ×𝑁 𝑓 . The operator ℒ∇· ◦ ∈ R𝑁𝑛 ×𝑁𝑒 where 𝑁𝑛 is the number of nodes. Consider the example mesh in Figure 2.7. The corresponding curl matrix definition would be h i𝑇 ℒ∇× ◦ = −1 1 −1 1 . (2.30) The divergence matrix would be 1 1 0 0     −1 0 1 0    ℒ∇· ◦ =     (2.31)  0 −1 0 1       0 0 −1 −1   It is apparent that each entry in the resulting ℒ∇· ◦ ℒ∇× ◦ is zero, satisfying (2.17). 2.4.2.2 Satisfaction of Condition 2 The basis function for the current density spans the gradient of the basis function for the charge density. The divergence operator ℒ∇· ◦ sums the edges connected to a node to 29 Figure 2.7: FEM 2D curl and divergence mesh. The curl is defined using the edges 𝑙ˆ𝑖 around face center 𝑛ˆ 𝑖 . The divergence is defined with edges that point into node 𝑔𝑖 . © 2021 IEEE satisfy (2.18). For PEC boundary conditions, the current density normal to the face is zero, eliminating the surface integral term. This can be seen by summing the rows of (2.31). For non-PEC boundaries, the domain is padded with ghost cells where the boundary nodes are not included. The non-terminating edges then make up the surface integral term when ℒ∇· ◦ is applied. Therefore, (2.18) is satisfied. 2.4.2.3 Satisfaction of Condition 3 For each cell a particle visits, the corresponding current density is projected onto the edges of that cell. Using (2.18), the discretized continuity equation that satisfies (2.19) is r𝑛+1 d𝑛+1 − d𝑛 𝑞 ∫ ∫ [ℒ∇· ] = − [ℒ∇· ] 𝑑r 𝑑r̃∇𝑊 0 · p(r̃). (2.32) Δ𝑡 𝛿𝑡 Ω r𝑛 This is equivalent to particle-in-cell area weighting in FDTD and the method used in [72, 51]. Due to the relationship between the 0-form on the primal grid and the 3-form on the dual grid, this is also equivalent to the approach defined in [65] for finite particle shapes. 30 Figure 2.8: FEM 2D particle current mapping with point particles. The particle current is mapped to each edge that encloses a cell along the particle path. © 2021 IEEE 2.4.3 DH-MFEM-PIC Formulation The EB-MFEM-PIC formulation defined in 2.4.2 can be reformulated to such that the primal grid has the degrees of freedom associated with the magnetic field and electric flux density. The dual grid contains the degrees of freedom for the electric field and magnetic flux density, still defined using the edges and faces, respectively. The sources are represented on the primal grid using 3-forms for the charge density and 2-forms for the current density. The arguments for this formulation satisfying the three conditions are similar to that presented in section 2.4.1. Note, the principal ramification of this approach is that it is a natural framework for modeling sources as both the currents and electric flux density are modeled using the same function space, and the fields, currents and charges lie in the primal grid. This si unlike EM-MFEM-PIC formulation where one needs transformations between primal and dual spaces. The ramification of these mappings become more consequential in higher order schemes–a topic that will be visited in a later paper. 31 Figure 2.9: Discrete curl mesh on the dual grid. The curl is defined with the outward pointing edge 𝑙ˆ𝑖 . © 2021 IEEE 2.4.3.1 Satisfaction of Condition 1 The discrete differential operators ℒ∇× ◦ and ℒ∇· ◦ are defined on the primal grid, but have the same definitions as (2.20) and (2.21). Consider the example mesh in Figure 2.9. Assume the edge on the corners of the cell all point out of the page and the normal to the edges point out of the cell. The corresponding curl matrix definition for lowest order would be −1 1 0      [ℒ∇× ◦] =  1 0 −1 , (2.33)        0 −1 1    The lowest order divergence matrix would be h i [ℒ∇· ◦] = 1 1 1 (2.34) Analyzing ℒ∇· ◦ ℒ∇× ◦, each entry in the resulting vector is zero, satisfying (2.17). 2.4.3.2 Satisfaction of Condition 2 For lowest order, the continuity equation reduces to (2.24). Consider the matrix form of ℒ∇· ◦ based on figure 2.10 32 Figure 2.10: FEM divergence mesh. The divergence is defined using face normals 𝑛ˆ 𝑖 that point out of cell 𝑐 𝑖 . © 2021 IEEE 0 1 1 0 1   ℒ∇· ◦ =  . (2.35)  1 0 −1 1 0   A volume integral would be a summation over the cells of a mesh, 𝑁𝑣 Õ h i [ℒ∇· ◦]𝑖𝑗 = 1 1 0 1 1 (2.36) 𝑖=1 This would effect the surface integral in (2.24). Therefore, the discrete continuity equation is satisfied exactly and consistency condition is fully satisfied. 2.4.3.3 Satisfaction of Condition 3 Again, the satisfaction of condition three relies on the correct definition of current density on the faces of each cell. For consistency of the change of charge density in cells, the amount of charge that passes through the face must be measured exactly. This is equivalent to what is effectively done in [81, 82] through different techniques. Consider the mapping of current defined by the particle in Figure 2.11. Here, though the particle moves at an angle with respect to the face, the current density as seen by the mesh moves 33 Figure 2.11: FEM particle mapping with finite shape particles. The particle is mapped to cell 𝑐4 and 𝑐3 at first and ends in cell 𝑐 1 with current mapped at on the marked faces. The particle shape does not affect the conservation properties of this method. © 2021 IEEE in the direction of the face normal. The solver in effect sees a charge density move from cell 𝑐4 to 𝑐3 and the total charge density in 𝑐 3 moves to 𝑐1 by time step 𝑛 + 1. 2.5 Consequences of Inconsistency In the previous sections, the three rule framework was defined and applied to different formulations of EM-PIC. The references to previous works show how these rules encom- passes well known EM-PIC formulations. Next, we will examine ramifications of violating these three rules. By examining potential pitfalls and how it manifests in results, one can more easily diagnose errors in new EM-PIC schemes. 2.5.1 Violation of Condition 1 First, consider violations of the first condition. Failing to satisfy (2.17) can arise from choosing the wrong basis sets to discretize (2.2). One way to do so is to choose basis sets that do not obey the de-Rahm complex. For example, if one solves both Ampere’s law and Faraday’s law on the same grid, staggered in time. Like in FDTD, we consider representing the electric field and magnetic field; however, this time it is on the same grid. Assume the 34 fields are represented by 1-forms and test (2.2a) and (2.2b) with the same basis functions. The discretized system is h𝑛+ 2 = h𝑛+ 2 − Δ𝑡 [ℒ∇× ◦]e𝑛 1 1 0 (2.37a) e𝑛+1 = e𝑛+1 + Δ𝑡 ([ℒ∇× ◦]h𝑛+ 2 − j𝑛+ 2 ), 1 1 0 (2.37b) 0 where ℒ∇× ◦ is a modified curl operator, as the previously defined operator has the wrong dimensions. If this scheme is done with a structured mesh in FDTD, one would have a higher error as the curl would be represented using forward or backward difference rather than a central difference. Instead of the recurrence scheme found in (A.10b), one would have 𝑛+1 𝑛 Δ𝑡 𝑛+ 1 𝑛+ 1 𝐸 𝑥,(𝑖+ = 𝐸 𝑥,(𝑖+ + (𝐻 2 1 −𝐻 2 1 ) 1 2 ,𝑗,𝑘) 1 2 ,𝑗,𝑘) 𝜀Δ 𝑦 𝑧,(𝑖+ 2 ,𝑗+1,𝑘) 𝑧,(𝑖+ 2 ,𝑗,𝑘) Δ𝑡 𝑛+ 1 𝑛+ 1 − (𝐻 2 1 −𝐻 2 1 ) (2.38) 𝜀Δ𝑧 𝑦,(𝑖+ 2 ,𝑗,𝑘+1) 𝑦,(𝑖+ 2 ,𝑗,𝑘) Δ𝑡 𝑛+1/2 − 𝐽 𝜀 𝑥,(𝑖+ 12 ,𝑗,𝑘) if a forward difference was used. The effects of this are more clear from an FEM perspective. Consider trying to find modes of a perfect electric conductor rectangular cavity using the generalized eigenvalue problem formed by the time harmonic form of (2.39) ∗ 0 [ℒ∇× ◦] [I] 0           b b     =Λ    (2.39) ∗ 𝑇         −[ℒ∇× ◦] [★𝜇−1 ] 0  e  0 [★𝜀 ] e         where the eigenvalues in Λ = 𝑗𝜔. The dimensions of the cavity are 0.6 m×0.4 m×1 m. In ∗ a consistent scheme, [ℒ∇× ◦] = [ℒ∇× ◦] and [★𝜇−1 ] is as defined in (??). However, in the inconsistent scheme we have defined, ∫ ∗ (1) (1) [ℒ∇× ◦]𝑖,𝑗 = 𝑑ΩW𝑖 · ∇ × W 𝑗 (2.40) and [★𝜇−1 ] = [★𝜀 ]. When the inconsistent scheme is used, spurious modes are found that do not correspond to the true modes of the cavity as seen in table 2.1. Using this inconsistent scheme could permit incorrect solutions into the field solver, invalidating results. 35 Table 2.1: First four rectangular cavity modes for consistent and inconsistent FEM solver in MHz 𝑓true 𝑓consistent 𝑓inconsistent 291.346 290.829 0.280 390.242 389.157 1.036 403.608 402.697 1.841 450.382 449.654 2.451 2.5.2 Violation of Condition 2 Next, we consider the repercussions of violating Condition 2. Consider the linear mapping scheme typically found in early PIC literature where the weighting function is defined as Ö |𝑟𝑛,𝑔 − 𝑟𝑛,𝑝 | 𝑊 r − r𝑝 (𝑡) = .  1− (2.41) Δ𝑛 𝑛∈{𝑥,𝑦,𝑧} When the particle is mapped to the nodes as the charge density, 𝑟𝑛,𝑔 is the 𝑥, 𝑦, 𝑧− coordinate of the node. When mapped to the edges as the current density, 𝑟𝑛,𝑔 is the 𝑥, 𝑦, 𝑧− coordinate of the center of the edge. It is well known that this mapping does not satisfy the continuity equation. For example, consider the particle motion in Figure 2.8 where the particle moves from (.3, .4) to (.7, .6). Table 2.2 shows the error in the continuity equation using this mapping. The result is an inconsistency in the solved Ampere’s law and measured Gauss’ law due to spurious charge seen by the field solver as the next mapping is done. This mapping scheme does not appropriately satisfy (2.18), as the combination of basis functions for the current density at the node is not exactly equal to the gradient of the charge density basis function. Despite being in the same direction, the basis function defined by the linear mapping for the current density is different than the basis functions used to measure the fields. Therefore, the de-Rham complex is not truly satisfied and the spatial link between Gauss’ law and Ampere’s law is no longer preserved. 36 Table 2.2: Continuity equation error for linear mapping 𝜕𝑡 𝜌 ∇·J 𝜕𝜌 + ∇ · J node 1 -0.100 -0.050 -0.150 node 2 0.300 -0.850 -0.550 node 3 -0.300 0.850 0.550 node 4 0.100 0.050 0.150 Table 2.3: Error of Analytic Integration vs Quadrature Rule for Linear Path Analytic 1 pt Quad Error Edge 1 1.248 1.248 -4.441e-16 Edge 2 0.505 0.505 0.000 Edge 3 0.248 0.248 0.000 Table 2.4: Continuity Equation Error for Linear Path 𝑄 𝑛+1 −𝑄 𝑛 𝑄 𝑛+1 −𝑄 𝑛 𝑄 𝑛+1 𝑄𝑛 Δ𝑡 ∇·J Δ𝑡 +∇·J Node 1 0.284 0.459 -1.753 1.753 0.000 Node 2 0.550 0.450 1.000 -1.000 -4.440e-16 Node 3 0.166 0.091 0.753 -0.753 0.000 2.5.3 Violation of Condition 3 Next, we examine how failing to correctly integrate the spatial integral in (2.19) results in a violation of Gauss’ law and the continuity equation due to spurious charge on the grid. Consider the movement of a particle on a linear path that remains in a single cell using a 0-form basis set to measure the charge density. For a linear path, the path integral can be exactly computed either using an analytic scheme [52] or an one point integration rule, as seen in Table 2.3. This translates to machine precision error in the continuity equation in Table 2.4. Similarly, consider the case when a particle crosses at least one face on its path. When the path is in multiple cells, an error occurs when the path integral is not calculated in each cell, as seen in Table 2.5. This is an inconsistency in (2.19) as though the time integral may be computed correctly, parts of the path integral are missing. 37 Table 2.5: Continuity Equation Error for Multiple Cells 𝑄 𝑛+1 −𝑄 𝑛 𝑛+1/2 𝑄 𝑛+1 −𝑄 𝑛+𝛼 𝑄 𝑛+𝛼 −𝑄 𝑛 𝑛+1/2 Δ𝑡 + ∇ · j𝑛𝑐 Δ𝑡2 + Δ𝑡1 + ∇ · j𝑐𝑐 Node 1 1.008 0.000 Node 2 0.394 -4.441e-16 Node 3 -0.008 2.220e-16 Node 4 -1.394 -4.441e-16 Table 2.6: Error of Analytic Integration vs Quadrature Rule for Quadratic Path Analytic 1 pt Quad 1 pt error 2 pt Quad 2 pt error Edge 1 1.253 1.250 -0.003 1.253 -6.661e-16 Edge 2 0.500 0.500 0.000 0.500 2.220e-16 Edge 3 0.253 0.250 0.003 0.253 0.000 2.5.4 Results for Consistent DH-MFEM-PIC Finally, we examine results when all three conditions are satisfied. We have chosen to illustrate these using DH-MFEM-PIC as results for EB-MFEM-PIC are available [51, 52, ?]. First, we consider a particle beam in a 0.1 m PEC cylindrical cavity with a radius of 0.02 m. The particle beam has a voltage of 7.107 kV and current of 0.25A, simulated by inserting ten particles per time step. The starting velocity of the particles entering the tube is 5×104 km/s with a beam width of 0.008 m. In this test, when the three conditions are satisfied, the beam will spread as it progresses through the tube with a uniform distribution of particles. Using a DH-based EM-PIC formulation and point particles for particle shape, by following section 2.4.3, a self-consistent charge conserving EM-PIC scheme was created. Point particles are used because without introducing errors, it removes the need to compute the volume of a particle passing through a face, which can be extremely difficult in three dimensions, The error in the continuity equation and Gauss law is found in Figure 2.13, showing that satisfying conditions two and three is possible with a nearest grid point algorithm in FEM. The trend of the spreading beam is compared to that of an EB-MFEM-PIC and a FDTD based EM-PIC formulation in Figure 2.14. Lastly, we examine the expansion of a plasma ball. A spherical plasma ball is given a Gaussian distribution of ions and electrons in a geometry with first order absorbing 38 Figure 2.12: Particle beam comparison of DH-MFEM-PIC and FDTD-PIC. © 2021 IEEE (a) (b) Figure 2.13: Error in the continuity equation (2.13a) and Gauss’ law (2.13b) for the particle beam case. © 2021 IEEE boundary conditions. The experiment was simulated with the DH-FEM formulation with a point particle shape. A combination of 𝑆𝑟 + ions at 1K and electrons at 100K are placed at a density of 5 × 108 particles per cubic meter such that the system is initially charge neutral. Further details on the experiment setup can be found in [?] and derivations of the analytic results in [83]. As is evident from Figure 2.15, we achieve excellent agreement with the analytic results. Again, because all the conditions are satisfied, the continuity equation and errors in Gauss’ law are kept to machine precision. 39 Figure 2.14: Particle Beam for FDTD, 𝐸𝐵, and 𝐷𝐻 formulation © 2021 IEEE Figure 2.15: Electron density of expanding adiabatic plasma ball using DH-MFEM-PIC compared to analytic expected densities. © 2021 IEEE 40 2.6 Summary and conclusions In this work, we have presented a framework to define self-consistent EM-PIC schemes. By considering the space of functions used to represent fields, fluxes, and sources, one can ensure that all the relationships between Maxwell’s equations are preserved in a consistent manner. This has been demonstrated for FDTD and two formulations of FEM, and shown how it generalizes and relates to well known previous work. The framework is easily related to other field solvers not mentioned directly provided one considers the relationship of the basis functions. Future work will apply this framework to higher order basis functions and particle paths. © 2021 IEEE. Reprinted, with permission, from Crawford, Zane D., Scott O’Connor, John Luginsland, and B. Shanker. "Rubrics for Charge Conserving Current Mapping in Finite Element Electromagnetic Particle in Cell Methods." IEEE Transactions on Plasma Science (2021). 41 CHAPTER 3 HIGHER ORDER CHARGE CONSERVING ELECTROMAGNETIC FINITE ELEMENT PARTICLE IN CELL METHOD 3.1 Introduction Modeling novel beam-wave interaction devices, such as accelerators, vacuum electronics, and directed energy devices [84, 85, 86] relies on robust numerical tools capable of self- consistent analysis of the interaction of a plasma with electromagnetic fields. This is typically done using an electromagnetic particle-in-cell (EM-PIC) method to evolve a given plasma distribution in time and space [11]. It consists of a method that discretizes both the Newton’s equations of motion and a Maxwell’s field solver. The coupling between the two is effected through the Lorentz force due to the electric field and magnetic flux density. Given the range of applications, there has been extensive interest in developing PIC solvers; a majority of those used in field are based on finite difference time domain (FDTD) methods [50]. The simplicity of the formulation, ease of particle position updates, and readily available parallelization algorithms make this approach an attractive workhorse for PIC. In what follows, we will use the a concatenation of abbreviations to denote regime and method to solve PIC. For instance, EM-FDTDPIC denotes an electromagnetic PIC using FDTD. While EM-FDTDPIC has a number of advantages, there has been significant recent effort to explore the advantages offered by finite element methods to solving PIC problems. To a large part, this is driven by success of this methodology in microwave and millimeter frequency regimes wherein the success of this method has been demonstrated in analysis and design of complex topologies and electrically large objects. The foray of FEM into PIC is not without challenges, the principal of which is charge conservation. To understand this, note that in evolving the fields, we only solve the two curl equations (Faraday’s and 42 Ampere’s laws) and need a framework wherein Gauss’ law are satisfied as well. This implies that discretization in space and time should be such that these laws are satisfied. Ref. [87] rigorously develops the conditions that should be satisfied, and demonstrates how current EM-PIC formulations satisfy these conditions. For instance, the spatial and temporal basis sets used in an explicit FDTD time-stepping scheme, together with an appropriate integration of the path, satisfies these constraints. Developing such a method for for FEM was a long standing challenge. This was solved recently; see pioneering papers by [88, 51, 89, 90]. The methods introduced here were based on explicit updates of field solution and particle position, and on the proper representation of quantities on the underlying discrete mesh. In the same vein, a Poisson bracket approach that utilizes Whitney forms defined by B-spline FEM formulations [91] to define a structure-preserving EM-PIC scheme [53, 92], with several of these methods using higher order basis sets. Note, in manner akin to FDTD, one solves all of Maxwell’s equations. In an explicit setting, this approach avoids exciting null spaces (and corruption of Gauss’ law due to these null spaces) [93, 94, 95]. But restriction of field solve to an explicit field update has challenges; it is only conditionally stable and the smallest time step is governed by the finest feature in the model and not the physics. Overcoming this bottleneck has a well known remedy. Indeed, Newmark-Beta methods are unconditionally stable and the constraint on time comes from the physics that one needs to captures. But implicit field solve implies the need to rethink PIC solves such that the Gauss’ laws are satisfied. A problem that is unstated is that implicit solves introduce a null space; for Maxwell solvers, this null space is of the form ∇𝜙(r) and for the wave equation, this null space is of the form 𝑡∇𝜙(r). It is apparent that the null space will corrupt the satisfaction Gauss’ laws in addition to other challenges. This problem was solved recently [96, 97]. Specifically, imposition of Coulomb Gauge using a quasi-Helmholtz decomposition (in simply connected systems) in [97] enables satisfying Gauss’ law to machine precision for both the Maxwell solver and the wave equations. This 43 implies that the infrastructure that is already in place to solve the vector wave equation can readily used for PIC analysis. As we build this line of progress, the next ingredient that is missing is higher order basis sets for field, current and particle representation within a PIC framework. Hierarchical and interpolatory basis functions are known for FEM field solvers [98, 99, 100]. For smoothly varying geometries, higher order bases provide more accurate fields while utilizing fewer degrees of freedom. Furthermore, these higher order formulations satisfy the relationships of the de-Rahm complex [101, 102, 103, 104]. Developing a higher order EM-FEMPIC framework will be the key contribution of this paper. Specifically, we will present an unconditionally stable, exact current mapping FEM EM-PIC scheme that uses higher order basis functions on tetrahedral meshes. The rest of this paper is organized as follows: In Section 4.2 provide a brief problem statement. Next, in Section 3.3, we define the spatial and temporal discretization of the problem. Section 3.3.2 describes the current mapping scheme used to conserve charge regardless of the time marching scheme. In Section 4.4, we present results that demonstrate the use of the higher order FEM-PIC scheme. Finally, we conclude this paper in Section 4.5 outlining future directions of research. 3.2 Problem Statement Consider a region of free space Ω containing charged species that is bounded by 𝜕Ω. For simplicity we consider only a single species. The permittivity and permeability of free √ space are denoted as 𝜀0 and 𝜇0 , and the speed of light denoted using 𝑐 = 1/ 𝜇0 𝜀0 . There also exists a time-varying electromagnetic field due to moving charges and potentially an impressed electromagnetic field. The distribution of the charge species can be represented by a phase space distribution function (PSDF) 𝑓 (𝑡, r, v) that satisfies the Vlasov equation 𝜕𝑡 𝑓 (𝑡, r, v) + v · ∇ 𝑓 (𝑡, r, v)+ (3.1) 𝑞 [E(𝑡, r) + v × B(𝑡, r)] · ∇𝑣 𝑓 (𝑡, r, v) = 0. 𝑚 44 Next, will describe how the equations to be solved and how the problem is discretized in space and time. 3.3 Overview of Discrete Solutions To solve the above equations, we start with defining charge and current density in ∫ ∫ terms of the PSDF, 𝜌(𝑡, r) = 𝑞 Ω 𝑓 (𝑡, r, v)𝑑v and J(𝑡, r) = 𝑞 Ω v 𝑓 (𝑡, r, v)𝑑v. Using a particle approximation using 𝑁𝑝 samples with shape functions 𝑆(r), one obtains 𝑁𝑝 Õ 𝜌(𝑡, r) = 𝑞 𝑆(r − r𝑝 (𝑡)) (3.2a) 𝑝=1 𝑁𝑝 Õ J(𝑡, r) = 𝑞 v(𝑡)𝑆(r − r𝑝 (𝑡)) (3.2b) 𝑝=1 where r𝑝 (𝑡) and v𝑝 (𝑡) are the position and velocity of particle 𝑝. In this work the shape functions are chose to eb Dirac delta functions, though generalization to other shape functions is possible [87]. The particular choice shape function is immaterial to the results of this paper. The electromagnetic fields, E(𝑡, r) and B(𝑡, r), are solutions to Maxwell’s curl equations 𝜕B(𝑡, r) − = ∇ × E(𝑡, r) (3.3a) 𝜕𝑡 𝜕D(𝑡, r) = ∇ × H(𝑡, r) − J(𝑡, r) (3.3b) 𝜕𝑡 and satisfy Gauss’ laws ∇ · B(𝑡, r) = 0 (3.4a) ∇ · D(𝑡, r) = 𝜌(𝑡, r). (3.4b) The fields are also subject to boundary conditions which are either Dirichlet, Neumann, or impedance boundary conditions on 𝜕Ω𝐷 , 𝜕Ω𝑁 , or 𝜕Ω𝐼 which bound the domain as 𝑛ˆ × E(𝑡, r) = Ψ𝐷 (𝑡, r) on 𝜕Ω𝐷 , (3.5a) 𝑛ˆ × H(𝑡, r) = Ψ𝑁 (𝑡, r) on 𝜕Ω𝑁 , (3.5b) 45 B(𝑡, r) 𝑛ˆ × − 𝑌 𝑛ˆ × 𝑛ˆ × E(𝑡, r) = Ψ𝐼 (𝑡, r) on 𝜕Ω𝐼 . (3.5c) 𝜇 The fields, which are defined with respect to the free space background medium, are related to the fluxes through constitutive relations D(𝑡, r) = 𝜀0 E(𝑡, r) (3.6a) H(𝑡, r) = 𝜇−10 B(𝑡, r). (3.6b) The source definitions are used to evolve the particle positions using Newton’s equations of motion and Lorentz force, viz., F(𝑡, r) = 𝑞(𝑡, r)(E(𝑡, r) + v(𝑡, r) × B(𝑡, r)). The simulation follows the usual PIC cycle: particles are mapped to a discretized space to solve for the electric field and magnetic flux density, which are in turn used to push the particles, defining a new current and particle positions, restarting the cycle. Though the PSDF is sampled with particles, the total description of the electromagnetic problem is continuous, and must be discretized in space and time. 3.3.1 Discretization in Space Assume that the domain Ω is represented using a collection of finite elements 𝒦 = {𝒩 , ℰ, ℱ , 𝒯 } defined using 𝑁𝑛 nodes, 𝑁𝑒 edges, 𝑁 𝑓 faces and 𝑁𝑡 tetrahedron. Each tetrahedron contains basis functions to represent fields, fluxes, and sources that follow the de-Rham sequence as seen in Fig. 3.1[99, 100]. This sequence preserves the differential relations between the quantities of interest, such that the curl of a field is a flux, the divergence of a flux is a charge. The Hodge star operator ★ maps a field to a flux on a dual mesh.It is well known that Whitney basis functions can be used to represent the electric field and magnetic flux [105, 106, 51]. For a 𝑘th order interpolatory basis functions, the electric field is represented using Í 𝑁1 (1) 𝑁1 higher order Whitney edge basis functions, E(𝑡, r) = 𝑒 (𝑡)W𝑖 (r), 𝑖=1 𝑖 where there are 6(𝑘 + 1) degrees of freedom associated with edges, 4𝑘(𝑘 + 1) degrees of freedom associated with faces, and 𝑘(𝑘 2 − 1)/2 associated with the cell volumes. The magnetic flux density is 46 Figure 3.1: Relationship of fields and fluxes with respect to the de-Rham complex. Í 𝑁2 (2) represented using 𝑁2 Whitney face basis function, B(𝑡, r) = 𝑖=1 𝑏 𝑖 (𝑡)W𝑖 (r) where there are (𝑘 + 1)(𝑘 + 2)/2 degrees of freedom associated with faces and 𝑘(𝑘 + 1)(𝑘 + 2)/2 associated with the cell volumes. A more complete description of these basis functions is provided in the Appendix. We have chosen to define the problem such that Faraday’s law (3.3a), E(𝑡, r), and B(𝑡, r) are defined on the primal grid. Therefore, Ampere’s law (3.3b), H(𝑡, r), D(𝑡, r), as well as J(𝑡, r),and 𝜌(𝑡, r) are defined in the dual space on the corresponding dual grid. This means that while the electric field and magnetic flux can be directly represented using a Whitney basis on the primal mesh, the sources 𝜌(𝑡, r), J(𝑡, r) cannot. On structured grids it is straightforward to define dual basis function spaces to represent dual quantities, however, on an unstructured FEM mesh the dual fields are only indirectly accessible via Hodge operators. Therefore the source distribution 𝜌(𝑡, r) cannot be directly represented, and is instead measured with higher order nodal basis functions on the primal mesh, which are defined in the Appendix. The current density J(𝑡, r), which lies in the same space as D(𝑡, r), is measured by the the electric field basis, the higher order Whitney edge functions. Before we proceed with prescribing the discrete framework, consider an auxillary function ∫ 𝑡 G(𝑡, r) = J(𝜏, r)𝑑𝜏 (3.7) 0 such that Ampere’s law is rewritten as 𝜕D(𝑡, r) 𝜕G(𝑡, r) = ∇ × H(𝑡, r) − . (3.8) 𝜕𝑡 𝜕𝑡 Using (3.8) and Faraday’s law and spatial basis functions defined earlier, one may write the 47 discrete system as  ¯  [★𝜇−1 ] 0  𝜕𝑡 𝐵(𝑡)      𝜕 ¯ 𝐸(𝑡)     0  [★ 𝜖 ] 0   𝑡   | {z } 𝑀¯¯ (3.9)  0  ˜  𝐵(𝑡) [∇×] ¯   0  +   =  ¯    𝜕𝑡 𝐺(𝑡) ˜ 𝑇 ¯ −[∇×] [★𝐼 ]  𝐸(𝑡) − 𝜖        | {z } | {z } 𝑆¯¯ 𝐹¯¯ where the degree of freedom vectors 𝐸(𝑡) ¯ = [𝑒1 (𝑡), 𝑒2 (𝑡), . . . , 𝑒 𝑁𝑒 (𝑡)], 𝐵(𝑡) ¯ = [𝑏 1 (𝑡), 𝑏2 (𝑡), . . . , 𝑏 𝑁 (𝑡)], 𝑓 ¯ (1) and 𝐺(𝑡) = [𝑔1 (𝑡), 𝑔2 (𝑡), ...𝑔𝑁𝑒 (𝑡)] with 𝑔𝑖 (𝑡) = hW𝑖 (r), G(𝑡, r)i.The coupled system matrix is composed of discrete Hodge matrix operators (1) (1) [★𝜖 ]𝑖,𝑗 = hW𝑖 (r), 𝜀 · W 𝑗 (r)i (3.10) (2) (2) [★𝜇−1 ]𝑖,𝑗 = hW𝑖 (r), 𝜇−1 · W 𝑗 (r)i, (3.11) the surface impedance matrix (1) (1) [★𝐼 ]𝑖,𝑗 = h𝑛ˆ 𝑖 × W𝑖 (r), 𝜇−1 𝑛ˆ 𝑗 × W 𝑗 (r)i𝜕Ω (3.12) and a discrete curl matrix (2) ˜ 𝑖,𝑗 = hW , ∇ × W (r)i. (1) [∇×] 𝑖 𝑗 (3.13) The operator h·, ·i and h·, ·i𝜕Ω define a volume and surface integral, respectively, over the support of the basis functions, which is either a tetrahedron or face. The curl matrix [∇×] ˜ includes the metric information, unlike the definition usually obtained through discrete exterior calculus [∇×], which has entries of only 0,+1, or -1. For lowest order 𝑘 = 1, this definition can be obtained as [∇×] = [★𝜇−1 ]−1 [𝑀 𝑐 ]; (3.14) however, the simplicity of defining the matrix by inspection is lost for 𝑘 > 1. 48 3.3.2 Evolution of Particle Path and Current Mapping The fields are evolved in time using an unconditionally stable Newmark-beta time marching scheme [107]. This allows larger time steps than would be afforded by a leapfrog method. In this framework, the current mapping (or evolution of charge) has to be consistent with that used for evolution of fields. Unfortunately, a naive approach to incorporate the particle current as the forcing function will violate conservation of charge. The method presented in [96] is agnostic to any time stepping method used for a field solve, and overcomes this bottleneck. In this paper, a similar method is used, but adapted to a higher order basis function in space. Using the definition of the particle current density in (3.7), the Newmark-beta time marching scheme is defined as (𝛾 𝑀¯¯ + 𝛽Δ 𝑆) ¯¯ ¯ 𝑛+1 − 𝛾Δ 𝑆¯¯ 𝑋¯ 𝑛 + (𝛾 𝑀 ¯¯ + 𝛽Δ 𝑆¯¯ )𝑋¯ 𝑛−1 𝑡 𝑋 𝑡 𝑡 𝑀 (3.15) + 𝛾Δ𝑡 𝐹¯¯ 𝑛+1 + 𝛾Δ𝑡 𝐹¯¯ 𝑛−1 = 0 where 𝑋¯ 𝑚 = [𝐵¯ 𝑇 (𝑡 𝑚 ) 𝐸¯ 𝑇 (𝑡 𝑚 )] and 𝐹¯¯ 𝑚 = [0 − 𝜀−1 𝐺¯ 𝑇 (𝑡 𝑚 )]. The degree of freedom vector (1) 𝐺¯ = [𝑔1 (𝑡), 𝑔2 (𝑡), · · · , 𝑔𝑁𝑒 (𝑡)] with 𝑔𝑖 (𝑡) = hW𝑖 , G(𝑡, r)i. The parameters 𝛾 and 𝛽 are used to change the stability of the scheme, often choosing 𝛾 = .5 and 𝛽 = .25. To define the forcing function 𝐹, ˜ it is necessary to use an integration rule appropriate for the product of the higher order edge basis function and the particle path, which may also be a higher order polynomial. A key point that should be noted is the de-linking of the time stepping algorithm used for particle push and field updates. Consistency is ensured by proper inclusion in the right hand side of (4.3). Note, that this assumes non-relavistic motion. In this work, a fourth order Adams-Bashforth push is used, making the particle path a fourth order Lagrange polynomial. 3.3.3 Satisfaction of Gauss’s Magnetic Law While the above equation can be used, as shown [97] one would need to impose a Coulomb gauge to ensure that the null spaces do not corrupt the solution. The level 49 of the null space that is excited depends on the accuracy of the solution at every time step; which in turn, is related to how well the initial conditions are remembered in a time marching system. While [97] used a topological framework, here we investigate an alternate approach. Our goal is more modest and targeted. We seek to examine particle motion due to impressed magnetic lenses, we want to ensure the discrete representation of impressed magnetic flux density is divergence free when represented using Whitney basis. To that end, consider reconstructing an impressed magnetic flux density B𝑖 (𝑡, r) which satisfies (3.4a). The usual approach to obtain the coefficients to approximate B𝑖 (𝑡, r) would be to use Galerkin testing with the divergence-conforming basis set [★𝜇−1 ]b = f (3.16) (2) where f𝑖 = hW𝑖 (r), B𝑖 (𝑡, r)i. However, this construction of B̃𝑖 (𝑡, r) will not satisfy Gauss’s law unless B̃𝑖 (𝑡, r) = B𝑖 (𝑡, r). This is accomplished by solving an optimization problem where Gauss’s magnetic law is the constraint. [★𝜇−1 ] [∇·]𝑇  b ℬ         =  (3.17)  𝜆  0        [∇·] 0      For the lowest order spatial basis functions, the discrete divergence operator [∇·] can be written by inspection. The discrete divergence operator for higher orders, like the discrete curl operator, cannot be written as easily. It can be written as ˜ [∇·] = [★3 ]−1 [∇·] (3.18) where (3) (3) [★3 ]𝑖𝑗 = h𝑊𝑖 (r), 𝑊𝑗 (r)i (3.19) and (3) ˜ 𝑖𝑗 = h𝑊 (r), ∇ · W (r)i (2) [∇·] 𝑖 𝑗 (3.20) 50 Table 3.1: Satisfaction of Gauss’ Magnetic Law for Linear Field 𝑚 err in B err in ∇ · B err in B err in ∇ · B 1 1.0604 × 10 −1 2.9681 × 10 −5 1.1367 × 10 4.5062 × 10−19 −1 2 1.1049 × 10−15 6.1606 × 10−19 1.2201 × 10−15 4.4173 × 10−19 3 3.8299 × 10−15 1.9566 × 10−18 9.8539 × 10−15 3.9481 × 10−18 Table 3.2: Satisfaction of Gauss’ Magnetic Law for Sinusoidal Field 𝑚 err in B err in ∇ · B err in B err in ∇ · B 1 1.3198 × 10 1.6310 × 10 1.5048 × 10 1.1404 × 10−18 −2 −5 −2 2 6.1965 × 10−4 3.1300 × 10−7 5.8211 × 10−4 1.2061 × 10−18 3 6.3122 × 10−6 2.7822 × 10−9 7.0793 × 10−6 1.0762 × 10−17 with 𝑊 ( 3)(r) as the higher order volumetric basis function, which is defined in the Appendix. This allows the reconstructed field, regardless of the error in the representation of the original function, to still satisfy (3.4a). Consider the tests done in Table 3.1 and 3.2. In this test, a divergence free function was reconstructed in a volume 50 cm × 15 cm × 20 cm, for several orders of spatial basis functions. In Table 3.1, the function B(r) = 𝑦 𝑥ˆ + 𝑥 𝑦ˆ can be reconstructed exactly with second order and higher basis functions. Therefore, once the function was modeled correctly to machine precision, the divergence free nature of the reconstructed field is seen in both the constrained and non-constrained formulation. However, in Table 3.2, the function was B(r) = sin 𝑦 𝑥ˆ + cos(𝑥) 𝑦),  ˆ which cannot be represented exactly by a finite set of polynomials. Despite having similar accuracy, only the constrained formulation satisfies (3.4a). It is evident that the formulation has the desired properties with respect to representation of fields. As alluded to earlier, we will use this only for representing only the impressed field. A topological approach, akin to [97], is being developed an will be presented in later papers. 3.4 Results In this section, we present several numerical tests using higher order FEM-PIC. Particle free results are provided here to demonstrate correctness of our implementation. We 51 demonstrate that the higher order basis functions presented satisfies the continuity equation and Gauss’ law for a number of cases using the stated particle mapping scheme. 3.4.1 Higher Order Cost Here we present the error in computed fields with respect to number of unknowns and basis function order. The relative error in a field propagating through a region of free space is shown in Fig. 3.2. The region is 0.25 m × 0.2 m × 0.5 m in size. The normally incident electric field is defined as ¯ 𝑟 , 𝑡) = 𝑦ˆ cos 2𝜋 𝑓0 𝜏 𝑒 −(𝜏−8𝜎)2 /2𝜎2 (V/m), 𝐸(¯  (3.21) where where 𝜏 = 𝑡 − 𝑟¯ · 𝑧ˆ /𝑐, 𝜎 = 3/[2( 𝑓max − 𝑓0 )] with the center 𝑓0 = and maximum frequency denoted 𝑓max = 195𝑀𝐻𝑧. The relative error is defined as (𝑟) (𝑎) ||E𝑖 − E𝑖 || 2 relative error = (3.22) ||E(𝑎) || 2 . It take several orders of magnitude more unknowns for a first order basis function to reach the same level of accuracy as a second order basis function. The trade off is the condition number of the system increases by roughly an order of magnitude as the order increases which effects the rate at which an iterative solver will converge. Therefore, consideration can be taken in balancing the size of the problem and simulation time for a given error. 3.4.2 Higher Order Particle Motion The first example with particles is the orbit of a single particle around a nucleus. This test demonstrates when higher order bases give more accurate particle trajectories. At a certain distance, a particle with initial velocity perpendicular to the radial electric field will result ˆ in a circular orbit in a plane. In this test, a particle with initial velocity v = 3 × 106 𝜙m/s is set .25m from the centroid of a cylindrical ring geometry. The geometry has an inner radius of .2m and outer radius of .3m. Five meshes were generated to compare refinement 52 Figure 3.2: Relative Error in E𝑡,r for field propagating through free space. in the edge length to basis function order. The electric field due to a nucleus with a charge of 𝑄 𝑛 = 1.423𝑛𝐶 is reconstructed using the higher order interpolatory Whitney edge basis set and used to push the electron in the geometry. The analytic definition of the electric field is 𝑄 𝑛 𝑟ˆ E(r) = (3.23) 4𝜋𝜀r · r which, unlike the constant magnetic field used for cyclotron motion experiments, can be improved by using higher order spatial basis functions. The experiment was run for 6000 time steps with Δ𝑡 =0.4ns, which corresponds to approximately 4 cycles. The effect of the higher order basis functions can be seen in 3.3 and 3.4 where the relative error is defined in (3.22). For both the 𝑧ˆ and 𝜌ˆ components of the electric field, the error in the fields converges. The improvement of the field error translates to an improvement of the particle trajectory. In 3.5, the relative error is shown for the 𝑧ˆ component of the particle position. A key takeaway is that using a higher order basis function leads to more accurate particle trajectories than simply refining the mesh. 53 Figure 3.3: Relative Error in E𝜌 for orbital motion Figure 3.4: Relative Error in E𝑧 for orbital motion. 3.4.3 Plasma Ball In this example, we simulate an adiabatic expansion of a plasma ball. This example has both approximate analytic solutions [108] as well as experimental data [109]. A Gaussian distribution of 12000 𝑆𝑟 + ions and electrons a placed at the center of a spherical geometry that enclosed with a first order absorbing boundary condition. The initial temperature of 𝑆𝑟 + ions is 1K and electrons are 100K placed at a density of 5 × 108 particles per cubic meter such that the particles are and will remain sufficiently away from the boundary. 54 Figure 3.5: Relative Error in r for orbital motion. Here, we use three geometries; the first with a radius of 6 cm with an avergae edge length of 1.02 cm, the second 12 cm with average edge length of 2.04 cm, and the third 18 cm with average edge length of 3.06 cm. The experiment was run with first and second order basis functions, with a comparison to the analytic solution in Table 3.3. Though there is good agreement between all of the experimental and analytic data, there is not a clear improvement as order increases. This is due to the fields being well behaved enough in this example that they are approximated well enough by the first order basis functions. 3.4.4 Expanding Particle Beam In this test, we demonstrate an expanding plasma beam in the PEC cavity. This test is a standard test to confirm charge conservation as errors will accumulate and cause striations in the beam. Additionally, quasi-analytic solutions and comparisons to this method and other discretization schemes can be found in [110]. Macroparticles are injected into the cavity at an initial velocity and are allowed to repel each other as they progress down the cavity. The parameters used in this experiment are included in Table 3.4. Key here is that the higher order basis functions are defined such that the differential relations between the basis functions for the electric field and magnetic flux are preserved. 55 Table 3.3: Plasma Ball expansion for different basis function orders and geometry radii 𝑘 1 2 Radius 6 cm 12 cm 18 cm Table 3.4: Expanding Particle Beam Parameters Parameter Value Cavity Radius 20 mm Cavity Length 100 mm Boundary Conditions PEC 𝑣0 1.02 · 107 m/s 𝑣0 /𝑐 0.16678 beam radius 𝑟𝑏 8.00 mm Number particles per time step 10 species electrons Turn on time 2 ns beam current 5 mA macro-particle size 103921.12 min edge length 3.89 mm max edge length 13.5 mm Δ𝑡 33.3 ps 56 Figure 3.6: Relative Error of 𝜌(𝑡, ˜ r) vs −[∇]G(𝑡, r) for particle beam. Table 3.5: Error in impressed Quadrupole Fields basis order 1 basis order 2 basis order 3 𝑁2 error 𝑁2 error 𝑁2 error 140629 9.835 × 10−2 89502 1.143 × 10−2 28488 1.726 × 10−2 597623 3.8299 × 10−2 343812 3.803 × 10−3 103876 3.9481 × 10−3 637050 2.424 × 10−3 240270 6.932 × 10−4 First, consider the error in the measured charge density as defined by (3.22) shown in 3.6. Here, the error is small and saturates as the total number of particles in the cavity stabilizes. The increase of error as the basis function order increases can be attributed to the increase of the condition number of the mass matrix that is inverted to compute the divergence of the integrated current. The error between the divergence of the electric flux and integrated current is shown in 3.7. The error is near machine precision, demonstrating the higher order divergence operator acting on the curl of the magnetic field does go to zero. Lastly, different time marching schemes can be taken by varying the values of 𝛾 and 𝛽 in (4.4) using first order spatial basis functions. The error in Gauss’ law for three different time marching schemes is shown in 3.8, where average acceleration has 𝛾 = .5, 𝛽 = .25, backward difference has 𝛾 = 1.5, 𝛽 = 1„ and galerkin has 𝛾 = 1.5, 𝛽 = .8. The error remains at machine precision as expected. In total, this shows that the representation of the particle in the simulation is correct and that Gauss’ law is satisfied. 57 Figure 3.7: Relative Error of [∇·]D(𝑡, r) vs−[∇]G(𝑡, r) for particle beam. Figure 3.8: Relative Error of [∇·]D(𝑡, r) vs−[∇]G(𝑡, r) for particle beam with different time marching schemes. 58 Figure 3.9: Schematic for perfect electrical conducting box with Pafonsky quadrapole (shaded). 3.4.5 Pafonsky Quadrupole In this example, the Panofsky quadrupole used in [111] is modeled. This example uses higher order basis functions to get more accurate impressed fields using fewer degrees of freedom than a lower order method. From Table 3.5, a mesh with more than several million degrees of freedom would be needed to obtain a result with similar error in the impressed fields. This would significantly effect the time needed to simulate the result. The quadrupole is 13.4 cm in width, 5.6cm in height, and .044 cm in height, shown in Fig. 3.9. We define a rectangular PEC cavity with the quadrupole .011cm from the xy-plane. The fields due to the quadrapole are calculated by evaluating Bio-Savart’s law using the defined current density. It would be too expensive to recalculate the fields as particles pass through the system, therefore the coeffiecients needed to reconstruct the fields are obtained through (3.17), which also enforces (3.4a). Particles are emitted from the xy-plane at 𝑧 = 0, passing through the fields generated by the quadrupole. The beam is emitted such that the beam will focus in the 𝑦-dimension and defocus in the 𝑥-dimension. The iterative solver tolerance for both the fields and the divergence enforcement was set to 1e-5. A snapshot of the particle positions at 35.7 ns is shown in Fig 3.10 and 3.11, where the background magnetic flux density B 𝑦 (𝑡, r) and B𝑥 (𝑡, r) are plotted, respectively. The particle beam smoothly expands in the 𝑥-dimension and compresses in the 𝑦-dimension which is a sign that charge is being conserved correctly. 59 Figure 3.10: Particle beam traveling in 𝑥 − 𝑧 plane in Pafonsky quadrapole with background B 𝑦 (𝑡, r). Figure 3.11: Particle beam traveling in 𝑦 − 𝑧 plane in Pafonsky quadrapole with background B𝑥 (𝑡, r). 3.5 Summary In this work, we have presented a higher order, exact current mapped FEM-based EM-PIC with unconditionally stable time marching. When higher order fields dominate the physics being modeled, higher order basis functions can reduce the number of degrees of freedom needed to model it while also getting more accurate fields. A method was also provided to ensure that Gauss’ magnetic law is satisfied, which can occur when time marching schemes other than leap frog are used. Future work will create quasi-Helmholtz decomposition for higher order basis sets which will further increase the efficiency of using higher order basis while also allowing fields obtained by the vector wave equation to be used. Even with the use of higher order basis, there is a limit to the size of problem that can fit on a single node. To handle this challenge, domain decomposition approaches for Maxwell solvers will be developed. 60 CHAPTER 4 DOMAIN DECOMPOSITION FRAMEWORK FOR MAXWELL FINITE ELEMENT SOLVERS AND APPLICATION TO PIC 4.1 Introduction The design of novel accelerator and vacuum devices require computational tools capable of simulating self-consistent interactions between plasma and electromagnetic fields [85, 112, 113, 114]. This is commonly done by evolving a particle distribution in time due to the Lorentz force defined by an electric field and magnetic flux density; this includes both external/impressed and those generated by the acceleration of charged species. The predominant approach is the finite difference time domain based particle in cell method (EM-FDTDPIC) [11]. This method has seen decades of use and advancement, primarily due its ease of use and efficiency, as the fields can be determined from a simple recurrence formula rather than matrix inversion. Furthermore, parallel implementations and fast techniques have enabled EM-FDTDPIC to simulate large problems [115, 116, 117, 118]. However, there has been interest in exploring alternatives to EM-FDTDPIC. To a large extent, this is motivated by the higher order fidelity that one can obtain in terms of both geometry as well as physics. In response to this need, the past decade has seen extensive effort to create an electro- magnetic particle in cell method based around the finite element method (EM-FEMPIC); see [50, 51] and papers therein. It has been shown that by choosing proper basis sets for the fields and currents, it is possible to obtain charge conserving EM-FEMPIC schemes through exact current mapping [119, 90, 51]. Unfortunately, these methods use explicit time stepping with first order basis functions to represent both the fields, paticles, currents, and time. More recent work has utilized higher order spatial basis better represent fields, but the error in representation is still constrained by the representation for time [120]; 61 note, the moniker “structure-preserving” is used when referring to a method that uses Whitney basis/ensures that the de-Rham relations hold. The downside of using explicit time stepping is that it is conditionally stable, and the time step to ensure stability is determined by the smallest feature size. Obviously, this is a challenge when analyzing plasma in geometrically complex in structures. More recent work provided a general framework for charge conservation to be satisfied, irrespective of time stepping scheme [87]. A practical demonstration of satisfaction of conservation laws for implicit unconditionally stable finite element scheme for Maxwell equation was shown for a set of test problems in [96]. It was also shown in this paper, that the proposed methodology was valid for a wave equation FEM solver as well. Unfortunately, it is well known that there exists a null space that grows as 𝑡∇𝜙(r) for the wave equation solvers and a null space that varies as ∇𝜙(r) for implicit Maxwell solvers, where 𝜙(r) is a scalar field. While it may be possible to control the magnitude of the null space excited in Maxwell solvers, the one for the wave equation unfortunately rears its head with time. The development of a discrete Coulomb gauge in [97] shows how one may obtain a method that satisfies Gauss’ laws for both wave equation and Maxwell equation solvers. Despite this progress, a fundamental bottleneck remains. All EM-FEMPIC methods are more expensive than EM-FDTDPIC methods. The reason for the efficiency of FDTD schemes fairly simple; there is no inversion thanks to the structure of the method [121]. Unfortunately, FEM methods typically involve 3 inversion of matrices. As is well known, a naive direct inversion would cost 𝒪(𝑁𝑑𝑜 𝑓 ) per time step where 𝑁𝑑𝑜 𝑓 is the number of degrees of freedom. This can be reduced to 𝜈 𝒪(𝑁𝑑𝑜 𝑓 log 𝑁𝑑𝑜 𝑓 ) if a multifrontal sparse solver is used [?, 122] with 𝒪(𝑁𝑑𝑜 𝑓 𝑝𝑜𝑙𝑦𝑙𝑜 𝑔(𝑁𝑑𝑜 𝑓 )) for factorization where 𝜈 < 2 [123]. Alternatively, one typically uses an iterative solver whose cost scales as 𝒪(𝑁𝑖𝑡𝑒𝑟 𝑁𝑑𝑜 𝑓 ), where 𝑁𝑖𝑡𝑒𝑟 is the number of iterations which depends on the error threshold. As a result, one takes recourse to preconditioners. But it has been shown that the most computationally effective approach has been to use domain 62 decomposition [124, 125]. A bulk of the effort has been applying this technique to amortize the cost of wave equations and largely in the Fourier domain [126, 127, 128, 129], with some scant papers in the time domain [130, 131]. The primary benefit of this approach is to reduce the complexity of the problem by effectively reducing the size of 𝑁𝑑𝑜 𝑓 that needs to be inverted. The domain decomposition approach in this work is the finite element tearing and interconnecting (FETI) framework. Here, the domain is partitioned into multiple regions/- subdomains. For each subdomain, the unknowns are divided into interior and boundary quantities. The quantities on the boundary are shared between multiple subdomains, and are related via appropriate transmission conditions. Via manipulations, one only solves for unknowns on the boundary of the subdomains which is far fewer that the overall number of unknowns. All interior unknowns can be related to the ones on the boundary. Note, this differs from the domain decomposition approaches currently used for FDTD, which takes a Schwarz approach which uses overlapping subdomains, requiring iteration of each subdomain to get a converged solution[132]. Thus, the main contributions of this paper are as follows: (a) we derive two FETI schemes for time domain Maxwell equation solvers, (b) these methods are then used within an unconditionally stable EM-FEMPIC framework using an exact current mapping scheme, and (c) we present a number of results that demonstrate the efficacy of the proposed technique. The rest of this paper is organized as follows. Section 4.2 will define the problem setup. Section 4.2.1 will briefly review particle mapping approach used followed by the field discretization in space and time in Section 4.2.2. Section 4.3 will define the two domain decomposition formulations. Results using both of these methods will be presented in Section 4.4 followed by conclusions and future work in Section 4.5. 63 4.2 Preliminaries Let a region Ω ∈ R3 be bounded by surface 𝜕Ω and contain at least one charged species. The charges are represented by a phase space distribution function (PSDF) that is subject to the Maxwell-Vlasov equation 𝜕𝑡 𝑓 (𝑡, r, v) + v · ∇ 𝑓 (𝑡, r, v)+ (4.1) 𝑞 [E(𝑡, r) + v × B(𝑡, r)] · ∇𝑣 𝑓 (𝑡, r, v) = 0. 𝑚 The electric field E(𝑡, r) and magnetic flux density B(𝑡, r) are obtained through the mixed finite element method, which discretizes Faraday’s law and Ampere’s law. The PSDF is evolved through the PIC cycle, which consists of: 1. Gathering the particles onto the mesh to be mapped into the field solve, 2. Solving for the fields due to the current from the particle motion, 3. Interpolating the fields at the particle positions, 4. Pushing the particles to new positions. The following subsections will describe how Maxwell’s equations are discretized in space and time, and how the particles are mapped into the field solver. 4.2.1 Current Mapping The charge and current density is defined using the PSDF conventional definition, ∫ ∫ 𝜌(𝑡, r) = 𝑞 Ω 𝑓 (𝑡, r, v)𝑑v and J𝜌 (𝑡, r) = 𝑞 Ω v 𝑓 (𝑡, r, v)𝑑v. The PSDF is approximated by using 𝑁𝑝 samples, or particles, with shape 𝑆(r) such that the charge and current density Í 𝑁𝑝 Í 𝑁𝑝 are defined as 𝜌(𝑡, r) = 𝑞 𝑝=1 𝑆(r − r𝑝 (𝑡)) and J𝜌 (𝑡, r) = 𝑞 𝑝=1 v(𝑡)𝑆(r − r𝑝 (𝑡)) where r𝑝 (𝑡) and v𝑝 (𝑡) are the position and velocity of particle 𝑝. In this work, the shape functions are delta functions, though other shape functions are also valid and can be used to improve noise quality. The proper function to measure the charge density is the 0-form and the 64 1-form for the current density. Different schemes to evolve Newton’s equations can be used to preserve certain quantities like phase. In this work, to maintain accuracy of the time evolution at larger time step sizes, the Adams-Bashforth scheme is used. Lastly, the time integral of the current density ∫ 𝑡 G(𝑡, r) = 𝑑𝜏J𝜌 (𝜏, r) (4.2) 0 is used to allow the particle current to evolve consistently with the fields. Otherwise, there is a disconnect between the continuity equation and Gauss’ law. 4.2.2 Field discretization 4.2.2.1 Spatial Discretization The FEM formulation used in this work is the mixed finite element method. The basis functions used are the 1-form W(1) (r) ∈ 𝐻(𝑐𝑢𝑟𝑙; Ω) which allows for tangential continuity of fields, and the 2-form W(2) (r) ∈ 𝐻(𝑑𝑖𝑣; Ω) which allows for normal continuity of fluxes. The Í 𝑁𝑒 reconstructed fields used to accelerate the particles are defined as E(𝑡, r) = 𝑒 (𝑡)W(1) (r) 𝑖=1 𝑖 Í𝑁 𝑓 and B(𝑡, r) = 𝑖=1 𝑏 𝑖 (𝑡)W(2) (r) where 𝑁𝑒 are the number of edges and 𝑁 𝑓 are the number of faces. Testing Faraday’s law with the 2-form and Ampere’s law with the 1-form leads to the semidiscrete Maxwell system [★𝜇−1 ] 0  𝜕𝑡 𝐵¯         𝜕𝑡 𝐸¯      0  [★𝜖 0 ]  | {z } 𝐴¯¯ 1 (4.3) 0 [𝑀 𝑐 ]  𝐵¯  0       +   =    𝐿¯  −𝑐 2 [𝑀 𝑐 ]𝑇 𝑐[★𝐼 ] 𝐸¯    𝜀      | {z } |{z} 𝐴¯¯ 0 𝐹¯¯ where the degree of freedom vectors 𝐸¯ = [𝑒1 (𝑡), 𝑒2 (𝑡), . . . , 𝑒 𝑁𝑒 (𝑡)], 𝐵¯ = [𝑏 1 (𝑡), 𝑏2 (𝑡), . . . , 𝑏 𝑁 𝑓 (𝑡)], (1) and 𝐿¯ = [𝑙1 (𝑡), 𝑙2 (𝑡), ...𝑙 𝑁𝑒 (𝑡)] with 𝑙 𝑖 (𝑡) = hW𝑖 (r), J𝑠 (𝑡, r) − 𝜕𝑡 G(𝑡, r)i. The surface currents 65 J𝑠 (𝑡, r) exist on 𝜕Ω due to a Neumann or Robin boundary condition. The integrated particle current is chosen to satisfy charge conservation for arbitrary time marching schemes. 4.2.2.2 Temporal Discretization This work utilizes the Newmark-𝛽 time marching scheme to evolve the fields in time with parameters 𝛾 and 𝛽, which creates an unconditionally stable time marching scheme when 𝛾 = .5 and 𝛽 = .25. Other choices can lead to different stability conditions. The parameter choices effect the definition of a temporal basis function 𝑁(𝑡) that is tested by a function 𝑊(𝑡) [133]. When (4.3) is discretized with Newmark-𝛽 with 𝛾 = 0.5 and 𝛽 = 0.25, the fully discretized system becomes (0.5𝐴¯¯ 1 + 0.25Δ𝑡 𝐴¯¯ 0 )𝑋¯ 𝑛+1 − .5Δ𝑡 𝐴¯¯ 0 𝑋¯ 𝑛 + (0.5𝐴¯¯ 1 + .25Δ𝑡 𝐴¯¯ 0 )𝑋¯ 𝑛−1 + 0.5Δ𝑡 𝐺˜ 𝑛+1 + 0.5Δ𝑡 𝐺˜ 𝑛−1 (4.4) + 0.25Δ𝑡 𝐹˜ 𝑛+1 − 0.5Δ𝑡 𝐹˜ 𝑛 + 0.25𝐹˜ 𝑛−1 = 0. Here, 𝑋¯ 𝑚 = [𝐵¯ 𝑚 𝐸¯ 𝑚 ]𝑇 , 𝐺˜ 𝑚 = [0 𝜀−1 𝐺] ¯ 𝑚,𝑇 , and 𝐹˜ 𝑚 = [0 𝜀−1 𝐽¯𝑠𝑚 ] at time step 𝑚. The boundary current 𝐽¯𝑠𝑚 is tangential to a surface. As it has no normal component, it does not contribute to Gauss’ law and therefore does violate charge conservation. 4.3 Domain Decomposition Consider the region Ω depicted in Fig. 4.1 that is divided into 𝑁𝑣 non-overlapping subdomains where subdomain Ω𝑖 and Ω 𝑗 are separated by boundary Γ𝑖𝑗 . The junction of more than two volumes is a "corner" denoted by Γ𝑐 . In each subdomain Ω𝑖 , (4.4) is defined and is assumed to be a self contained “primal" problem that has fictitious excitations from "dual" unknowns. The dual unknowns, the Lagrange multipliers denoted by Λ, effect either Neumann or Robin boundary conditions such that the equivalence theorem can be applied for each Ω𝑖 with respect to Ω using the external currents on 𝜕Ω𝑖 , the fictitious boundary current from Λ, and the particle current in Ω𝑖 . By solving for the dual unknowns 66 Figure 4.1: Region Ω partitioned into 𝑁𝑣 subdomains. and a small subset of the primal unknowns, the overall computational cost of solving the original monolithic system in (4.4) is reduced. 4.3.1 MFEM-DP 1 The first formulation presented uses a Neumann-like boundary condition to enforce the continuity conditions of the fields at the interior boundaries, 𝑛ˆ × E𝑖 (𝑡, r) = 𝑛ˆ × E 𝑗 (𝑡, r) r ∈ Γ𝑖𝑗 (4.5a) 𝑛ˆ · B𝑖 (𝑡, r) = 𝑛ˆ · B 𝑗 (𝑡, r) r ∈ Γ𝑖𝑗 . (4.5b) This is accomplished by ensuring x𝑛𝑖 = x𝑛𝑗 | Γ𝑖𝑗 and defining 𝑛ˆ × H𝑖 (𝑡, r) = Λ(𝑡, r) r ∈ Γ𝑖𝑗 . (4.6) The system of equations to be solved in each sub-region   0.5𝐴¯¯ 1 + 0.25Δ𝑡 𝐴¯¯ 0 𝑋¯ (𝑖),𝑛+1 = ℒ (𝑖) − Λ(𝑖) (𝑖) (𝑖) (4.7) 67 where ℒ (𝑖) = 0.5Δ𝑡 𝐴¯¯ 0 𝑋¯ (𝑖),𝑛   − 0.5𝐴¯¯ 1 + .25Δ𝑡 𝐴¯¯ 0 𝑋¯ (𝑖),𝑛−1 − 0.5Δ𝑡 𝐺˜ (𝑖),𝑛+1 (𝑖) (𝑖) (4.8) − 0.5Δ𝑡 𝐺˜ (𝑖),𝑛−1 − 0.25Δ𝑡 𝐹˜ (𝑖),𝑛+1 + 0.5Δ𝑡 𝐹˜ (𝑖),𝑛 − 0.25𝐹˜ (𝑖),𝑛−1 . It is noted that the time history of the Lagrange multiplier is not explicitly stored, unlike the other currents in the system. Let the unknowns defined on Γ𝑐 be denoted by the subscript 𝑐 and all other unknowns in Ω𝑖 denoted by subscript 𝑟, such that (4.7) can be written as  ¯¯ (𝑖) ¯¯ (𝑖)   ¯ (𝑖),𝑛+1   (𝑖)   (𝑖)  𝐴𝑟𝑟 𝐴𝑟𝑐  𝑋𝑟  ℒ 𝑟  Λ̄𝑟   ¯ (𝑖) ¯ (𝑖)   (𝑖),𝑛+1  =  (𝑖)  −  (𝑖)  (4.9)        ¯ ¯ 𝐴 𝑐𝑟 𝐴 𝑐𝑐  𝑋𝑐 ¯  ℒ 𝑐  Λ̄𝑐         where 𝐴¯¯ = .5𝐴¯¯ 1 + .25Δ𝑡 𝐴¯¯ 0 . The first set of equations can be rewritten as = 𝐴¯¯ 𝑟𝑟 (ℒ 𝑟 − 𝐴¯¯ 𝑟𝑐 𝑋¯ 𝑐 (𝑖),𝑛+1 (𝑖),−1 (𝑖) (𝑖),−1 (𝑖),𝑛+1 (𝑖) 𝑋¯ 𝑟 − Λ̄𝑟 ) (4.10) which can be substituted into the second set of equations   𝐴¯¯ 𝑐𝑐 − 𝐴¯¯ 𝑐𝑟 𝐴¯¯ 𝑟𝑟 𝐴¯¯ 𝑟𝑐 𝑋¯ 𝑐 (𝑖) (𝑖),−1 (𝑖),−1 (𝑖) (𝑖),𝑛+1 (𝑖) = ℒ𝑐   (4.11) − 𝐴¯¯ 𝑐𝑟 𝐴¯¯ 𝑟𝑟 (𝑖) (𝑖),−1 (𝑖) (𝑖) (𝑖) ℒ 𝑟 − Λ̄𝑟 − Λ̄𝑐 . Define a signed Boolean matrix 𝐵¯¯ 𝑟 that maps unknowns in Ω𝑖 , excluding those on (𝑖) Γ𝑐 , to a set of unique global unknowns defined on ∪Γ𝑖𝑗 . The definition of 𝐵¯¯ 𝑟 is such that (𝑖) Í𝑁𝑣 ¯¯ (𝑖) (𝑖) 𝐵 𝑥¯ = 0, which enforces (4.5a) and (4.5b). Furthermore, the transpose 𝐵¯¯ 𝑟 forms (𝑖),𝑇 𝑖=1 𝑟 a map from the global set of boundary unknowns to those of a particular volume. A similar unsigned Boolean matrix can be defined for unknowns that lie on Γ𝑐 within Ω𝑖 , 𝐵¯¯ 𝑐 , which maps to a set of unique global unknowns. Using 𝐵¯¯ 𝑟 on (4.10) and summing over all volumes results in 𝐾¯¯ 𝑟𝑟 Λ𝑟 = ℒ̄ 𝑟 − 𝐾¯¯ 𝑟𝑐 𝑋¯ 𝑐 (4.12) 68 where 𝑁𝑣 𝐾¯¯ 𝑟𝑟 = 𝐵¯¯ 𝑟 𝐴¯¯ 𝑟𝑟 𝐵¯¯ 𝑟 Õ (𝑖) (𝑖),𝑇 (4.13a) 𝑖=1 𝑁𝑣 𝐾¯¯ 𝑟𝑟 = 𝐵¯¯ 𝑟 𝐴¯¯ 𝑟𝑟 𝐴¯¯ 𝑟𝑐 𝐵¯¯ 𝑐 Õ (𝑖) (𝑖) (𝑖),𝑇 (4.13b) 𝑖=1 𝑁𝑣 𝐵¯¯ 𝑟 𝐴¯¯ 𝑟𝑟 ℒ 𝑟 Õ (𝑖) (𝑖),−1 (𝑖) ℒ̄ 𝑟 = (4.13c) 𝑖=1 and 𝑋¯ 𝑐 the degree of freedom vector for corner unknowns. The boundary condition (4.5a) is enforced at the corners of Ω𝑖 as 𝐵¯¯ 𝑐 is defined so that (𝑖) Í𝑁𝑣 ¯¯ (𝑖),𝑇 (𝑖) 𝐵 𝑖=1 𝑐 𝜆 𝑐 = 0. Using that definition, (4.11) becomes 𝐾¯¯ 𝑐𝑐 𝑋¯ 𝑐𝑛+1 = ℒ̄ 𝑐 + 𝐾¯¯ 𝑐𝑟 Λ𝑟 . (4.14) where 𝑁𝑣 𝐾¯¯ 𝑐𝑐 = 𝐵¯¯ 𝑐 (𝐴¯¯ 𝑐𝑐 − 𝐴¯¯ 𝑐𝑟 𝐴¯¯ 𝑟𝑟 𝐴¯¯ 𝑟𝑐 )𝐵¯¯ 𝑐 Õ (𝑖) (𝑖) (𝑖),−1 (𝑖),−1 (𝑖),−1 (𝑖),𝑇 (4.15a) 𝑖=1 𝑁𝑣 𝐾¯¯ 𝑐𝑟 = 𝐵¯¯ 𝑐 𝐴¯¯ 𝑐𝑟 𝐴¯¯ 𝑟𝑟 𝐵¯¯ 𝑟 Õ (𝑖) (𝑖) (𝑖),−1 (𝑖),𝑇 (4.15b) 𝑖=1 𝑁𝑣 𝐵¯¯ 𝑐 (ℒ 𝑐 − 𝐴¯¯ 𝑐𝑟 𝐴¯¯ 𝑟𝑟 ℒ 𝑟 ). Õ (𝑖) (𝑖) (𝑖) (𝑖),−1 (𝑖) ℒ̄ 𝑐 = (4.15c) 𝑖=1 As in the vector wave equation domain decomposition formulations, 𝐾¯¯ 𝑐𝑐 is sparse and generally much smaller than 𝐾¯¯ 𝑟𝑟 . Therefore, if solving the boundary problem directly, it is possible to define   𝑋¯ 𝑐𝑛+1 = 𝐾¯¯ −1 𝑐𝑐 ℒ̄ 𝑐 + 𝐾 𝑐𝑟 Λ𝑟 ¯¯ (4.16) and then use substitution into (4.12) to yield 𝒦¯¯𝑟𝑟 Λ𝑟 = ℒ̄ 𝑟 − 𝐾¯¯ 𝑟𝑐 𝐾¯¯ −1 𝑐𝑐 ℒ̄ 𝑐 (4.17)   where 𝒦¯¯𝑟𝑟 = 𝐾¯¯ 𝑟𝑟 + 𝐾¯¯ 𝑟𝑐 𝐾¯¯ −1 ¯¯ 𝑐𝑐 𝐾 𝑐𝑟 . If using an iterative solver, (4.12) and (4.14) can be solved as a coupled system instead. 69 4.3.2 MFETI-DP2 In the previous formulation, a unique Lagrange multiplier was defined for each edge and face on Γ. This formulation enforces field continuity though a Neumann-like boundary condition. A more robust enforcement of field continuity is to use a Robin boundary condition [106]. The boundaries between Ω𝑖 and Ω 𝑗 , Γ𝑖𝑗 and Γ 𝑗𝑖 , are distinct and have a unique Lagrange multiplier defined for each subdomain. Instead of (4.6), the boundary condition on the Lagrange multipliers becomes Λ(𝑖) (𝑡, r) = 𝑛ˆ 𝑖 × H(𝑡, r) − 𝜂−1 𝑛ˆ 𝑖 × 𝑛ˆ 𝑖 × E(𝑡, r) ∈ Γ𝑖𝑗 . (4.18) The edges on Γ𝑖𝑗 , as well as the interior unknowns and unknowns on 𝜕Ω𝑖 are included (𝑖) in a vector 𝑋¯ 𝐼 . The Lagrange multipliers for the corner edges remain the same, with a unique Lagrange multiplier for each feature. The faces on an Γ𝑖𝑗 are grouped together with the corner edges, which we will denote as 𝑋¯ 𝑁 . The system of equations solved in each subdomain is  ¯¯ (𝑖) ¯¯ (𝑖)   ¯ (𝑖),𝑛+1   (𝑖)    𝐴𝐼𝐼 𝐴𝐼𝑁   𝑋𝐼  ℛ 𝐼  .5Δ𝑡 Λ̄(𝑖),𝑛+1  𝐼   (𝑖),𝑛+1  =  (𝑖)  −          ∈ Ω𝑖 (4.19)  ¯ (𝑖) ¯ (𝑖)  (𝑖) 𝐴¯ 𝑁 𝐼 𝐴¯ 𝑁 𝑁  𝑋¯ 𝑁   ℛ 𝑁   Λ̄ 𝑁          where given ℛ = .5Δ𝑡 𝐴¯¯ 0 𝑋¯ 𝑛 − (.5𝐴¯¯ 1 + .25Δ𝑡 𝐴¯¯ 0 )𝑋¯ 𝑛−1 − .5Δ𝑡 𝐺˜ 𝑛+1 − .5Δ𝑡 𝐺˜ 𝑛−1 (4.20) − .25Δ𝑡 𝐹˜ 𝑛+1 + .5Δ𝑡 𝐹˜ 𝑛 − .25𝐹˜ 𝑛−1 − .5Δ𝑡 Λ𝑛 + .25Λ𝑛−1 , (𝑖) ℛ 𝐼 contains ℛ associated with interior unknowns, external unknowns, and unknowns of (𝑖) edges on Γ𝑖𝑗 and ℛ 𝑁 contain the remaining unknowns in subdomain Ω𝑖 . To derive the set of equations to solve for the dual unknowns, first consider the summation of (4.18) for two volumes Λ(𝑖) (𝑡, r) + Λ(𝑗) (𝑡, r) + (𝑌𝑖 + 𝑌𝑗 )𝑛ˆ 𝑖 × 𝑛ˆ 𝑖 × E(𝑡, r). (4.21) Define an unsigned Boolean matrix such that 𝐵¯¯ 𝐼 𝑋¯ 𝐼 = 𝑒¯𝐼 , the unknowns for the electric (𝑖) (𝑖) field on Γ𝑖𝑗 . Also, define the map from the boundary edges on Γ𝑖𝑗 to its counterpart on Γ 𝑗𝑖 , 70 𝑇¯𝑖𝑗 . Testing (4.21) with curl-conforming basis functions W𝑒 , generalizing the case to any number of subdomains, and using 𝐵¯¯ 𝐼 and 𝑇¯¯𝑗𝑖 restrict and map the unknowns of Ω 𝑗 to Ω𝑖 (𝑖) yields 𝑇¯¯𝑖𝑗 (Λ(𝑗) − (𝑌𝑖 + 𝑌𝑗 ) 𝑀 ¯¯ (𝑗) 𝐵¯¯ (𝑗) 𝑋¯ (𝑗) ) = 0 Õ Λ(𝑖) + 𝑇 𝐼 𝐼 (4.22) 𝑗∈𝜎(𝑖) where 𝜎(𝑖) are all subdomains Ω 𝑗 which share a boundary with Ω𝑖 and the matrix 𝑀 ¯¯ (𝑖) h𝑛ˆ × W𝑒 , 𝑛ˆ × W𝑒 i| . 𝑇 Γ𝑖𝑗 Rewrite the first set of equations in (4.19) as = 𝐴¯¯ 𝐼𝐼 (ℛ 𝐼 − .5Δ𝑡 Λ̄𝐼 − 𝐴¯¯ 𝐼𝑁 𝑋¯ 𝑁 (𝑖),𝑛+1 (𝑖),−1 (𝑖) (𝑖),𝑛+1 (𝑖) (𝑖),𝑛+1 𝑋¯ 𝐼 ) (4.23) Now, (4.23) can be plugged into (4.22) to yield 𝑇𝑖𝑗 (𝐼¯¯(𝑗) − (𝑌𝑖 + 𝑌𝑗 ) 𝑀 ¯¯ (𝑗) 𝐵¯¯ (𝑗) 𝐴¯¯ (𝑗),−1 𝐵¯¯ (𝑗),𝑇 )Λ(𝑗) Õ Λ(𝑖)+ 𝑇 𝐼 𝐼𝐼 𝐼 𝑗∈𝜎(𝑖) (4.24) 𝑇𝑖𝑗 𝐵¯¯ 𝐼 𝐴¯¯ 𝐼𝐼 (𝐴¯¯ 𝐼𝑁 𝐵¯¯ 𝑁 𝑋¯ 𝑁 − ℛ 𝐼 ). Õ (𝑗) (𝑗),−1 (𝑗) (𝑗) (𝑗) = 𝑗∈𝜎(𝑖) In order to complete the definition of the equation to solve for Λ𝐼 over ∪Γ𝑖𝑗 , we first define the equation for 𝑋¯ 𝑁 , the unknowns associated with faces on Γ𝑖𝑗 and edges on Γ𝑐 . The derivation of the equation to solve for 𝑋¯ 𝑁 follows a similar path as 𝑋¯ 𝑐 from MFET-DP1. First, substitute (4.23) into the second set of equations in (4.22) to yield   𝐴¯¯ 𝑁 𝑁 − 𝐴¯¯ 𝑁 𝐼 𝐴¯¯ 𝐼𝐼 𝐴¯¯ 𝐼𝑁 𝑋¯ 𝑁 (𝑖) (𝑖),−1 (𝑖),−1 (𝑖) (𝑖),𝑛+1 (𝑖) = ℛ𝑁   (4.25) − 𝐴¯¯ 𝑁 𝐼 𝐴¯¯ 𝐼𝐼 (𝑖) (𝑖),−1 (𝑖) (𝑖) (𝑖) ℛ 𝐼 − Λ̄𝐼 − Λ̄𝑁 . The Boolean matrix 𝐵¯¯ 𝑁 is defined as (𝑖)  ¯¯ (𝑖) 𝐵 0   ¯ 𝐵¯ 𝑁 =  𝑟, 𝑓 (𝑖)  (4.26)  0 ¯ ¯ (𝑖)  𝐵 𝑐   to enforce (4.5b) for the faces on Ω𝑖 and (4.5a) edges on corners where 𝐵¯¯ 𝑟, 𝑓 is the portion of (𝑖) 𝐵¯¯ 𝑟 that acts on the face degrees of freedom. Using a similar process as (4.14), the equation (𝑖) for 𝐾¯¯ 𝑁 𝑁 𝑋¯ 𝑁 𝑛+1 = ℛ̄ 𝑁 + 𝐾¯¯ 𝑁 𝐼 Λ𝐼 . (4.27) 71 where 𝑁𝑣 𝐾¯¯ 𝑁 𝑁 = 𝐵¯¯ 𝑁 (𝐴¯¯ 𝑁 𝑁 − 𝐴¯¯ 𝑁 𝐼 𝐴¯¯ 𝐼𝐼 𝐴¯¯ 𝐼𝑁 )𝐵¯¯ 𝑁 Õ (𝑖) (𝑖) (𝑖), (𝑖),−1 (𝑖),−1 (𝑖),𝑇 (4.28a) 𝑖=1 𝑁𝑣 𝐾¯¯ 𝑁 𝐼 = 𝐵¯¯ 𝑁 𝐴¯¯ 𝑁 𝐼 𝐴¯¯ 𝐼𝐼 𝐵¯¯ 𝐼 𝐵¯¯ 𝐿 Õ (𝑖) (𝑖) (𝑖),−1 (𝑖),𝑇 (𝑖),𝑇 (4.28b) 𝑖=1 𝑁𝑣 𝐵¯¯ 𝑁 (ℛ 𝑁 − 𝐴¯¯ 𝑁 𝐼 𝐴¯¯ 𝐼𝐼 ℛ 𝐼 ). Õ (𝑖) (𝑖) (𝑖) (𝑖),−1 (𝑖) ℛ̄ 𝑁 = (4.28c) 𝑖=1 The global system of equations for the edge Lagrange multipliers Define a Boolean matrix 𝐵¯¯ 𝐿 that maps the local Lagrange multipliers to a global vector. The Boolean matrix (𝑖) is applied to (4.24) and discretized in time to define the equation 𝐾¯¯ 𝐼𝐼 Λ𝐼 = ℛ̄ 𝐼 − 𝐾¯¯ 𝐼𝑁 𝑁¯ 𝑁𝑛+1 (4.29) where 𝑁𝑣 𝐾¯¯ 𝐼𝐼 = 𝐼 + Õ (𝑖) 𝐵𝐿 𝑖=1 (4.30a) 𝑇¯¯𝑖𝑗 (𝐼¯¯(𝑗) − 𝑀 𝑇 𝐵¯¯ 𝐼 𝐴¯¯ 𝐼𝐼 )𝐵¯¯ 𝐼 𝐵 𝐿 Õ (𝑗) (𝑗) (𝑗),−1 (𝑗),𝑇 (𝑗),𝑇 𝑗∈𝜎(𝑖) 𝑁𝑣 𝐾¯¯ 𝐼𝑁 = − 𝑇 𝑗𝑖 𝑀 𝑇 𝐵 𝐼 𝐴¯¯ 𝐼𝐼 𝐴¯¯ 𝐼𝑁 𝐵 𝑁 Õ Õ (𝑗) (𝑗) (𝑗),−1 (𝑗) (𝑗),𝑇 (𝑖) 𝐵𝐿 (4.30b) 𝑖=1 𝑗∈𝜎(𝑖) 𝑁𝑣 ¯¯ (𝑗) 𝐵(𝑗) 𝐴¯¯ (𝑗),−1 ℛ (𝑖) ) Õ (𝑗) Õ ℛ̄ 𝐼 = − 𝐵𝐿 𝑇 𝑗𝑖 𝑀 𝑇 𝐼 𝐼𝐼 𝐼 (4.30c) 𝑖=1 𝑗∈𝜎(𝑖) It can be solved as a coupled system with (4.27) or in a Schur complement system similar to (4.17) using 𝒦¯¯𝐼𝐼 Λ𝐼 = ℒ̄ 𝐼 − 𝐾¯¯ 𝐼𝑁 𝐾¯¯ −1 𝑁 𝑁 ℒ̄ 𝑁 . (4.31)   where 𝒦¯¯𝐼𝐼 = 𝐾¯¯ 𝐼𝐼 − 𝐾¯¯ 𝐼𝑁 𝐾¯¯ −1𝑁𝑁 𝑁𝐼 𝐾¯¯ . However, 𝐾¯¯ 𝑁 𝑁 is much larger than its MFETI-DP1 counterpart 𝐾¯¯ 𝑐𝑐 and will take more time and memory to compute and store its inverse. 72 Table 4.1: Average number of unknowns per volume and ratio of interior to boundary unknowns for different decompositions of Ω𝑟𝑒 𝑐𝑡 . MFETI-DP1 MFETI-DP2 𝑁𝑣 avg unk. 𝛼 avg unk. 𝛼 2 1815 N/A 1755 6.75 5 776 1.71 705 1.00 10 404 0.60 351 0.34 25 170 0.15 133 0.08 50 89 0.05 60 0.03 4.3.3 Discussion The proposed domain decomposition formulations reduces the cost of solving for the fields by defining a set of equations on fictitious interior boundaries of the true problem. The definitions of these equations uses small matrix inverses, and therefore are generally not sparse. Furthermore, those matrix inverses have to be stored in memory. At the very least, one matrix inverse for each subdomain is needed, making the memory storage more than solving for the original system. However, it is expected that for direct inverses, inverting either the coupled equations of (4.12) and (4.14) or (4.29) and (4.27) would be more efficient than solving the original problem. If storing the inverse of the matrix associated with Γ𝑐 is possible, (4.17) and (4.31) are even smaller in size. The following discussion uses a sample problem to further probe properties of MFETI-DP1 and MFETI-DP2. In particular, we examine the performance of the method with a sample geometry with several configurations of subdivision. Define Ω𝑟𝑒 𝑐𝑡 as a rectangular region in free space with dimension .25 m × .2 m × .3 m discretized with average edge length of 5.5 cm. The region is the split into subdomains in several configurations listed in Table 4.1 so that MFETI-DP1 and MFETI-DP2 can be applied, where 𝛼 is the ratio of interior unknowns to the number of boundary unknowns. This table is provided to provide the connection between the number of subdomains and 𝛼 for the test example, as 𝛼 will be more extensible to other geometries than a set number of subdomains. Figure 4.2 shows the condition number of the matrix formed by the coupled 73 equations (4.12) and (4.14) for MFETI-DP1, (4.29) and (4.27) for MFETI-DP2, as well as the Schur complement forms 𝒦¯¯𝑟𝑟 and 𝒦¯¯𝐼𝐼 with different numbers of subdivisions compared to the condition number of of the monolithic system matrix 𝐴. ¯¯ Solving the coupled equations results in condition number that are generally at or worse than the monolithic solve. Using the Schur complement results in condition numbers comparable to the monolithic solve for MFETI-DP1 and several orders of magnitude smaller for MFETI-DP2. Typically, MFETI-DP1 and MFETI-DP2 when used with iterative solvers converge faster than the monolithic solve. Consider a domain Ω𝑟𝑒 𝑐𝑡 is illuminated by an incident field defined by ¯ 𝑟 , 𝑡) = 𝑦ˆ cos 2𝜋 𝑓0 𝜏 𝑒 −(𝜏−8𝜎)2 /2𝜎2 (V/m). 𝐸(¯  (4.32) The bandwidth of the modulated Gaussian is 𝜎 = 3(𝜋 𝑓𝑏𝑤 , where the bandwidth 𝑓𝑏𝑤 = 𝑓𝑚𝑎𝑥 − 𝑓0 with maximum frequency 𝑓𝑚𝑎𝑥 = 500MHz and center frequency 𝑓0 = 300MHz. The time step size is 0.333 ns which is ten times larger than what could be taken by a similar FDTD solver or leapfrog method for MFEM. Figure 4.3 shows the iteration counts for the coupled solves for MFETI-DP1 and MFETI-DP2 compared to the monolithic solve using GMRES with a tolerance of 10−6 and restarted after 20 iterations when Ω𝑟𝑒 𝑐𝑡 is illuminated by an incident field. With a diagonal preconditioner, both domain decomposition approaches perform better than the monolithic solve, with MFETI-DP2 performing best overall. Next, as alluded to earlier, the principal advantage of Newmark-𝛽 schemes is the unconditional stability. As a result, it is important to verify that the both MFETI-DP1 and MFETI-DP2 retain this feature. To ascertain this, we examine the eigenspectra of these formulations to understand the stability of the time marching scheme. For the time marching scheme used in this work to be stable, the matrix used to evolve the system in time must have eigenvalues Re(𝜆) >= 0 for the system to be unconditionally stable. Figure 4.4 shows that for solving both the coupled system and using the Schur complement approach, the domain decomposition methods will be stable. From the insert in Figure 4.4, the effect of using the schur complement shifts the eigenvalues of 𝒦¯¯𝑟𝑟 and 𝒦¯¯𝐼𝐼 away from 74 Figure 4.2: Condition number of system matrix. Figure 4.3: Convergence of history of GMRES for monolithic, MFETI-DP1, and MFETI-DP2. the imaginary axis, though MFETI-DP2 is shifted further away than MFETI-DP1. 4.4 Results In this section, we present numerical results that show how the proposed algorithms perform compared to a non-DDM approach. Additionally, we will demonstrate that as the 75 (b) MFETI-DP matrix eigenvalues using Schur (a) MFETI-DP matrix eigenvalues. complement. Figure 4.4: System matrix eigenvalues. fields are almost identical to machine precision, the overall trajectory of the particles are the identical to machine precision as well, i.e., the approach does not induce additional error. 4.4.1 Field Comparison The first example shows performance of the methods and difference in the fields. In this example, a 0.25𝑚 × 0.2𝑚 × 1𝑚 box is illuminated an incident field defined by ¯ 𝑟 , 𝑡) = 𝑦ˆ cos 2𝜋 𝑓0 𝜏 𝑒 −(𝜏−8𝜎)2 /2𝜎2 (V/m). 𝐸(¯  (4.33) The bandwidth of the modulated Gaussian is 𝜎 = 3(𝜋 𝑓𝑏𝑤 , where the bandwidth 𝑓𝑏𝑤 = 𝑓𝑚𝑎𝑥 − 𝑓0 with maximum frequency 𝑓𝑚𝑎𝑥 = 595MHz and center frequency 𝑓0 = 295MHz. The time step size 𝛿 𝑡 = (150 𝑓0 )−1 and the experiment is run for 1500 time steps. The boundary conditions on the external boundaries are robin boundary conditions. Three approaches are compared: monolithic (no DDM), MFETI-DP1, and MFETI-DP2 using three levels of discretization with average edge lengths of ℎ1 = 5.20cm, ℎ2 = 3.63cm, and ℎ 3 = 2.59cm. The DDM meshes were split using METIS, such that no volumes have a regular shape. Figure 4.5, shows that the fields obtained through domain decomposition for edge lengths ℎ1 , ℎ2 , and ℎ 3 respectively. Regardless of the number of subdivisions, 76 Table 4.2: Ratio of interior to boundary unknowns for different decompositions for particle beam test. MFETI-DP1 MFETI-DP2 𝑁𝑣 𝛼 runtime (s) 𝛼 runtime (s) 5 3.000 571.4 1.893 694.1 10 0.955 614.8 0.652 711.9 25 0.257 466.7 0.181 488.6 50 0.094 562.5 0.077 507.6 the methods remain accurate to near machine precision with respect to the monolithic solve. This implies that using either domain decomposition formulation would cause the same acceleration on particles as the fields that are similar to each other within machine precision. The next result demonstrates the effect of MFETI-DP on simulation time. As seen in Figure 4.6, both formulations are faster than the original monolithic solve for any number of subdivisions. However, there is an optimal ratio of interior volume unknowns to boundary unknowns for the fastest runtime. For each case, there is a crossover point at which further subdividing Ω incurs a cost in solving the boundary problem that is greater than the savings by reducing the unknowns in each volume. 4.4.2 Particle Beam The next two examples demonstrate MFETI-DP1 and MFETI-DP2 with particles. In these experiments, the fields that are obtained through domain decomposition are used to push the particles, which are not subject to any decomposition scheme. The first example is a particle beam traveling through a cylindrical PEC cavity. This case is a standard test to show charge conservation as errors readily appear as striations in the particle distribution. Furthermore, approximate solutions can be used to further verify the result. The first test is done on a cavity that is 10 cm in length and 0.02 cm in radius disctretized with an average edge length of 5 mm resulting in roughly 21,000 unknowns. The ten particles at each time step enter at the bottom face of the cavity with a radius of 8 mm. The beam voltage and 77 (a) Geometry with average edge length ℎ 1 . (b) Geometry with average edge length ℎ2 . (c) Geometry with average edge length ℎ 3 . Figure 4.5: Relative error in electric field for MFETI-DP. 78 (a) Geometry with average edge length ℎ 1 . (b) Geometry with average edge length ℎ2 . (c) Geometry with average edge length ℎ 3 . Figure 4.6: Timing data for MFETI-DP. 79 current are 500 V and 50 mA respectively, causing the an initial velocity 1.326 × 104 km/s. The geometry was subdivided in four configurations, the details of which are listed in Table 4.2. The experiment was run for 1000 time steps at a time step size of 1 ns, which for the monolithic case took 1,342 seconds. As shown in figure 4.7, all the simulations maintained charge conservation over the course of the run. There is no notable difference in the simulations, regardless of the number of subdivisions used. Seven particles were tracked over the course of the run to track the particle trajectories with the domain decomposition. Because the fields are virtually equivalent, the particle trajectories are consistent across the different geometries. Figure 4.8 shows the error in the tracked particles’ path with rep sect to the path in the monolithic solve. The next result was a particle beam in a larger PEC cavity. In this case, the PEC cylinder was 40 cm in length and 10 cm in radius. The geometry has an average edge length of 7.38 mm and a total of 218,059 unknowns. The beam voltage and current are 3.2 kV and 50 mA respectively, resulting in an initial velocity of 3.35 × 104 km/s. The experiment was run for both MFETI-DP1 and MFETI-DP2 using the iterative solver GMRES using a tolerance of 10−6 . The particle distribution is compared to a well-validated axialsymmetric EM-FDTDPIC code, XOOPIC [134] and shows good agreement in Figure 4.9. 4.4.3 Plasma Ball The last example presented is the case of an adiabatic expansion of a plasma ball using MFETI-DP2. This example takes a considerable number or degrees of freedom to resolve the debye length and allow the particles to expand for an extended period of time. However, this can be made more tractable by using either domain decomposition schemes and unconditionally stable time marching to evolve the system with larger time steps than a similar EM-FDTDPIC formulation would allow. Additionally, this experiment has analytic solutions [108] that can be used for validation. A charge neutral, Gaussian distribution of 24000 particles is set at the center of a spherical region 24 cm in diameter. The 𝑆𝑟 + ions are 80 (a) MFETI-DP1 (b) MFETI-DP2 Figure 4.7: Error in Gauss’ law for particle beam. at 1K and the electrons are at 100K with a density of 5 × 108 particles per cubic meter. The outer boundary has an impedance boundary condition and the boundary is sufficiently far from the particles that none leave during the simulation. The average edge length is 1 cm with a total of 681,214 unknowns. Figure 4.10 shows good agreement of the particle distribution with the analytic result at 150 ns, 300 ns, and 800 ns. 4.5 Summary In this work, we presented two formulations for domain decomposition using non- overlapping subdomains for the mixed finite element method. The methods were used to 81 (a) Difference in particle path for MFETI-DP1. (b) Difference in particle path for MFETI-DP2. Figure 4.8: Error in tracked particle positions. create a MFETI-PIC scheme which allows for computing the fields due to particle motion for larger problems. For sufficiently large problems, this results in time savings without loss of accuracy. All the results shown were computed on a single node; however, the method can be made parallel to take full advantage of high performance computing resources. Doing so suggest creating a domain decomposition framework for particle mapping and pushing consistent with the field solver. Future work will focus on creating a qausi-helmholtz decomposition for the FETI- framework. This ensures that the nullspace of the field solver does not corrupt the satisfaction of charge conservation. This is a key issue for using the time domain vector 82 (a) XOOPIC vs MFETI-DP1. (b) XOOPIC vs MFETI-DP2. Figure 4.9: Particle beam expansion for enlarged PEC cavity. 83 Figure 4.10: Adiabatic plasma ball expansion after 150, 300, and 800 ns. 84 wave equation approach which has a nullspace that grows linearly in time. The combination of domain decomposition and qausi-Helmholtz decomposition will allow for solving even larger problems as the number of degrees of freedom is reduced to those associated with edges, rather than have both edges and faces as presented here. 85 CHAPTER 5 CONCLUSION AND FUTURE WORK The focus of this work is to define what is needed to create charge conserving EM- FEMPIC schemes that are accurate and efficiently simulate large problems. These capabili- ties are needed to simulate real problems in a reasonable time frame. The mathematics needed to conserve charge in an electromagnetic finite element particle-in-cell method were examined first. The rubrics provide a rationale behind the particular choice of basis functions in this work and the literature, as well as the needed insight to use other time marching schemes while satisfying all of Maxwell’s equations. A connection was drawn between three Maxwell solver approaches: FDTD, EB-MFEM, and DH-MFEM. The insight provided by the rubrics has been applied throughout this thesis, as well as in related works. Counterexamples were provided to show the consequences of failing to follow the rubrics. Examples also showed that the DH-based method yields results comparable to FDTD and EB-based EM-FEMPIC schemes. Next, a charge conserving higher order EM-FEMPIC scheme was demonstrated. Whit- ney basis functions for nodes, edges, faces, and volumes were used to ensure both of Gauss’ laws were satisfied. This is possible because they satisfy the de-Rham complex. The number of degrees of freedom needed to obtain a certain level of accuracy was shown to be reduced for higher order approaches. The results showed that the increased accuracy of the fields improved the accuracy of the particle trajectories. Finally, two domain decomposition methods for Maxwell solvers were presented and applied to EM-FEMPIC. The two methods followed the finite element tearing and interconnect approach, utilizing both Neumann and Robin boundary conditions for the interior boundaries. The fields obtained by both approaches were shown to give the same fields as solving the system without the fictitious divisions. It was also demonstrated that both approaches decreased simulation time compared to the standard monolithic 86 approach. Results with particles showed that introducing the domain decomposition did not affect the satisfaction of Gauss’ laws. Furthermore, as the fields obtained were the same as the original method, the particle trajectories remained unchanged while obtained in a shorter period of time. Several challenges and advances remain with the work presented within this thesis. First, the techniques applied in the latter chapters for EB-based Maxwell solvers can be similarly applied for DH-based solvers. More work can be done to better understand the benefits and drawbacks between the two approaches. The impact of different particle shape functions on the methods presented in this work also remains an open topic. Though it does not affect charge conservation, other factors such as noise would be impacted which may influence the performance of designed devices. However, the bulk of future work is built upon expanding on the domain decomposition techniques presented. The nullspaces of the proposed approaches should be studied more in depth. Furthermore, it inherently creates a multiply connected system that would benefit from an appropriate qausi-Helmholtz decomposition. Not only would it ensure that Gauss’ law is satisfied exactly and not corrupted by any nullspace, it also enables using the vector wave equation as a field solver. This would further reduce the computational cost as the system of equations would no longer need to solve for the magnetic flux and it could be obtained through a forward calculation instead. Other preconditioning schemes and direct methods for domain decomposition are also highly desirable. 87 APPENDICES 88 APPENDIX A SUPPLEMENTARY INFORMATION FOR "RUBRICS FOR CHARGE CONSERVING CURRENT MAPPING IN FINITE ELEMENT PARTICLE IN CELL METHODS" A.0.1 Calculus on Distribution functions The following must be understood in the sense of distributions. Specifically, we consider derivatives and integral of the delta function. These appear in (2.8) onwards where one represents the spatial variation of PDF using 𝑆(r − r𝑝 (𝑡)) = 𝑆(r) ★𝑠 𝛿 r − r𝑝 (𝑡) . Calculus on  delta function is critical for developing insight into implementing equations of continuity and Gauss’ law on the grid. To do so, we start with a transform to k-space to prove (2.11) and (2.13). A.0.1.1 Time derivatives ∫𝑡 𝑑𝜏v𝑝 (𝜏), the time derivative of 𝛿 r − r𝑝 (𝑡) can be written as  Using r𝑝 (𝑡) = r0𝑝 + 0  ∫ 𝑡  𝜕𝑡 𝛿 r − r𝑝 (𝑡) = 𝜕𝑡 𝛿(r − r0𝑝 ) ★ 𝛿 r − 𝑑𝜏v𝑝 (𝜏)  0  ∫ 𝑡  (A.1) = 𝛿(r − r0𝑝 ) ★ 𝜕𝑡 𝛿 r − 𝑑𝜏v𝑝 (𝜏) 0 Next, we expliot ∫ ∞ 1 𝛿(r − a) = 𝑑3 k𝑒 𝑗k·(r−a) (A.2) (2𝜋)3 −∞ to yield  ∫ 𝑡  ∫ ∞  ∫𝑡  1 𝑗k· r− 𝑑𝜏v𝑝 (𝜏) 𝜕𝑡 𝛿 r − 𝑑𝜏v𝑝 (𝜏) = 𝑑 3 k𝜕𝑡 𝑒 0 0 (2𝜋)3 −∞  ∫ 𝑡  (A.3) = −∇ · v𝑝 (𝑡)𝛿 r − 𝑑𝜏v𝑝 (𝜏) . 0 Therefore  ∫ 𝑡  𝜕𝑡 𝛿(r − r𝑝 (𝑡)) = −𝛿(r − r𝑜𝑖 ) ★ ∇ · v𝑝 (𝑡)𝛿 r − 𝑑𝜏v𝑝 (𝑡) 0 (A.4)   = −∇ · v𝑝 (𝑡)𝛿 r − r𝑝 (𝑡) 89 A.0.2 Time integrals Next, starting from the time integral of (2.11) and using the definition of current density in equation (2.9b), ∫ 𝑡 ∫ 𝑡  ∫ 𝜏  0 0 𝑑𝜏∇ · v𝑝 (𝜏)𝑆(r − r𝑝 (𝜏)) = 𝑆(r) ★𝑠 𝑑𝜏∇ · v𝑝 (𝜏)𝛿 r − 𝑑𝜏 v𝑝 (𝜏 ) (A.5) 0 0 0 Again, using (A.2) and 𝜕𝑡 r̃𝑝 (𝑡) = v𝑝 (𝑡), ∫ 𝑡  ∫ 𝜏  𝑑𝜏∇ · v𝑝 (𝜏)𝛿 r − 𝑑𝜏0v𝑝 (𝜏0) 0 0 𝑡 ∞ 𝑑𝜏 ∫ ∫ ∫𝜏 𝑗k·(r− 0 𝑑𝜏0 v𝑝 (𝜏0 )) = 3 ∇· 𝑑 kv𝑝 (𝜏)𝑒 3 0 (2𝜋) −∞ ∫ ∞ ∫ 𝑡 ∫𝜏 (A.6) 1 𝑗k·r −𝑗k· 0 𝑑𝜏0 v𝑝 (𝜏0 ) =∇· 𝑑 3 k𝑒 𝑑𝜏v 𝑝 (𝜏)𝑒 (2𝜋)3 −∞ 0 ∫ ∞ ∫ r𝑝 (𝑡) 1 𝑗k·r =∇· 𝑑 k𝑒 3 𝑑r̃𝑝 𝑒 −𝑗k·r̃𝑝 (2𝜋)3 −∞ r𝑝 (0) Interchanging the order of the integrals in the above equation results in ∫ 𝑡 ∫ r𝑝 (𝑡) 𝑑𝜏∇ · v𝑝 (𝜏)𝑆(r − r𝑝 (𝜏)) = ∇ · 𝑑r̃𝑆 (r − r̃) (A.7) 0 r𝑝 (0) This equation can be used to define the charge density in terms of the particle’s path; it can also be used to define the time integral of the current density. A.0.3 Maxwell Solvers The following sections are presented purely for completeness and to conform with material presented earlier. We present details of the three Maxwell solvers used here; FDTD, EB-MFEM, and DH-MFEM. A.0.3.1 FDTD Define the rooftop function Λ𝑖 (r) ˆ | 𝑖·r| 0 < | 𝑖ˆ · r| ≤ Δ𝑖  1 −    Δ𝑖 Λ𝑖 (r) = (A.8)  0  o.w.  90 where 𝑖 = 𝑥, 𝑦, 𝑧. The basis functions to discretize Maxwell’s equations for FDTD are 1-form := u𝑖 Λ𝑖+1 (r − l𝑚 )Λ𝑖+2 (r − l𝑚 ) (A.9a) 2-form := u𝑖 Λ𝑖 (r − f𝑚 ) (A.9b) with u𝑖 = 𝑥, ˆ 𝑦, ˆ 𝑧ˆ , l𝑚 is the midpoint of each edge, and f𝑚 is the midpoint of each face. Testing Faraday’s law with the normal to cell faces and Ampere’s law is with the edge unit vectors yields ˆ 𝜕𝑡 𝜇H(𝑡, r)i = −hu𝑖 𝛿((r − f𝑚 ) · 𝑖), hu𝑖 𝛿((r − f𝑚 ) · 𝑖), ˆ ∇ × E(𝑡, r)i 𝜕𝑡 H𝑖 (𝑡, f𝑚 ) = −𝜇−1 0 u 𝑖 · ∇ × E(𝑡, f𝑚 ) (A.10a) = 𝜇−1 −1 0 𝜕𝑖+1 E 𝑖+2 (𝑡, f𝑚 ) − 𝜇0 𝜕𝑖+2 E 𝑖+1 (𝑡, f𝑚 ) ˆ 𝜕𝑡 𝜀E(𝑡, r)i = −hu𝑖 𝛿((r − l𝑚 ) · 𝑖), hu𝑖 𝛿((r − f𝑚 ) · 𝑖), ˆ ∇ × H(𝑡, r)i − hu𝑖 𝛿((r − l𝑚 ) · 𝑖), ˆ J(𝑡, r)i 𝜕𝑡 E𝑖 (l𝑚 ) = −𝜀0−1 u𝑖 · ∇ × H(𝑡, l𝑚 ) − 𝜀0−1 u𝑖 · J(𝑡, l𝑚 ) = 𝜀0−1 𝜕𝑖+1 H𝑖+2 (𝑡, l𝑚 ) − 𝜀0−1 𝜕𝑖+2 H𝑖+1 (𝑡, l𝑚 ) − 𝜀0−1 J𝑖 (𝑡, l𝑚 ). (A.10b) The addition on index 𝑖 follows cycle {𝑥, 𝑦, 𝑧}. This testing is equivalent to the usual presentation of the FDTD recurrence relation. For example, consider 𝑖 = 𝑥 in (A.10a) would yield 𝑛+ 12 𝑛− 12 Δ𝑡 𝑛 𝑛 𝐻 =𝐻 − (𝐸 𝑧,(𝑖,𝑗+1,𝑘+ 1 − 𝐸 1 ) 𝑥,(𝑖,𝑗+ 21 ,𝑘+ 21 ) 𝑥,(𝑖,𝑗+ 12 ,𝑘+ 12 ) 𝜇Δ 𝑦 2 ) 𝑧,(𝑖,𝑗,𝑘+ 2) (A.11) Δ𝑡 + (𝐸 𝑛 − 𝐸 𝑛𝑦,(𝑖,𝑗+ 1 ,𝑘) ). 𝜇Δ𝑧 𝑦,(𝑖,𝑗+ 12 ,𝑘+1) 2 Likewise, 𝑖 = 𝑥 in (A.10b) would yield 𝑛+1 𝑛 Δ𝑡 𝑛+ 1 𝑛+ 1 𝐸 𝑥,(𝑖+ = 𝐸 𝑥,(𝑖+ + (𝐻 2 1 1 − 𝐻 2 1 1 ) 1 2 ,𝑗,𝑘) 1 2 ,𝑗,𝑘) 𝜀Δ 𝑦 𝑧,(𝑖+ 2 ,𝑗+ 2 ,𝑘) 𝑧,(𝑖+ 2 ,𝑗− 2 ,𝑘) Δ𝑡 𝑛+ 1 𝑛+ 12 − (𝐻 2 1 − 𝐻 ) (A.12) 𝜀Δ𝑧 𝑦,(𝑖+ 2 ,𝑗,𝑘+ 12 ) 𝑦,(𝑖+ 12 ,𝑗,𝑘− 12 ) Δ𝑡 𝑛+1/2 − 𝐽 . 𝜀 𝑥,(𝑖+ 21 ,𝑗,𝑘) 91 A.0.3.2 EB-MFEM Formulation Í 𝑁𝑒 (1) Í𝑁 𝑓 (2) In this formulation, E(𝑡, r) ≈ 𝑖=1 𝑖 𝑒 (𝑡)W𝑖 (r) and and B(𝑡, r) ≈ 𝑏 (𝑡)W𝑖 (r). 𝑖=1 𝑖 Fara- day’s law is measured with the 2-form basis function and Ampere’s law with the 1-form basis function to yield 𝜕𝑡 b = [ℒ∇× ◦]e (A.13a) 𝜕𝑡 e = 𝑐 2 [ℒ∇× ◦]𝑇 [★𝜇−1 ]b − 𝜀−1 j (A.13b) √ assuming the domain is freespace with the speed of light 𝑐 = 𝜀0 𝜇0 −1 . The Hodge operator matrix is defined as ∫ (2) (2) [★𝜇−1 ]𝑖𝑗 = 𝑑ΩW𝑖 (r) · W 𝑗 (r). (A.14) Ω The degrees of freedom for the current are defined as 𝑁𝑝 r𝑛+1 𝑞𝑛 Õ ∫ (1) j𝑖 = 𝑑p𝑛 · W𝑖 (r) (A.15) Δ𝑡 r𝑛 𝑛=1 where p𝑛 is the particle path for the 𝑛-th particle. A.0.3.3 DH-MFEM Formulation Í 𝑁𝑒 (1) Í𝑁 𝑓 (2) In this formulation, H(𝑡, r) ≈ 𝑖=1 ℎ 𝑖 (𝑡)W𝑖 (r) and D(𝑡, r) ≈ 𝑖=1 𝑑 𝑖 (𝑡)W𝑖 (r). Ampere’s law is measured with the 2-form basis function and Faraday’s law with the 1-form basis function to yield 𝜕𝑡 h = 𝑐 2 [ℒ∇× ◦][★𝜀−1 ]d (A.16a) 𝜕𝑡 d = [ℒ∇× ◦]h − j. (A.16b) The Hodge operator matrix is defined as ∫ (2) (2) [★𝜀−1 ]𝑖𝑗 = 𝑑ΩW𝑖 (r) · W 𝑗 (r). (A.17) Ω The degrees of freedom for the current are defined as 𝑁𝑝 r𝑛+1 𝑞𝑛 Õ ∫ j𝑖 = 𝑑p𝑛 𝑛ˆ 𝑖 · J𝑛 (r) (A.18) Δ𝑡 r𝑛 𝑛=1 92 where the particle path for the 𝑛-th particle and is only evaluated on the surface. © 2021 IEEE. Reprinted, with permission, from Crawford, Zane D., Scott O’Connor, John Luginsland, and B. Shanker. "Rubrics for Charge Conserving Current Mapping in Finite Element Electromagnetic Particle in Cell Methods." IEEE Transactions on Plasma Science (2021). 93 APPENDIX B HIGHER ORDER BASIS DEFINITIONS In this appendix, we define the higher order interpolatory basis functions. The basis functions are formed by multiplying the lowest order Whitney elements by an interpolatory Lagrange polynomial. We use the definition by Silvester [?]  Î𝑖−1  𝑖!1 𝑚=0 (𝑘𝜆 − 𝑚), 1 ≤ 𝑖 ≤ 𝑘    𝑃𝑖 (𝜆) = (B.1) 𝑖=0   1,   and the shifted Silvester Polynomial 1 Î𝑖−1  𝑚=1 (𝑘𝜆 − 𝑚), 2 ≤ 𝑖 ≤ 𝑘 + 1    (𝑖−1)!  𝑃ˆ 𝑖 (𝜆) = (B.2) 𝑖=1   1,   B.0.1 Node Basis Function (0) The k-th order nodal basis function 𝑊𝑠 (r) ∈ 𝐻 1 is defined as (0) ˜ 𝑖 (r) 𝑊𝑖 (r) = 𝑄 (0)𝑊 (B.3) where 𝑊˜ 𝑠(0) (r) = 𝜆 𝑠 (B.4) and 𝑄 (0) = 𝑃ˆ 𝑠 (𝜆 𝑠 )𝑃𝑡 (𝜆𝑡 )𝑃𝑢 (𝜆𝑢 )𝑃𝑣 (𝜆𝑣 ). (B.5) There are a total of (𝑘 + 1)(𝑘 + 2)(𝑘 + 3)/6 degrees of freedom in each tetrahedron with 4 associated with nodes, 𝑘 − 1 for each edge, (𝑘 − 1)(𝑘 − 2) for each face, and 𝑘(𝑘 − 1)(𝑘 − 2)/6 internal to each tetrahedron. 94 B.0.2 Whitney Edge Basis Function (1) The k-th order edge basis function W𝑠 (r) ∈ 𝐻(∇×) is defined as (1) (1) W𝑠 (r) = 𝑄 (1) W̃𝑠𝑡 (r) (B.6) where (1) W̃𝑠 (r) = 𝜆 𝑠 ∇𝜆𝑡 − 𝜆𝑡 ∇𝜆 𝑠 (B.7) and 𝑄 (1) = 𝑃ˆ 𝑠 (𝜆 𝑠 )𝑃ˆ 𝑡 (𝜆𝑡 )𝑃𝑢 (𝜆𝑢 )𝑃𝑣 (𝜆𝑣 ). (B.8) There are a total of 𝑘(𝑘 +2)(𝑘 +3)/2 degrees of freedom in each tetrahedron with 𝑘 associated with each edge, 𝑘(𝑘 − 1)/2 with each face, and 𝑘(𝑘 − 1)(𝑘 − 2)/6 internal to each tetrahedron. B.0.3 Whitney Face Basis Function (2) The k-th order face basis function W𝑠 (r) ∈ 𝐻(∇·) is defined as (2) (2) W𝑠 (r) = 𝑄 (2) W̃𝑠 (r) (B.9) where (2) W̃𝑠 (r) = 𝜆 𝑠 ∇𝜆𝑡 × ∇𝜆𝑢 + 𝜆𝑡 ∇𝜆𝑢 × ∇𝜆𝑣 + 𝜆𝑢 ∇𝜆 𝑠 × ∇𝜆𝑡 (B.10) and 𝑄 (2) = 𝑃ˆ 𝑠 (𝜆 𝑠 )𝑃ˆ 𝑡 (𝜆𝑡 )𝑃ˆ 𝑢 (𝜆𝑢 )𝑃𝑣 (𝜆𝑣 ). (B.11) There are a total of 𝑘(𝑘 + 1)(𝑘 + 3)/2 degrees of freedom in each tetrahedron with 𝑘(𝑘 + 1)/2 with each face and 𝑘(𝑘 − 1)(𝑘 + 1)/6 internal to each tetrahedron. B.0.4 Volume Basis Function The k-th order volume basis function 𝑊 ˜ (3) (r) ∈ 𝐿2 is defined as (3) ˜ (3) (r) 𝑊𝑠 (r) = 𝑄 (3)𝑊 (B.12) 95 where 𝑊˜ (3) (r) = 𝜆 𝑠 ∇𝜆𝑡 · ∇𝜆𝑢 × 𝜆𝑣 + 𝜆𝑡 ∇𝜆𝑢 · ∇𝜆𝑣 × 𝜆 𝑠 (B.13) + 𝜆𝑢 ∇𝜆𝑣 · ∇𝜆 𝑠 × 𝜆𝑡 + 𝜆𝑣 ∇𝜆 𝑠 · ∇𝜆𝑡 × 𝜆𝑢 and 𝑄 (3) = 𝑃ˆ 𝑠 (𝜆 𝑠 )𝑃ˆ 𝑡 (𝜆𝑡 )𝑃ˆ 𝑢 (𝜆𝑢 )𝑃ˆ 𝑣 (𝜆𝑣 ). (B.14) The function is associated with tetrahedra with 𝑘(𝑘 + 1)(𝑘 + 2)/6 unknowns in each tetrahedron, all of which are internal to the cell. 96 APPENDIX C UNCONDITIONALLY STABLE TIME STEPPING METHOD FOR MIXED FINITE ELEMENT MAXWELL SOLVERS C.1 Introduction Efficient solutions to Maxwell’s equations is of significant interest to problems in fields ranging from microwave engineering to optics to radar to remote sensing to plasma physics [135, 136, 137, 138, 139]. As a result, it has seen extensive development over the past four decades. So much so that there exist reliable commercial simulators that are trusted for design and analysis; to wit, most microwave engineering and antenna design is carried out on a laptop before being fabricated and tested. Despite these advances, there are a number of open problems and challenges, largely brought on by our desire to improve the fidelity of simulation for increasingly complex geometric and material layouts. Adding temporal dimension to the analysis increases the challenges as well. The methods used for solution to Maxwell’s equations are based either on differential or integral equations; this paper will focus on the former. Over the years, differential equation based methods have become the mainstay of computational electromagnetic tools; starting with the classical Yee-algorithm for direct discretization of Maxwell’s equations [18] to finite element methods (FEMs) [32, 38]. A number of challenges had to be overcome to make FEMs robust and viable for use in design and analysis of systems [77, 106]. Over the years, research on FEM has focused on developing methods so as to increase the accuracy and reduce computational costs; these include higher order interpolatory basis functions [99], higher order hierarchical basis sets [?], and domain decomposition methods [140]. These techniques have largely been applied to solving the vector wave equation in the frequency domain. In the time domain, the 97 temporal discretization follows a Newmark-Beta prescription [141, 142]. With the right choice of parameters, it has been shown that this prescription is unconditionally stable. The time step size is not intimately tied to the smallest mesh feature, but chosen to capture the physics to a prescribed order. However, the solutions to the vector wave equation are prone to inaccuracies arising from a non-trivial null space that corresponds to a function of the form 𝑡∇𝜙(¯𝑟 ), where 𝜙(¯𝑟 ) is a scalar function of space. There has been an effort to use finite elements to directly solve Maxwell’s equations [46, 143, 144], albeit less so than vector wave formulations. These methods follow the prescription of the Yee-algorithm. Faraday’s and Ampere’s law are staggered in time leading to a leap frog time stepping. As a result, they do not have a time-growing null space, are conditionally stable, and energy conserving. Unfortunately, the time step size depends on the size of the features of the mesh and materials defined within it. The resulting time step sizes are very small, significantly more so than needed to represent content at the highest frequency. This is in contrast to the state of the art of solving the vector wave equation. It would be ideal if one could render TD-MFEM provably unconditionally stable (it would enable reaping the advances made for vector wave equation FEM solvers). Earlier work has analyzed the stability for leapfrog time stepping in TD-MFEM [145]. Likewise, unconditionally stable time stepping schemes have been used in TD-MFEM, but without any rigorous analysis to prove stability [146]. In this work, we will prove that using Newmark- Beta method for TD-MFEM can produce unconditionally stable time stepping schemes that are up to second order accurate in time. With the correct choice in parameters, the time stepping schemes are shown to be energy conserving and non-dissipative. Furthermore, these properties hold for higher-order spatial basis functions used to discretize TD-MFEM. Thus, the specific contributions of this work are (i) a general framework for creating time stepping schemes of up to second order for the TD-MFEM of arbitrary spatial order, (ii) proof of stability for different boundary conditions and error analysis for the time stepping 98 scheme, (iii) analysis of the null space of the discrete system, and (iv) numerical results that verifying the analytic proofs. The rest of the paper is organized as follows: Section C.2 describes the general problem. Section C.3 presents the spatial and temporal discretization of Maxwell’s equations. The section contrasts the Newmark-Beta method to the leapfrog method and provides stability and error analysis of the Newmark-Beta method applied to the MFEM. Section C.4 presents several numerical results verifying properties discussed in previous sections and the overall efficacy of the method. Finally, conclusion and other remarks are provided in Section C.5. C.2 Problem Statement Consider the volume Ω ∈ R3 bounded by Γ shown in figure C.1. In this domain, assume the region has homogeneous, isotropic, and time-independent constitutive parameters, denoted by 𝜀 (F/m) for the permittivity and 𝜇 (H/m) for the permeability. Assume that y x z Figure C.1: Sample Volume Ω with Different Boundary Conditions. Reproduced courtesy of The Electromagnetics Academy the volume is source free, being free of impressed charges or currents. The magnetic ¯ 𝑟 , 𝑡) (A/m) and electric flux density 𝐷(¯ field 𝐻(¯ ¯ 𝑟 , 𝑡) (C/m2 ) are defined in terms of the ¯ 𝑟 , 𝑡) (Wb/m2 ) and electric field 𝐸(¯ magnetic flux density 𝐵(¯ ¯ 𝑟 , 𝑡) (V/m), respectively, through 99 the constitutive parameters such that ¯ 𝑟 , 𝑡) = 𝜇−1 𝐵(¯ 𝐻(¯ ¯ 𝑟 , 𝑡) (C.1a) ¯ 𝑟 , 𝑡) = 𝜀𝐸(¯ 𝐷(¯ ¯ 𝑟 , 𝑡). (C.1b) The electric field and magnetic flux density are governed by Maxwell’s equations ¯ 𝑟 , 𝑡) 𝜕 𝐵(¯ = −∇ × 𝐸(¯ ¯ 𝑟 , 𝑡) 𝜕𝑡 ¯ 𝑟 , 𝑡) 𝜕𝜀𝐸(¯ ¯ 𝑟 , 𝑡) 𝐵(¯ =∇× 𝜕𝑡 𝜇 (C.2) ¯ 𝑟 , 𝑡)) = 0 ∇·(𝜀𝐸(¯ ¯ 𝑟 , 𝑡) = 0. ∇· 𝐵(¯ The boundary Γ may be partitioned into several subdomains, Γ𝑖 , such that Γ := ∪𝑖 Γ𝑖 . Each Γ𝑖 has a boundary condition, either Dirichlet (Γ𝐷 ), Neumann (Γ𝑁 ), or impedance (Γ𝐼 ) [106]. The boundary conditions are defined as ¯ 𝑟 , 𝑡) = 𝚿𝐷 (¯𝑟 , 𝑡) on Γ𝐷 𝑛ˆ × 𝐸(¯ (C.3a) ¯ 𝑟 , 𝑡) 𝐵(¯ 𝑛ˆ × = 𝚿𝑁 (¯𝑟 , 𝑡) on Γ𝑁 (C.3b) 𝜇 ¯ 𝑟 , 𝑡) 𝐵(¯ 𝑛ˆ × ¯ 𝑟 , 𝑡) = 𝚿𝐼 (¯𝑟 , 𝑡) on Γ𝐼 − 𝑌 𝑛ˆ × 𝑛ˆ × 𝐸(¯ (C.3c) 𝜇 where 𝑌 = 𝜀/𝜇 (0) is the surface admittance. The functions Ψ𝐷 (¯𝑟 , 𝑡), Ψ𝑁 (¯𝑟 , 𝑡), and p Ψ𝐼 (¯𝑟 , 𝑡) represent some arbitrary function that is defined by the user. Additionally, we assume initial conditions on the fields, 𝐸(r, ¯ 0) and 𝐵(r, ¯ 0), satisfy (C.2). C.3 MFEM Formulation C.3.1 Variational Formulation Solutions to (C.2) are found by creating variational equations with appropriate function spaces. The suitable function space for the electric field is the Hilbert space 𝐻(𝑐𝑢𝑟𝑙; Ω) = { 𝑢¯ ∈ 𝐿¯2 (Ω); ∇ × 𝑢¯ ∈ 𝐿¯2 (Ω)}, (C.4) 100 which ensures that the tangential component of the field are continuous across tetrahedral faces. Likewise, the Hilbert space 𝐻(𝑑𝑖𝑣; Ω) = { 𝑢¯ ∈ 𝐿¯2 (Ω); ∇ · 𝑢¯ ∈ 𝐿2 (Ω)} (C.5) is used for the magnetic flux density, which permits functions that have normal continuity across tetrahedral faces. The Whitney edge basis function and Whitney face basis function reside in 𝐻(𝑐𝑢𝑟𝑙; Ω) and 𝐻(𝑑𝑖𝑣; Ω), respectively, and can be used to discretize the fields [46]. The variational equations are formed by taking inner products of Faraday’s law with the Whitney face basis function 𝐵¯ ∗ and Ampere’s law with the Whitney edge basis function 𝐸¯ ∗ such that 1 𝜕 ¯ ¯∗ ∫ ∫ 1 𝑑Ω 𝐵·𝐵 =− 𝑑Ω ∇ × 𝐸¯ · 𝐵¯ ∗ , (C.6a) 𝜇 𝜕𝑡 𝜇 𝜕 ∫ ∫ 1 ¯ ¯∗ 𝑑Ω𝜀 𝐸¯ · 𝐸¯ ∗ = 𝑑Ω∇ × 𝐵·𝐸 . (C.6b) 𝜕𝑡 𝜇 C.3.2 Spatial Discretization The spatial discretization is obtained by applying Galerkin’s method to the variational equations. As previously stated, Whitney edge and face elements as seen in [46] define the lowest order representation of 𝐸¯ ∗ and 𝐵¯ ∗ . The electric field is represented using Whitney edge elements and the magnetic flux using Whitney face elements such that 𝑁𝑒 𝑁𝑓 Õ Õ ¯ 𝑟 , 𝑡) = 𝐸(¯ ¯ 1 (¯𝑟 ), 𝐵(¯ 𝑒 𝑖 (𝑡)𝑊 ¯ 𝑟 , 𝑡) = ¯ 2 (¯𝑟 ), 𝑏 𝑖 (𝑡)𝑊 (C.7) 𝑖 𝑖 𝑖=1 𝑖=1 where 𝑁𝑒 is the number of edges and 𝑁 𝑓 is the number of faces. It is natural to use higher order basis function for more accuracy. We use the interpolatory higher order basis functions defined in [99] for Whitney edge elements and for Whitney face elements defined in [?]. We define vectors 𝑒¯ = [𝑒1 (𝑡), 𝑒2 (𝑡), ..., 𝑒 𝑁𝑒 (𝑡)]𝑇 and 𝑏¯ = [𝑏 1 (𝑡), 𝑏2 (𝑡), ..., 𝑏 𝑁 𝑓 (𝑡)]𝑇 that have the degrees of freedom 𝑒 𝑖 (𝑡) and 𝑏 𝑖 (𝑡) as elements. Using the defined expansion for the field and flux quantities, applying the Galerkin method to (C.6) yields the semi-discrete 101 Maxwell’s equations 𝜕M𝜇−1 𝑏¯ = −M𝑐 𝑒¯ 𝑑𝑡 (C.8) 𝜕M𝜀 𝑒¯ = M𝑇𝑐 𝑏¯ − M𝐼 𝑒¯ + 𝑗¯𝑁 + 𝑗¯𝐼 . 𝑑𝑡 Here matrices are defined as ∫ [M𝜀 ]𝑖,𝑗 = 𝑑Ω𝜀𝑊 ¯ 1 (¯𝑟 ) · 𝑊¯ 1 (¯𝑟 ), (C.9a) 𝑖 𝑗 Ω ∫ 1 ¯2 ¯ 2 (¯𝑟 ), [M𝜇−1 ]𝑖,𝑗 = 𝑑Ω 𝑊 𝑖 (¯𝑟) · 𝑊 𝑗 (C.9b) Ω 𝜇 ∫ 1 ¯2 ¯ 1 (¯𝑟 ), [M𝑐 ]𝑖,𝑗 = 𝑑Ω 𝑊 𝑖 (¯ 𝑟) · ∇ × 𝑊 𝑗 (C.9c) Ω 𝜇 ∫ 1 ¯2 [M𝑇𝑐 ]𝑖,𝑗 = 𝑑Ω 𝑊 𝑗 (¯ 𝑟) · ∇ × 𝑊 ¯ 1 (¯𝑟 ), 𝑖 (C.9d) Ω 𝜇 ∫ [M𝐼 ]𝑖,𝑗 = 𝑑Γ𝑌 𝑛ˆ × 𝑊 ¯ 1 (¯𝑟 ) · 𝑛ˆ × 𝑊 ¯ 1 (¯𝑟 ). (C.9e) 𝑖 𝑗 Γ𝐼 The elements of the coefficient vectors 𝑗¯𝑁 and 𝑗¯𝐼 defined as ∫ [ 𝑗¯𝑁 ]𝑖 (𝑡) = 𝑑Γ𝑁 𝑊 ¯ 1 (¯𝑟 ) · Ψ𝑁 (¯𝑟 , 𝑡) 𝑖 ∫Γ𝑁 (C.10) [ 𝑗¯𝐼 ]𝑖 (𝑡) = ¯ 1 (¯𝑟 ) · Ψ𝐼 (¯𝑟 , 𝑡) 𝑑Γ𝐼 𝑊 𝑖 Γ𝐼 and represent the boundary current due to Neumann and impedance boundary conditions. Another interpretation of the TD-MFEM is through discrete exterior calculus. The Galerkin discrete Hodge matrices M𝜀 and M𝜇−1 are defined as [80]. The discrete curl operator D𝑐 is defined as 1 −1 D𝑐 = M M𝑐 . (C.11) 𝜇 𝜇−1 The discrete curl operator and its transpose for the lowest order maps edges of a tetrahedron to the faces it bounds, as well as preserves the direction of the curl on that face [79]. This follows from the concept of the exterior derivative in the generalized Stokes theorem from exterior calculus. Equation (C.11) holds for all orders, however only has a simple mapping of faces to edges for the lowest order. 102 C.3.3 Temporal Discretization The semi-discrete Maxwell’s equations, derived in the previous section, are fully discretized by choosing an appropriate temporal discretization. Discretization schemes exist that solve the electric field and magnetic flux density concurrently or in a staggered fashion. In a staggered scheme, the solutions to Faraday’s law and Ampere’s law lie at different time steps, which allows a different time stepping method for each equation. A concurrent time stepping scheme couples the equations into a single system, written in the form  ¯  I 0  𝜕𝑡 𝑏   0 D𝑐  𝑏¯   0           +  = , (C.12)    𝑇  𝜕𝑡 𝑒¯  −D𝑐 M𝜇−1 M𝐼   𝑒¯   𝑗¯𝑁 + 𝑗¯𝐼        0 M 𝜀           where I is the identity matrix. The single time derivative allows for the choice of a single time stepping scheme and implies the solved quantities lie at the same time step. The following subsections discuss the leapfrog method (for sake of completeness) and the Newmark-Beta method and their applications to finding solutions to the semi-discrete Maxwell’s equations. Note, that time step sizes for leapfrog (Δ𝑡,𝐿𝐹 ) and Newmark-Beta (Δ𝑡,𝑁𝛽 ) will be different. C.3.3.1 Leapfrog Method In the leapfrog method, the time derivatives in Faraday’s law and Ampere’s law are discretized using central finite differences staggered at half-time step intervals. Applying this to (C.8) yields 𝑏¯ 𝑛+ 2 = 𝑏¯ 𝑛− 2 − Δ𝑡,𝐿𝐹 D𝑐 𝑒¯ 𝑛 1 1 (C.13a) 𝑛+ 12 𝑛+ 12 (M𝜀 + M𝐼 )¯𝑒 𝑛+1 = M𝜀 𝑒¯ 𝑛 + Δ𝑡,𝐿𝐹 D𝑇𝑐 M𝜇−1 𝑏¯ 𝑛+ 2 + Δ𝑡,𝐿𝐹 ( 𝑗¯𝑁 1 + 𝑗¯𝐼 ). (C.13b) The time stepping scheme consists of an explicit update in Faraday’s law and a matrix inversion in Ampere’s law as the Galerkin Hodge matrix is sparse, but not diagonal. It is 103 conditionally stable [147, 148] provided time step size satisfies 2 Δ𝑡,𝐿𝐹 ≤ q . (C.14) 𝑇 𝜌(M−1𝜀 D 𝑐 M𝜇−1 D 𝑐 ) When the time step size satisfies the stability requirement, the time stepping scheme also conserves energy. The mesh that discretizes the volume Ω determines the maximum time step size. The spectral radius of the matrices tends to grow as elements in the mesh shrink, limiting the time step size in relation to the maximum frequency the mesh could resolve. C.3.3.2 Spectral radius One drawback of this time-stepping scheme is that its stability is tied to the discretization of the mesh and the dielectric and magnetic material within Ω. As the electrical size of the elements shrink, so does the maximum time step size, regardless of coarser features of the mesh. Table C.1: Leapfrog maximum time-step vs mesh edge lengths 𝜆𝑚𝑎𝑥 𝑐Δ𝑡,𝐿𝐹 ℎ 𝑚𝑎𝑥 ℎ 𝑎𝑣 𝑔 ℎ 𝑚𝑖𝑛 1.0587E+03 6.1468E-02 5.0000E-01 3.2400E-01 1.9376E-01 2.5661E+03 3.9481E-02 3.2773E-01 2.1477E-01 1.3689E-01 6.4293E+03 2.4943E-02 2.2165E-01 1.4365E-01 7.8437E-02 9.8165E+03 2.0186E-02 1.6865E-01 1.0925E-01 5.7826E-02 1.2769E+02 1.7699E-01 1.4142E+00 1.0424E+00 8.5782E-01 2.3979E+02 1.2916E-01 1.5326E+00 8.8836E-01 4.2409E-01 1.4249E+03 5.2983E-02 1.4142E+00 5.6001E-01 1.8142E-01 1.1535E+04 1.8622E-02 1.4142E+00 2.2648E-01 6.2023E-02 In Table C.1, the relationship between the edge lengths ℎ, spectral radius 𝜆𝑚𝑎𝑥 , and time step sizes Δ𝑡,𝐿𝐹 for the leapfrog scheme are presented. The spectral radius increases as the edge lengths shrink, resulting in a smaller maximum time step. Instead, it is preferred that time step is dictated by the bandwidth of the signals present in the simulation. In the results, it will show that the spatial errors tend to dominate error far more than temporal error, especially for the lowest order spatial basis functions. 104 C.3.3.3 Newmark-Beta The Newmark-Beta method is a concurrent time stepping method to solve second order differential equations [149, 150], and is a popular method of choice for solving TD-FEM. To begin, the first order differential equation from (C.12) is recast as a second order differential equation, A2 𝜕𝑡2 𝑢¯ + A1 𝜕𝑡 𝑢¯ + A0 𝑢¯ + 𝑓¯ = 0, (C.15) where  ¯ 𝑏  0   𝑢¯ =   , 𝑓¯ =    ,  𝑒¯  −( 𝑗¯𝐼 + 𝑗¯𝑁 )  (C.16)     I 0  0      D𝑐  A2 = 0, A1 =   , A0 =  .   𝑇  0 M 𝜀  −D𝑐 M𝜇−1 M𝐼      Applying a 1-dimensional finite element method in time and identifying particular weights in the formulation results in a three-step recurrence relation that uses two parameters to yield different accuracy and stability properties [151]. The three-step recurrence relation for (C.16) is defined as [𝛾A1 + 𝛽Δ𝑡,𝑁𝛽 A0 ]𝑢¯ 𝑛+1 1 + [(1 − 2𝛾)A1 + ( + 𝛾 − 2𝛽)Δ𝑡,𝑁𝛽 A0 ]𝑢¯ 𝑛 2 (C.17) 1 𝑛−1 + [(𝛾 − 1)A1 + ( − 𝛾 + 𝛽)Δ𝑡,𝑁𝛽 A0 ]𝑢¯ 2 1 1 + (𝛽Δ𝑡,𝑁𝛽 ) 𝑓¯ 𝑛+1 + ( + 𝛾 − 2𝛽)Δ𝑡,𝑁𝛽 𝑓¯ 𝑛 + ( − 𝛾 + 𝛽)Δ𝑡,𝑁𝛽 𝑓¯ 𝑛−1 = 0 2 2 which has the degrees of freedom for the electric field and magnetic flux density at the same timestep. The parameters 𝛾 and 𝛽 affect the accuracy and stability properties. The following sections prove and demonstrate that with appropriate choice of parameters, the TD-MFEM system with Dirichlet, Neumann, and impedance boundary conditions is not only unconditionally stable, but also energy conserving. 105 C.3.4 Stability Analysis In order to determine the stability of the Newmark-Beta time-stepping scheme, the zero-input response of the system is examined using a z-transform. For the system to be stable, the roots of the z-transformed scheme must lie on or within the unit circle in the complex plane. The spectral theorem is used on (C.17) such that the matrix A−1 1 A0 is replaced by its eigenvalue 𝜆 [106, 152]. Taking the z-transform yields the quadratic equation  1 [𝛾 + 𝛽Δ𝑡,𝑁𝛽 𝜆]𝑧 2 + [(1 − 2𝛾) + ( + 𝛾 − 2𝛽)Δ𝑡,𝑁𝛽 𝜆]𝑧 2  (C.18) 1 + [(𝛾 − 1) + ( − 𝛾 + 𝛽)Δ𝑡,𝑁𝛽 𝜆] 𝑢(𝑧)¯ =0 2 with zeros q 1 − 2𝛾 + (1/2 + 𝛾 − 2𝛽)Δ𝑡,𝑁𝛽 𝜆 1 + (1 − 2𝛾)Δ𝑡,𝑁𝛽 𝜆 + [(1/2 + 𝛾)2 − 4𝛽](Δ𝑡,𝑁𝛽 𝜆)2 𝑧= ± . 2𝛾 + 2𝛽Δ𝑡,𝑁𝛽 𝜆 2𝛾 + 2𝛽Δ𝑡,𝑁𝛽 𝜆 (C.19) The parameters 𝛾 and 𝛽 must be chosen to have the desired accuracy and stability characteristics for a given time step of size Δ𝑡,𝑁𝛽 . For a time stepping scheme to be unconditionally stable, |𝑧| ≤ 1 for all Δ𝑡,𝑁𝛽 . This leaves the characterization of the eigenvalue 𝜆 of the coupled Maxwell system, which will be shown to have a nonnegative real part. C.3.4.1 Eigenanalysis of System with Dirichlet and Impedance Boundary Conditions The zeros of the stability equation values depend on the parameters 𝛾 and 𝛽, the time step size Δ𝑡,𝑁𝛽 , and the eigenvalues of the matrix system. First, consider a system with no impedance boundary conditions. The coupled Maxwell system can be analyzed by forming the generalized eigenvalue system D𝑐  𝑏¯   ¯   ¯ 𝑏   0 𝑏    A−1 1 A0  =    =𝜆 , (C.20)  𝑒¯  −M−1 𝑇 0   𝑒¯   𝑒¯            𝜀 D 𝑐 M𝜇−1   106 which leads to ¯ D𝑐 𝑒¯ = 𝜆𝑏, (C.21a) 𝑇 ¯ − M−1 𝜀 D𝑐 M𝜇−1 𝑏 = 𝜆 𝑒¯ . (C.21b) It follows that 𝑇 ¯ 2¯ ¯ D𝑐 M−1 𝜀 D𝑐 M𝜇−1 𝑏 = −𝜆 𝑏 = Λ𝑏 (C.22a) 𝑇 M−1 𝜀 D 𝑐 M𝜇−1 D 𝑐 𝑒¯ = −𝜆 𝑒¯ = Λ¯ 2 𝑒. (C.22b) The eigenvalues 𝜆 of the system matrix A−1 1 A0 are purely imaginary when there are no impedance boundary conditions. This fact can be used to simplify (C.19) to find values 𝛾 and 𝛽 that will result in a stable time stepping scheme. The eigenvalues of the system change when impedance boundary conditions are imposed. The eigenvalues of A−1 1 A0 have nonnegative real parts with impedance boundary conditions. First, the matrices in (C.16) are redefined in a more general form, M𝜇−1 0  0 𝜇−1      0 M 𝑐  A1 =   , A0 =   (C.23)   −1 𝑇   0 M 𝜀  −𝜇0 M𝑐 M 𝐼      where the basis functions can be of any order. To determine the behaviour of A−1 1 A0 , we create a similar matrix A𝑆 = FA𝑆 F−1 , where F = A−1 1 , leading to −1 −1 0 −𝜇−1 M 2 M M 2    0 𝜇 −1 𝑐 𝜀  A𝑠 =  −1 −1 −1 −1 . (C.24) 𝜇−1 𝑇 2 2 2 2   0 M𝜀 M𝑐 M𝜇−1 M 𝜀 M𝐼 M 𝜀   This matrix is a non-symmetric saddle point matrix, whose properties are discussed in Theorem 3.6 in [153]. The matrix A𝑠 can be proven to be positive semi-definite because the matrix M𝐼 is positive semi-definite. Therefore, the eigenvalues of A𝑠 , and by the properties of similar matrices A−1 1 A0 , have nonnegative real parts. This is demonstrated for the zeroth and second order vector basis function for both Dirichlet and impedance boundary conditions in figure C.2. The eigenvalues for the Dirichlet case lie on the imaginary axis, while the impedance boundary condition case is symmetric about the real axis in right half plane. Knowing that the eigenvalues have nonnegative real parts, it is possible to find 107 ×10 10 Eigenvalues of A−1 1 A0 , p=0 ×10 10 Eigenvalues of A−1 1 A0 , p=2 1 4 Impedance B.C. 0.8 Imepdance B.C. 3 Dirichlet B.C. Dirichlet B.C. 0.6 2 0.4 1 0.2 Im( λ) 0 Im( λ) 0 -0.2 -1 -0.4 -2 -0.6 -3 -0.8 -1 -4 0 5 10 15 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 9 Re( λ) ×10 Re( λ) ×10 10 (a) Zeroth Order Basis Functions(p=0) (b) Second Order Basis Functions(p=2) Figure C.2: Eigenvalues of A−1 1 A0 . Reproduced courtesy of The Electromagnetics Academy general regions of stability and conditions for which the Newmark-Beta time-stepping scheme is nondissipative. The scheme is unconditionally stable for 𝛾 ≥ 0.5 and 𝛽 ≥ 𝛾/2. When 𝛾 = 0.5, 𝛽 ≥ 0.25, |𝑧| = 1 for all the possible eigenvalues of A−1 1 A0 , (C.17) forms a nondissipative time-stepping scheme, which is essential in many applications. When 𝛾 = 0.5 and 𝛽 = 0.25 or 𝛽 = 0.5, we obtain a Crank-Nicholson scheme, which is known to be unconditionally stable and second order. Other choices of 𝛾 and 𝛽 lead to disippative, yet still stable time stepping schemes. C.4 Results In this Section, we demonstrate numerically the utility of using the Newmark-Beta time stepping scheme with TD-MFEM. In particular, the resulting time stepping scheme is unconditionally stable while maintaining the convergence and accuracy properties of the more common leapfrog method. C.4.1 Field Convergence To numerically demonstrate convergence, a TEM plane wave is simulated propagating in free space normally incident through a 0.5 m x 0.5 m x 0.5 m cube. The electric field is 108 defined as ¯ 𝑟 , 𝑡) = 𝑦ˆ cos 2𝜋 𝑓0 (𝑡 − 𝑟¯ · 𝑧ˆ /𝑐 − 8𝜎) 𝑒 −(𝑡−¯𝑟 ·𝑧ˆ /𝑐−8𝜎)2 /2𝜎2 (V/m), 𝐸(¯  (C.25) where 𝜎 = 3/[2( 𝑓max − 𝑓0 )] with the center and maximum frequency denoted as 𝑓0 and 𝑓max , respectively. The relative error is defined as ∫ || 𝑑Ω(𝜙¯ 𝑒 𝑥𝑝 (¯𝑟 , 𝑡) − 𝜙¯ 𝑎𝑛𝑎 (¯𝑟 , 𝑡))|| 2 Relative Error = ∫ (C.26) || 𝑑Ω𝜙¯ 𝑎𝑛𝑎 (¯𝑟 , 𝑡)|| 2 ¯ 𝑟 , 𝑡) can be either the electric field or magnetic flux density. The speed of where 𝜙(¯ √ √ propagation is the speed of light 𝑐 = 1/ 𝜀𝜇 = 1/ 𝜀0 𝜇0 (m/s). Dirichlet boundary ¯ 𝑟 , 𝑡) is defined conditions as defined in (C.3) are used on each face of the volume, where 𝐸(¯ ¯ 𝑟 , 𝑡) follows from Maxwell’s equations. in (C.25) and 𝐵(¯ For the spatial convergence tests, the center frequency is 11 MHz, maximum frequency 21 MHz. Figure C.3 shows the convergence properties of the fields with respect to space and basis function order. The convergence rate is the same between the two time stepping schemes, with similar errors as well. This is expected as the spatial discretization is the same. It is also noted that increasing the order of the basis function yields a greater improvement in accuracy than simply decreasing the element sizes. -1.5 -0.5 log10(Relative Error) log10(Relative Error) -3.5 -2.5 Newmark- , p=0 Newmark- , p=0 Newmark- , p=1 Newmark- , p=1 Newmark- , p=2 Newmark- , p=2 leapfrog, p=0 Leapfrog, p=0 -5.5 -4.5 Leapfrog, p=1 Leapfrog, p=1 leapfrog, p=2 Leapfrog, p=2 -7.5 -6.5 -1 -0.75 -0.5 -0.25 -1 -0.75 -0.5 -0.25 log10(hmin) log10(hmin) (a) Electric Field (b) Magnetic Flux Density Figure C.3: Space Convergence. Reproduced courtesy of The Electromagnetics Academy Figure C.4 shows the relative error in the electric field and magnetic flux as the time step size is changed using the Newmark-Beta time stepping scheme with 𝛾 = .5 and 𝛽 = .25 109 -1.5 -0.5 log10(Relative Error) log10(Relative Error) -3.5 -2.5 -5.5 -4.5 Leapfrog max( t ) Leapfrog max( t) Newmark- , p=0 Newmark- , p=0 Newmark- , p=1 Newmark- , p=1 Newmark- , p=2 Newmark- , p=2 -7.5 -6.5 -11 -10 -9 -8 -7 -11 -10 -9 -8 -7 log 10( t ) log 10( t ) (a) Electric Field (b) Magnetic Flux Density Figure C.4: Time Convergence. Reproduced courtesy of The Electromagnetics Academy for different orders of vector basis functions. The domain Ω is a 0.5 m x 0.5 m x 0.5 m cube partitioned into 439 tetrahedra with an average edge length ℎavg = 143.7 mm. As the order of the Newmark-Beta time stepping scheme is the same as the Leapfrog scheme, it behaves similarly as the Newmark-Beta , but only for time step sizes smaller than the maximum time step size that is marked. Increasing the order of basis functions results in a smaller maximum time step size, whereas the unconditionally stable method allows for much larger time step sizes while maintaining the advantages of using higher order basis functions. C.4.2 Computational Cost Comparison In the previous section, the spatial and temporal convergence of the Newmark-Beta time stepping scheme was compared to the leapfrog method. Here we compare the cost of using the Newmark-Beta time marching scheme. The leapfrog method requires requires matrix inversion for (C.13b). The Newmark-Beta method in (C.17) requires matrix inversion as well. The matrices in both (C.13b) and (C.17) are sparse, having 𝒪(𝑁𝑒 ) and 𝒪(𝑁𝑒 + 𝑁 𝑓 ) nonzero elements, respectively. For practical problems, direct inversion is expensive and iterative solvers are used instead. We assume that a Krylov subspace based method like GMRES converges in 𝑁𝑖𝑡𝑒𝑟 steps for both leapfrog and Newmark-Beta. Assuming that the 110 duration of analysis is 𝑇 gives 𝑁𝑡,𝑙 𝑓 and 𝑁𝑡,𝑁𝛽 as the number of time steps for leapfrog and Newmark-Beta, respectively. It follows that runtime scales as 𝒪(𝑁𝑡,𝑙 𝑓 𝑁𝑖𝑡𝑒𝑟 𝑁𝑒 ) for leapfrog and 𝒪 𝑁𝑡,𝑁𝛽 𝑁𝑖𝑡𝑒𝑟 (𝑁𝑒 + 𝑁 𝑓 ) for Newmark-Beta. It is apparent that the increased degrees  of freedom required for a Newmark system could lead to higher costs; however, the time step sizes can be vastly different. Indeed, depending on the geometry (especially when discretization is driven by the need to capture feature sizes), the cost of the increased degrees of freedom can be counterbalanced by the number of time steps taken when there is flexibility in the time step size. Thus, as the ceiling of Δ𝑡,𝐿𝐹 is fixed by the mesh, so is the number of time steps for a fixed duration of time. However, this is not the case for Newmark-Beta. The numerical experiment conducted next illustrates the cost tradeoff. This experiment is shown in figure C.5 and is conducted using three 1 m x 1 m x 1 m cubes with Dirichlet boundary conditions. A normally-incident TEM plane wave is simulated propagating in free space. The electric field is defined in (C.25) with a center frequency of 110 MHz and a maximum frequency of 210 MHz. The Newmark-Beta time stepping scheme was run three times for the same amount of simulated time. Once at the maximum time step size for leapfrog, once at half that time step size, and once with a time step size of 1/20 𝑓𝑚𝑎𝑥 ≈ 238 ps. GMRES was used with a diagonal preconditioner and a tolerance of 10−4 . When Δ𝑡,𝑁𝛽 = 1/20 𝑓𝑚𝑎𝑥 , the runtime for Newmark-Beta was equal to 4.5 4 3.5 log10(Runtime) 3 2.5 LF @ max t 2 Newmark- h avg 8.7 10-2 m h avg 4.9 10-2 m 1.5 h avg 3.0 10-2 m 1 -11 -10.5 -10 -9.5 log10( t) Figure C.5: Runtime Comparison. Reproduced courtesy of The Electromagnetics Academy or faster than leapfrog. While this was a simple example, we expect to reap significant 111 benefits as the model representation drives discretization. C.4.3 Eigenvalues of Space-time matrix In this test, we further verify the stability of the proposed Newmark-Beta time-marching scheme. To determine the stability of a time stepping schemes, we analyze the eigenvalues of a space-time matrix A that satisfies 𝑢¯ 𝑛+1 = A𝑢¯ 𝑛 (C.27) for some input vector 𝑢¯ 𝑛 . The space-time matrix is an amplification matrix for some input and for it to have a bounded output, the absolute value of the eigenvalues must lie on or within the unit circle in the complex domain. For an energy-conserving system, all the eigenvalues lie on the unit circle, signifying that there is no numerical loss in the system. In this example, two different meshes discretize the same 1 m cube. The first mesh is nearly uniformly discretized, with a ratio of ℎmax /ℎmin ≈ 1.5, where ℎ is edge length. The second mesh is not uniformly discretized, in which the edges become smaller near one face. The ratio of the largest edge to smallest edge in the nonuniform mesh is ≈ 8. Figures C.6, C.7, and C.8 compare the eigenvalues of the space-time matrices for both the leapfrog and Newmark-Beta time stepping schemes for varying time step size. The reference time step size Δ𝑡 = 1/(30 𝑓max ). When the time step size of the leapfrog method satisfies the stability criterion, the eigenvalues all lie on the unit circle. However, when the time step size no longer satisfies the stability criterion, the eigenvalues close to −1 leave the unit circle and head toward 0 and −∞. Figure C.7 and C.8 show that the space-time eigenvalues of the Newmark-Beta time stepping scheme with 𝛾 = .5 and 𝛽 = .25 always lie on the unit circle for a wide range of time step sizes, lending credence to it unconditionally stable nature. The eigenvalues of the system when these particular values are chosen are analytically ±1. Figure C.7 uses the uniformly discretized volume and figure C.8 uses the nonuniformly discretized volume, 112 Leapfrog Eigenvalues Leapfrog Eigenvalues 1 1 0.8 0.8 0.25* ∆t 0.6 0.25* ∆t 0.6 0.0625* ∆t 0.4 0.0625* ∆t 0.4 0.2 0.2 Im( λ) Im( λ) 0 0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1 -1 -12 -10 -8 -6 -4 -2 0 2 -20 -15 -10 -5 0 5 Re(λ) Re( λ) (a) uniform airbox (b) nonuniform airbox Figure C.6: Leapfrog Eigenvalues. Reproduced courtesy of The Electromagnetics Academy Figure C.7: Newmark-Beta Eigenvalues for Uniform Airbox. Reproduced courtesy of The Electromagnetics Academy with both showing the behavior with Dirichlet boundary conditions and impedance boundary conditions. As the overall eigenvalue behavior holds as the basis function order increases, similar results are seen in figures C.9 and C.10, which shows the space-time eigenvalues of the same uniform airbox with Dirichlet boundary conditions, but with first and second order basis functions. Closer inspection of the eigenvalues show that there is an error, as the eigenvalues lie outside the unit circle, but it is strictly numerical in nature and related to the condition number of the matrix A−1 1 A0 . Finally, to numerically confirm late-time stability, figure C.11 shows the electric field and magnetic flux density run for 100,000 time steps at a time step size of Δ𝑡,𝑁𝛽 = 1/10 𝑓𝑚𝑎𝑥 113 Figure C.8: Newmark-Beta Eigenvalues for Nonuniform Airbox. Reproduced courtesy of The Electromagnetics Academy or 4.762 ns. There is a noise floor that corresponds to a DC null space. However, the noise does not grow over time, indicating the nullspace does not contain 𝑡∇𝜙(¯𝑟 ) present in TD-FEM. C.5 Summary In this work, it was proven that the Newmark-Beta algorithm can be used to create an unconditionally stable time stepping framework for solving TD-MFEM. We proved stability for systems arising from different boundary conditions, as well as showed convergence and energy conservation. The benefit of being able to use different time step sizes is apparent. In future work, this method will be applied to simulating plasmas. Utilizing higher order basis sets with optimal particle pushing strategies will be investigated. Preconditioning strategies will also be developed to improve the use of iterative solvers for TD-MFEMs. This work was supported by the Department of Energy Computational Science Graduate fellowship under grant DE-FG02-97ER25308. The authors would also supported in part by the HPCC Facility at Michigan State University under Grant NSF ECCS-1408115, and AFOSR Grant RC107859. © 2020 . Reprinted, with permission, from Crawford, Zane, Jie Li, Andrew Christlieb, and Balasubramaniam Shanker. "Unconditionally stable time stepping method for mixed 114 Newmark Beta Eigenvalues, γ = .5, β = .25 1 0.8 0.6 0.4 0.2 Im( λ) 0 -0.2 0.5* δ t 0.1* δ t -0.4 0.033333* δt -0.6 -0.8 -1 -1.5 -1 -0.5 0 0.5 1 1.5 Re( λ) (a) First Order ×10 -5 Newmark Beta Eigenvalues, γ = .5, β = .25 ×10 -4 Newmark Beta Eigenvalues, γ = .5, β = .25 1 4 2 Im( λ) Im( λ) 0 0 -2 -4 0.5* δt 0.5* δt 0.1* δt 0.1* δt 0.033333* δt 0.033333* δt -6 -1 -1.00006 -1.00004 -1.00002 -1 -0.99998 -0.99996 -0.99994 0.9999 0.99995 1 1.00005 1.0001 Re( λ) Re( λ) (b) p=1, Centered at Re(𝜆) = -1 (c) p=1, Centered at Re(𝜆) = 1 Figure C.9: Newmark-Beta Eigenvalues for Uniform Airbox with First Order Basis functions. Reproduced courtesy of The Electromagnetics Academy finite element maxwell solvers." Progress In Electromagnetics Research C 103 (2020): 17-30. 115 Newmark Beta Eigenvalues, γ = .5, β = .25 1 0.8 0.6 0.4 0.2 Im( λ) 0 -0.2 0.5* δ t -0.4 0.1* δ t 0.033333* δt -0.6 -0.8 -1 -1.5 -1 -0.5 0 0.5 1 1.5 Re( λ) (a) Second order ×10 -4 Newmark Beta Eigenvalues, γ = .5, β = .25 ×10 -3 Newmark Beta Eigenvalues, γ = .5, β = .25 1 6 0.8 4 0.6 0.4 2 0.2 Im( λ) 0 Im( λ) 0 -0.2 -2 -0.4 -0.6 -4 0.5* δt 0.5* δt 0.1* δt 0.1* δt -0.8 0.033333* δt 0.033333* δt -6 -1 -1.0006 -1.0004 -1.0002 -1 -0.9998 -0.9996 -0.9994 0.999 0.9995 1 1.0005 1.001 Re( λ) Re( λ) (b) p=2, centered at Re(𝜆) = -1 (c) p=2, centered at Re(𝜆) =1 Figure C.10: Newmark-Beta Eigenvalues for Uniform Airbox with Second Order Basis Functions. Reproduced courtesy of The Electromagnetics Academy 0 0 Experimental Experimental Analytic Analytic -3 -3 0 0 -5 -5 -6 -10 -6 log 10 ( E y ) log 10 ( B x ) -10 -15 -15 -20 -20 -9 -25 -9 -25 -30 -30 -35 0 0.5 1 1.5 -12 -12 -35 0 0.5 1 1.5 -15 -15 0 50 100 150 200 250 300 350 400 450 0 50 100 150 200 250 300 350 400 450 time ( µ s) time ( µ s) (a) Electric Field (b) Magnetic Flux Density Figure C.11: Plane Wave Simulation for 100,000 Time Steps. Reproduced courtesy of The Electromagnetics Academy 116 BIBLIOGRAPHY 117 BIBLIOGRAPHY [1] Gwynne Evans, Jonathan Blackledge, and Peter Yardley. Numerical methods for partial differential equations. Springer Science & Business Media, 2012. [2] George B Arfken and Hans J Weber. Mathematical methods for physicists, 1999. [3] BA Finlayson and LE Scriven. The method of weighted residuals—a review. Appl. Mech. Rev, 19(9):735–748, 1966. [4] Jay P Boris. Relativistic plasma simulation-optimization of a hybrid code. In Proc. Fourth Conf. Num. Sim. Plasmas, pages 3–67, 1970. [5] J-L Vay. Simulation of beams or plasmas crossing at relativistic velocity. Physics of Plasmas, 15(5):056701, 2008. [6] JP Quintenz. Nonuniform mesh diode simulation code. Journal of Applied Physics, 49(8):4377–4382, 1978. [7] Richard Brownell True. Space-charge-limited beam forming systems analyzed by the method of self-consistent fields with solution of Poisson’s equation on a deformable relaxation mesh. PhD thesis, University of Connecticut, 1972. [8] Richard True. A general purpose relativistic beam dynamics code. In AIP Conference Proceedings, volume 297, pages 493–499. American Institute of Physics, 1993. [9] Nikolaos A Gatsonis and Anton Spirkin. A three-dimensional electrostatic particle-in- cell methodology on unstructured delaunay–voronoi grids. Journal of Computational Physics, 228(10):3742–3761, 2009. [10] John P Verboncoeur. Particle simulation of plasmas: review and advances. Plasma Physics and Controlled Fusion, 47(5A):A231, 2005. [11] Charles K Birdsall and A Bruce Langdon. Plasma physics via computer simulation. CRC press, 2018. [12] Roger W Hockney and James W Eastwood. Computer simulation using particles. crc Press, 1988. [13] JJ Watrous, JW Lugisland, and GE Sasser III. An improved space-charge-limited emission algorithm for use in particle-in-cell codes. Physics of Plasmas, 8(1):289–296, 2001. [14] Andrew D. Greenwood, Keith L. Cartwright, John W. Luginsland, and Ernest A. Baca. On the elimination of numerical cerenkov radiation in pic simulations. Journal of Computational Physics, 201(2):665 – 684, 2004. 118 [15] David Smithe, Peter Stoltz, John Loverich, Chet Nieter, and Seth Veitzer. Development and application of particle emission algorithms from cut-cell boundaries in the vorpal em-fdtd-pic simulation tool. In 2008 IEEE International Vacuum Electronics Conference, pages 217–218. IEEE, 2008. [16] J Loverich, C Nieter, D Smithe, S Mahalingam, and P Stoltz. Charge conserving emission from conformal boundaries in electromagnetic pic simulations. Comput Phys Commun, 2009. [17] Body-fitted electromagnetic pic software for use on parallel computers. Computer Physics Communications, 87(1):155 – 178, 1995. Particle Simulation Methods. [18] Kane Yee. Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media. 14(3):302–307. [19] Panayiotis A Tirkas and Constantine A Balanis. Finite-difference time-domain method for antenna radiation. IEEE Transactions on Antennas and Propagation, 40(3):334–340, 1992. [20] Raymond J Luebbers, Karl S Kunz, Michael Schneider, and Forrest Hunsberger. A finite-difference time-domain near zone to far zone transformation (electromagnetic scattering). IEEE Transactions on Antennas and Propagation, 39(4):429–433, 1991. [21] Min Qiu. Analysis of guided modes in photonic crystal fibers using the finite- difference time-domain method. Microwave and Optical Technology Letters, 30(5):327– 330, 2001. [22] Jian-Ming Jin. Theory and computation of electromagnetic fields. John Wiley & Sons, 2011. [23] Ardavan Farjadpour, David Roundy, Alejandro Rodriguez, Mihai Ibanescu, Peter Bermel, John D Joannopoulos, Steven G Johnson, and Geoffrey W Burr. Improving accuracy by subpixel smoothing in the finite-difference time domain. Optics letters, 31(20):2972–2974, 2006. [24] Richard Courant et al. Variational methods for the solution of problems of equilibrium and vibrations. Lecture notes in pure and applied mathematics, pages 1–1, 1994. [25] PL Arlett, AK Bahrani, and OC Zienkiewicz. Application of finite elements to the solution of helmholtz’s equation. In Proceedings of the Institution of Electrical Engineers, volume 115, pages 1762–1766. IET, 1968. [26] P Silvester. A general high-order finite-element analysis program waveguide. IEEE Transactions on Microwave Theory and Techniques, 17(4):204–210, 1969. [27] Bo-nan Jiang, Jie Wu, and Louis A Povinelli. The origin of spurious solutions in computational electromagnetics. Journal of computational physics, 125(1):104–123, 1996. 119 [28] M Hara, T Wada, T Fukasawa, and F Kikuchi. A three dimensional analysis of rf electromagnetic fields by the finite element method. IEEE Transactions on Magnetics, 19(6):2417–2420, 1983. [29] Jie Wu and Bo-nan Jiang. A least-squares finite element method for electromagnetic scattering problems. 1996. [30] BM Azizur Rahman and J Brian Davies. Penalty function improvement of waveguide solution by finite elements. IEEE Transactions on Microwave Theory and Techniques, 32(8):922–928, 1984. [31] Jean-Claude Nédélec. Mixed finite elements in R3 . Numerische Mathematik, 35(3):315– 341, 1980. [32] Zoltan J Cendes. Vector finite elements for electromagnetic field computation. IEEE Transactions on Magnetics, 27(5):3958–3966, 1991. [33] Hassler Whitney. Geometric integration theory. Courier Corporation, 2012. [34] 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. [35] Zheng Lou and Jian-Ming Jin. Modeling and simulation of broad-band antennas using the time-domain finite element method. IEEE transactions on antennas and propagation, 53(12):4099–4110, 2005. [36] Dan Jiao, A Arif Ergin, Balasubramaniam Shanker, Eric Michielssen, and Jian-Ming Jin. A fast higher-order time-domain finite element-boundary integral method for 3-d electromagnetic scattering analysis. IEEE Transactions on Antennas and Propagation, 50(9):1192–1202, 2002. [37] Su Yan and Jian-Ming Jin. A fully coupled nonlinear scheme for time-domain modeling of high-power microwave air breakdown. IEEE Transactions on Microwave Theory and Techniques, 64(9):2718–2729, 2016. [38] Jin-Fa Lee, R. Lee, and A. Cangellaris. Time-domain finite-element methods. 45(3):430– 442. [39] GB Jacobs and Jan S Hesthaven. Implicit–explicit time integration of a high-order particle-in-cell method with hyperbolic divergence cleaning. Computer Physics Communications, 180(10):1760–1767, 2009. [40] Jiefu Chen and Qing Huo Liu. Discontinuous galerkin time-domain methods for multiscale electromagnetic simulations: A review. Proceedings of the IEEE, 101(2):242– 254, 2012. [41] Stephen D Gedney, John C Young, Tyler C Kramer, and J Alan Roden. A discontinuous galerkin finite element time-domain method modeling of dispersive media. IEEE Transactions on Antennas and Propagation, 60(4):1969–1977, 2012. 120 [42] Annalisa Buffa and Ilaria Perugia. Discontinuous galerkin approximation of the maxwell eigenproblem. SIAM Journal on Numerical Analysis, 44(5):2198–2226, 2006. [43] Zheng Lou and Jian-Ming Jin. A novel dual-field time-domain finite-element domain- decomposition method for computational electromagnetics. IEEE transactions on antennas and propagation, 54(6):1850–1862, 2006. [44] Marinos N Vouvakis, Zoltan Cendes, and Jin-Fa Lee. A fem domain decomposition method for photonic and electromagnetic band gap structures. IEEE Transactions on Antennas and Propagation, 54(2):721–733, 2006. [45] Charbel Farhat and Francois-Xavier Roux. A method of finite element tearing and interconnecting and its parallel solution algorithm. International journal for numerical methods in engineering, 32(6):1205–1227, 1991. [46] A. Bossavit. Whitney forms: a class of finite elements for three-dimensional compu- tations in electromagnetism. 135(8):493–500. [47] Alain Bossavit and Lauri Kettunen. Yee-like schemes on staggered cellular grids: A synthesis between fit and fem approaches. IEEE Transactions on Magnetics, 36(4):861– 867, 2000. [48] R Hiptmair. Finite elements in computational. Acta Numerica 2002: Volume 11, (11):237, 2002. [49] RN Rieben, GH Rodrigue, and DA White. A high order mixed vector finite element method for solving the time dependent maxwell equations on unstructured grids. Journal of Computational Physics, 204(2):490–519, 2005. [50] 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. [51] 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. [52] 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. 194:43–53. [53] Michael Kraus, Katharina Kormann, Philip J Morrison, and Eric Sonnendrücker. Gempic: Geometric electromagnetic particle-in-cell methods. Journal of Plasma Physics, 83(4), 2017. [54] Barry Marder. A method for incorporating gauss’ law into electromagnetic pic codes. Journal of Computational Physics, 68(1):48 – 55, 1987. 121 [55] Jan Deca, Andrey Divin, Giovanni Lapenta, Bertrand Lembège, Stefano Markidis, and Mihály Horányi. Electromagnetic particle-in-cell simulations of the solar wind interaction with lunar magnetic anomalies. 112(15):151102. [56] Roland Lichters, Robert EW Pfund, and Jürgen Meyer-ter Vehn. Lpic++. a parallel one-dimensional relativistic electromagnetic particle-in-cell code for simulating laser-plasma-interaction. Technical report, Max-Planck-Institut fuer Quantenoptik. [57] S Peter Gary, Shinji Saito, and Hui Li. Cascade of whistler turbulence: Particle-in-cell simulations. 35(2). [58] Yang Chen and Scott E Parker. Electromagnetic gyrokinetic 𝛿f particle-in-cell turbulence simulation with realistic equilibrium profiles and geometry. 220(2):839– 855. [59] Peter Mardahl, Hae-June Lee, Gregg Penn, JS Wurtele, and Nathaniel J Fisch. Intense laser pulse amplification using raman backscatter in plasma channels. 296(2-3):109– 116. [60] S Chakrabarti, JJ Martin, JB Pearson, and RA Lewis. Developing antimatter contain- ment technology: Modeling charged particle oscillations in a penning-malmberg trap. [61] Ricardo A Fonseca, Luis O Silva, Frank S Tsung, Viktor K Decyk, Wei Lu, Chuang Ren, Warren B Mori, S Deng, S Lee, T Katsouleas, et al. Osiris: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators. In International Conference on Computational Science, pages 342–351. Springer. [62] Guangye Chen, Luis Chacón, and Daniel C Barnes. An energy-and charge-conserving, implicit, electrostatic particle-in-cell algorithm. 230(18):7018–7036. [63] Guangye Chen and Luis Chacón. An energy-and charge-conserving, nonlinearly im- plicit, electromagnetic 1d-3v vlasov–darwin particle-in-cell algorithm. 185(10):2391– 2402. [64] Igor V Sokolov. Alternating-order interpolation in a charge-conserving scheme for particle-in-cell simulations. 184(2):320–328. [65] John Villasenor and Oscar Buneman. Rigorous charge conservation for local electro- magnetic field solvers. 69(2-3):306–316. [66] A Bruce Langdon. On enforcing gauss’ law in electromagnetic particle-in-cell codes. 70(3):447–450. [67] Dong-Yeop Na, Haksu Moon, Yuri A Omelchenko, and Fernando L Teixeira. Lo- cal, explicit, and charge-conserving electromagnetic particle-in-cell algorithm on unstructured grids. 44(8):1353–1362. 122 [68] GB Jacobs and Jan S Hesthaven. High-order nodal discontinuous galerkin particle- in-cell method on unstructured grids. Journal of Computational Physics, 214(1):96–121, 2006. [69] Bradley A Shadwick, John C Bowman, and PJ Morrison. Exactly conservative integrators. 59(3):1112–1133. [70] Oscar Buneman. Dissipation of currents in ionized media. 115(3):503. [71] John Dawson. One-dimensional plasma model. 5(4):445–459. [72] James W Eastwood. The virtual particle electromagnetic particle-mesh method. Computer Physics Communications, 64(2):252–266, 1991. [73] Enzo Tonti. Finite formulation of the electromagnetic field. 32:1–44. [74] G. A. Deschamps. Electromagnetics and differential forms. 69(6):676–696. [75] Leszek Demkowicz, Peter Monk, Leon Vardapetyan, and Waldemar Rachowicz. De rham diagram for hp finite element spaces. Computers & Mathematics with Applications, 39(7-8):29–38, 2000. [76] Douglas Arnold, Richard Falk, and Ragnar Winther. Finite element exterior calculus: from hodge theory to numerical stability. 47(2):281–354. [77] Peter Monk et al. Finite element methods for Maxwell’s equations. Oxford University Press. [78] Douglas N Arnold, Richard S Falk, and Ragnar Winther. Finite element exterior calculus, homological techniques, and applications. [79] F. L. Teixeira and W. C. Chew. Lattice electromagnetic theory from a topological viewpoint. 40(1):169–187. [80] Timo Tarhasaari, Lauri Kettunen, and Alain Bossavit. Some realizations of a discrete hodge operator: a reinterpretation of finite element techniques [for em field analysis]. IEEE Transactions on magnetics, 35(3):1494–1497, 1999. [81] Mark L Stowell and Daniel A White. Discretizing transient current densities in the maxwell equations. Technical report, Lawrence Livermore National Lab.(LLNL), Livermore, CA (United States). [82] Dong-Yeop Na, Fernando L Teixeira, and Weng C Chew. Polynomial finite-size shape functions for electromagnetic particle-in-cell algorithms based on unstructured meshes. 4:317–328. [83] VF Kovalev and V Yu Bychenkov. Analytic solutions to the vlasov equations for expanding plasmas. 90(18):185004. [84] Robert E Peterkin and John W Luginsland. A virtual prototyping environment for directed-energy concepts. Computing in Science & Engineering, 4(2):42–49, 2002. 123 [85] Young-Min Shin, Jian-Xun Wang, Larry R Barnett, and Neville C Luhmann. Particle- in-cell simulation analysis of a multicavity w-band sheet beam klystron. IEEE transactions on electron devices, 58(1):251–258, 2010. [86] Simon J. Cooke and George M. Stantchev. Conformal time-domain particle-in-cell simulation of vacuum electronic devices with accurate surface loss. In 2013 IEEE 14th International Vacuum Electronics Conference (IVEC), pages 1–2, 2013. [87] 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. [88] 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. [89] 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. [90] Scott O’Connor, Zane D Crawford, John P Verboncoeur, John Luginsland, and B Shanker. A set of benchmark tests for validation of 3-d particle in cell methods. IEEE Transactions on Plasma Science, 49(5):1724–1731, 2021. [91] Annalisa Buffa, Giancarlo Sangalli, and Rafael Vázquez. Isogeometric analysis in electromagnetics: B-splines approximation. Computer Methods in Applied Mechanics and Engineering, 199(17-20):1143–1152, 2010. [92] 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. [93] Martin Campos Pinto, Katharina Kormann, and Eric Sonnendrücker. Variational framework for structure-preserving electromagnetic particle-in-cell methods. arXiv e-prints, pages arXiv–2101, 2021. [94] FAN Peifeng, QIN Hong, and XIAO Jianyuan. Discovering exact, gauge-invariant, local energy–momentum conservation laws for the electromagnetic gyrokinetic system by high-order field theory on heterogeneous manifolds. Plasma Science and Technology, 23(10):105103, 2021. [95] Benedikt Perse, Katharina Kormann, and Eric Sonnendrücker. Geometric particle- in-cell simulations of the vlasov–maxwell system in curvilinear coordinates. SIAM Journal on Scientific Computing, 43(1):B194–B218, 2021. [96] Scott O’Connor, Zane D Crawford, OH Ramachandran, John Luginsland, and B Shanker. Time integrator agnostic charge conserving finite element pic. arXiv preprint arXiv:2102.06248, 2021. 124 [97] Scott O’Connor, Zane D Crawford, OH Ramachandran, John Luginsland, and B Shanker. Quasi-helmholtz decomposition, gauss’ laws and charge conservation for finite element particle-in-cell. arXiv preprint arXiv:2103.06737, 2021. [98] Ralf Hiptmair. Higher order whitney forms. Progress in Electromagnetics Research, 32:271–299, 2001. [99] Roberto D Graglia, Donald R Wilton, and Andrew F Peterson. Higher order interpolatory vector bases for computational electromagnetics. IEEE transactions on antennas and propagation, 45(3):329–342, 1997. [100] R. D. Graglia and A. F. Peterson. Hierarchical divergence-conforming Nédélec elements for volumetric cells. 60(11):5215–5227. [101] Georges A Deschamps. Electromagnetics and differential forms. Proceedings of the IEEE, 69(6):676–696, 1981. [102] Alain Bossavit. Whitney forms: A class of finite elements for three-dimensional computations in electromagnetism. IEE Proceedings A (Physical Science, Measurement and Instrumentation, Management and Education, Reviews), 135(8):493–500, 1988. [103] Douglas Arnold, Richard Falk, and Ragnar Winther. Finite element exterior calculus: from hodge theory to numerical stability. Bulletin of the American mathematical society, 47(2):281–354, 2010. [104] Karl F Warnick and Peter H Russer. Differential forms and electromagnetic field theory. Progress In Electromagnetics Research, 148:83–112, 2014. [105] Peter Monk. Finite element methods for Maxwell’s equations. Oxford University Press, 2003. [106] Jian-Ming Jin. The finite element method in electromagnetics. John Wiley & Sons, 2015. [107] 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. [108] VF Kovalev and V Yu Bychenkov. Analytic solutions to the vlasov equations for expanding plasmas. Physical review letters, 90(18):185004, 2003. [109] 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. [110] Martin Reiser and Patrick O’Shea. Theory and design of charged particle beams, volume 312. Wiley Online Library, 1994. [111] Y Li, P Chin, R Kishek, M Reiser, M Venturini, JG Wang, Y Zou, and TF Godlove. Design, simulation and test of pulsed panofsky quadrupoles. In Proceedings of the 1999 Particle Accelerator Conference (Cat. No. 99CH36366), volume 5, pages 3369–3371. IEEE, 1999. 125 [112] RH Abrams, B Levush, AA Mondelli, and RK Parker. Vacuum electronics for the 21st century. IEEE Microwave magazine, 2(3):61–72, 2001. [113] Samuel F Martins, RA Fonseca, Wei Lu, Warren B Mori, and LO Silva. Exploring laser- wakefield-accelerator regimes for near-term lasers using particle-in-cell simulation in lorentz-boosted frames. Nature Physics, 6(4):311–316, 2010. [114] Valery A Dolgashev and Sami G Tantawi. Simulations of currents in x-band accelerator structures using 2d and 3d particle-in-cell code. In PACS2001. Proceedings of the 2001 Particle Accelerator Conference (Cat. No. 01CH37268), volume 5, pages 3807–3809. IEEE, 2001. [115] Gerald Sasser, Joseph Blahovec, Lester Bowers, John Luginsland, John Watrous, and Shari Colella. Current emission, resistive losses, and other challenging problems in the simulation of high power microwave components. In 30th Plasmadynamic and Lasers Conference, page 3730, 1999. [116] Jean-Luc Vay, DP Grote, RH Cohen, and Alex Friedman. Novel methods in the particle-in-cell accelerator code-framework warp. Computational Science & Discovery, 5(1):014019, 2012. [117] Ricardo A Fonseca, Luis O Silva, Frank S Tsung, Viktor K Decyk, Wei Lu, Chuang Ren, Warren B Mori, S Deng, S Lee, T Katsouleas, et al. Osiris: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators. In International Conference on Computational Science, pages 342–351. Springer, 2002. [118] 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. [119] 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. [120] Alexander S Glasser and Hong Qin. The geometric theory of charge conservation in particle-in-cell simulations. arXiv preprint arXiv:1910.12395, 2019. [121] Enzo Tonti. Finite formulation of the electromagnetic field. Progress in electromagnetics research, 32:1–44, 2001. [122] Yang Liu, Pieter Ghysels, Lisa Claus, and Xiaoye Sherry Li. Sparse approximate multifrontal factorization with butterfly compression for high-frequency wave equations. SIAM Journal on Scientific Computing, 43(5):S367–S391, 2021. [123] X. Sherry Li and P. Ghysels. Direct sparse linear solvers, preconditioners. [124] Charbel Farhat and Francois-Xavier Roux. A method of finite element tearing and interconnecting and its parallel solution algorithm. International journal for numerical methods in engineering, 32(6):1205–1227, 1991. 126 [125] Charbel Farhat, Michel Lesoinne, Patrick LeTallec, Kendall Pierson, and Daniel Rixen. Feti-dp: a dual–primal unified feti method—part i: A faster alternative to the two-level feti method. International journal for numerical methods in engineering, 50(7):1523–1544, 2001. [126] Yujia Li and J-M Jin. A vector dual-primal finite element tearing and interconnecting method for solving 3-d large-scale electromagnetic problems. IEEE Transactions on Antennas and Propagation, 54(10):3000–3009, 2006. [127] Yu-Jia Li and Jian-Ming Jin. Implementation of the second-order abc in the feti- dpem method for 3d em problems. IEEE Transactions on Antennas and Propagation, 56(8):2765–2769, 2008. [128] Ming-Feng Xue and Jian-Ming Jin. Nonconformal feti-dp methods for large-scale electromagnetic simulation. IEEE transactions on antennas and propagation, 60(9):4291– 4305, 2012. [129] CT Wolfe, U Navsariwala, and Stephen D Gedney. A parallel finite-element tearing and interconnecting algorithm for solution of the vector wave equation with pml absorbing medium. IEEE Transactions on Antennas and Propagation, 48(2):278–284, 2000. [130] Ali Akbarzadeh-Sharbaf, Vahid Mohtashami, and Dennis D Giannacopoulos. An unconditionally stable and energy-preserving domain-decomposition method for transient modeling of large-scale electromagnetic problems. IEEE Transactions on Antennas and Propagation, 67(11):6989–7000, 2019. [131] Umesh D Navsariwala and Stephen D Gedney. An efficient implementation of the finite-element time-domain algorithm on parallel computers using a finite-element tearing and interconnecting algorithm. Microwave and Optical Technology Letters, 16(4):204–208, 1997. [132] Zhi-Hong Lai, Jean-Fu Kiang, and Raj Mittra. A domain decomposition finite difference time domain (fdtd) method for scattering problem from very large rough surfaces. IEEE Transactions on Antennas and Propagation, 63(10):4468–4476, 2015. [133] 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. [134] 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. [135] Jian-Ming Jin and Douglas J Riley. Finite element analysis of antennas and arrays. Wiley Online Library. [136] Kunimasa Saitoh and Masanori Koshiba. Full-vectorial imaginary-distance beam propagation method based on a finite element scheme: application to photonic crystal fibers. 38(7):927–933. 127 [137] M. Fivaz, S. Brunner, G. de Ridder, O. Sauter, T.M. Tran, J. Vaclavik, L. Villard, and K. Appert. Finite element approach to global gyrokinetic particle-in-cell simulations using magnetic coordinates. 111(1):27 – 47. [138] J-F Lee, D-K Sun, and Zoltan J Cendes. Full-wave analysis of dielectric waveguides using tangential vector finite elements. 39(8):1262–1271. [139] BM Azizur Rahman and J Brian Davies. Finite-element analysis of optical and microwave waveguide problems. 32(1):20 – 28. [140] Seung-Cheol Lee, Marinos N Vouvakis, and Jin-Fa Lee. A non-overlapping domain decomposition method with non-matching grids for modeling large finite antenna arrays. 203(1):1–21. [141] Stephen D Gedney and Umesh Navsariwala. An unconditionally stable finite element time-domain solution of the vector wave equation. 5(10):332–334. [142] Thomas Rylander and Anders Bondeson. Stability of explicit–implicit hybrid time- stepping schemes for maxwell’s equations. 179(2):426–438. [143] Lauri Kettunen, Kimmo Forsman, and Alain Bossavit. Gauging in whitney spaces. IEEE transactions on magnetics, 35(3):1466–1469, 1999. [144] Man-Fai Wong, O. Picon, and V. Fouad Hanna. A finite element method based on whitney forms to solve maxwell equations in the time domain. 31(3):1618–1621. [145] Garry Rodrigue and Daniel White. A vector finite element time-domain method for solving maxwell’s equations on unstructured hexahedral grids. 23(3):683–706. [146] Ali Akbarzadeh-Sharbaf and Dennis D Giannacopoulos. Finite-element time-domain solution of the vector wave equation in doubly dispersive media using möbius transformation technique. 61(8):4158–4166. [147] S. D. Gedney and J. A. Roden. Numerical stability of nonorthogonal fdtd methods. 48(2):231–239. [148] Shumin Wang and F. L. Teixeira. Some remarks on the stability of time-domain electromagnetic simulations. 52(3):895–898. [149] Nathan Mortimore Newmark. Computation of dynamic structural response in the range approaching failure. Department of Civil Engineering, University of Illinois. [150] W. L Wood. Practical time-stepping schemes. Oxford [England] : Clarendon Press ; New York : Oxford University Press, 1990. [151] O. C. Zienkiewicz. A new look at the newmark, houbolt and other time stepping formulas. a weighted residual approach. 5(4):413–418. [152] O. C Zienkiewicz and 1945 Morgan, K. (Kenneth). Finite Elements and Approximation. New York Wiley, 1983. 128 [153] Michele Benzi, Gene H Golub, and Jörg Liesen. Numerical solution of saddle point problems. 14:1–137. [154] José A Bittencourt. Fundamentals of plasma physics. Springer Science & Business Media, 2013. [155] 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. [156] Jay P Boris. Relativistic plasma simulation-optimization of a hybrid code. In Proc. 4th Conf. Num. Sim. Plasmas, pages 3–67, 1970. [157] 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. [158] 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. [159] VH Rumsey. Reaction concept in electromagnetic theory. Physical Review, 94(6):1483, 1954. [160] G. Chen, L. Chacón, and D.C. Barnes. An energy- and charge-conserving, implicit, electrostatic particle-in-cell algorithm. Journal of Computational Physics, 230(18):7018 – 7036, 2011. [161] Rodney J Mason. Implicit moment particle simulation of plasmas. Journal of Computational Physics, 41(2):233 – 244, 1981. [162] D.W. Hewett. A global method of solving the electron-field equations in a zero- inertia-electron-hybrid plasma simulation code. Journal of Computational Physics, 38(3):378 – 395, 1980. [163] Charles K Birdsall and A Bruce Langdon. Plasma physics via computer simulation. CRC press, 2004. [164] Burkay Donderici and Fernando L Teixeira. Mixed finite-element time-domain method for transient maxwell equations in doubly dispersive media. IEEE transactions on microwave theory and techniques, 56(1):113–120, 2008. [165] Alain Bossavit. ’generalized finite differences’ in computational electromagnetics. Progress In Electromagnetics Research, 32:45–64, 2001. [166] John Leonidas Volakis, John L Volakis, Arindam Chatterjee, and Leo C Kempel. Finite element method for electromagnetics. Universities Press, 1998. 129 [167] Douglas N Arnold. Mixed finite element methods for elliptic problems. Computer methods in applied mechanics and engineering, 82(1-3):281–300, 1990. [168] Peter Monk. A finite element method for approximating the time-harmonic maxwell equations. Numerische mathematik, 63(1):243–261, 1992. [169] Peter Monk. Analysis of a finite element method for maxwell’s equations. SIAM Journal on Numerical Analysis, 29(3):714–729, 1992. [170] Peter Monk. A comparison of three mixed methods for the time-dependent maxwell’s equations. SIAM journal on scientific and statistical computing, 13(5):1097–1122, 1992. [171] Bo He and F.L. Teixeira. Geometric finite element discretization of maxwell equations in primal and dual spaces. 349(1–4):1 – 14. [172] S. P. Walker, M. J. Bluck, and I. Chatzis. The stability of integral equation time-domain scattering computations for three-dimensional scattering; similarities and differences between electrodynamic and elastodynamic computations. 15(5-6):459–474. [173] A. J. Pray, N. 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. 60(8):3772–3781. [174] A. Bossavit and L. Kettunen. Yee-like schemes on a tetrahedral mesh, with diagonal lumping. 12(1-2):129–142. [175] K. S. Yee and J. S. Chen. The finite-difference time-domain (fdtd) and the finite-volume time-domain (fvtd) methods in solving maxwell’s equations. 45(3):354–363. [176] Ru-Shan Chen, Lei Du, Zhenbao Ye, and Yang Yang. An efficient algorithm for implementing the crank–nicolson scheme in the mixed finite-element time-domain method. 57(10):3216–3222. [177] Richard Cousin, Jean Larour, Jacques Gardelle, Bruno Cassany, Patrick Modin, Philippe Gouard, and Pierre Raymond. Gigawatt emission from a 2.4-ghz compact magnetically insulated line oscillator (milo). 35(5):1467–1475. [178] Douglas N Arnold. Spaces of finite element differential forms. In Analysis and numerics of partial differential equations, pages 117–140. Springer. [179] T Zh Esirkepov. Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor. 135(2):144–153. [180] C-D Munz, Pascal Omnes, Rudolf Schneider, Eric Sonnendrücker, and Ursula Voss. Divergence correction techniques for maxwell solvers based on a hyperbolic model. 161(2):484–511. [181] Karl F Warnick and Peter H Russer. Differential forms and electromagnetic field theory. 148:83–112. 130 [182] PJ Mardahl, JP Verboncoeur, and Ck Birdsall. Progress on a 3d particle-in-cell model of a w-band klystron. In IEEE Conference Record-Abstracts. 1999 IEEE International Conference on Plasma Science. 26th IEEE International Conference (Cat. No. 99CH36297), page 161. IEEE. [183] Guangye Chen, Luis Chacón, Lin Yin, Brian J Albright, David J Stark, and Robert F Bird. A semi-implicit, energy-and charge-conserving particle-in-cell algorithm for the relativistic vlasov-maxwell equations. 407:109228. [184] DB Seidel, ML Kiefer, RS Coats, TD Pointon, JP Quintenz, and WA Johnson. The 3d, electromagnetic, particle-in-cell code, quicksilver. Technical report, Sandia National Labs., Albuquerque, NM (USA). [185] TD Arber, Keith Bennett, CS Brady, A Lawrence-Douglas, MG Ramsay, NJ Sircombe, P Gillies, RG Evans, Holger Schmitz, AR Bell, et al. Contemporary particle-in-cell approach to laser-plasma modelling. 57(11):113001. [186] Andrew J Christlieb, Robert Krasny, John P Verboncoeur, Jerold W Emhoff, and Iain D Boyd. Grid-free plasma simulation techniques. 34(2):149–165. [187] 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. [188] M Surendra. Radiofrequency discharge benchmark model comparison. Plasma Sources Science and Technology, 4(1):56, 1995. [189] 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. [190] Miles M Turner. Verification of particle-in-cell simulations with monte carlo collisions. Plasma Sources Science and Technology, 25(5):054007, 2016. [191] 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. [192] 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. [193] 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. [194] Nathan M Newmark. A method of computation for structural dynamics. Journal of the engineering mechanics division, 85(3):67–94, 1959. [195] Charles K Birdsall and A Bruce Langdon. Plasma physics via computer simulation. CRC press, 2004. 131 [196] David Smithe, Peter Stoltz, John Loverich, Chet Nieter, and Seth Veitzer. Development and application of particle emission algorithms from cut-cell boundaries in the vorpal em-fdtd-pic simulation tool. In 2008 IEEE International Vacuum Electronics Conference, pages 217–218. IEEE, 2008. [197] 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. [198] 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. [199] Richard Marchand. Ptetra, a tool to simulate low orbit satellite–plasma interaction. IEEE Transactions on Plasma Science, 40(2):217–229, 2011. [200] 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. [201] 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. [202] Guangye Chen, Luis Chacón, and Daniel C Barnes. An energy-and charge-conserving, implicit, electrostatic particle-in-cell algorithm. Journal of Computational Physics, 230(18):7018–7036, 2011. [203] Jean-Luc Vay, Irving Haber, and Brendan B Godfrey. A domain decomposition method for pseudo-spectral electromagnetic simulations of plasmas. Journal of Computational Physics, 243:260–268, 2013. [204] Martin Reiser and Patrick O’Shea. Theory and design of charged particle beams, volume 312. Wiley Online Library, 1994. 132