SUPERCONVERGENCE AND ACCURACY ENHANCEMENT OF DISCONTINUOUS GALERKIN SOLUTIONS FOR VLASOV-MAXWELL EQUATIONS AND NUMERICAL ANALYSIS OF A HYBRID METHOD FOR RADIATION TRANSPORT By Andrés Felipe Galindo Olarte A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Applied Mathematics—Doctor of Philosophy 2023 ABSTRACT In this thesis we will analyze and enhance two schemes for kinetic equations. Namely the discontinous Galerkin (DG) methods for solving the the Vlasov-Maxwell (VM) system and a hybrid method for solving the time-dependent radiation transport equation (RTE). In Chapter 2 we will consider the DG methods for solving the VM system, a funda- mental model for collisionless magnetized plasma. The DG methods provide accurate numerical description with conservation and stability properties. However, to resolve the high dimensional probability distribution function, the computational cost is the main bottleneck even for modern-day supercomputers. The first part of this thesis studies the applicability of a post-processing technique to the DG solution to enhance its accuracy and resolution for the VM system. This postprocessor is applied at the final time of the simulation, and its cost is negligible, it succeeds by producing a high-resolution solution with the same cost of computing a low-resolution one, thus saving computational time in the process. In particular, we prove the superconvergence of order (2k + 12 ) in the neg- ative order norm for the probability distribution function and the electromagnetic fields when piecewise polynomial degree k is used. Numerical tests including Landau damp- ing, two-stream instability and streaming Weibel instabilities are considered showing the performance of the post-processor. This is based on joint work with Yingda Cheng, Juntao Huang and Jennyfer Ryan [1]. In Chapter 3, we prove rigorous error estimates for a hybrid method introduced in [2] for solving the time-dependent RTE. The method relies on a splitting of the kinetic distribution function for the radiation into uncollided and collided components. A high- resolution method (in angle) is used to approximate the uncollided components and a low-resolution method is used to approximate the the collided component. After each time step, the kinetic distribution is reinitialized to be entirely uncollided. For this anal- ysis, we consider a mono-energetic problem on a periodic domains, with constant ma- terial cross-sections of arbitrary size. We assume the uncollided equation is solved ex- actly and the collided part is approximated in angle via a spherical harmonic expansion (PN method). Using a non-standard set of semi-norms, we obtain estimates of the form C(ε, σ, ∆t)N −s where s ≥ 1 denotes the regularity of the solution in angle, ε and σ are scattering parameters, ∆t is the time-step before reinitialization, and C is a complicated function of ε, σ, and ∆t. These estimates involve analysis of the multiscale RTE that includes, but necessarily goes beyond, usual spectral analysis. We also compute error estimates for the monolithic PN method with the same resolution as the collided part in the hybrid. Our results highlight the benefits of the hybrid approach over the monolithic discretization in both highly scattering and streaming regimes. This is based in a joint work with Cory D. Hauck and Victor Decaria [3]. Copyright by ANDRÉS FELIPE GALINDO OLARTE 2023 “Las cosas solo son puras si uno las mira desde lejos, es muy importante conocer nuestras raíces , saber de donde venimos, conocer nuestra historia , pero al mismo tiempo tan importante como saber de donde somos, es entender que todos en el fondo somos de ningún lado del todo y de todos lados un poco.” Jorge Drexler v ACKNOWLEDGEMENTS The last five years at Michigan State University have been the most fulfilling and profes- sionally rewarding of my life thus far. This journey would not have been possible without Professor Yingda Cheng’s guidance. I thank her for taking me as her student, for her pa- tience, valuable feedback, and for helping me navigate the job market. I will always be grateful for the opportunity to work with her. I also want to thank Irene M Gamba and the Oden institute at University of Texas at Austin for giving me a job opportunity for the next two years. At MSU, there are many people to thank. First, I would like to acknowledge my com- mittee members, Professor Daniel Appelö, Professor Andrew Christlieb, and Professor Zhengfang Zhou, for their guidance and for taking the time to read and evaluate my work. I would also like to thank Professor Peter Bates for continuously believing in me and recruiting me to this special place called Michigan State University (MSU). Thank you for always giving me advice and helping me prepare my job application material and in- terviews. I would also like to thank my fellow group members, Amit Rotem, Kai Huang, Zhichao Peng, Yann-Meing Law, and Spencer Lee, for their support. I would also like to thank my friends from the math department who have made my stay more pleasant: Luis Suarez, Danika Van Niel, Cullen Haselby, Shih-Fang Yeh, Yuta Hozumi, Chloe Lewis, Gokul Bhusal, Bowen Su, Azzam Alfarraj, Joe Melby, Leonardo Abbrescia, Rodrigo Be- cerra Matos, Ben Jones, Zhixin Wang, Chen Zhang, Max Throm, Shikha Bhutani, Rithwik Vidyarthi, Owen Ekblad, Rolando Ramos, Aldo Garcia, Valerie Jean Pierre, Christopher Potvin, and Albert Chua. I want to thank the math department staff, particularly Laura Willoughby, Estrella Starn, Taylor Alvarado, and Sabrina Walton. Thanks for always help- ing me! I want to acknowledge my collaborators in all the projects in this dissertation; it would have been impossible to carry out my research without their help, Jennifer Ryan, Juntao Huang, Cory D Hauck, and Victor Decaria. Thank you also for guiding me in my profes- vi sional career. My bachelor’s degree was fundamental to my formation. I want to thank the Univer- sidad Distrital Francisco José de Caldas for giving this kid from the south of Bogotá the opportunity to afford high-level education and, thus, the first step to becoming a mathe- matician. The first person I would like to acknowledge is Professor Arturo Sanjuán. He was my first mentor, and I am grateful to call him my friend. Thanks for all the advice on life and mathematics throughout the years and for trusting me to become his first student. He showed me that I could be a person and a mathematician. I also want to thank pro- fessors Alejandro Masmela, Deccy Trejos, Luis Fernando Villarraga, Herbert Sarmiento, Oriol Mora, and Carlos Ochoa. My time at the Universidad de los Andes was terrific. Those three years and a half were fundamental for me. I learned more about math and the profession. On top of that, I had the best friends/colleagues to support me at every step. Thank you Jerson Caro, Daniel Avila, Rodolfo Quintero, Sebastián Osorio, Nicolás Escobar, Santiago Pin- zon, Hernán Garcia, Edison Leguizamon, Gustavo Chaparro, David Moreno Paris, An- drés Patiño, Mariana Vicaria, Duvan Cardona, Weimar Astaiza, Juanita Duque and Laura Gamboa. Your friendship means a lot to me! I also want to thank my advisor, Jean Cortissoz, for helping me through my thesis process and navigating my doctoral pro- gram search. I also want to thank professors Monika Winklemeir, Florent Schaffhauser, Alexander Getmanenko, and Andrei Giniatoulline. Thank you for all your lessons and mentoring. I would also like to recognize all my friends from the Comunidad Latinoamericana at MSU. I appreciate their vote of confidence in me as president of this organization for three years. Thank you for helping me create a space where international students from Latin America and the Caribbean could find a place where they feel at home, correctly represented, respected, and understood. I would like to especially mention Francisco Flores, Marisol Masso, Angelica Herreño, Guillermo Huanes, Vanessa Maldonado, Diego vii Sierra, Diego Granados, Akash Saxena, Jorge Nevares, Erik Amezquita, Skyy Pineda, Jennifer Mojica Santana, Viridiana Garcia, Giovanny Salazar, Henry Gonzalez, Gerardo Melgar, Martina Borges, Luisa Parrado, Marcela Tabares, Paulo Izquierdo, Viviana Or- tiz, Eloy Moreno, Astrid Olave, Ian Fisher, Tayna Carrasquillo, Joon Chung, Joelis Lama, Santiago Hernandez, Santiago Rodriguez, Lina Gomez, Jose Quintero, Francisco Santos, Maria Buitrago, Angelica Bernabe, Julian Bello, William Torres, Juan Sebastian Hernan- dez, Dennise Celis, Carolina Vargas, Catalina Barraza, Daniel Hoffman, Juan Sandoval, Angie Vega, Marco Lopez, Mayra Florez, Pepe Montero, Andrea Bernardes, Christian Gonzalez, Daniel Maldonado, Rafael Castro, Laura Castro, María Alejandra Garcia, Nico- las Scamardi, Martina Gimenez, Omar Posos, Francisco Campos, Pedro Chu, Gianna Mendez, Andres Lanzas, Julian Quiroga, Valeria Obando, and Jean-Paul Ortiz. Thanks to the Center for Latin American and Caribbean Studies (CLACS) for supporting us. Thank you, Laura Medina, for serving as our faculty advisor for three years. I cannot finish this paragraph without thanking the Office for International Students and Scholars (OISS) for helping us with resources and support at different events—thanks to Leidy Egan, Clara Graucob, and Krista McCallum. I cannot express my gratitude to my parents, Martha Olarte and José Galindo. Thank you for always supporting me and my brother through difficult times. You always gave me the utmost care and filled my life with love and support. Thanks for constantly be- lieving in me. Last but not least, I want to give all my love and thanks to my partner Melina Jimenez. Thank you for the gift of your love and company. You make my life better. This Ph.D. journey would not have been the same without you. TE AMO. viii TABLE OF CONTENTS CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 CHAPTER 2 SUPERCONVERGENCE AND ACCURACY ENHANCEMENT OF DISCONTINUOUS GALERKIN SOLUTIONS FOR VLASOV- MAXWELL EQUATIONS . . . . . . . . . . . . . . . . . . . . . . . . 8 CHAPTER 3 NUMERICAL ANALYSIS OF A HYBRID METHOD FOR RADIATION TRANSPORT . . . . . . . . . . . . . . . . . . . . . . . 40 CHAPTER 4 SUMMARY AND CONCLUSION . . . . . . . . . . . . . . . . . . . 66 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 ix CHAPTER 1 INTRODUCTION The evolution of a large particle system such as gases or plasmas at the microscopic level is described by systems of ODE’s, however in general solving these systems numerically is extremelly costly and brings little to no insight into the macroscopic behaviour of the system. Thus we seek a reduced model of particle dynamics, that bridges the gap be- tween the microscopic and macroscopic description of the physical phenomena and that is also accurate. One of such models are the kinetic equations [4],[5]. Kinetic equations intend to describe these particles systems by means of a distribution functions f in phase space. This phase space includes macroscopic variables i.e the position in physical space, but also microscopic variables which describe the state of the particle, in this thesis the only microscopic variables that will be considered are the velocity components, other ex- amples of microscopic variables are, internal energy and spin variables, etc. This object f represents a particle density in phase space i.e. f dx dv is the number of particles in a small volume [4]. Kinetic equations have applications in different fields such as: gas dynamics, plasma physics, biology, socieconomics, nuclear engineering, etc. In the classical kinetic theory of rarified gases, the variation of a non-negative function, f = f (x, v, t), characterizing the particle densities having velocity v ∈ R3 in position x ∈ R3 at time t, is obtained via the equation ∂f + v · ∇x f = Q(f ) ∂t The operator Q(f ), on the right-hand side of equation (1), describes the effects of internal forces due to particle interactions, and its form depends on the details of the microscopic dynamic [5]. The most well-known examples are Boltzmann’s equation and the Vlasov mean-field equation. There are numerous challenges for numerical solvers for kinetic equations. The first one coming from the intrinsic high-dimensionality of the problem, in general (x, v, t) ∈ R7 . 1 There are additional difficulties and requirements specific to kinetic equations: • Conservation properties We would like to preserve physicial laws, such as conserva- tion. [5]. • Computational cost. Aside from the computational cost that comes from high-dimensions. Computations using the operator Q(f ) involve the computation of high-dimensional integrals in velocity space at each point in physical space [5]. • Presence of multiple scales. Usually one have to deal with multiple space-time scales, which span different regimes. [5]. In this work we will consider deterministic numerical methods for kinetic equations. For a survey on these methods see [5]. In this thesis we intend to analyze two different schemes for kinetic equations. In the rest of this chapter we would like introduce the reader to two schemes: Discontinuous Galerkin method for the Vlasov-Maxwell equations and a hybrid method for radiation transport. 1.1 Discontinuous Galerkin method for the Vlasov-Maxwell equations In the first part of this thesis, we consider numerical solutions of the Vlasov-Maxwell (VM) system, a fundamental model for collisionless magnetized plasma. The dimension- less form of the equations that describes the evolution of a single species of non-relativistic electrons under the self-consistent electromagnetic field while the ions are treated as uni- form fixed background is given by ∂t f + v · ∇x f + (E + v × B) · ∇v f = 0, (1.1a) ∂E ∂B = ∇x × B − J, = −∇x × E, (1.1b) ∂t ∂t ∇x · E = ρ − ρi , ∇x · B = 0, (1.1c) with Z Z ρ(x, t) = f (x, v, t) dv, J(x, t) = f (x, v, t)v dv, Ωv Ωv 2 where the equations are defined on Ω = Ωx × Ωv , x ∈ Ωx denotes the position in physical space, and v ∈ Ωv in velocity space. Here f (x, v, t) ≥ 0 is the distribution function of electrons at position x with velocity v at time t, E(x, t) is the electric field, B(x, t) in the magnetic field, ρ(x, t) is the electron charge density, and J(x, t) is the current density. The charge density of background ions is denoted by ρi , which is chosen to satisfy total charge R neutrality, Ωx (ρ(x, t) − ρi ) dx = 0. Periodic boundary conditions in Ωx and compact sup- port in Ωv are assumed. The VM system has wide applications in plasma physics for describing space and laboratory plasmas, with application to fusion devices, high-power microwave generators, and large scale particle accelerators. Much work has been carried out in the literature aiming at accurate deterministic de- scription of the probability density function for nonlinear behavior of charged particles in plasma. Califano et al. used a semi-Lagrangian approach to compute the streaming Weibel instability [6], current filamentation instability [7], magnetic vortices [8], magnetic reconnection [9]. Also, various methods have been proposed for the relativistic VM sys- tem [10, 11, 12, 13]. This work concerns the discontinuous Galerkin (DG) method for solv- ing the VM system. The DG method is a class of finite element method that uses discon- tinuous polynomial spaces, and they have desirable properties for convection-dominated problems [14]. In particular, DG methods have been used to simulate the Vlasov-Poisson system in plasmas [15, 16, 17] and for a gravitational infinite homogeneous stellar system [18]. They have been also used to solve VM system [19, 20] and the relativistic VM system [21]. The DG methods have nice properties such as stability, charge and energy conserva- tion and high order accuracy, which are highly desirable for long time simulations. 1.2 Accuracy enhanement and Superconvergence techniques for DG methods The main computational challenge for any grid based solver for the VM system is the high-dimensionality of the Vlasov equation. This makes the computation extremely expensive even on modern-day exa-scale supercomputers. Post-processing techniques, which can greatly enhance the resolution of the numerical solution at any given time, 3 are therefore desirable because it is only applied once at the end of the simulation with negligible computational cost. Post-processing for finite element methods is a mature technology. The post-processing technique presented here takes advantage of the infor- mation contained in the negative-order norm and was originally developed by Bramble and Schatz [22] in the context of continuous finite element methods for elliptic prob- lems. It consists of a convolution of the finite element solution with a local averaging operator. We can then establish the convergence in the negative order norm which is higher than that one obtained in the usual L2 -norm. In [23], Cockburn, Luskin, Shu and Süli applied this technique to the DG methods for solving linear hyperbolic equations. This technique was further extended to the DG methods for solving nonlinear conser- vational laws [24, 25] and nonlinear symmetric systems of hyperbolic conservation laws [26]. This method is currently part of a filtering family known as a Smoothness-Increasing Accuracy-Conserving (SIAC) filters [27]. This chapter will demonstrate the performance of post-processing by the SIAC filter for DG solutions to the VM system. In particular, we consider benchmark numerical tests for Vlasov-Ampére (VA) and VM systems, and study the numerical error for short and long time simulations with varying polynomial order. In order to validate the enhanced accuracy of the post-processed solution, an impor- tant step is to establish the superconvergence of the negative order norm of the error and its divided differences. In [23], Cockburn, Luskin, Shu and Süli established a framework to prove negative-order estimates for the DG solutions to linear conservational laws of order 2k + 1 using polynomials of degree k. After this, there have been important exten- sions. L2 and L∞ superconvergence estimates were established for DG solutions for lin- ear constant coefficient hyperbolic systems with the position-dependent SIAC filter [28]. Ji, Meng et al [24, 25, 26] proved superconvergence for non-linear conservation laws and nonlinear symmetric hyperbolic systems of the DG solutions of order at least ( 23 k +1). It is highly nontrivial to establish superconvergence for nonlinear problems because a suitable 4 dual problem has to be identified, and additionally the divided difference of the solution does not satisfy the PDE, which makes the proof highly technical [25, 26]. In Chapter 2, we aim to prove negative-order estimates of DG solutions to the VM system. Since the VM system is nonlinear, it is nontrivial to extend the proof in [23]. We identify a proper dual problem, which aids the estimates of the consistency term. In the end, we proved superconvergence of order (2k + 21 ) in the negative norm for the probability distribution function and the electromagnetic fields. 1.3 A hybrid method for the radiation transport equation The second part of this thesis will deal with the radiation transport equation (RTE) [29, 30, 31]. This equation describes the movement of particles through a material medium by means of a kinetic distribution function that gives the density of particles with respect to the local phase space measure. In a general setting, the phase space is six dimen- sional: three dimensions for particle position and three for particle momentum, the latter of which is typically decomposed into energy and direction (or angle) of flight. Thus in the time-dependent setting, the RTE is defined over a seven-dimensional domain. The RTE describes two basic processes: particle advection and interactions with the material medium. These interactions can be of various types and include scattering and emission/absorption processes. The rate at which these processes occur is determined by the properties of the material, expressed via cross-sections. Material cross-sections may vary in space and depend on the particle energy and, in situations that the material evolves, the cross-sections may evolve as well. When cross-sections vary significantly, the RTE may exhibit multiscale behavior. It is the combination of this multiscale behavior with the high-dimensional phase space that makes simulating the RTE a challenging task. A well-known multi-scale feature of the RTE is the diffusion limit. In regions where the scattering cross-section is large, the solution of the RTE can be accurately approxi- mated by its angular average [32, 33]. Moreover this average is well-approximated by the solution of a diffusion equation. This solution to diffusion equation has long been used 5 as a cheap approximation the solution to the RTE in scattering dominated regimes. Another common limit is the absorption limit, which is characterized by a complete lack of scattering. In this case, the RTE does not have a simple asymptotic approximation. However, due the abscence of scattering, there is no coupling between the angles and energies of the kinetic distribution. Thus, with a proper discretiation, the RTE solution can be easily parallelized. In problems for which the scattering cross-section varies dramatically, both of the lim- its above can exist simultaneously, along with a range of transition regimes in between. A consequence of this fact is that a monolithic numerical treat of the RTE will require many degrees of freedom that are strongly coupled. In practice, the time-dependent RTE is of- ten updated in time with an implicit scheme. In such cases, designing the linear solvers can be a challenge. A variety of approaches have been proposed for addressing the multiscale challenges posed of the RTE. These include micro-macro decompositions [34], high order-low order (HOLO) methods [35], diffusion-based acceleration [30, 36], and preconditioned Krylov approaches [37]. In this thesis, we consider a hybrid formulation [2] that is based on the notion of first-collision source [38]. In this hybrid formulation, the RTE is split into two components: an uncollided component that tracks the particles up to point of their first material interaction and a collided component that track the particles that remain. The resulting system is then approximated with two different angular discretizations: a high-resolution discretization for the uncollided equation and a low-resolution for the collided equation. The intuition that drives this strategy is that scattering produces a smoother solution; hence the collided equation should require less resolution to recover an accurate solution. The uncollided equation, on the other, requires higher resolution; however it takes the form of a purely absorbing RTE and can therefore be solved much more efficiently the original RTE using the same number of degrees of freedom. The efficiency of the hybrid approach for the RTE has been demonstrated in several papers 6 [39, 40, 41], including generalizations to hybrid energy discretizations [42] and hybrid spatial discretizations [43]. A key component of the hybrid implementation for the time-dependent RTE is a rela- beling procedure that, after a given time step, maps the collided numerical solution into the space of the uncollided numerical solution and then uses the sum to re-initialize the uncollided equation. Meanwhile, the collided equation is re-initialized to zero. This re- labeling step is critical, since otherwise the hybrid numerical solution would eventually convergence to a low-resolution numerical solution of the collided equation. Despite the intuitive motivation of the hybrid and the success of the hybrid approach in numerical simulations, the method still lacks rigorous justification. This is due in part to complications introduced by the multiscale behavior of the RTE. For example, spec- tral approximations of the RTE in angle are fairly straightforward to analyze [44], but a multiscale analysis that takes into account the degree of scattering is significantly more complicated [45]. The relabeling step of the hybrid formulation complicates the situation even further. In Chapter 3, we take a first step in analyzing the hybrid method for the time-dependent mono-energetic version of the RTE with isotropic scattering. We focus only on the angu- lar discretization of the RTE, comparing the standard spectral approximation (PN ) for the full system with a discretization of the hybrid that features a spectral approximation of the same resolution for the collided equation but assumes an exact solution for the un- collided equation. Clearly, the hybrid formulated in this way is more expensive than the monolithic approach. Thus the goal of the analysis is determine what is gained from the extra work involves in a high-resolution simulation of the uncollided equation, which in practice is computed with a high-fidelity collocation method or with a Monte-Carlo method. In Chapter 4, we provide a brief conclusion and future work. 7 CHAPTER 2 SUPERCONVERGENCE AND ACCURACY ENHANCEMENT OF DISCONTINUOUS GALERKIN SOLUTIONS FOR VLASOV-MAXWELL EQUATIONS In this chapter, we will provide the first step towards proving rigorously accuracy en- hancement of the postprocessed DG solution to the VM system (1.1), by proving super- convergence of the negative norm of the DG solution. The remainder of the chapter is organized as follows. In Section 2.1, we introduce the DG method for the VM system as well as relevant notations that will be required for the negative order estimates. In Section 2.2 we introduce SIAC filtering. In Section 2.3 we prove the negative-order norm estimates of the DG solutions to the VM system. The superconvergence results are con- firmed numerically in Section 2.4. 2.1 Discontinuous Galerkin Numerical Scheme 2.1.1 Notations, Definitions and Projections We begin by introducing the necessary notation used in the chapter. Without loss of generality, we assume the spatial and velocity domain to be Ωx = [−Lx , Lx ]dx and Ωv = [−Lv , Lv ]dv , where Lv is chosen large enough so that f = 0 at ∂Ωv . Through out the chapter, standard notations will be used for the Sobolev spaces. Given a bounded domain D ∈ R⋆ (with ⋆ = dx ,dv , or dx + dv ) and any nonnegative integer m, H m (D) denotes the L2 -Sobolev space of order m with the standard Sobolev norm ∥·∥m,D , W m,∞ denotes the L∞ -Sobolev space of order m with the standard Sobolev norm ∥·∥m,∞,D and the semi-norm | · |m,∞,D . When m = 0, we also use H 0 (D) = L2 (D) and W 0,∞ (D) = L∞ (D). Let Thx = {Kx } and Thv = {Kv } be partitions of Ωx and Ωv , respectively, with Kx and Kv being Cartesian elements or simplices; then Th = {K : K = Kx × Kv , ∀Kx ∈ Thx , ∀Kv ∈ Thv } defines a partition of Ω. Let Ex be the set of the edges of Thx and Ev the set of the edges of Thv ; then the edges of Th will be E = {Kx × ev : ∀Kx ∈ Thx , ∀ev ∈ Ev } ∪ {ex × Kv : ∀ex ∈ Ex , ∀Kx ∈ Thx }. 8 Furthermore, Ev = Evi ∪ Evb with Evi and Evb being the set of interior and boundary edges of Thv respectively. In addition, we denote the mesh size of Th as h = max(hx , hv ) = maxK∈Th hK , where hx = maxKx ∈Thx hKx with hKx = diam(Kx ), hv = maxKv ∈Thv hKv with hKv = diam(Kv ), and hK = max(hKx , hKv ) for K = Kx × Kv . When the mesh is refined, we hx hv assume both hx,min and hv,min are uniformly bounded from above by a positive constant σ0 . Here hx,min = minKx hKx ∈Thx and hv,min = minKv ∈Thv hKv . It is further assumed that {Th⋆ }h is shape-regular with ⋆ = x or v. That is, if ρK⋆ denotes the diameter of the largest sphere included in K⋆ , there is hK⋆ ≤ σ⋆ , ∀K⋆ ∈ Th⋆ ρK⋆ for a positive constant σ⋆ independent of h⋆ . Furthermore the inner products are defined as Z X Z (g, h)Ω = gh dx dv = gh dx dv, (2.1) Ω K∈Th K Z X Z (U, W)Ωx = U · W dx = U · W dx. (2.2) Ωx Kx ∈Thx Kx Now for g ∈ L2 (Ω), U, W ∈ (L2 (Ωx ))dx , we define the L2 -norm of (g, U, W) as q ∥(g, U, W)∥0,Ω = ∥g∥20,Ω + ∥U∥20,Ωx + ∥W∥20,Ωx (2.3) This will be helpful in the error analysis of the negative-order norm. The negative order norm of order l is defined as: given l > 0 and domain Ω, (g, ϕ)Ω + (U, U)Ωx + (W, W)Ωx ∥(g, U, W)∥−l,Ω = sup q ϕ∈C0∞ (Ω),U ,W∈[C0∞ (Ωx )]dx ∥ϕ∥2l,Ω + ∥U∥2l,Ωx + ∥W∥2l,Ωx Next we define the discrete spaces Ghk = g ∈ L2 (Ω) : g|K=Kx ×Kv ∈ P k (Kx × Kv ), ∀Kx ∈ Thx , ∀Kx ∈ Thx , ∀Kv ∈ Thv ,  (2.4) = g ∈ L2 (Ω) : g|K ∈ P k (K), ∀K ∈ Th ,  n dx o Uhr = U ∈ L2 (Ωx ) : U|Kx ∈ [P r (Kx )]dx , ∀Kx ∈ Thx ,  (2.5) 9 where P r (D) denotes the set of polynomials of total degree at most r on D, and k and r are nonnegative integers. For piecewise functions defined with respect to Thx or Thv , we further introduce the jumps and averages as follows. For any edge e = {Kx+ ∩ Kx− } ∈ Ex , with n± x as the outward unit normal to ∂Kx± , g ± = g|Kx± and U± = U|Kx± , the jump across e are defined as − − − − − − [g]x = g + n+ x + g nx , [U]x = U+ · n+ x + U · nx , [U]τ = U+ × n+ x + U × nx and the averages are 1 1 {g}x = (g + + g − ), {U}x = (U+ + U− ). 2 2 By replacing the subscript x with v, one can define [g]v , [U]v , {g}v , and {U}v for an interior edge of Thv in Evi . For a boundary edge e ∈ Evb with nv being the outward unit normal we use 1 1 [g]v = gnv , {g}v = g, {U}v = U. (2.6) 2 2 This is consistent with the fact that the exact solution f is compactly supported in v. R R P R For convenience, we introduce some shorthand notations, Ω⋆ = T ⋆ = K⋆ ∈T ⋆ K⋆ , h h R R P R R P R Ω = Th = K∈Th K , E⋆ = e∈E⋆ e , where again ⋆ is x or v. In addition, ∥g∥0,E = R R 1/2 2 2 1/2 2 (∥g∥0,Ex ×Thv + ∥g∥0,Thx ×Ev ) with ∥g∥0,Ex ×Th =v Ex Thv g dv dsx , and we have that R R 1/2 ∥g∥0,Thx ×Ev = T x Ev g 2 dsv dx . We will make use of the following equality, which h can be easily verified using the definition of averages and jumps. 1 2 [g ]⋆ = g⋆ [g]⋆ , with ⋆ = x or v. (2.7) 2 2.1.2 The DG method for the Vlasov-Maxwell system Now we review the DG method for the VM system proposed in [19]. The scheme seeks a numerical solution fh ∈ Ghk and (Eh , Bh ) ∈ Uhk × Uhk such that for any g ∈ Ghk , 10 U, W ∈ Uhk , Z Z Z ∂t fh g dxdv − fh v · ∇x g dxdv − fh (Eh + v × Bh ) · ∇v g dxdv K K K Z Z Z Z + fh\v · nx g dsx dv + (fh (Eh +\ v × Bh ) · nv )g dsv dx = 0, Kv ∂Kx Kx ∂Kv (2.8a) Z Z Z Z ∂t Eh · U dx = Bh · ∇x × U dx + n\x × Bh · U dsx − Jh · U dx, (2.8b) Kx Kx ∂Kx Kx Z Z Z ∂t Bh · W dx = − Eh · ∇x × W dx − n\x × Eh · W dsx (2.8c) Kx Kx ∂Kx with Z Jh (x, t) = fh (x, v, t)v dv. (2.9) Thv Here nx and nv are outward unit normals of ∂Kx and ∂Kv , respectively. All “hat” func- tions are numerical fluxes that are determined by upwinding, i.e.,   |v · n x | fh\v · nx := fg h v · nx = {fh v}x + [fh ]x · nx (2.10a) 2 fh (Eh +\ v × Bh ) · nv := fh (Eh + ^ v × Bh ) · nv   |(Eh + v × Bh ) · nv | = {fh (Eh + v × Bh )}v + [fh ]v · nv , (2.10b) 2   1 nx × Eh := nx × Eh = nx × {Eh }x + [Bh ]τ \ f (2.10c) 2   1 nx × Bh := nx × Bh = nx × {Bh }x − [Eh ]τ \ f (2.10d) 2 where these relations define the meaning of “tilde”. In [19], alternating and central fluxes for the Maxwell’s equation are also considered. The discussions will be similar to what will be presented in this chapter for the upwind flux, and thus are omitted. Upon summing up (2.8a) with respect to K ∈ Th and similarly summing (2.8b) and (2.8c) with respect to Kx ∈ Thx , the scheme (2.8) becomes the following: look for fh ∈ Ghk , Eh , Bh ∈ Uhk , such that ((fh )t , g)Ω + ah (fh , Eh , Bh ; g) = 0 (2.11a) ((Eh )t , U)Ωx + ((Bh )t , W)Ωx + bh (Eh , Bh ; U, W) = lh (Jh ; U), (2.11b) 11 for any g ∈ Ghk , U, W ∈ Uhk , where Z ah (fh , Eh , Bh ; g) = ah,1 (fh ; g) + ah,2 (fh , Eh , Bh ; g), lh (Jh ; U) = − Jh · U dx, Thx Z Z bh (Eh , Bh ; U, W) = − Bh · ∇x × U dx − fh · [U] dsx B τ Thx Ex Z Z + Eh · ∇x × W dx + Efh · [W] dsx , τ Thx Ex and Z Z Z ah,1 (fh ; g) = − fh v · ∇x g dxdv + h v · [g]x dsx dv fg Th Thv Ex Z ah,2 (fh , Eh , Bh ; g) = − fh (Eh + v × Bh ) · ∇v g dxdv Th Z Z + fh (Eh^+ v × Bh ) · [g]v dsv dx Thx Ev The semi-discrete formulation (2.8) can then be solved by a numerical ODE solver, see the description in [19]. The L2 and energy stability of (2.8) are established in [19]. The main result in [19] for the semi-discrete L2 error estimates of the approximations fh , Eh , Bh , is as follows. Theorem 2.1.1 ([19]). For k ≥ 2 when dx = 3 and k ≥ 1 when dx = 1, 2, the semi-discrete DG method of (2.11a)-(2.11b), for the Vlasov-Maxwell equations with the upwind fluxes of (2.10a)- (2.10d), has the following error estimate ∥(f − fh )(t)∥20,Ω + ∥(E − Eh )(t)∥20,Ωx + ∥(B − Bh )(t)∥20,Ωx ≤ Ch2k+1 , ∀t ∈ [0, T ]. (2.12) Here the constant C is independent of h, but depends on the upper bounds of ∥∂t f ∥k+1,Ω ,∥f ∥k+1,Ω , |f |1,∞,Ω , ∥E∥1,∞,Ωx , ∥B∥1,∞,Ωx , ∥E∥k+1,Ωx , ∥B∥k+1,Ωx over the time interval [0, T ], and it also depends on the polynomial degree k, mesh parameters σ0 , σx and σv , and domain parameters Lx and Lv . In this work, we also consider (1.1) when there is no magnetic field (i.e. when B = 0). This reduced problem is called the VA system, and the DG discretizations would follow a similar discussion by setting Bh = 0 in (2.8) at all times. 12 2.2 Smoothness-Increasing Accuracy-Conserving Filters We extract the higher-order accuracy of the DG method solved over a uniform mesh contained in the negative-order norm by using the SIAC filter. This technique could also be applied over nonuniform meshes, however this would force us to compute the post- processing coefficients in each element in the mesh, increasing the computational com- plexity of the implementation [46]. This filter improves the order of accuracy by reducing the spurious oscillations in the error. This is done by convolving the numerical approxi- mation with a specially chosen kernel, 2(k+1),k+1 (fh∗ (x, v), E∗h (x), B∗h (x)) = Kh ⋆ (fh , Eh , Bh )(x, v), (2.13) where (fh∗ , E∗h , B∗h ) is the filtered solution, (fh , Eh , Bh ) is an approximated solution com- 2(k+1),k+1 puted at the final time, and Kh is the convolution kernel. The kernel is translation- invariant and composed of a linear combination of B-splines of order k + 1 obtained by convolving the characteristic function over the interval (− 12 , 12 ) with itself k times and scaled by the uniform mesh size. Using B-splines makes this kernel computationally ef- ficient, provided the mesh is uniform, as the kernel is translation invariant and is locally supported in at most 2k + 2 elements. The one-dimensional convolution kernel is of the form: k 2(k+1),k+1 1 X 2(k+1),k+1 (k+1)  x  Kh (x) = cγ ψ −γ . (2.14) h γ=−k h 2(k+1),k+1 The weights of the B-splines, cγ , are chosen so that accuracy is not destroyed (the 2(k+1),k+1 kernel can reproduce polynomials of degree up to 2k), i.e. Kh ⋆ p = p for p = 1, x, · · · , x2k , see [23] for details. For the general case, assume the mesh size is uniform in each direction, given arbitrary (x, v) = (x1 , · · · , xdx , v1 , · · · , vdv ) ∈ Rdx +dv , we set Ydx Ydv (k+1) (k+1) ψ (x, v) = ψ (xi ) ψ (k+1) (vj ) (2.15) i=1 j=1 13 The kernel for our case is of the form 2(k+1),k+1 1 Kh (x, v) = Q  Q  dx dv i=1 hxi j=1 hvj    X x1 xd v 1 vd × cγ2(k+1),k+1 ψ (k+1) ,··· , x, ,··· , v −γ (2.16) dx +dv hx1 hdx hv1 hdv γ∈{−k,...,k} where hxi and hvi denote the mesh size in xi and vi direction, resp. The success of the filter relies on the following results. Theorem 2.2.1. (Bramble and Schatz [22]) For T > 0, let u = (f, E, B) be the exact solution 2(k+1),k+1 of the problem (1.1). Let Ω0 + 2supp(Kh (x, v)) ⊂⊂ Ω and U = (fh , Eh , Bh ) is any approximation to u, then 2(k+1),k+1 h2k+2 X ∥u(T ) − Kh ⋆ U ∥0,Ω0 ≤ |u|2k+2,Ω + CP ∥∂hλ (u − U )∥−(k+1),Ω . (2.17) (2k + 2)! |λ|≤k+1 2(k+1),k+1 where CP depends solely on Ω0 , Ω, dx , dv , k, cγ , and it is independent of h. In (2.17), we used the notation of the divided differences. We define      1 1 1 ∂hxi w(x, v) = w x + hxi ei , v − w x − hxi ei , v , (2.18) hxi 2 2 here ei is the unit multi-index whose i-th component is 1 and all others 0. Analogously for velocity space variables vj , the difference quotients are defined as      1 1 1 ∂hvj w(x, v) = w x, v + hvj ej − w x, v − hvj ej , (2.19) hvj 2 2 For any multi-index λ = (αx1 , · · · , αdx , βv1 , · · · , βdv ) we set λ-th order difference quo- tient to be α β ∂hλ w(x, v) = (∂hαx1 · · · ∂hxdx )(∂hβv1 · · · ∂hvdv )w(x, v). (2.20) 1 dx 1 dv 2.3 Superconvergent Error Estimates for the DG method In this section, we prove the superconvergence error estimate in the negative norm of the DG solution for the VM system. In Section 2.3.1, we review basic approximation and regularity properties. Section 2.3.2 will construct the dual problem which is the key to our estimates. The main result and the proof will be given in Section 2.3.3. 14 2.3.1 Preliminaries We summarize some of the standard approximation properties of the above discrete spaces, as well as some inverse inequalities [47]. For any nonnegative integer k, Let Πk be g k the L2 projection onto Ghk , and Πm 2 m x be the L projection onto Uh . We define ζh = Π g − g and ζhU = Πkx U − U, as the Projection errors of g and U respectively. Lemma 2.3.1. (Approximation properties) There exist a constant C > 0, such that for any g ∈ H k+1 (Ω) and U ∈ [H k+1 (Ω)]dx , the following hold: 1/2 ∥ζhg ∥0,K + hK ∥∇⋆ ζhg ∥0,K + hK ∥ζhg ∥0,∂K ≤ Chk+1 K ∥g∥k+1,K , ∀K ∈ Th 1/2 ∥ζhU ∥0,Kx + hKx ∥∇x × ζhU ∥0,Kx + hKx ∥ζhU ∥0,∂Kx ≤ Chk+1 Kx ∥U∥k+1,Kx , ∀Kx ∈ Thx ∥ζhU ∥0,∞,Kx ≤ Chk+1 Kx ∥U∥k+1,∞,Kx , ∀Kx ∈ Thx where the constant C is independent of the mesh sizes hK and hKx , but depends on k and the shape regularity parameters σx and σv of the mesh. Here ⋆ = x or v. Lemma 2.3.2 (Inverse inequality). There exists a constant C > 0, such that for any g ∈ P k (K) or P k (Kx ) × P k (Kv ) with K = (Kx × Kv ) ∈ Th , and for any U ∈ [P k (Kx )]dx , the following hold: ∥∇x g∥0,K ≤ Ch−1 Kx ∥g∥0,K , ∥∇v g∥0,K ≤ Ch−1 Kv ∥g∥0,K , −d /2 −1/2 ∥U∥0,∞,Kx ≤ ChKxx ∥U∥0,Kx , ∥U∥0,∂Kx ≤ ChKx ∥U∥0,Kx , where the constant C is independent of the mesh sizes hKx , hKv , but depends on k and the shape regularity parameters σx and σv of the mesh. To assist the proof, we also need a regularity result for a linear PDE system. 15 Lemma 2.3.3. Consider the following system of equations with periodic boundary conditions in x and zero boundary condition in v for all t ∈ [0, T ]: ∂t φ + A1 (x, v, t) · ∇x φ + A2 (x, v, t) · ∇v φ + A3 (x, v, t) · F = 0, (2.21a) Z ∂t F = ∇x × D + g∇v φ dv, (2.21b) Ωv Z ∂t D = −∇x × F − g(v × ∇v φ) dv, (2.21c) Ωv where the given functions A1 , A2 ∈ W l+1,∞ (Ω) satisfy the divergence free constraint ∇x ·A1 = 0 and ∇v · A2 = 0. For any l ≥ 0 and the fixed time t, the solution to (2.21) satisfy the following estimate ∥φ(·, ·, t)∥2l,Ω + ∥F(·, t)∥2l,Ωx + ∥D(·, t)∥2l,Ωx ≤ C ∥φ(·, ·, 0)∥2l,Ω + ∥F(·, 0)∥2l,Ωx + ∥D(·, 0)∥2l,Ωx .   (2.22) Here C depends on ∥A3 ∥L∞ ((0,T );W l,∞ (Ω)) and ∥g∥L∞ ((0,T );W l+1,∞ (Ω)) . Proof. See the appendix. 2.3.2 The dual problem In order to prove negative-order estimates for the system, the key is to find the dual problem associated to (1.1). We note that, for the nonlinear problem, the dual problem is not unique, see [48]. We construct the dual problem as follows: find functions φ(·, ·, t), F(·, t) and D(·, t) such that φ(·, v, t) is periodic in all dimensions in space and φ(x, ·, t) vanishes in the boundary of the velocity region for all t ∈ [0, T ] and ∂t φ + v · ∇x φ + (E + v × B) · ∇v φ − v · F = 0 (2.23a) Z ∂t F = ∇x × D − f ∇v φ dv, (2.23b) Ωv Z ∂t D = −∇x × F + f (v × ∇v φ) dv (2.23c) Ωv with final time conditions φ(x, v, T ) = Φ(x), F(x, T ) = F(x) and D(x, T ) = D(x), Φ ∈ C0∞ (Ω) and D, F ∈ [C0∞ (Ωx )]dx . 16 Notice that by multiplying (φ, F, D) on both sides of (1.1a)-(1.1b), and multiplying by (f, E, B) on both sides of (2.23a)-(2.23c), and then summing up and integrating over velocity and physical space, we obtain Z Z Z ∂t (f φ) dx dv + ∇x · (f φv) dx dv + ∇v · (f φ(E + v × B)) dx dv Ω Ω Ω Z − f v · F dx dv = 0, Z Ω Z ∂t (E · F + B · D) dx + f v · F dx dv Ωx Z Z Ω Z = ∇x · (B × F) dx + ∇x · (E × D) dx − f (E + v × B) · ∇v φ dxdv, Ωx Ωx Ω where we used the identities ∇⋆ (ϕU) = ϕ∇⋆ · U + U · ∇⋆ ϕ, ∇⋆ · (U × W) = W · (∇⋆ × U) − U · (∇⋆ × W), for scalar functions ϕ and vector functions U and W and the fact that ∇v · (E + v × B) = 0. By adding all equations above and using boundary conditions, we arrive at d [(f, φ)Ω + (E, F)Ωx + (B, D)Ωx ] + F(f, E, B; φ) = 0, (2.24) dt where Z F(f, E, B; φ) = f (E + v × B) · ∇v φ dxdv. (2.25) Ω 2.3.3 The main result In this part, we give our main theorem on the negative-norm of the error for the DG solutions. Note that superconvergence of the negative norm of the solution itself is not sufficient in proving high order convergence of the post-processed solution according to Theorem 2.2.1. However, it is a necessary first step. As shown in [25], it is highly non- trivial to prove superconvergence of the divided difference of the solution for nonlinear problems, we will leave this to explore in our future work. 17 Theorem 2.3.1. If (fh , Eh , Bh ) is a solution to (2.11a)-(2.11b) with the numerical initial condi- tion fh = Πk f and Eh = Πkx E, Bh = Πkx B and k ≥ (dx + dv )/2, then ∥(f − fh , E − Eh , B − Bh )∥−(k+1),Ω ≤ Ch2k+1/2 , where C is a constant independent of h and depends on the upper bounds of ∥∂t f ∥k+2,Ω ,∥f ∥k+2,Ω , |f |1,∞,Ω , ∥E∥1,∞,Ωx , ∥B∥1,∞,Ωx , ∥E∥k+2,Ωx , ∥B∥k+2,Ωx over the time interval [0, T ], and it also depends on the polynomial degree k, mesh parameters σ0 , σx and σv , and domain parameters Lx and Lv . Proof. We define efh = f − fh = ϵfh − ζhf , where ϵfh = Πk f − fh and ζhf is defined just as in B k E E E Section 2.3.1. Analogously ϵE k h = Πx E−Eh , ϵh = Πx B−Bh , then eh = E−Eh = ϵh −ζh and ∞ ∞ dx eB B B h = B−Bh = ϵh −ζh . We follow the ideas in [23]. For any Φ ∈ C0 (Ω), F, D ∈ [C0 (Ωx )] , we estimate the term (efh (T ), Φ)Ω + (eE B h (T ), F)Ωx + (eh (T ), D)Ωx =(efh (T ), φ(T ))Ω + (eE B h (T ), F(T ))Ωx + (eh (T ), D(T ))Ωx =(f (T ), φ(T ))Ω + (E(T ), F(T ))Ωx + (B(T ), D(T ))Ωx − [(fh (T ), φ(T ))Ω + (Eh (T ), F(T ))Ωx + (Bh (T ), D(T ))Ωx ] Z T =(f0 , φ(0))Ω + (E0 , F(0))Ωx + (B0 , D(0))Ωx − F(f, E, B; φ) dτ 0 − (fh (0), φ(0))Ω − (Eh (0), F(0))Ωx − (Bh (0), D(0))Ωx Z T d − [(fh , φ)Ω + (Eh , F)Ωx + (Bh , D)Ωx ] dτ 0 dt h i = − (ζhf0 , φ(0))Ω + (ζhE0 , F(0))Ωx + (ζhB , D(0))Ωx Z T − ((fh )t , φ)Ω + ((Eh )t , F)Ωx + ((Bh )t , D)Ωx dτ 0 Z T − (fh , φt )Ω + (Eh , Ft )Ωx + (Bh , Dt )Ωx + F(f, E, B; φ) dτ, 0 where for the first equality we used (2.24), and the numerical initial condition is used in the last equality. 18 Notice that for any χ ∈ Ghk , ξ, η ∈ Uhk Z T ((fh )t , φ)Ω + ((Eh )t , F)Ωx + ((Bh )t , D)Ωx dτ 0 Z T Z T Z T = ((fh )t , φ − χ)Ω dτ + ((fh )t , χ)Ω dτ + ((Eh )t , F − ξ)Ωx + ((Bh )t , D − η)Ωx dτ 0 0 0 Z T + ((Eh )t , ξ)Ωx + ((Bh )t , η)Ωx dτ 0 Z T Z T = ((fh )t , φ − χ)Ω dτ − ah (fh , Eh , Bh ; χ) dτ 0 0 Z T + ((Eh )t , F − ξ)Ωx + ((Bh )t , D − η)Ωx dτ 0 Z T − bh (Eh , Bh ; ξ, η) − lh (Jh , ξ) dτ 0 Z T Z T = ((fh )t , φ − χ)Ω + ah (fh , Eh , Bh ; φ − χ) dτ + ((Eh )t , F − ξ)Ωx + ((Bh )t , D − η)Ωx dτ 0 0 Z T Z T + bh (Eh , Bh ; F − ξ, D − η) − lh (Jh , F − ξ) dτ − ah (fh , Eh , Bh ; φ) dτ 0 0 Z T − bh (Eh , Bh ; F, D) − lh (Jh , F) dτ. 0 After this calculation we can conclude that (efh (T ), Φ)Ω + (eE B h (T ), F)Ωx + (eh (T ), D)Ωx = ΘM + ΘN + ΘC , (2.26) where h i ΘM = − (ζhf0 , φ(0))Ω + (ζhE0 , F(0))Ωx + (ζhB0 , D(0))Ωx , Z T ΘN = − ((fh )t , φ − χ)Ω + ah (fh , Eh , Bh ; φ − χ) dτ 0 Z T − ((Eh )t , F − ξ)Ωx + ((Bh )t , D − η)Ωx + bh (Eh , Bh ; F − ξ, D − η) − lh (Jh , F − ξ) dτ, 0 Z T ΘC = − (fh , φt )Ω − ah (fh , Eh , Bh ; φ) dτ 0 Z T Z T − (Eh , Ft )Ωx + (Bh , Dt )Ωx − bh (Eh , Bh ; F, D) + lh (Jh , F) dτ − F(f, E, B; φ) dτ. 0 0 In the following we will estimate ΘM , ΘN and ΘC . 19 Lemma 2.3.4 (Projection Estimate). ΘM satisfies q |ΘM | ≤ Ch2k+2 ∥φ(0)∥2k+1,Ω + ∥F(0)∥2k+1,Ωx + ∥D(0)∥2k+1,Ωx (2.28) where C depends on ∥f0 ∥k+1,Ω , ∥E0 ∥k+1,Ωx and ∥B0 ∥k+1,Ωx . Proof. By the definition of Πk , (f0 − Πk f0 , φ(0))Ω = (f0 − Πk f0 , φ(0) − Πk φ(0))Ω ≤ ∥f0 − Πk f0 ∥∥φ(0) − Πk φ(0)∥ ≤ Chk+1 ∥f0 ∥k+1,Ω hk+1 ∥φ(0)∥k+1,Ω . The last line was an application of the first part of Lemma 2.3.1. By the same lines we obtain analogous results for the E and B parts. The conclusion follows by grouping them all together and an application of Cauchy-Schwarz inequality. For the second term, we have the following result: Lemma 2.3.5 (Residual). Let χ = Πk f, ξ = Πkx F, η = Πkx D, we have Z T 1/2 2k+1/2 |ΘN | ≤ Ch ∥φ∥2k+1,Ω + ∥F∥2k+1,Ωx + ∥D∥2k+1,Ωx dt 0 where C depends on the upper bounds of ∥f ∥k+2,Ω , ∥f ∥1,∞,Ω , ∥E∥0,∞,Ωx , ∥B∥0,∞,Ωx , ∥E∥k+2,Ωx , ∥B∥k+2,Ωx over the time interval [0, T ], and it also depends on the polynomial degree k, mesh parameters σ0 , σx and σv , and domain parameters Lx and Lv . Proof. Due to the definition of the projection operators, ((fh )t , φ − χ)Ω = 0, ((Eh )t , F − ξ)Ωx = 0, and ((Bh )t , D − η)Ωx = 0, and lh (Jh ; F − ξ) = −(Jh , F − ξ)Ωx = 0, we have Z T ΘN = −ah (fh , Eh , Bh ; ζhφ ) − bh (Eh , Bh ; ζhF , ζhD ) dτ. 0 20 From its definition, Z Z bh (Eh , Bh ; ζhF , ζhD ) = Eh · ∇ x × ζhD dx − Bh · ∇x × ζhF dx Thx Thx Z Z + Efh · [ζ D ]τ dsx − fh · [ζ F ]τ dsx B h h Ex Ex Z Z F =− eE D h · ∇x × ζh dx + eBh · ∇x × ζh dx Thx Thx Z Z E D B F − efh · [ζh ]τ dsx + efh · [ζh ]τ dsx Ex ZEx Z D + (∇x × E) · ζh dx − (∇x × B) · ζhF dx. Thx Thx By Lemma 2.3.1, Z (eE D k E h ) · ∇x × ζh dx ≤ Ch ∥eh ∥0,Ωx ∥D∥k+1,Ωx , Thx Z k B eB F h · ∇x × ζh dx ≤ Ch ∥eh ∥0,Ωx ∥F∥k+1,Ωx , Thx Z E D B F k+1/2 h ) · [ζh ]τ − (eh ) · [ζh ]τ dsx ≤ Ch (ef f (∥D∥k+1,Ωx + ∥F∥k+1,Ωx ) Ex × ∥eE B  h ∥0,Ex + ∥eh ∥0,Ex . Now notice that E E ∥eEh ∥0,Ex ≤ ∥ϵh ∥0,Ex + ∥ζh ∥0,Ex ≤ C[h−1/2 ∥ϵE h ∥0,Ωx + h k+1/2 ] ≤ Ch−1/2 [∥eE h ∥0,Ωx + h k+1 ]. Analogously −1/2 ∥eBh ∥0,Ex ≤ Ch [∥eBh ∥0,Ωx + h k+1 ]. Therefore, Z E D B F k h ) · [ζh ]τ − (eh ) · [ζh ]τ dsx ≤ Ch (∥D∥k+1,Ωx + ∥F∥k+1,Ωx ) (ef f Ex × ∥eE B k+1  h ∥0,Ωx + ∥eh ∥0,Ωx + h . 21 Now by the properties of the orthogonal projection Πkx Z Z D (∇x × E) · ζh dx = (∇x × E − Πkx (∇x × E)) · ζhD dx ≤ Ch2k+2 ∥D∥k+1,Ωx , Thx Thx where C depends on ∥E∥k+2,Ωx . By an analogous procedure Z (∇x × B) · ζhF dx ≤ Ch2k+2 ∥F∥k+1,Ωx , Thx where C depends on ∥B∥k+2,Ωx . Putting all the above calculations together, we arrive at, |bh (Eh , Bh ; ζhF , ζhD )| ≤ Chk (∥D∥k+1,Ωx + ∥F∥k+1,Ωx ) ∥eE B k+1  h ∥0,Ωx + ∥eh ∥0,Ωx + h , (2.29) where C depends on ∥E∥k+2,Ωx , ∥B∥k+2,Ωx . We will deal now with the term ah , which is ah (fh , Eh , Bh , ζhφ ) = ah,1 (fh , ζhφ ) + ah,2 (fh , Eh , Bh ; ζhφ ). (2.30) First, we have Z Z Z Z ah,1 (fh ; ζhφ ) = efh v · ∇x ζhφ dxdv + efh v[ζhφ ]x dsx dv − ∇x f · vζhφ dxdv g Th Thv Ex Th The first term can be easily bounded, by using Lemma 2.3.1. Z efh v · ∇x ζhφ dxdv ≤ Chk ∥efh ∥0,Ω ∥φ∥k+1,Ω . Th Similarly, Z Z efh v[ζhφ ]x dsx dv ≤ C∥efh ∥Thv ×Ex ∥ζhφ ∥Thv ×Ex g Thv Ex ≤ Chk+1/2 ∥efh ∥Thv ×Ex ∥φ∥k+1,Ω ≤ Chk+1/2 (∥ϵfh ∥Thv ×Ex + ∥ζhf ∥Thv ×Ex )∥φ∥k+1,Ω ≤ Chk (∥efh ∥0,Ω + hk+1 )∥φ∥k+1,Ω . For the last term notice that by the properties of the projection Πk and the fact that Πk (∇x f · v) is a polynomial of degree k, Z Z φ ∇x f · vζh dxdv = (∇x f · v − Πk (∇x f · v))ζhφ dxdv Th Th ≤ Ch2k+2 ∥φ∥k+1,Ω , 22 where C depends on ∥f ∥k+2,Ω . By using all the calculations above, we can conclude that |ah,1 (fh ; ζhφ )| ≤ Chk ∥efh ∥0,Ω ∥φ∥k+1,Ω + Ch2k+1 ∥φ∥k+1,Ω , (2.31) where C depends on ∥f ∥k+2,Ω . To conclude our proof, we only need to bound ah,2 , this time we will do things a little bit different, notice that ah,2 (fh , Eh , Bh , ζhφ ) = ah,2 (f, Eh , Bh , ζhφ ) − ah,2 (efh , Eh , Bh , ζhφ ). We will get started by noting that f (Eh ^ + v × Bh ) = f {Eh + v × Bh }v = f (Eh + v × Bh ), then Z ah,2 (f, Eh , Bh , ζhφ ) =− f (Eh + v × Bh ) · ∇v ζhφ dxdv Z TZh + f (Eh + v × Bh ) · [ζhφ ]v dxdv Tx Ev Zh Z Z = f (eEh +v× eB h) · ∇v ζhφ dxdv − f (eE B φ h + v × eh ) · [ζh ]v dxdv Th Thx Ev Z + ∇v f · (E + v × B)ζhφ dxdv. Th f (E + v × B) · ∇v ζhφ dxdv , R We obtained the last inequality by adding and subtracting Th integration by parts, and the fact that ∇v · (E + v × B) = 0. in this way Z φ f (eE B k E B h + v × eh ) · ∇v ζh dxdv ≤ Ch (∥eh ∥0,Ωx + ∥eh ∥0,Ωx )∥φ∥k+1,Ω , Th and Z Z φ f (eE B h + v × eh ) · [ζh ]v dvdx ≤ Ch k+1/2 (∥eE B h ∥0,Ωx + ∥eh ∥0,Ωx )∥φ∥k+1,Ω . Thx Ev Last but not least by the same arguments as previous estimates Z Z φ ∇v f · (E + v × B)ζh dxdv = (∇v f · (E + v × B) − Πk ∇v f · (E + v × B))ζhφ dxdv Th Th ≤ Ch2k+2 ∥φ∥k+1,Ω , where C depends on ∥f ∥k+2,Ω , ∥E∥k+1,Ωx , ∥B∥k+1,Ωx . We can conclude that |ah,2 (f, Eh , Bh ; ζhφ )| ≤ Chk (∥eE B h ∥0,Ωx + ∥eh ∥0,Ωx )∥φ∥k+1,Ω + Ch 2k+2 ∥φ∥k+1,Ω . (2.32) 23 Finally we just need to estimate Z ah,2 (efh , Eh , Bh ; ζhφ ) =− efh (Eh + v × Bh ) · ∇v ζhφ dxdv Z TZh + efh (Eh^+ v × Bh ) · [ζhφ ]v dsv dx Thx Ev We have Z efh (Eh + v × Bh ) · ∇v ζhφ dxdv ≤ Chk ∥efh ∥0,Ω (∥Eh ∥0,∞,Ωx + ∥Bh ∥0,∞,Ωx )∥φ∥k+1,Ω Th ≤ Chk ∥efh ∥0,Ω (∥ϵE B k k h ∥0,∞,Ωx + ∥ϵh ∥0,∞,Ωx + ∥Πx E∥0,∞,Ωx + ∥Πx B∥0,∞,Ωx )∥φ∥k+1,Ω ≤ Chk−dx /2 ∥efh ∥0,Ω (∥ϵE B h ∥0,Ωx + ∥ϵh ∥0,Ωx )∥φ∥k+1,Ω + Chk ∥efh ∥0,Ω (∥E∥0,∞,Ωx + ∥B∥0,∞,Ωx )∥φ∥k+1,Ω ≤ Chk−dx /2 ∥efh ∥0,Ω (∥eE B h ∥0,Ωx + ∥eh ∥0,Ωx + h k+1 )∥φ∥k+1,Ω + Chk ∥efh ∥0,Ω ∥φ∥k+1,Ω ≤ Chk ∥efh ∥0,Ω (h−dx /2 ∥eE h ∥0,Ωx + h −dx /2 B ∥eh ∥0,Ωx + 1)∥φ∥k+1,Ω , Here we used the fact that whenever dx = 1, 2, 3, k + 1 − dx /2 > 0, Lemma 2.3.2 and the fact that Πx is bounded in any Lp -norm (1 ≤ p ≤ ∞) [49, 50], ∥Πx E∥0,∞,Ωx ≤ C∥E∥0,∞,Ωx , ∥Πx B∥0,∞,Ωx ≤ C∥B∥0,∞,Ωx . Finally Z Z efh (Eh^ + v × Bh ) · [ζhφ ]v dsv dx Thx Ev ≤ Chk+1/2 (∥Eh ∥0,∞,Ωx + ∥Bh ∥0,∞,Ωx )∥efh ∥0,Thx ×Ev ∥φ∥k+1,Ω ≤ Chk+1/2 (∥Eh ∥0,∞,Ωx + ∥Bh ∥0,∞,Ωx )h−1/2 (∥efh ∥0,Th + hk+1 ∥f ∥k+1,Ω )∥φ∥k+1,Ω ≤ Chk (∥Eh ∥0,∞,Ωx + ∥Bh ∥0,∞,Ωx )(∥efh ∥0,Th + hk+1 ∥f ∥k+1,Ω )∥φ∥k+1,Ω ≤ Chk (∥efh ∥0,Th + hk+1 )(h−dx /2 ∥eE h ∥0,Ωx + h −dx /2 B ∥eh ∥0,Ωx + 1)∥φ∥k+1,Ω . In this way we conclude that |ah,2 (eh , Eh , Bh ; ζhφ )| ≤ Chk (∥efh ∥0,Th + hk+1 )(h−dx /2 ∥eE h ∥0,Ωx + h −dx /2 B ∥eh ∥0,Ωx + 1)∥φ∥k+1,Ω (2.33) 24 Then by putting together (2.29), (2.31), (2.32), (2.33), and using Theorem 2.1.1, we have |ah (fh , Eh , Bh ; ζhφ ) + bh (Eh , Bh ; ζhF , ζhD )| ≤ Chk (∥D∥k+1,Ωx + ∥F∥k+1,Ωx + ∥φ∥k+1,Ω ) ∥eE B k+1  h ∥0,Ωx + ∥eh ∥0,Ωx + h + Chk (∥efh ∥0,Th + hk+1 )(h−dx /2 ∥eE h ∥0,Ωx + h −dx /2 B ∥eh ∥0,Ωx + 1)∥φ∥k+1,Ω ≤ Ch2k+1/2 (∥D∥k+1,Ωx + ∥F∥k+1,Ωx + ∥φ∥k+1,Ω ) . where we have used k + 1/2 − dx /2 > 0. An application of Cauchy-Schwarz inequality concludes the proof. Lastly, we need to estimate the third term, ΘC . Lemma 2.3.6 (Consistency). We have Z T 1/2 2k+1 |ΘC | ≤ Ch ∥φ∥2k+1,Ω dt (2.34) 0 where C depends on the upper bounds of ∥∂t f ∥k+1,Ω ,∥f ∥k+1,Ω , |f |1,∞,Ω , ∥E∥1,∞,Ωx , ∥B∥1,∞,Ωx , ∥E∥k+1,Ωx , ∥B∥k+1,Ωx over the time interval [0, T ], and it also depends on the polynomial degree k, mesh parameters σ0 , σx and σv , and domain parameters Lx and Lv . Proof. The terms inside the integral of ΘC can be split in I + II, where I = (fh , φt )Ω − ah (fh , Eh , Bh ; φ) + lh (Jh , F) II = (Eh , Ft )Ωx + (Bh , Dt )Ωx − bh (Eh , Bh ; F, D)+F(f, E, B; φ) since φ is a smooth function, [φ]x = 0 and [φ]v = 0, in this way, by using (2.23a), and the definition of lh , we conclude that, I = (fh , −v · ∇x φ − (E + v × B) · ∇v φ + v · F)Ω − ah (fh , Eh , Bh ; φ) + lh (Jh ; F) Z Z =− fh v · ∇x φ dxdv − fh (E + v × B) · ∇v φ dxdv − lh (Jh ; F) Th Ω Z Z + fh v · ∇x φ dx + fh (Eh + v × Bh ) · ∇v φ dxdv + lh (Jh ; F) Th Ω Z = − fh (eE B h + v × eh ) · ∇v φ dx dv. Ω 25 On the other hand, by using (2.23b) and (2.23c), since F and D are smooth functions [F]τ = [D]τ = 0,we have that II = (Eh , ∇x × D)Thx − (Bh , ∇x × F)Thx − bh (Eh , Bh ; F, D) + F(f, E, B; φ) Z Z − f Eh · ∇v φ dxdv + f Bh · (v × ∇v φ) dxdv Ω Ω = (Eh , ∇x × D)Thx − (Bh , ∇x × F)Thx − (Eh , ∇x × D)Thx + (Bh , ∇x × F)Thx Z Z − f (Eh + v × Bh ) · ∇v φ dxdv + f (E + v × B) · ∇v φ dxdv ZΩ Ω = f (eE B h + v × eh ) · ∇v φ dxdv. Ω We obtain Z I + II = efh (eE B h + v × eh ) · ∇v φ dxdv Ω ≤ C∥efh ∥Ω (∥eE B h ∥Ωx + ∥eh ∥Ωx )∥∇v φ∥∞,Ω ≤ C∥efh ∥Ω (∥eE B h ∥Ωx + ∥eh ∥Ωx )∥φ∥k+1,Ω where we used the Sobolev inequality [51], ∥∇v φ∥∞,Ω ≤ C∥φ∥k+1,Ω , which requires k > (dx + dv )/2. Using Theorem 2.1.1, we conclude the proof. It is easy to transform the dual problem (2.23) to an initial value problem (2.21) by changing time t′ = T − t. Then using Lemma 2.3.3, where A1 (x, v, t) = −v, A2 (x, v, t) = −(E + v × B), A3 (x, v, t) = v, g = f and l = k + 1, ∥φ∥2k+1,Ω + ∥F∥2k+1,Ωx + ∥D∥2k+1,Ωx ≤ C[∥Φ∥2k+1,Ω + ∥F∥2k+1,Ωx + ∥D∥2k+1,Ωx ] (2.35) where C depends on ∥f ∥L∞ ((0,T );W k+2,∞ (Ω)) . Then an application of Theorem 2.1.1 gives us q |(efh (T ), Φ)Ω + (eE h (T ), F) Ωx + (e B h (T ), D) Ωx | ≤ Ch 2k+1/2 ∥Φ∥2k+1,Ω + ∥F∥2k+1,Ωx + ∥D∥2k+1,Ωx (2.36) Therefore the estimate for the zero-divided difference negative-order norm is given by ∥(f − fh , E − Eh , B − Bh )∥−(k+1),Ω (f − fh , Φ)Ω + (E − Eh , F)Ωx + (B − Bh , D)Ωx = sup q ≤ Ch2k+1/2 . ϕ∈C0∞ (Ω),F,D∈[C ∞ (Ωx )]dx ∥Φ∥2k+1,Ω + ∥F∥2k+1,Ωx + ∥D∥2k+1,Ωx 26 2.4 Numerical Experiments In this section, we validate our theoretical results using several numerical tests. In particular, we want to demonstrate the performance of the post-processing technique for the VA system and the VM system. We heavily use the fact that the VM (VA) system is time reversible to provide quantitative measurements of the errors. In particular, let f (x, v, 0), E(x, 0), B(x, 0) denote the initial conditions and f (x, v, T ), E(x, T ), B(x, T ) be the solution of the VM system at t = T . If we choose f (x, −v, T ), E(x, T ), −B(x, T ) as the initial condition at t = 0, then evolving the VM system to t = T, we will recover f (x, −v, 0), E(x, 0), −B(x, 0). 2.4.1 Vlasov-Ampére examples We consider two classical benchmark examples. • Landau damping: f (x, v, 0) = fM (v)(1 + A cos(kx)), x ∈ [0, L], v ∈ [−Vc , Vc ], (2.37) 2 where A = 0.5, k = 0.5, L = 4π, Vc = 6π, and fM (v) = √1 e−v /2 . 2π • Two-stream instability: f (x, v, 0) = fT S (v)(1 + A cos(kx)), x ∈ [0, L], v ∈ [−Vc , Vc ], (2.38) 2 where A = 0.05, k = 0.5, L = 4π, Vc = 6π, and fT S (v) = √1 v 2 e−v /2 . 2π Notice that in both examples we have taken Vc to be larger than the usual values in the literature in order to completely eliminate the boundary effects and accurately reflect the accuracy enhancement property. In Tables 2.1, we run the VA system with initial condition from Landau damping to T = 1 and then back to T = 0 and then we apply the SIAC filter, and compare it with the initial conditions. We use the third order TVD-RK method as the time integrator [52]. To make sure the spatial error dominates, we take ∆t = CFL/(Vc /∆x + Emax /∆v) for P1 , 27 Emax denotes the maximum value of E(·, T ) in Ωx , for P2 we take ∆t = CFL/(Vc /(∆x)5/3 + Emax /(∆v)5/3 ), and ∆t = CFL/(Vc /(∆x)7/3 + Emax /(∆v)7/3 ) for P3 . For P1 and P3 we take the CFL = 0.1, and we take the CFL = 0.2 for P2 . From the table, we observe (k + 1)-th order of convergence for the DG solution before post-processing for both f and E. We can clearly see that we improve the order of the error to at least O(h2k+1/2 ) after post- processing. In Figure 2.1 we plot the errors of the numerical solution before and after post-processing for P1 and using 128 × 128 elements. We can see that the errors before post-processing are highly oscillatory, and that the post-processing smooths out the error and greatly reduces its magnitude. In Figure 2.2, we plot the errors of the approximations for E obtained when solving using a 128 × 128 mesh with P1 and 32 × 32 mesh with P3 . We can clearly see that the errors before post-processing are highly oscillatory, and the post-processing gets rid of the oscillations and dramatically reduces the magnitude of the error. Another point that we want to make is the following: if we look at Table 2.1, for k = 2 and a mesh of 64 × 64, the L2 -errors before and after post-processing are similar in magni- tude. However, if we look at Figure 2.3 which plots the absolute value of the error in f in this case, we can clearly see that the L∞ -norm of the error of the filtered solution is much smaller than the unfiltered solution. Therefore, by removing the spurious oscillations, even if the L2 -error is comparable, the L∞ error is further reduced by the post-processor. This is probably due to the high oscillatory nature of the solution. In Tables 2.2, we run the VA system with initial condition from two stream instability to T = 1 and then back to T = 0 and then we apply the SIAC filter, and compare it with the initial conditions. To integrate in time we used Fourth order Runge-Kutta for P1 and P3 and third order TVD-RK method for P2 . For P1 and P2 we take ∆t just as in the Landau damping example, and ∆t = CFL/(Vc /(∆x)7/4 + Emax /(∆v)7/4 ) for P3 . We use CFL = 0.2 for all cases. We can observe (k + 1/2)-order of convergence for the DG solution before post-processing for both f and E. Just as in the Landau damping example we can see the 28 Before post-processing After post-processing mesh error f order error E order error f ∗ order error E∗ order P1 16 × 16 1.42E-02 - 1.19E-02 - 2.28E-02 - 1.04E-02 - 32 × 32 6.22E-03 1.19 3.16E-03 1.91 6.16E-03 1.89 2.84E-03 1.88 64 × 64 1.59E-03 1.97 5.65E-04 2.48 8.74E-04 2.82 4.36E-04 2.70 128 × 128 4.08E-04 1.96 1.12E-04 2.33 1.10E-04 2.99 6.31E-05 2.79 256 × 256 1.03E-04 1.98 2.51E-05 2.16 1.37E-05 3.00 9.01E-06 2.81 512 × 512 2.60E-05 1.99 6.14E-06 2.03 1.71E-06 3.00 1.71E-06 2.39 P2 16 × 16 7.08E-03 - 1.97E-03 - 2.09E-02 1.88E-03 - 32 × 32 1.08E-03 2.71 1.13E-04 4.12 2.87E-03 2.87 1.08E-04 4.12 64 × 64 1.35E-04 3.00 6.62E-06 4.10 1.20E-04 4.58 5.15E-06 4.39 128 × 128 1.63E-05 3.04 5.59E-07 3.57 2.70E-06 5.47 2.04E-07 4.66 256 × 256 2.01E-06 3.03 6.57E-08 3.09 5.29E-08 5.67 5.75E-09 5.15 P3 16 × 16 1.73E-03 - 2.19E-04 - 2.16E-02 - 9.71E-05 - 32 × 32 1.52E-04 3.51 7.18E-06 4.93 2.60E-03 3.05 3.09E-06 4.97 64 × 64 1.06E-05 3.84 1.30E-07 5.79 5.65E-05 5.52 7.52E-08 5.36 128 × 128 6.45E-07 4.04 3.42E-09 5.25 3.95E-07 7.16 8.24E-10 6.51 Table 2.1 L2 errors for the numerical solution and the post-processed solution for Landau Damping. Figure 2.1 Errors for f before (on the left) and after post-processing (on the right) for 128 × 128 elements and P1 . Landau damping. 29 (a) 128 × 128 and P1 (b) 32 × 32 and P3 Figure 2.2 Errors before (solid line) and after post-processing (dashed line) for E for dif- ferent mesh sizes and Pk . Landau damping. T = 2. Figure 2.3 Absolute value of errors for f before (on the left) and after post-processing (on the right) for 64 × 64 elements and P2 . Landau damping. improvement of the error to O(h2k+1/2 ) after post-processing. In Figure 2.4 we plot the errors of the numerical solution before and after post-processing for P1 and using 128×128 elements. We can arrrive to similar conclusions as in the Landau damping case. Figure 2.5, we plot the errors of the approximations for E obtained when solving using a 16 × 16 mesh with P2 and 32 × 32 mesh with P3 . We can clearly see how the SIAC filter, gets rid of the spurious oscillations and dramatically reduce the magnitude of the error. Now we provide plots comparing the solution profile before and after post-processing 30 Before post-processing After post-processing mesh error f order error E order error f ∗ order error E∗ order P1 16 × 16 2.63E-02 - 3.24E-03 - 3.46E-02 - 3.12E-03 - 32 × 32 5.58E-03 2.23 2.13E-04 3.93 9.83E-03 1.82 1.50E-04 4.38 64 × 64 2.23E-03 1.32 4.54E-05 2.23 1.19E-03 3.05 2.45E-05 2.62 128 × 128 6.08E-04 1.88 1.04E-05 2.13 8.91E-05 3.74 3.66E-06 2.74 256 × 256 1.78E-04 1.77 2.49E-06 2.06 7.06E-06 3.66 5.15E-07 2.83 512 × 512 4.98E-05 1.84 6.13E-07 2.02 6.80E-07 3.38 6.30E-08 3.03 P2 16 × 16 7.60E-03 - 1.07E-04 3.54E-02 3.91E-05 - 32 × 32 2.36E-03 1.69 5.45E-06 4.30 9.81E-03 1.85 4.88E-06 3.00 64 × 64 2.96E-04 2.99 4.18E-07 3.70 5.48E-04 4.16 3.22E-07 3.92 128 × 128 4.69E-05 2.66 3.67E-08 3.51 1.28E-05 5.42 1.14E-08 4.81 256 × 256 7.15E-06 2.71 4.72E-09 2.96 2.23E-07 5.84 5.16E-10 4.47 P3 16 × 16 5.06E-03 - 2.31E-05 - 3.68E-02 - 1.88E-05 - 32 × 32 1.09E-04 5.54 1.18E-06 4.29 1.02E-02 1.85 8.97E-08 7.71 64 × 64 2.69E-05 2.01 1.99E-08 5.89 3.47E-04 4.88 6.56E-09 3.77 128 × 128 2.01E-06 3.75 2.74E-10 6.18 2.83E-06 6.94 9.68E-11 6.08 Table 2.2 L2 errors for the numerical solution and the post-processed solution. Two stream instability. Figure 2.4 Errors for f before (on the left) and after post-processing (on the right) for 128 elements and P1 , using two-stream instability. 31 (a) 16 × 16 and P2 (b) 32 × 32 and P3 Figure 2.5 Errors before (solid line) and after post-processing (dashed line) for E for dif- ferent mesh sizes and Pk . Two stream instability. T = 2. for a longer computational time. To compute those plots, we use a third-order Runge- Kutta method with ∆t = CFL/(Vc /∆x + Emax /∆v)) and CFL = 0.1. In Figures 2.6 to 2.9, we show a comparison of contour plots of the numerical solution for f before and after post-processing with different mesh size and k = 1, 2. There is visible improvement of the resolution of the solution, particularly for k = 1. 2.4.2 Vlasov-Maxwell example In this part, we will test our post-processor for the VM system. Specifically we will use the streaming Weibel (SW) instability as an example. This is a reduced version of the VM equations with one spatial variable, x2 , and two velocity variables v1 and v2 . The variables under consideration are the distribution function f (x2 , v1 , v2 , t), a 2D electric field E = (E1 (x2 , t), E2 (x2 , t), 0) and a 1D magnetic field B = (0, 0, B3 (x2 , t)) and the reduced VM system reads as ∂t f + v2 fx2 + (E1 + v2 B3 )fv1 + (E2 − v1 B3 )fv2 = 0, (2.39a) ∂B3 ∂E1 ∂E1 ∂B3 ∂E2 = , = − j1 , = −j2 , (2.39b) ∂t ∂x2 ∂t ∂x2 ∂t 32 (a) 16 × 16 (b) 32 × 32 (c) 64 × 64 Figure 2.6 Comparison of contour plots before (left) and after post-processing (right) for different mesh-sizes. Landau damping, k = 1 and T = 10. 33 (a) 16 × 16 (b) 32 × 32 (c) 64 × 64 Figure 2.7 Comparison of the definition of the contour plots before (left) and after post- processing (right) for different mesh-sizes. Landau damping, k = 2 and T = 10. 34 (a) 16 × 16 (b) 32 × 32 (c) 64 × 64 Figure 2.8 Comparison of contour plots before (left) and after post-processing (right) for different mesh-sizes. Two stream instability, k = 1 and T = 20. 35 (a) 16 × 16 (b) 32 × 32 (c) 64 × 64 Figure 2.9 Comparison of the definition of the contour plots before (left) and after post- processing (right) for different mesh-sizes. Two stream instability, k = 2 and T = 20. 36 where Z ∞ Z ∞ Z ∞ Z ∞ j1 = f (x2 , v1 , v2 , t)v1 dv1 dv2 , j2 = f (x2 , v1 , v2 , t)v2 dv1 dv2 . (2.40) −∞ −∞ −∞ −∞ The initial conditions are given by 1 −v22 /β −(v1 −ω0,1 )2 /β 2 f (x2 , v1 , v2 , 0) = e [δe + (1 − δ)e−(v1 +ω0,2 ) /β ], (2.41a) πβ E1 (x2 , v1 , v2 , 0) = E2 (x2 , v1 , v2 , 0) = 0, B3 (x2 , v1 , v2 , 0) = b sin(κ0 x2 ), (2.41b) which for b = 0 is an equilibrium state composed of counter-streaming beams propagat- ing perpendicular to the direction of inhomogeneity. Following [6, 19], we trigger the instability by taking β = 0.01 , b = 0.001 (the amplitude of the initial perturbation of the magnetic field). Here, Ωx = [0, Ly ] , where Ly = 2π/κ0 , and we set Ωv = [−1.8, 1.8]2 . We consider the following set of parameters, δ = 0.5, ω0,1 = ω0,2 = 0.3, κ0 = 0.2. In Table 2.3, we run the VM system with initial condition from SW instability to T = 5 and then back to T = 0, we then apply the SIAC filter and compare it with the initial conditions. We use a third order TVD-RK method as the time integrator. To make sure the spatial error dominates, we take ∆t = O(∆x) for P1 and ∆t = O(∆x5/3 ) for P2 , in both cases we used CFL = 0.1. From the table we can observe (k + 1)-th order of convergence for the DG solution before post-processing for f , E1 , E2 and B3 . After post-processing we can see overall the order of convergence improves to O(h2k+1/2 ). In Figure 2.10 we plot a cross-section of the errors of the numerical solution at x2 ≈ 0.15π before and after post-processing for P1 using 80 × 80 × 80 elements. We can see that before post-processing that the errors are highly oscillatory, and after post-processing the error surface is smooth out and the error is much smaller in magnitude. In Figure 2.11 we plot the errors of E1 , E2 and B3 , we used the same number of elements as in Figure 2.10, We can clearly see similar conclusions. 37 Before post-processing Mesh Error f Order Error B3 Order Error E1 Order Error E2 Order P1 20 × 20 × 20 2.20E-01 - 2.61E-06 - 2.12E-06 - 5.31E-06 - 40 × 40 × 40 7.17E-02 1.61 6.54E-07 2.00 7.06E-07 1.58 5.46E-07 3.28 80 × 80 × 80 1.92E-02 1.90 1.63E-07 2.00 1.96E-07 1.85 7.05E-08 2.95 160 × 160 × 160 4.89E-03 1.98 4.07E-08 2.00 5.13E-08 1.94 6.40E-09 3.46 P2 20 × 20 × 20 1.07E-01 - 2.56E-07 - 2.49E-07 - 1.02E-06 40 × 40 × 40 1.64E-02 2.70 3.14E-08 3.03 2.93E-08 3.09 9.72E-08 3.40 80 × 80 × 80 2.23E-03 2.88 1.63E-09 4.27 1.90E-09 3.95 6.93E-09 3.81 160 × 160 × 160 2.92E-04 2.93 1.41E-10 3.52 1.72E-10 3.46 2.46E-10 4.81 After post-processing Mesh Error f ∗ Order Error B∗3 Order Error E∗1 Order Error E∗2 Order P1 20 × 20 × 20 2.95E-01 - 3.17E-07 - 1.08E-07 - 5.08E-06 - 40 × 40 × 40 6.13E-02 2.27 7.16E-08 2.14 1.49E-08 2.87 4.38E-07 3.54 80 × 80 × 80 5.87E-03 3.38 1.12E-08 2.68 3.11E-09 2.26 6.33E-08 2.79 160 × 160 × 160 4.19E-04 3.81 2.01E-09 2.48 7.47E-10 2.06 6.22E-09 3.35 P2 20 × 20 × 20 2.89E-01 - 1.24E-08 - 9.06E-09 - 4.41E-07 - 40 × 40 × 40 4.58E-02 2.66 5.61E-10 4.46 2.97E-10 4.93 2.63E-08 4.07 80 × 80 × 80 2.03E-03 4.49 2.94E-11 4.25 1.31E-11 4.50 2.57E-09 3.36 160 × 160 × 160 4.43E-05 5.52 1.65E-12 4.15 5.55E-13 4.56 1.12E-10 4.53 Table 2.3 L2 errors for the numerical solution (Above) and the post-processed solution (Below). SW instability. Figure 2.10 Cross-sectional plot of the error for f at x2 ≈ 0.15π, before (on the left) and after post-processing (on the right) for 803 elements and P1 . SW instability. 38 (a) Error for B3 . (b) Error for E1 . (c) Error for E2 . Figure 2.11 Errors before (solid line) and after post-processing (dashed line) for the differ- ent fields using mesh size of 80 × 80 × 80 and P1 . T = 10. SW instability. 39 CHAPTER 3 NUMERICAL ANALYSIS OF A HYBRID METHOD FOR RADIATION TRANSPORT In this Chapter we will analyze a Hybrid method for Radiation transport introduced in [2]. The remaining of the chapter is organized as follows: In Section 3.1, we introduce the RTE, reduce it to the purely scattering problem, recall the PN method, and then describe the setup of the hybrid. Having established the setting of of the problem, we then sum- marize the main results of this chapter . In Section 3.2, we derive error estimates for the PN equations. In Section 3.3, we analyze the hybrid problem. In Section 3.4, we generalize results back to the original RTE with non-zero absorption. The appendix contains some generic results used for the estimates of this chapter. 3.1 Background 3.1.1 The radiation transport equation We consider a time-dependent transport equation with periodic boundaries, isotropic scattering, unit-speed particles, and diffusion scaling: Z σt ε  σt  1 ε ε ε∂t Ψ + Ω · ∇x Ψ + Ψ = − εσa Ψε + εQ, Ψε = Ψε dΩ, (3.1a) ε ε 4π S2 Ψε |t=0 = g. (3.1b) Here Ψε = Ψε (x, Ω, t) is a function of position x ∈ X = [0, 2π)3 , direction of flight Ω ∈ S2 , and time t > 0. It can be interpreted physically as the density of particles at time t with respect to the measure dΩdx. Particles interact with a material background characterized by an absorption cross-section σa ≥ 0, total cross-section σt ≥ σa (which accounts for scat- tering and absorption), and a known source Q = Q(x, Ω, t). The quantity σt − ε2 σa is the scaled scattering cross-section, where the non-dimensional parameter ϵ > 0 characterizes the strength of the scattering as well as the relevant time scale. Indeed, it is well-known [32] that in the limit ε → 0, Ψε → Ψ0 where Ψ0 is independent of angle and satisfies the 40 diffusion equation of the form   Z 1 1 0 ∂t Ψ − ∇x · ∇x Ψ + σa Ψ0 = 0, 0 Ψ0 t=0 = g dΩ. (3.2a) 3σt 4π S2 For g ∈ L2 (X × S2 ), σt , σa ∈ L∞ (X), Q ∈ L2 (X × S2 × [0, ∞)) (3.1), is known to have a semi-group solution Ψε ∈ C([0, ∞); L2 (X × S2 )) [53, Theorem XXI.2.3]. If in addition, Ω · ∇x g ∈ L2 (X × S2 ), then Ψε ∈ C 1 ([0, ∞); L2 (X × S2 )). We assume this is the case for remainder of the chapter. In order to facilitate a clear stability and error analysis, we assume that the cross- sections σa and σt are constant in space. This assumption on σa allows us to convert (3.1) to a purely scattering system for the function ψ = eσa t Ψε :1 Z σ σ 1 ε∂t ψ + Ω · ∇x ψ + ψ = ψ + εq, ψ= ψ dΩ, (3.3a) ε ε 4π S2 ψ|t=0 = g, (3.3b) where q = eσa t Q and σ := σt . Henceforth, we focus our analysis on (3.3). The results can then be translated back to the case of non-zero absorption by undoing the transformation, which gives exponential decay if σa > 0. This assumption is made for simplicity, but it does introduce a measure of regularity into the solution that is not typical in applications. Indeed, a more reasonable assumption is that the cross-sections are piece-wise smooth and that the boundaries are equipped with inflow data. Hence the analysis here can be viewed as a localized proxy for a more realistic scenario. A more sophisticated analysis to include boundary and interior layers is the subject of future work. 3.1.2 The PN approximation Given N ∈ N≥0 , the PN method is a spectral discretization of the transport equation with respect to the angular variable Ω. Let {mℓ,k }ℓ,k be the real-valued, orthonormal basis of spherical harmonics, where ℓ ≥ 0 denotes the degree and k ∈ {−ℓ, . . . , ℓ} denotes the 1 To reduce notation, we suppress the dependence of ψ on ε. 41 order. For any u ∈ L2 (S2 ), the angular moment uℓ,k is given by Z uℓ,k = mℓ,k u dΩ. (3.4) S2 For convenience, we collect the basis elements of degree ℓ into vectors mℓ = (mℓ,−ℓ , . . . , mℓ,ℓ )T , and we denote by uℓ = (uℓ,−ℓ , . . . , uℓ,ℓ )T the vector of corresponding moments. Let PN (S2 ) ⊂ L2 (S2 ) to be the span of all spherical harmonics with degree at most N . Then the orthog- onal projections PN : L2 (S2 ) → PN (S2 ) and P eN : L2 (S2 ) → PN (S2 ) are given by X N XN X ℓ X ∞ X ℓ PN u = mTℓ uℓ = mℓ,k uℓ,k and P eN u = (I − PN )u = mℓ,k , (3.5) ℓ=0 ℓ=0 k=−ℓ ℓ=N +1 k=−ℓ where I is the identity operator. The PN approximation of (3.3) seeks a function X N X N X ℓ N ψ (x, Ω, t) = mTℓ (Ω)ψN ℓ (x, t) = mℓ,k (Ω)ψNℓ,k (x, t) (3.6) ℓ=0 ℓ=0 k=−ℓ such that Z σ σ 1 ε∂t ψ + PN (Ω · ∇x ψ ) + ψ N = ψ N + εPN q, N N ψN = ψ N dΩ, (3.7a) ε ε 4π S2 ψ N |t=0 = PN g. (3.7b) When expressed in terms of the moments ψN ℓ , the PN method yields the following linear, symmetric hyperbolic system: X3 (i) ε∂t ψN 0 + a1 ∂xi ψN 1 = εq0 , for ℓ = 0, (3.8a) i=1 3 3 X (i) X (i) σ N ε∂t ψN ℓ + (aℓ )T ∂xi ψN ℓ−1 + aℓ+1 ∂xi ψN ℓ+1 + ψℓ = εqℓ , for 1 ≤ ℓ ≤ N − 1, (3.8b) i=1 i=1 ε 3 X (i) σ ε∂t ψN N + (aN )T ∂xi ψN N −1 + ψN = εqN , for ℓ = N. (3.8c) i=1 ε N (i) Formulas for the elements in the matrices aℓ ∈ R(2ℓ−1)×(2ℓ+1) can be found in the appendix of [44]. In the current work, we rely only on the fact the they are bounded in the operator (i) norm, specifically that ∥aℓ ∥2 ≤ 4. 42 R The exact moments ψℓ = S2 mℓ ψdΩ satisfy an infinite system of equations with a structure similar to (3.8): X 3 (i) ε∂t ψ0 + a1 ∂ x i ψ 1 = εq0 , for ℓ = 0, (3.9) i=1 3 3 X (i) X (i) σ ε∂t ψℓ + (aℓ )T ∂xi ψℓ−1 + aℓ+1 ∂xi ψℓ+1 + ψℓ = εqℓ , for ℓ ≥ 1. (3.10) i=1 i=1 ε In particular, the PN equations for ψN can be obtained by truncating (3.9) at ℓ = N and then neglecting the moment ψN N +1 that would otherwise appear in (3.8c). 3.1.3 The hybrid method The hybrid method is based on a separation of ψ into a collided component ψc and an uncollided component ψu . These components satisfy the coupled system σ ε∂t ψu + Ω · ∇x ψu + ψu = εq, (3.11a) ε σ σ σ ε∂t ψc + Ω · ∇x ψc + ψc = ψc + ψu , (3.11b) ε ε ε where, as before, a bar denotes the angular average of S2 . The idea of the hybrid is to solve (3.11) using a high-resolution angular discretization for ψu and a low-resolution angular discretization for ψc over a time step ∆t = T /M , where M ∈ N>0 and then perform a reconstruction to reinitialize ψu and ψc for the next time step. To formalize this procedure, define a set of temporal grid points 0 = t0 < t1 < ... < tM = T , and for m ∈ {1, 2, . . . , M }, let f (t− m ) = limδ→0+ f (tm − δ) for any function f of t that is continuous on [tm−1 , tm ). Then for m ∈ {1, 2, . . . , M }, ψu,m and ψc,m satisfy the following system of equations over the interval [tm−1 , tm ) σ ε∂t ψu,m + Ω · ∇x ψu,m + ψu,m = εq, (3.12a) ε σ σ ε∂t ψc,m + Ω · ∇x ψc,m + ψc,m = (ψu,m + ψc,m ), (3.12b)  ε ε  g,  m = 1, ψu,m |t=tm−1 = (3.12c) ψu,m−1 (t− −  m−1 ) + ψc,m−1 (tm−1 ) m > 1.  ψc,m |t=tm−1 = 0. (3.12d) 43 The intuition behind this splitting is that (3.12a) can be discretized with a high-resolution angular discretization but solved more efficiently than (3.3) since angular unknowns are no longer coupled. Although (3.12b) features the same type of angular coupling as (3.3), it can be solved with fewer degrees of freedom because the source ψu,m is, in general, more regular than q. However, because ψu,m decays exponentially whenever σ > 0, the hybrid is only solved for a time step ∆t before the relabeling in (3.12d)- (3.12c) is implemented. The enables the hybrid to capture more high-resolution features than (3.12b) can do alone. 3.1.4 Angular discretization of the hybrid We focus now on the angular discretization of (3.12). The strategy of the hybrid is to discretize (3.12a) in angle with a high-resolution method and (3.12b) in angle with a low-resolution method. In practice, there are a variety of strategies and combinations available to do so. For the purposes of analysis, we assume that (3.12a) is solved exactly and that (3.12b) is discretized with a PN method. That is, we seek ψ N = ψu,m N N + ψc,m where for each m ∈ {1, 2, . . . M } N N ∈ C([tm−1 , tm ); X × L2 (S2 )) × C([tm−1 , tm ); X × PN (S2 ))  ψu,m , ψc,m (3.13) satisfies N N σ N ε∂t ψu,m + Ω · ∇x ψu,m + ψ = εq, (3.14a) ε u,m N N  σ N σ N ε∂t ψc,m + PN Ω · ∇x ψc,m + ψc,m = (ψu,m N ), + ψc,m (3.14b) ε ε (3.14c)   g,  m = 1, N N ψc,m t=tm−1 = 0, ψu,m t=tm−1 = (3.14d) N (t− N −  ψu,m−1  m−1 ) + ψc,m−1 (tm−1 ) m > 1. We compare the accuracy of the solution defined in (3.14) with the monolithic PN method (3.7a), using the same value of N . To simplify the numerical analysis of these two models, we keep the time and space variables continuous. 44 Since (3.7a) and (3.14b) have the same computational complexity, the goal is to asses the additional benefit of solving (3.14a). Clearly the additional cost of solving (3.14a) involves both memory and run-time; both are fairly easy to quantify. However, assessing the gains in accuracy is not as simple. Thus is important to understand these gains in order to better quantify observed improvements in run-time efficiency provided by the hybrid. 3.1.5 Preview of main results. Let s ≥ 1 be the number of angular L2 derivatives in the solution ψ and let N be an integer such that N ≥ s − 1. Let eN = ψ − ψ N be the error in the PN approximation N N (3.7a), and let eN M = ψ − (ψu,M + ψc,M ) be the error in the hybrid at the M -th time step. The main results of the chapter are the PN error estimate in Theorem 3.2.1 and the hybrid error estimate in Theorem 3.3.2. We compare these errors for two different regimes: first, when σ ≍ 12 and ε → 0 (the diffusion regime) and second, when ε ≍ (1) and σ → 0 (the purely absorbing regime).3 When ε ≪ 1 and σ ≍ 1, Theorem 3.2.1 implies that " s−1 # i s−1 T 2 X T ε ε ∥eN ∥L2 (X×S2 ) (T ) ≲ e−σT /ε + s +O . (3.15) (N + 1)s i=0 εi+1 σ σ Thus with sufficient regularity, the PN approximation is spectrally accurate and grows linearly in time for large T . The first and third term in brackets depend on whether or not q and g are isotropic (independent of angle). The first term is due to initial layers when g is non-isotropic, and the third terms arises when q is non-isotropic. When g and q are isotropic, the PN error reduces to (see Corollary 3.2.2) εs−1 T ∥eN ∥L2 (X×S2 ) (T ) ≲ . (3.16) σ s (N + 1)s For the hybrid method the initial condition and source are always isotropic. Thus 2 Recall that a ≍ b if and only if a = O(b) and b = O(a) 3 Technically, this is the streaming regime for Ψ, since there is no absorption. 45 Theorem 3.3.2 gives the following compact bound εs−1 T ∥eNM ∥L2 (X×S2 ) (T ) ≲ . (3.17) σ s (N + 1)s Thus the hybrid method comes equipped with a better error estimate than the PN method, and the estimates agree when the data is isotropic. This suggests that the hybrid method is at least as accurate as the PN approximation when ε ≪ 1 and σ ≍ 1. In addition, the hybrid estimate is independent of the time step ∆t used for re-initialization in this regime. For both the hybrid method and PN approximation, the errors converge to zero as ε → 0 whenever s > 1 (modulo the initial layer in (3.15)). This fact is consistent with the fact that the PN method recovers the diffusion limit (3.2a) whenever N ≥ 1 [54]. For the second regime of interest, σ ≪ 1 and ε ≍ 1, Theorem 3.2.1 gives the PN error estimate ps (T /ε) ∥eN ∥L2 (X×S2 ) (T ) ≲ , (3.18) (N + 1)s Ps+1 where ps (ω) = i=0 ci (T )ω i is a polynomial of degree s + 1 with non-negative coefficients ci (T ) = ai T + bi , ai , bi ≥ 0. Here the hybrid estimate provides a significant improvement: ∆ts T ∆tσ ∥eN M ∥L2 (X×S2 ) (T ) ≲ s+1 s min(1, 2 ) (3.19) ε (N + 1) ε In particular, ∥eN M ∥L2 (X×S2 ) (T ) = 0, when σ = 0. This result is expected since in that case the uncollided solution and the transport solution agree. As expected, the error is mono- tonic in ∆t; however, small time steps require more evaluations of the uncollided equation and more reintializions. In practice, this additional cost must be taken into account. 3.2 PN Analysis 3.2.1 Spherical harmonics preliminaries A natural space to analyze the transport equation and the PN approximation is the Sobolev space H◦s (S2 ). To describe this space, we recall some elementary facts about spherical harmonics which can be found, for example, in [55, 56]. For u, v ∈ L2 (S2 ), let X∞ X∞ (u, v)L2 (S2 ) = u⊤ℓ vℓ and ∥u∥2L2 (S2 ) = ∥uℓ ∥2 , (3.20) ℓ=0 ℓ=0 46 where Xℓ Xℓ u⊤ ℓ vℓ = uℓ,k vℓ,k and ∥uℓ ∥2 = u⊤ ℓ uℓ = |uℓ,k |2 . (3.21) k=−ℓ k=−ℓ A standard way to define Sobolev spaces on the sphere is via the Laplace-Beltrami op- erator ∆◦ , which is the spherical component of the Laplacian and for which the spherical harmonics are eigenfunctions: −∆◦ Yℓ,j = ℓ(ℓ + 1)Yℓ,j . For even integers s, the usual norm is !1/2  s/2 ∞  2s 1 X 2 1 ∥u∥H◦s (S2 ) := + ∆◦ u = bℓ ∥uℓ ∥ , bℓ = +ℓ . (3.22) 4 ℓ=0 2 The definition of this inner product extends naturally to all s ∈ R, and the space H◦s (S2 ) is then the completion of smooth functions under the H◦s (S2 ) norm [55]. Rather than working directly with the H◦s (S2 ) norm in (3.22), it is convenient in the analysis below to use an equivalent norm. For s ≥ 0, define the H s (S2 ) semi-norm and norm by ∞ !1/2 X  1/2 |u|H s (S2 ) := bℓ ∥uℓ ∥ 2 and ∥u∥H s (S2 ) = s∥u∥2L2 (S2 ) + |u|2H s (S2 ) , (3.23) ℓ=s respectively, where the sum in the semi-norm definition in (3.23) begins at s for technical arguments that are used in the proof of Lemma (3.2.5) below. When s = 0, the norms coincide: ∥u∥H 0 (S2 ) = ∥u∥H◦0 (S2 ) = ∥u∥L2 (S2 ) . More generally, the following equivalence holds. Lemma 3.2.1 (Norm equivalence). For any s ≥ 0, c1 (s)∥u∥H s (S2 ) ≤ ∥u∥H◦s (S2 ) ≤ c2 (s)∥u∥H s (S2 ) . (3.24) where     1  if s = 0, 1  if s = 0, c1 (s) = and c2 (s) = q . (3.25)  √1  5 s− 1 s     3s if s ≥ 1  s 2 if s ≥ 1 47 Proof. When s = 0, the norms are equal, so (3.24) holds trivially. Thus assume that s ≥ 1. The first inequality in (3.24) follows from the fact that ∞ ∞  2s ∞  2s 2 X 2 X 1 2 X 1 1 ∥u∥H s (S2 ) = s ∥uℓ ∥ + +ℓ ∥uℓ ∥ ≤ 3s +ℓ ∥uℓ ∥2 = 2 ∥u∥2H◦s (S2 ) . ℓ=0 ℓ=s 2 ℓ=0 2 [c1 (s)] (3.26) To prove the second inequality in (3.24), we use the elementary inequality  2s s 1 ≤ s− , s ≥ 1, (3.27) 4 2 to conclude that s−1  2s  2s X 1 1 1 ∥u∥2H◦s (S2 ) = +ℓ 2 2 ∥uℓ ∥ + |u|H s (S2 ) ≤ s− s∥u∥2L2 (S2 ) + |u|2H s (S2 ) ℓ=0 2 s 2  2s  2s 1 1 2 4 1 (3.28) ≤ s− s∥u∥L2 (S2 ) + s− |u|2H s (S2 ) s 2 s 2  2s  5 1  ≤ s− s∥u∥2L2 (S2 ) + |u|2H s (S2 ) = [c2 (s)]2 ∥u∥2H s (S2 ) . s 2 Lemma 3.2.2 (Approximation property). For s ≥ 0 and N ≥ max{0, s − 1}, 1 1 ∥(I − PN )u∥L2 (S2 ) ≤ s |(I − PN )u|H s (S2 ) ≤ |u|H s (S2 ) . (3.29) (N + 1) (N + 1)s Proof. From the definition of the projection PN in (3.5), ∞ ∞  2s ∞  2s X 2 1 X 1 2 1 X 1 ∥uℓ ∥ ≤ 2s +ℓ ∥uℓ ∥ ≤ 2s +ℓ ∥uℓ ∥2 ℓ=N +1 (N + 1) ℓ=N +1 2 (N + 1) ℓ=s 2 | {z } | {z } | {z } =∥(I−PN )u∥2 2 =|(I−PN )u|2 =|u|2 L (S2 ) H s (S2 ) H s (S2 ) (3.30) Taking square roots of each term above yields the desired result. For vector-valued functions of space, we define the usual L2 (X) inner-product and norm by Z Z ⊤ (v, w)L2 (X) = v(x) w(x) dx and ∥v∥2L2 (X) = (v, v)L2 (X) = ∥v(x)∥2 dx. (3.31) X X 48 The space H 0,s (X × S2 ) = L2 (X; H s (S2 )) is the space of measurable functions u : X × S2 → R with the semi-inner product Z X ∞  2s ∞  2s 1 ⊤ X 1 (u, v)H 0,s (X×S2 ) = +ℓ uℓ (x) vℓ (x) dx = +ℓ (uℓ , vℓ )L2 (X) (3.32) X ℓ=s 2 ℓ=s 2 such that the semi-norm Z X ∞  2s ∞  2s 1 X 1 |u|2H 0,s (X×S2 ) = +ℓ |uℓ (x)| dx =2 +ℓ ∥uℓ ∥2L2 (X) (3.33) X ℓ=s 2 ℓ=s 2 is bounded. For r ∈ N≥0 and s ≥ 0 we define H r,s (X ×S2 ) to be the space of functions u : X ×S2 → R such that the semi-norms. X 3 |u|H ς,s (X×S2 ) := |∂xi1 xi2 ...xiς u|H 0,s (X×S2 ) . (3.34) i1 ,i2 ,...,iς =1 are bounded for all positive integers ς ≤ r. These semi-norms above are equivalent to the standard semi-norms, but are more convenient in the context of the RTE since X 3 |Ω · ∇x u|H r,s (X×S2 ) ≤ |∂xi u|H r,s (X×S2 ) = |u|H r+1,s (X×S2 ) , (3.35) i=1 which will be used in the analysis below. Henceforth, the domain of integration for H s and H r,s will be left off when there is no ambiguity, i.e., |u|H s := |u|H s (S2 ) and |u|H r,s := |u|H r,s (X×S2 ) (3.36) Finally, for p ≥ 1 and measurable functions u : X × S2 × [α, β] → R, we denote the space- time semi-norms by Z β 1/p |u|Lp ([α,β];H r,s ) = |u|pH r,s dτ and |u|L∞ ([α,β]);H r,s ) = ess sup |u|H r,s . (3.37) α t∈[α,β] 3.2.2 Stability of the PN system In this section, we derive estimates on high-order semi-norms that arise in the sub- sequent error analysis. The analysis requires iterated inequalities of Grönwall type (see 49 Lemma .0.1 in the Appendix), so we define the following notation to simplify integrals that arise. Let 0 ≤ α ≤ t ≤ β < ∞, and for any function f ∈ L1 ([α, β]), define the bounded linear operator Aα : L1 ([α, β]) → C 0 ([α, β]) and its powers Akα , k ∈ N≥0 , by   1 Z t I,  k = 0, −σ(t−τ )/ε2 k Aα [f ](t) := e f (τ )dτ and Aα [f ](t) = (3.38) ε α Aα [Ak−1   α [f ]](t), k ≥ 1. In addition, let 1 ∈ L1 ([α, β]) be the function that is identically one, and let 2 Fα (t) = e−σ(t−α)/ε , t ∈ [α, β]. (3.39) It is clear from the definition in (3.38) Aα is monotonic; that is, if 0 ≤ f (t) ≤ g(t) for a.e. t ∈ [α, β], then Aα [f ](t) ≤ Aα [g](t) for all t ∈ [α, β]. Lemma 3.2.3. Let Fα be given as in (3.39). Then for all t ≥ α and every k ∈ N≥0 , k ! 2 εk 1 t − α (t − α)k e−σ(t−α)/ε  Aα [1](t) ≤ min k , k and Aα [Fα ](t) = (3.40) σ k k! ε k!εk 2 Proof. We first prove the bound in (3.40). Since e−σ(t−τ )/ε ≤ 1, a direct calculation gives k 1 t τk−1 Z τ1  1 t−α Z Z Aα [1](t) ≤ k k ··· 1dτ0 · · · dτk−1 = . (3.41) ε α α α k! ε On the other hand, it follows directly from the definition of Aα that ε ε εk Aα [1](t) = (1 − Fα (t)) ≤ =⇒ Akα [1](t) ≤ (3.42) σ σ σk Together (3.41) and (3.42) yield (3.40). We prove the second statement in (3.40) by induction on k. When k = 0, the statement follows trivially, 2 (t − α)0 e−σ(t−α)/ε A0α [Fα ](t) = Fα (t) = , (3.43) 0!e0 Now let us assume the statement is true for arbitrary k, and 2 1 t −σ(t−τ )/ε2 (τ − α)k e−σ(τ −α)/ε Z k+1  k  Aα [Fα ](t) = Aα Aα [Fα ] (t) = e dτ, (3.44) ε α k!εk direct integration leads to, 2 (t − α)k+1 e−σ(t−α)/ε Ak+1 α [Fα ](t) = (3.45) (k + 1)!εk+1 50 3.2.2.1 Estimates for the PN system We derive evolution equations for the H r,s semi-norms, defined in (3.34). For r = 0, we test each equation in (3.8) by bℓ ψN ℓ , integrate by parts, and then sum over ℓ ∈ {s, s + 1, ..., N }. This gives ε σ ∂t |ψ N |2H 0,s + |ψ N |2H 0,s = ε(PN q, ψ N )H 0,s 2 ε 3 X N  2s  2s !   X 1 1 N  T (i) N − ℓ+ − γs,ℓ ℓ − ψ ℓ , aℓ ∂xi ψℓ−1 , i=1 ℓ=s 2 2 L2 (X) (3.46) where γs,ℓ = (1 − δs,ℓ ) is used to handle the first non-zero term in the sum over ℓ. For the special case s = 0, (3.46) recovers the usual L2 energy equation: εd N 2 σ ∥ψ ∥L2 (X×S2 ) + ∥ψ N − ψ N ∥2L2 (X×S2 ) = ε(PN q, ψ N )2L2 (X×S2 ) . (3.47) 2 dt ε To find a closed estimate with respect to the H 0,s semi-norms, we focus on the summation in (3.46). Lemma 3.2.4. Let s ≥ 1 and ℓ ≥ s. Then  2s  2s  s  s−1 1 1 1 1 ℓ+ − γs,ℓ ℓ − ≤ 2es ℓ + ℓ− . (3.48) 2 2 2 2 Proof. We first establish an elementary inequality. Since ℓ ≥ s,         1 1 1 1 1 1 1/(s−1/2) ℓ+ = ℓ− 1+ ≤ ℓ− 1+ ≤ ℓ− e (3.49) 2 2 ℓ − 21 2 s − 12 2 Therefore  s−1  s−1  s−1 1 s−1 1 1 ℓ+ ≤e s−1/2 ℓ− ≤e ℓ− . (3.50) 2 2 2 We use (3.50) to show (3.48). There are two cases: Case 1 (ℓ = s): In this case, γs,ℓ = 0. Since 2s > s + 1/2 and by (3.50),  2s   s  s−1  s  s−1 1 1 1 1 1 1 ℓ+ = s+ ℓ+ ℓ+ ≤ 2es ℓ + ℓ− (3.51) 2 2 2 2 2 2 51 Case 2 (ℓ > s): In this case, γs,ℓ = 1. Applying a binomial expansion and then (3.50) gives  2s  2s X 2s    k  2s 2s−1 X 2s  k 1 1 2s 1 1 1 ℓ+ − ℓ− = ℓ− − ℓ− = ℓ− 2 2 k=0 k 2 2 k=0 k 2 2s−1 X 2s 2s − 1  k 2s−1 X 2s − 1  k 1 1 = ℓ− ≤ 2s ℓ− k=0 2s − k k 2 k=0 k 2  2s−1  s  s−1  s  s−1 1 1 1 1 1 = 2s ℓ + = 2s ℓ + ℓ+ ≤ 2es ℓ + ℓ− 2 2 2 2 2 (3.52) Lemma 3.2.5 (Semi-norm recurrence). Let s ≥ 1, q ∈ L1 ([0, T ]; H r,s ) and g ∈ H r,s . Then for all t ∈ [0, T ], |ψ N |H r,s (t) ≤ CsA0 [|ψ N |H r+1,s−1 ](t) + |PN g|H r,s F0 (t) + εA0 [|PN q|H r,s ](t), (3.53) where F0 is defined in (3.39) and C is a constant independent of the data. Proof. We assume first that r = 0 and focus on the last term in (3.46). It follows from (i) (i) the induced norm bound ∥aℓ ∥2 ≤ 4 [44], (ii) the bounds in Lemma 3.2.4, (iii) the Cauchy- Schwarz inequality, and (iv) the H r,s semi-norm definitions in (3.33) and (3.34) that N  2s  2s !   X 1 1 N  T (i) N ℓ+ − γs,ℓ ℓ − ψ ℓ , aℓ ∂xi ψℓ−1 ℓ=s 2 2 L2 (X) N   s   s−1 X 1 1 (i) ≤ 2es ℓ+ ℓ− ∥aℓ ∥2 ∥ψN N ℓ ∥L2 (X) ∥∂xi ψℓ−1 ∥L2 (X) 2 2 ℓ=s (3.54) N ! 1/2 N ! 1/2   2s  2(s−1) X 1 X 1 ≤ Cs ℓ+ |ψN | 2 ℓ L2 (X) ℓ − |∂xi ψN 2 ℓ−1 |L2 (X) ℓ=s 2 ℓ=s 2 ≤ Cs|ψ N |H 0,s |∂xi ψ N |H 0,s−1 . Applying the bound above to the right-hand side of (3.46) and applying Lemma .0.1 gives (3.53). The case r > 1 can be handled by differentiating the PN equations in space and then repeating the arguments above. 52 For a general function ϕ ∈ Lp ([α, β]; H r,s ), let |ϕ|Lp ([α,•];H 0,r+s ) : [α, β] → R be the map defined |ϕ|Lp ([α,•];H r,s ) (τ ) = |ϕ|Lp ([α,τ ];H r,s ) . Lemma 3.2.6 (Stability of higher-order semi-norms). Let q ∈ L1 ([0, T ]; H r,0 ) and g ∈ H r,0 . If t ∈ [0, T ], then |ψ N |H r,0 (t) ≤ |PN g|H r,0 + |PN q|L1 ([0,t];H r,0 ) . (3.55) If, in addition, s ≥ 1, q ∈ L1 ([0, T ]; H i,j ) and g ∈ H i,j for each i, j such that 0 ≤ i ≤ r, 0 ≤ j ≤ s, and i + j = r + s, |ψ N |H r,s (t) ≤ Cs s!As0 |PN g|H r+s,0 + |PN q|L1 ([0,•];H r+s,0 ) (t)   s−1 s−1 X s! X s! + Cs Ai0 [F0 |PN g|H r+i,s−i ](t) + Cs ε Ai+1 0 [|PN q|H r+i,s−i ](t), i=0 (s − i)! i=0 (s − i)! (3.56) where F0 is given in (3.39) where Cs is a constant depending only on s. Proof. First we will prove (3.55) for s = 0, in which case H r,0 = L2 (X × S2 ). From (3.47) and the Cauchy-Schwarz inequality, εd N 2 ∥ψ ∥L2 (X×S2 ) ≤ ε(PN q, ψ N )L2 (X×S2 ) ≤ ε∥PN q∥L2 (X×S2 ) ∥ψ N ∥L2 (X×S2 ) (3.57) 2 dt an application of Lemma .0.1, gives ∥ψ N ∥L2 (X×S2 ) (t) ≤ ∥PN g∥L2 (X×S2 ) + ∥PN q∥L1 ([0,t];L2 (X×S2 )) . (3.58) For r > 0, ϕN = ∂xi1 xi2 ...xir ψ N satisfies (3.7a), with initial condition given by the following ϕN t=0 = PN (∂xi1 xi2 ...xir g) = ∂xi1 xi2 ...xir PN g and source PN (∂xi1 xi2 ...xir q) = ∂xi1 xi2 ...xir PN q. Repeating the argument above gives, in analogy with (3.58), ∥∂xi1 xi2 ...xir ψ N ∥L2 (X×S2 ) (t) ≤ ∥∂xi1 xi2 ...xir PN g∥L2 (X×S2 ) + ∥∂xi1 xi2 ...xir PN q∥L1 ([0,t];L2 (X×S2 )) , (3.59) 53 Summing (3.59) over each permutation of the r derivatives ∂xi1 ∂xi2 · · · ∂xir and using the definitions in (3.34) recovers (3.55). We next prove (3.56). Let br,s (t) = |PN g|H r,s F0 (t) + εA0 [|PN q|H r,s ](t) and cr,s (t) = |ψ N |H r,s (t), (3.60) where F0 is defined in (3.39). Then Lemma 3.2.5 gives the following recursion relation: cr,s (t) ≤ CsA0 [cr+1,s−1 ](t) + br,s (t), s ≥ 1. (3.61) Unrolling this recursion in s gives s−1 s X s! cr,s (t) ≤ s!C As0 [cr+s,0 ](t) + C i Ai0 [br+i,s−i ](t), (3.62) i=0 (s − i)! and translating back to the semi-norms with (3.60) gives |ψ N |H r,s (t) ≤ s!C s As0 [|ψ N |H r+s,0 ](t) s−1 s−1 X s! X s! + C i |PN g|H r+i,s−i Ai0 [F0 ](t) + ε C i+1 Ai+1 0 [|PN q|H r+i,s−i ](t). i=0 (s − i)! i=0 (s − i)! (3.63) We then apply (3.55) (with r replaced with r + s) to the first term on the right hand side of (3.63). 3.2.2.2 Estimates for continuous system In this section we extend the stability results for |ψ N |H r,s to |ψ|H r,s , using the infinite moment hierarchy in (3.9). These estimates will be useful in deriving consistency esti- mates. For the case r = 0, as in the previous section, we test each equation in (3.9) by bℓ ψℓ , integrate by parts, and sum over ℓ ≥ s. The result is analogous to (3.46), namely ε σ ∂t |ψ|2H 0,s + |ψ|2H 0,s = ε(q, ψ)H 0,s 2 ε 3 X ∞  2s  2s !   (3.64) X 1 1  T (i) − ℓ+ − γs,ℓ ℓ − ψ ℓ , aℓ ∂xi ψℓ−1 , i=1 ℓ=s 2 2 L2 (X) 54 As before, when s = 0, (3.64) recovers the usual L2 stability result: εd σ ∥ψ∥2L2 (X×S2 ) + ∥ψ − ψ∥2L2 (X×S2 ) = ε(q, ψ)2L2 (X×S2 ) (3.65) 2 dt ε The following results are continuous analogues to Lemmas 3.2.5 and 3.2.6. Their proofs are nearly identical, so we only give a brief summary. Lemma 3.2.7 (Semi-norm recurrence for the continuous system). Let s ≥ 1, q ∈ L1 ([0, t]; H r,s ) and g ∈ H r,s . Then for all t ∈ [0, T ], |ψ|H r,s (t) ≤ CsA0 [|ψ|H r+1,s−1 ](t) + |g|H r,s F0 (t) + εA0 [|q|H r,s ](t), (3.66) where C is a constant independent of the data. Summary of Proof. The proof follows the same lines and with same constants as in Lemma 3.2.5 after changing ψN ℓ by ψℓ and then taking all the sums to infinity. To generalize the proof for r > 1, we differentiate the system (3.9) in space and repeat the process. Lemma 3.2.8 (Stability of higher order semi-norms for the continuous system). If q ∈ L1 ([0, T ]; H r,0 ) and g ∈ H r,0 , then for all t ∈ [0, T ], |ψ|H r,0 (t) ≤ |g|H r,0 + |q|L1 ([0,t];H r,0 ) . (3.67) If, in addition, q ∈ L1 ([0, T ]; H i,j ) and g ∈ H i,j for each i, j such that 0 ≤ i ≤ r, 0 ≤ j ≤ s, and i + j = r + s, |ψ|H r,s (t) ≤ Cs s!As0 |g|H r+s,0 + |q|L1 ([0,•];H r+s,0 ) (t)   s−1 s−1 (3.68) X s! X s! + Cs Ai0 [F0 |g|H r+i,s−i ](t) + Cs ε Ai+1 [|q|H r+i,s−i ](t), i=0 (s − i)! i=0 (s − i)! 0 where Cs is a constant depending only on s. Corollary 3.2.1 (Isotropic data and zero initial condition for the continuous systems). Let s ≥ 1. In the special case that g = 0 and q is isotropic, then if t ∈ [0, T ], |ψ|H r,s (t) ≤ Cs s!As0 [|q|L1 ([0,•];H r+s,0 ) ](t). (3.69) 55 Proof of Lemma 3.2.8. With the obvious changes, the proof is the same line by line as proof in Lemma 3.2.6, with some key steps replace by their continuous counterparts, namely, in the initial step we use (3.65) instead of (3.47) and we invoke Lemma 3.2.7 instead of Lemma 3.2.5. 3.2.3 PN error analysis In this section we will analyze the error produced by the solution of (3.3) when the PN approximation is used. Definition 3.2.1 (PN Error). The PN error is eN (t) = ψ(t) − ψ N (t) = η N (t) + ξ N (t), (3.70) where η N = ψ − PN ψ is the consistency error and ξ N = PN ψ − ψ N is the stability error. Lemma 3.2.9. For all t ∈ [0, T ], ξ N is controlled by η N via the following estimate: 1 t Z N ∥ξ ∥L2 (X×S2 ) (t) ≤ ∥PN (Ω · ∇x η N ∥L2 (X×S2 ) (τ ) dτ. (3.71) ε 0 Proof. Applying the projection PN to (3.3) and subtracting (3.7a) yields a PN equation for ξ N with a source that depends on η N : σ N σ ε∂t ξ N + PN (Ω · ∇x ξ N ) + ξ = ξ N − PN (Ω · ∇x η N ), ξN t=0 = 0. (3.72) ε ε Thus (3.72) follows immediately from the bound (3.58), replacing g by zero and q by ε−1 PN (Ω · ∇x η N ). Lemma 3.2.10. Let t ∈ [α, β) ⊆ [0, T ], then σ(t−α) ∥η N ∥L2 (X×S2 ) (t) ≤ e− ε2 ∥η N ∥L2 (X×S2 ) (α) + εAα [∥P eN q∥L2 (X×S2 ) ](t) (3.73) +Aα [∥P eN (Ω · ∇x ψ)∥L2 (X×S2 ) ](t) Proof. Applying the projection PN to (3.3), and subtracting it from (3.3), we see that η N satisfies 1 σ eN q − 1 [Ω · ∇x ψ − PN (Ω · ∇x PN ψ)] ∂t η N − PN (Ω · ∇x η N ) + 2 η N = P (3.74) ε ε ε 56 Since for any ϕ ∈ L2 (S2 ), (PN ϕ, η N )L2 (S2 ) = 0, testing the equation above against η N gives 1d N 2 σ 1 ∥η ∥L2 (X×S2 ) + 2 ∥η N ∥2L2 (X×S2 ) = (P eN q, η N )L2 (X×S2 ) − (Ω · ∇x ψ, η N )L2 (X×S2 ) 2 dt ε ε 1 e = (P eN q, η N )L2 (X×S2 ) − (PN (Ω · ∇x ψ), η N )L2 (X×S2 ) ε ≤ (∥P eN q∥L2 (X×S2 ) + 1 ∥P eN (Ω · ∇x ψ)∥L2 (X×S2 ) )∥η N ∥L2 (X×S2 ) ε (3.75) The conclusion then follows from Lemma .0.1. An immediate corollary of Lemma 3.2.10 is the following: Lemma 3.2.11. Let s ≥ 0 and N ≥ max{0, s − 1}. If g ∈ H 0,s and q ∈ L∞ ([0, T ]; H 0,s ). Then we have " # 1 T e−σT /ε |g|H 0,s + ε|q|L∞ ([0,T ];H 0,s ) A0 [1](T ) + 2 ∥η N ∥L2 (X×S2 ) (T ) ≤ s sup |ψ|H 1,s (τ ) . (N + 1) ε τ ∈[0,T ] (3.76) Proof. We apply Lemma 3.2.10 with α = 0. In this case ∥η N ∥L2 (X×S2 ) (α) = ∥P eN g∥. Mean- while, by Lemma 3.2.2, 1 1 ∥P eN g∥L2 (X×S2 ) ≤ |g|H 0,s and ∥PeN q∥L2 (X×S2 ) ≤ |q|H 0,s . (3.77) (N + 1)s (N + 1)s Thus since Aα [f ](t) ≤ ε−1 (t − α) supτ ∈[α,t] f (τ ) (cf. (3.40) with k = 1), another application of Lemma 3.2.2 gives A0 [∥PeN (Ω · ∇x ψ)∥L2 (X×S2 ) ](T ) ≤ T sup ∥P eN (Ω · ∇x ψ)∥L2 (τ ) ε τ ∈[0,T ] (3.78) T T ≤ s sup |Ω · ∇x ψ|H 0,s (τ ) ≤ sup |ψ|H 1,s (τ ) ε(N + 1) τ ∈[0,T ] ε(N + 1)s τ ∈[0,T ] Plugging the preceding bound into Lemma 3.2.10 gives the result. Lemma 3.2.12 (a priori estimate). Let s ≥ 1 and N ≥ s − 1. Let q ∈ L∞ ([0, T ]; H i,j ) and g ∈ H i,j for each i, j such that 0 ≤ i ≤ r, 0 ≤ j ≤ s, and i + j = r + s. Then for all t ∈ [0, T ], 57 (  s  εs s!    t |ψ|H r,s (t) ≤ Cs |g|H r+s, + t|q|L∞ ([0,t);H r+s,0 ) min , σs ε s−1   i −σt/ε2 X s t +e |g|H r+i,s−i (3.79) i=0 i εi s−1  i+1 ! ) X s! εi+1 1 t +ε |q|L∞ ([0,t];H r+i,s−i ) min i+1 , . i=0 (s − i)! σ (i + 1)! ε Proof. Recall the stability estimate from Lemma 3.2.8: |ψ|H r,s (t) ≤ Cs s!As0 |g|H r+s,0 + |q|L1 ([0,•];H r+s,0 ) (t)   s−1 s−1 (3.80) X s! X s! + Cs Ai0 [F0 |g|H r+i,s−i ](t) + Cs ε Ai+1 [|q|H r+i,s−i ](t). i=0 (s − i)! i=0 (s − i)! 0 Substituting the bounds for A0 [1] and the formula for A0 [Fα ] from Lemma 3.2.3 into the above estimate yields the stated result. Theorem 3.2.1 (PN error). Let s ≥ 1 and N ≥ s − 1. Let q ∈ L∞ ([0, T ]; H i,j ) and g ∈ H i,j for each i, j such that 0 ≤ j ≤ s, i + j ≤ s + 1. Then 2 e−σT /ε  2  N 1 ε ∥e ∥L2 (X×S2 ) (T ) ≤ s |g|H 0,s + s |q|L∞ ([0,T ];H 0,s ) min ,T (N + 1) (N + 1) σ (  s+1 ! 2Cs   εs−1 s!T T + s |g|H s+1,0 + T |q|L∞ ([0,T ];H s+1,0 ) min s , (N + 1) σ ε s−1   i+1 −σT /ε2 X s T +e |g|H 1+i,s−i i=0 i εi+1 s−1  i+1 ) X s! ε T 1 T i+2 + |q|L∞ ([0,T ];H 1+i,s−i ) min , i=0 (s − i)! σ i+1 (i + 1)! εi+1 (3.81) 58 Proof. By the triangle inequality, Lemma 3.2.9, ∥eN ∥L2 (X×S2 ) (T ) ≤ ∥η N ∥L2 (X×S2 ) (T ) + ∥ξ N ∥L2 (X×S2 ) (T ) 1 T Z N ≤ ∥η ∥L2 (X×S2 ) (T ) + ∥PN (Ω · ∇x η N )∥L2 (X×S2 ) (τ ) dτ ε 0 1 T Z N ≤ ∥η ∥L2 (X×S2 ) (T ) + ∥∇x η N ∥L2 (X×S2 ) (τ ) dτ ε 0 1 T (3.82) ≤ ∥η N ∥L2 (X×S2 ) (T ) + sup |ψ|H 1,s (τ ) (N + 1)s ε τ ∈[0,T ] 2 e−σT /ε ε ≤ s |g|H 0,s + |q|L∞ ([0,T ];H 0,s ) A0 [1](T ) (N + 1) (N + 1)s 2T + sup |ψ|H 1,s (τ ) ε(N + 1)s τ ∈[0,T ] In the last two lines, we applied spectral estimate in Lemma 3.2.2. Inthe last line we used Lemma 3.2.11. Applying Lemma 3.2.12 with r = 1 yields the result. Corollary 3.2.2 (PN error for Isotropic data). Let s ≥ 1 and N ≥ s−1. Let q ∈ L∞ ([0, T ]; H s+1,0 ) and g ∈ H s+1,0 . If g and q are isotropic, then  s+1 ! N 2Cs   εs−1 s!T T ∥e ∥L2 (X×S2 ) (T ) ≤ |g|H s+1,0 + T |q|L∞ ([0,T ];H s+1,0 ) min , (N + 1)s σs ε (3.83) 3.3 Hybrid error analysis 3.3.1 A priori estimates of the uncollided component Since our goal is to derive error estimates which only depend on data, and ψu,m ap- pears as a source in the collided equation, we require the following a priori estimates on ψu,m to bound |ψc,m |H 1,s in the proof of Theorem 3.3.2. Lemma 3.3.1 (Stability of the uncollided component). Let 1 ≤ m ≤ M , q ∈ L∞ ([0, T ]; H r,0 ), and g ∈ H r,0 for some r ≥ 0. Then for all t ∈ [tm−1 , tm ), 2 2 |ψu,m |H r,0 (t) ≤ e−σ(t−tm−1 )/ε |g|H r,0 + e−σ(t−tm−1 )/ε |q|L1 ([0,tm−1 ];H r,0 ) + εAtm−1 [|q|H r,0 ](t), (3.84a) σ |ψu,m |H r,0 (tm ) + |ψu,m |L1 ([tm−1 ,tm ];H r,0 ) ≤ |g|H r,0 + |q|L1 ([0,tm ];H r,0 ) , (3.84b) ε2 |ψu,m |L1 ([tm−1 ,t];H r,0 ) (t) ≤ ε(|g|H r,0 + t|q|L∞ ([0,t];H r,0 ) )Atm−1 [1](t). (3.84c) 59 Proof. We prove the result only for r = 0 since the other cases are obtained by applying the same techniques to (3.12a) differentiated r times in space. Testing (3.12a) with ψu,m and applying Cauchy-Schwarz inequality, we obtain the following differential inequality 1d σ ∥ψu,m ∥2 + 2 ∥ψu,m ∥2 ≤ ∥q∥∥ψu,m ∥. (3.85) 2 dt ε An application of (.4) in Lemma .0.1, and the fact that ψu,m (tm−1 ) = ψ(tm−1 ) gives 2 ∥ψu,m ∥L2 (X×S2 ) (t) ≤ e−σ(t−tm−1 )/ε ∥ψ∥L2 (X×S2 ) (tm−1 ) + εAtm−1 [∥q∥L2 (X×S2 ) ](t). (3.86) Using Lemma 3.2.8 on ∥ψ∥L2 (X×S2 ) (tm−1 ) gives the first result (3.84a). An application of (.3) in Lemma .0.1 over [tm−1 , tm ) gives σ ∥ψu,m ∥L2 (X×S2 ) (tm )+ ∥ψu,m ∥L1 ([tm−1 ,tm ],L2 (X×S2 )) ≤ ∥q∥L1 ([tm−1 ,tm ],L2 (X×S2 )) +∥ψ∥L2 (X×S2 ) (tm−1 ), ε2 (3.87) and then another application of Lemma 3.2.8 to bound ∥ψ∥L2 (X×S2 ) (tm−1 ) gives (3.84b). Estimate (3.84c) is obtained from integrating the first estimate (3.84a) from tm−1 to t. 3.3.2 Hybrid error analysis In this section we will analyze the error in the hybrid method using the formulation in (3.14). Definition 3.3.1 (Hybrid errors). Let 1 ≤ m ≤ M and t ∈ [tm−1 , tm ). The m-th hybrid error is N N eN m (t) = eu,m (t) + ec,m (t), (3.88) where eN N N N u,m (t) = ψu,m (t) − ψu,m (t) and ec,m (t) = ψc,m (t) − ψc,m (t) are the m-th errors in the uncollided and collided components. The collided error can be further decomposed as N N N ηc,m (t) = ψc,m (t) − PN ψc,m (t) and ξc,m = PN ψc,m (t) − ψc,m (t) (3.89) so that eN N N N c,m (t) = ηc,m (t) + ξc,m (t). Here ηc,m is the m-th collided consistency error and ξc,m is N the m-th collided stability error. The error eN M is simply called the hybrid error. 60 This next lemma gives a one-step analysis of the growth of the error in the uncollided and collided components from tm−1 to tm . Lemma 3.3.2. Let 1 ≤ m ≤ M , then if t ∈ [tm−1 , tm ), the m-th uncollided and collided errors satisfy, respectively, −σ(t−tm−1) /ε 2 − ∥eNu,m ∥L2 (X×S2 ) (t) ≤ e ∥eNm−1 ∥L2 (X×S2 ) (tm−1 ), (3.90a)   −σ(t−tm−1) /ε2 − N ∥ec,m ∥L2 (X×S2 ) (t) ≤ 1 − e ∥eNm−1 ∥L2 (X×S2 ) (tm−1 ) N 1 N + ∥ηc,m ∥L2 (X×S2 ) (t) + ∥Ω · ∇x ηc,m ∥L1 ([tm−1 ,t];L2 (X×S2 )) . (3.90b) ε Proof. Subtracting (3.14a) from (3.12a), yields the following evolution equation for eN u,m , σ N − ε∂t eN N u,m + Ω · ∇x eu,m + e = 0, eN = eNm−1 (tm−1 ), (3.91) ε u,m u,m t=tm−1 − where eN 0 (t0 ) = 0. Thus applying (3.84a) from Lemma 3.3.1, with r = 0 and a zero source term yields (3.90a). To prove (3.90b), we subtract from (3.14b) the projection applied to N (3.12b). This gives the following PN equation for ξc,m N N σ N σ N N ε∂t ξc,m + PN (Ω · ∇x ξc,m )+ ξc,m = (ξc,m + eN u,m ) − PN (Ω · ∇x ηc,m ), (3.92a) ε ε N ξc,m t=tm−1 = 0. (3.92b) −1 We apply Lemma 3.2.6 to (3.92a) with zero initial data and source ε−2 σeN u,m − ε PN (Ω · N N ∇x ηc,m ). Combined with bound (3.90a), the estimate on ξc,m becomes Z t Z t N σ 1 ∥ξc,m ∥L2 (X×S2 ) (t) ≤ 2 ∥eN u,m ∥L2 (X×S2 ) (τ ) dτ + ∥Ω · ∇x ηc,m N ∥L2 (X×S2 ) (τ ) dτ ε tm−1 ε tm−1 Z t σ N − 2 ≤ 2 ∥em−1 ∥(tm−1 ) e−σ(τ −tm−1 )/ε dτ ε tm−1 Z t 1 N + ∥Ω · ∇x ηc,m ∥L2 (X×S2 ) (τ ) dτ ε tm−1  2  = 1 − e−σ(t−tm−1) /ε ∥eN − m−1 ∥L2 (X×S2 ) (tm−1 ) 1 t Z N + ∥Ω · ∇x ηc,m ∥L2 (X×S2 ) (τ ) dτ. ε tm−1 (3.93) N Adding ∥ηc,m ∥L2 (X×S2 ) (t) to the both sides recovers (3.90b). 61 Now we can state an error for all time that only depends on the approximation prop- erties of the spherical harmonic discretization on the solution. Theorem 3.3.1. The hybrid error eN M satisfies M   − X − 1 ∥eN M ∥L2 (X×S2 ) (tM ) ≤ N N ∥ηc,m ∥L2 (X×S2 ) (tm ) + ∥Ω · ∇x ηc,m ∥L1 ([tm−1 ,tm ];L2 (X×S2 )) . (3.94) m=1 ε Proof. Adding the inequalities in Lemma 3.3.2 and taking the limit t → t− m+1 gives − N − N − ∥eNM ∥L2 (X×S2 ) (tM ) ≤ ∥eM −1 ∥L2 (X×S2 ) (tM −1 ) + ∥ηc,M ∥L2 (X×S2 ) (tM ) (3.95) 1 N + ∥Ω · ∇x ηc,M ∥L1 ([tM −1 ,tM ];L2 (X×S2 )) ε − Exhausting this recursion until eN 0 (t0 ) = 0 yields the result. Lemma 3.3.3. Let s ≥ 0, N ≥ max{0, s − 1}, and 1 ≤ m ≤ M . Then the m-th projection error N ηc,m satisfies, ∆t N ∥ηc,m ∥L2 (X×S2 ) (t− m) ≤ sup |ψc,m |H 1,s (τ ). (3.96) ε(N + 1)s τ ∈[tm−1 ,tm ] N Proof. An application of Lemma 3.2.10 with α = tm−1 , β = tm , ηc,m (tm−1 ) = 0, an isotropic source q = σ ψ , ε2 u,m along with the fact that Aα [f ](t) ≤ ε−1 (t − α) supτ ∈[α,t] f (τ ) (cf. (3.40) with k = 1)) and the spectral estimate in Lemma 3.2.2, gives ∆t N ∥ηc,m ∥L2 (X×S2 ) (t−m) ≤ sup ∥P eN (Ω · ∇x ψc,m ))∥L2 (X×S2 ) (τ ) dτ ε τ ∈[tm−1 ,tm ] ∆t ∆t ≤ s sup |Ω · ∇x ψc,m |H 0,s (τ ) dτ ≤ sup |ψc,m |H 1,s (τ ). ε(N + 1) τ ∈[tm−1 ,tm ] ε(N + 1)s τ ∈[tm−1 ,tm ] (3.97) 3.3.3 Estimating hybrid error in terms of the data Finally, we will apply the approximation properties and stability estimates to the es- timate in Theorem 3.3.1 to obtain an estimate that depends only on the regularity of the data. 62 Theorem 3.3.2. Let s ≥ 1 and N ≥ s − 1. If q ∈ L1 ([0, T ]; H s+1,0 ) and g ∈ H s+1,0 , then − 2Cs   ∥eN ∥ 2 M L (X×S ) 2 (T ) ≤ |g| H s+1,0 + T |q|L∞ ([0,T ];H s+1,0 ) (3.98) (N + 1)s  s−1 ε s!T ∆ts T   ∆tσ × min s , s+1 min 1, 2 . (3.99) σ ε ε Proof. Using ∥Ω · ∇x ψc,m ∥L2 (X×S2 ) ≤ ∥∇x ψc,m ∥L2 (X×S2 ) , and applying Lemma 3.2.2 and Lemma 3.3.3 to Theorem 3.3.1 yields, for N ≥ s − 1, M ! N − 1 X ∆t 1 ∥eM ∥L2 (X×S2 ) (tM ) ≤ sup |ψc,m |H 1,s (τ ) + |ψc,m |L1 ([tm−1 ,tm ];H 1,s ) (N + 1)s m=1 ε τ ∈[tm−1 ,tm ] ε M 2∆t X ≤ sup |ψc,m |H 1,s (τ ) ε(N + 1)s m=1 τ ∈[tm−1 ,tm ] σ Applying Corollary 3.2.1 with q = ψ ε2 u,m and r = 1, M − σ ∆t X ∥eN s   M (tM )∥L2 (X×S2 ) ≤ 2Cs s! sup A t |ψ u,m |L 1 ([t m−1 ,•];H s+1,0 ) (τ ). ε3 (N + 1)s m=1 τ ∈[tm−1 ,tm ] m−1 (3.100) For the summand above, it follows from (3.84b), the monotonicity of Aα , and Lemma 3.2.3 that Astm−1 |ψu,m |L1 ([tm−1 ,•];H s+1,0 ) (τ ) ≤ |ψu,m |L1 ([tm−1 ,tm ];H s+1,0 ) Astm−1 [1] (τ )   sup sup τ ∈[tm−1 ,tm ] τ ∈[tm−1 ,tm ] ε2   ≤ |g|H s+1,0 + |q|L1 ([0,T ];H s+1,0 ) Astm−1 [1] (tm ) σ  s  ε2  s   ε 1 ∆t ≤ |g|H s+1,0 + |q|L1 ([0,T ];H s+1,0 ) min , . σ σ s s! ε (3.101) Pluggin the above bound into (3.100) yields (since T = M ∆t)  s ε ∆ts  − ∆tM   ∥eN M ∥L2 (X×S2 ) (tM ) ≤ 2Cs s! |g|H s+1,0 + |q|L1 ([0,T ];H s+1,0 ) min , ε(N + 1)s σ s s!εs  s−1  (3.102) 2Cs   ε s!T ∆ts T = |g|H s+1,0 + |q|L1 ([0,T ];H s+1,0 ) min , s+1 . (N + 1)s σs ε 63 On the other hand, it follows from (3.84c) that sup Astm−1 |ψu,m |L1 ([tm−1 ,•];H s+1,0 ) (τ )   τ ∈[tm−1 ,tm ] ≤ ε(|g|H s+1,0 + T |q|L∞ ([0,T ];H s+1,0 ) ) sup tm−1 [1](τ ) As+1 (3.103) τ ∈[tm−1 ,tm ] s+1 ! εs+1  1 ∆t ≤ ε(|g|H s+1,0 + T |q|L∞ ([0,T ];H s+1,0 ) ) min , . σ s+1 (s + 1)! ε Plugging this into (3.100) yields, ε s!T ∆ts+1 σT  s−1  N − 2Cs   ∥eM ∥L2 (X×S2 ) (tM ) ≤ |g|H s+1,0 + T |q|L∞ ([0,T );H s+1,0 ) min , . (N + 1)s σs εs+3 (3.104) Taking a minimum of the right hand sides of (3.102) and (3.104) yields the result. 3.4 Return to the original transport model In this section we will show error estimates for the model (3.1). The analogous dis- cretizations for the models are the following. For the non-splitting PN discretization we seek Ψε,N ∈ C([0, T ); X × PN (S2 )), satisfying σt ε,N  σt  ε∂t Ψε,N + PN (Ω · ∇x Ψε,N ) + Ψ = − εσa Ψε,N + εPN Q, (3.105a) ε ε Ψε,N |t=0 = PN g (3.105b) ε,N ε,N and for the hybrid, we seek Ψε,N m = Ψu,m + Ψc,m where for each m ∈ {1, 2, . . . , M } Ψε,N ε,N 2 2 2  u,m , Ψc,m ∈ C([tm−1 , tm ); X × L (S )) × C([tm−1 , tm ); X × PN (S )) (3.106) satifies σt ε,N ε∂t Ψε,N ε,N u,m + Ω · ∇x Ψu,m + Ψ = εQ, (3.107a) ε u,m  σt ε,N  σt  ε∂t Ψε,N c,m + PN Ω · ∇x Ψε,Nc,m + Ψ = − εσ a Ψε,N ε,N u,m + Ψc,m , (3.107b) ε c,m ε (3.107c)   g,  m = 1, Ψε,N c,m t=tm−1 = 0, Ψε,N u,m t=tm−1 = (3.107d) Ψε,N − ε,N −  u,m−1 (tm−1 ) + Ψc,m−1 (tm−1 ) m > 1.  64 We define the correspondent PN error and m-th hybrid error respectively as, eε,N = Ψε − Ψε,N and eε,N m = Ψm − Ψm ε ε,N Since ψ N = eσa t Ψε,N and ψm N = eσa t Ψε,Nm , applying Theorems 3.2.1 and 3.3.2 gives the following estimate for the PN error for the original transport model (3.1). Theorem 3.4.1. Let s ≥ 1 and N ≥ s − 1, Q ∈ L∞ ([0, T ]; H i,j ) and g ∈ H i,j for each i, j such that 0 ≤ j ≤ s, i + j ≤ s + 1. Then 2 2 e−(σt +ε σa )T /ε  2  ε,N 1 ε ∥e ∥L2 (X×S2 ) (T ) ≤ |g|H 0,s + |Q| L∞ ([0,T ];H 0,s ) min ,T (N + 1)s (N + 1)s σt  s+1 ! 1  −σa T  εs−1 s!T T +2Cs e |g| H s+1,0 + T |Q|L∞ ([0,T ];H s+1,0 ) min s , (N + 1)s σt ε s−1   i+1 (3.108) −(σt +ε2 σa )T X s T +e |g|H 1+i,s−i i=0 i εi+1 s−1    i+1 i+2 ! X s ε T 1 T + |Q|L∞ ([0,T ];H 1+i,s−i ) min i+1 , i=0 i σt (i + 1)! εi+1 Meanwhile for the hybrid approximation. Theorem 3.4.2. Let s ≥ 1 and N ≥ s − 1, Q ∈ L1 ([0, T ]; H s+1,0 ) and g ∈ H s+1,0 , we have 2Cs  −σa T  ∥eε,N − M ∥L2 (X×S2 ) (T ) ≤ e |g|H s+1,0 + T |Q|L∞ ([0,T ];H s+1,0 ) (N + 1)s  s−1 ε s!T ∆ts T   ∆tσt × min , s+1 min 1, 2 . (3.109) σts ε ε 65 CHAPTER 4 SUMMARY AND CONCLUSION In this thesis, we studied two numerical schemes for kinetic equations. In Chapter 2, we proved theoretically and demonstrated computationally the effectiveness of the SIAC fil- ter to the DG solutions of the nonlinear VM system. We proved the superconvergence of order (2k + 12 ) in the negative norm of the DG solutions. This is nontrivial for non- linear systems, and is achieved by identifying a suitable dual problem. The numerical experiments verify the performance of the filter in reducing spurious oscillations in the numerical errors. For low order k, the resolution of the numerical solution is greatly en- hanced, which is highly desirable for long time kinetic simulations. In the future, we plan to prove superconvergence for the divided difference of the numerical solution to fully justify the enhanced resolution of the post-processed solution. Another interesting project will be to apply this post-processing technique to different kinetic equations that use the DG method. In Chapter 3, we derived multiscale error estimates for the PN approximation of the RTE and for a hybrid approximation for the RTE that is built using the PN approximation. By construction, the hybrid is is more expensive; we use these error estimates to under- stand the benefits of the additional expense for different parameter regimes. At each time step in the hybrid approximation, the collided equation is equipped with isotropic ini- tial conditions and zero initial condition. In scattering dominating regimes, this property is key to improved estimates over the monolithic PN approach. Meanwhile, in purely absorbing regimes, the hybrid captures the RTE solution exactly. In the future, we in- tend to revisit the current analysis for more general problems on non-periodic domains, with non-constant cross-sections and inflow boundary conditions. In addition, we intend to explicitly examine the effects of angular discretization errors in the treatment of the uncollided equation, which for the purposes of the current chapter was assumed to be solved exactly. 66 BIBLIOGRAPHY [1] Andrés Galindo-Olarte, Juntao Huang, Jennifer K Ryan, and Yingda Cheng. Su- perconvergence and accuracy enhancement of discontinuous galerkin solutions for vlasov-maxwell equations. arXiv preprint arXiv:2210.07908, 2022. [2] Cory D Hauck and Ryan G McClarren. A collision-based hybrid method for time- dependent, linear, kinetic transport equations. Multiscale Modeling & Simulation, 11(4):1197–1227, 2013. [3] Andrés Galindo-Olarte, Victor P DeCaria, and Cory D Hauck. Numerical analysis of a hybrid method for radiation transport. arXiv:2306.04714, 2023. [4] Lorenzo Pareschi. Kinetic equations: computation. arXiv:1311.7230v1, 2013. [5] Giacomo Dimarco and Lorenzo Pareschi. Numerical methods for kinetic equations. Acta Numerica, 23:369–520, 2014. [6] F. Califano, F. Pegoraro, S. V. Bulanov, and A. Mangeney. Kinetic saturation of the Weibel instability in a collisionless plasma. Phys. Rev. E, 57(6):7048–7059, 1998. [7] A. Mangeney, F. Califano, C. Cavazzoni, and P. Travnicek. A numerical scheme for the integration of the Vlasov-Maxwell system of equations. J. Comp. Phys., 179(2):495–538, 2002. [8] F. Califano, F. Pegoraro, and S. V. Bulanov. Impact of kinetic processes on the macro- scopic nonlinear evolution of the electromagnetic-beam-plasma instability. Phys. Rev. Lett., 84:3602–3605, 2000. [9] F. Califano, N. Attico, F. Pegoraro, G. Bertin, and S. V. Bulanov. Fast formation of magnetic islands in a plasma in the presence of counterstreaming electrons. Phys. Rev. Lett., 86(23):5293–5296, 2001. [10] N.J. Sircombe and T.D. Arber. VALIS: A split-conservative scheme for the relativistic 2d Vlasov-Maxwell system. J. Comp. Phys., 228(13):4773 – 4788, 2009. [11] Nicolas Besse, Guillaume Latu, Alain Ghizzo, Eric Sonnendrüker, and Pierre Bertrand. A wavelet-MRA-based adaptive semi-Lagrangian method for the rela- tivistic Vlasov-Maxwell system. J. Comp. Phys., 227(16):7889 – 7916, 2008. [12] Akihiro Suzuki and Toshikazu Shigeyama. A conservative scheme for the relativistic Vlasov-Maxwell system. J. Comp. Phys., 229(5):1643 – 1660, 2010. [13] F. Huot, A. Ghizzo, P. Bertrand, E. Sonnendrüker, and O. Coulaud. Instability of the time splitting scheme for the one-dimensional and relativistic Vlasov-Maxwell system. J. Comp. Phys., 185(2):512 – 531, 2003. [14] Bernardo Cockburn and Chi-Wang Shu. Runge-Kutta discontinuous Galerkin meth- ods for convection-dominated problems. J. Sci. Comput., 16:173–261, 2001. 67 [15] R. E. Heath, I. M. Gamba, P. J. Morrison, and C. Michler. A discontinuous Galerkin method for the Vlasov-Poisson system. J. Comp. Phys., 231:1140–1174, 2012. [16] R. E. Heath. Numerical analysis of the discontinuous Galerkin method applied to plasma physics. 2007. Ph. D. dissertation, the University of Texas at Austin. [17] Y. Cheng, I. M. Gamba, and P. J. Morrison. Study of conservation and recurrence of Runge-Kutta discontinuous Galerkin schemes for Vlasov-Poisson systems. J. Sci. Comp. accepted, 2012. preprint arXiv:1209.6413v2 [math.NA]. [18] Y. Cheng and I. M. Gamba. Numerical study of Vlasov-Poisson equations for infinite homogeneous stellar systems. Comm. Nonlin. Sci. Num. Sim., 17, 2012. [19] Yingda Cheng, Irene M Gamba, Fengyan Li, and Philip J Morrison. Discontinuous Galerkin methods for the Vlasov-Maxwell equations. SIAM Journal on Numerical Analysis, 52(2):1017–1049, 2014. [20] Yingda Cheng, Andrew J Christlieb, and Xinghui Zhong. Energy-conserving discon- tinuous Galerkin methods for the Vlasov–Ampere system. Journal of Computational Physics, 256:630–655, 2014. [21] He Yang and Fengyan Li. Discontinuous galerkin methods for relativistic Vlasov– Maxwell system. Journal of Scientific Computing, 73(2):1216–1248, 2017. [22] James H Bramble and Alfred H Schatz. Higher order local accuracy by averaging in the finite element method. Mathematics of Computation, 31(137):94–111, 1977. [23] Bernardo Cockburn, Mitchell Luskin, Chi-Wang Shu, and Endre Süli. Enhanced accuracy by post-processing for finite element methods for hyperbolic equations. Mathematics of Computation, 72(242):577–606, 2003. [24] Liangyue Ji, Yan Xu, and Jennifer K Ryan. Negative-order norm estimates for non- linear hyperbolic conservation laws. Journal of Scientific Computing, 54(2):531–548, 2013. [25] Xiong Meng and Jennifer K Ryan. Discontinuous Galerkin methods for nonlinear scalar hyperbolic conservation laws: divided difference estimates and accuracy en- hancement. Numerische mathematik, 136(1):27–73, 2017. [26] Xiong Meng and Jennifer K Ryan. Divided difference estimates and accuracy en- hancement of discontinuous Galerkin methods for nonlinear symmetric systems of hyperbolic conservation laws. IMA Journal of Numerical Analysis, 38(1):125–155, 2018. [27] Michael Steffan, Sean Curtis, Robert M Kirby, and Jennifer Ryan. Investigation of smoothness enhancing accuracy-conserving filters for improving streamline integra- tion through discontinuous fields. IEEE Transactions on Visualization and Computer Graphics, 14(3):680–692, 2008. 68 [28] Liangyue Ji, Paulien Van Slingerland, Jennifer K Ryan, and Kees Vuik. Super- convergent error estimates for position-dependent smoothness-increasing accuracy- conserving (SIAC) post-processing of discontinuous Galerkin solutions. Mathematics of computation, pages 2239–2262, 2014. [29] Gerald C Pomraning. The equations of radiation hydrodynamics. Courier Corporation, 2005. [30] Elmer Eugene Lewis and Warren F Miller. Computational methods of neutron transport. John Wiley and Sons, Inc., New York, NY, 1984. [31] K.M. Case and P.F. Zweifel. Linear Transport Theory. Addison-Wesley series in nuclear engineering. Addison-Wesley Publishing Company, 1967. [32] Edward W Larsen and Joseph B Keller. Asymptotic solution of neutron transport problems for small mean free paths. Journal of Mathematical Physics, 15(1):75–81, 1974. [33] Alain Bensoussan, Jacques L Lions, and George C Papanicolaou. Boundary layers and homogenization of transport processes. Publications of the Research Institute for Mathematical Sciences, 15(1):53–157, 1979. [34] Mohammed Lemou and Luc Mieussens. A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit. SIAM Journal on Scientific Computing, 31(1):334–368, 2008. [35] Luis Chacon, Guangye Chen, Dana A Knoll, C Newman, H Park, William Taitano, Jeff A Willert, and Geoffrey Womeldorff. Multiscale high-order/low-order (holo) algorithms and applications. Journal of Computational Physics, 330:21–45, 2017. [36] Marvin L Adams and Edward W Larsen. Fast iterative methods for discrete- ordinates particle transport calculations. Progress in nuclear energy, 40(1):3–159, 2002. [37] James S Warsa, Todd A Wareing, and Jim E Morel. Krylov iterative methods and the degraded effectiveness of diffusion synthetic acceleration for multidimensional sn calculations in problems with material discontinuities. Nuclear science and engineer- ing, 147(3):218–248, 2004. [38] Raymond E Alcouffe. A first collision source method for coupling monte carlo and discrete ordinates for localized source problems. In Monte-Carlo Methods and Applica- tions in Neutronics, Photonics and Statistical Physics: Proceedings of the Joint Los Alamos National Laboratory-Commissariat à l’Energie Atomique Meeting Held at Cadarache Castle, Provence, France April 22–26, 1985, pages 352–366. Springer, 2006. [39] Michael M Crockatt, Andrew J Christlieb, C Kristopher Garrett, and Cory D Hauck. An arbitrary-order, fully implicit, hybrid kinetic solver for linear radiative transport using integral deferred correction. Journal of Computational Physics, 346:212–241, 2017. 69 [40] Michael M Crockatt, Andrew J Christlieb, C Kristopher Garrett, and Cory D Hauck. Hybrid methods for radiation transport using diagonally implicit runge–kutta and space–time discontinuous galerkin time integration. Journal of Computational Physics, 376:455–477, 2019. [41] Michael M Crockatt, Andrew J Christlieb, and Cory D Hauck. Improvements to a class of hybrid methods for radiation transport: Nyström reconstruction and defect correction methods. Journal of Computational Physics, 422:109765, 2020. [42] Ben Whewell, Ryan G McClarren, Cory D Hauck, and Minwoo Shin. Multigroup Neutron Transport Using a Collision-Based Hybrid Method. Nuclear science and en- gineering, 2023. [43] Vincent Heningburg and Cory D Hauck. Hybrid solver for the radiative trans- port equation using finite volume and discontinuous galerkin. arXiv preprint arXiv:2002.02517, 2020. [44] Martin Frank, Cory Hauck, and Kerstin Küpper. Convergence of filtered spherical harmonic equations for radiation transport. Communications in Mathematical Sciences, 14(5):1443–1465, 2016. [45] Zheng Chen and Cory Hauck. Multiscale convergence properties for spectral ap- proximations of a model kinetic equation. Mathematics of Computation, 88(319):2257– 22 93, 2019. [46] Sean Curtis, Robert M Kirby, Jennifer K Ryan, and Chi-Wang Shu. Postprocessing for the discontinuous Galerkin method over nonuniform meshes. SIAM Journal on Scientific Computing, 30(1):272–289, 2008. [47] Philippe G Ciarlet. The finite element method for elliptic problems. Bull. Amer. Math. Soc, 1:800–802, 1979. [48] Guri Ivanovich Marchuk. Construction of adjoint operators in non-linear problems of mathematical physics. Sbornik: Mathematics, 189(10):1505, 1998. [49] Michel Crouzeix and Vidar Thomée. The stability in Lp and Wp1 of the L2 -projection onto finite element function spaces. Mathematics of Computation, 48(178):521–532, 1987. [50] Blanca Ayuso de Dios, José Antonio Carrillo de la Plata, and C-W Shu. Discontinu- ous Galerkin methods for the one-dimensional Vlasov-Poisson system. 2009. [51] Susanne C Brenner, L Ridgway Scott, and L Ridgway Scott. The mathematical theory of finite element methods, volume 3. Springer, 2008. [52] Sigal Gottlieb and Chi-Wang Shu. Total variation diminishing Runge-Kutta schemes. Mathematics of computation, 67(221):73–85, 1998. 70 [53] Robert Dautray and Jacques-Louis Lions. Mathematical analysis and numerical methods for science and technology: volume 6 evolution problems II, volume 6. Springer Science & Business Media, 1999. [54] Cory D Hauck and Robert B Lowrie. Temporal regularization of the PN equations. Multiscale Modeling & Simulation, 7(4):1497–1524, 2009. [55] Kendall Atkinson and Weimin Han. Spherical harmonics and approximations on the unit sphere: an introduction, volume 2044. Springer Science & Business Media, 2012. [56] Feng Dai and Yuan Xu. Approximation Theory and Harmonic Analysis on Spheres and Balls. Springer, 01 2013. 71 APPENDIX Proof of Lemma 2.3.3 By using equation (2.21a), the divergence free properties of A1 , A2 and the boundary conditions, we have the following Z 1d ∥φ∥2 = − (A3 · F)φ dxdv ≤ C(∥φ∥2 + ∥F∥2 ), 2 dt Ω where C depends on ∥A3 ∥L∞ ((0,T );L∞ (Ω)) . On the other hand using equations (2.21b) and (2.21c), Gauss theorem on the physical space integrals and integration by parts on the velocity space variables, Z Z 1d 1d ∥F∥2 + ∥D∥2 = (∇x × D · F − ∇x × F · D) dx − φ∇v g · F dxdv 2 dt 2 dt Ωx Ωv Z + φ(v × ∇v g)D dxdv Ωv Z Z = − φ∇v g · F dxdv + φ(v × ∇v g)D dxdv Ω Ω ≤ C ∥F∥2 + ∥D∥2 + ∥φ∥ 2  , where C depends on ∥g∥L∞ ((0,T );W 1,∞ (Ω)) . Now we add the tow inequalities above, to obtain 1d 1d 1d ∥φ∥2 + ∥F∥2 + ∥D∥2 ≤ C ∥F∥2 + ∥D∥2 + ∥φ∥2 ,  (.1) 2 dt 2 dt 2 dt where C depends on ∥A3 ∥L∞ ((0,T );L∞ (Ω)) and ∥g∥L∞ ((0,T );W 1,∞ (Ω)) . An application of Gron- wall’s inequality allow us to conclude. Now since we are considering the full Sobolev norm, we still need to estimate the L2 norms of the higher order derivatives ∂xβ ∂vγ , to do so we apply ∂xβ ∂vγ to the system (2.21) and then we repeat the same steps that we took above. Fundamental ODE result If | · | denotes either one of the semi-norms and norms defined throughout the chap- ter, one of the usual procedures is finding a bound for |φ|, by analyzing a differential 72 inequality, that looks like 1d 2 |φ| + κ|φ|2 ≤ χ(t)|φ|, 2 dt where κ ≥ 0 is a constant and χ(t) ≥ 0 for all times t. Formally one would divide by |φ| and integrate in time. However this computation is not rigorous if there exists some time t∗ , for which |φ|(t∗ ) = 0. The following lemma makes the formal calculation rigorous. Lemma .0.1. Assume χ is a non-negative continuous function on [α, β]. Assume ϕ ∈ C 1 ([α, β]), ϕ ≥ 0, and satisfies the following differential inequality 1 (ϕ(t)2 )′ + κϕ(t)2 ≤ χ(t)ϕ(t), ϕ(α) = ϕα ≥ 0. (.2) 2 Then for all t ∈ [α, β], Z t Z t ϕ(t) + κ ϕ(τ ) dτ ≤ ϕα + χ(τ ) dτ. (.3) α α Furthermore Z t ϕ(t) ≤ e −κ(t−α) ϕα + e−κ(t−τ ) χ(τ ) dτ. (.4) α Proof. We prove first (.3). Since ϕ and χ are non-negative functions, it follows that for any arbitrary δ > 0, the following differential inequality holds ϕ(t)ϕ′ (t) + κϕ(t)2 ≤ χ(t)(ϕ(t) + δ), (.5) dividing both sides of the inequality by ϕ + δ, and integrating in time, we arrive at Z t Z t Z t ϕ(t) + δ ϕ(τ ) ϕ(t) + κ ϕ(τ ) dτ ≤ ϕα + χ(τ ) dτ + δ ln + κδ dτ, α α ϕα + δ α ϕ(τ ) + δ Z t ϕ(t) + δ ≤ ϕα + χ(τ ) dτ + δ ln + κδ(t − α), α ϕα + δ the conclusion follows taking δ → 0+ . We next prove (.4). When κ = 0 the result follows immediately from (.3): Z t ϕ(t) ≤ ϕα + χ(τ ) dτ. (.6) α For the general case we multiply (.2) by e2κt , obtaining 1 ′ Φ(t)2 ≤ eκt χ(t)Φ(t), (.7) 2 73 where Φ(t) = eκt ϕ(t). Applying (.6) to Φ and undoing the transformation yields (.4). 74