Section: Scientific Foundations
Keywords : Numerical methods, Vlasov equation, unstructured grids, adaptivity, numerical analysis, convergence, semiLagrangian method, PIC method.
Development of simulation tools
Introduction
The numerical integration of the Vlasov equation is one of the key challenges of computational plasma physics. Indeed, kinetic models are posed in phase space and thus the number of dimensions is doubled. Since the early days of this discipline, an intensive work on this subject has produced many different numerical schemes. We have been interested in two of the most promising. The first and by far the most widely used is the Particle In Cell (PIC) method which approximates the distribution function, solution of the Vlasov equation, by Dirac measures which act like macroparticles. It is independent of dimension and thus becomes very efficient when dimension increases which is interesting for the Vlasov equation posed in phase space. However its convergence rate is slow. The second uses a phasespace grid for the solution of the Vlasov equation. In both cases, the Vlasov solver is coupled with a physical grid based field solver. At the beginning of the project we were mainly involved in the latter approach. In order to make such methods efficient, it is essential to consider means for optimizing the number of mesh points. This is done through different adaptive strategies. In order to understand the methods, it is also important to perform their mathematical analysis. However it is still out of reach to address full threedimensional problems with a phasespace grid method. For these kind of large problems we have also started recently to work on efficient PIC methods for threedimensional problems with complex geometries.
Convergence analysis of numerical schemes
Exploring grid based methods for the Vlasov equation, it becomes obvious that they have different stability and accuracy properties. In order to fully understand what are the important features of a given scheme and how to derive schemes with the desired properties, it is essential to perform a thorough mathematical analysis of this scheme, investigating in particular its stability and convergence towards the exact solution.
The semiLagrangian method
The semiLagrangian method consists in computing a numerical approximation of the solution of the Vlasov equation on a phase space grid by using the property of the equation that the distribution function f is conserved along characteristics. More precisely, for any times s and t , we have
where are the characteristics of the Vlasov equation which are solution of the system of ordinary differential equations
with initial conditions , .
From this property, f^{n} being known one can induce a numerical method for computing the distribution function f^{n+ 1} at the grid points consisting in the following two steps:

For all i, j , compute the origin of the characteristic ending at , i.e. an approximation of , .

As f^{n+ 1} can be computed by interpolating f^{n} which is known at the grid points at the origin of the characteristics .
This method can be simplified using a timesplitting procedure separating the advection phases in physical space and velocity space, as in this case the characteristics can be solved explicitly.
Adaptive semiLagrangian methods
Uniform meshes are most of the time not efficient to solve a problem in plasma physics or beam physics as the distribution of particles is evolving a lot as well in space as in time during the simulation. In order to get optimal complexity, it is essential to use meshes that are fitted to the actual distribution of particles. If the global distribution is not uniform in space but remains locally mostly the same in time, one possible approach could be to use an unstructured mesh of phase space which allows to put the grid points as desired. Another idea, if the distribution evolves a lot in time is to use a different grid at each time step which is easily feasible with a semiLagrangian method. And finally, the most complex and powerful method is to use a fully adaptive mesh which evolves locally according to variations of the distribution function in time. The evolution can be based on a posteriori estimates or on multiresolution techniques.
ParticleInCell codes
The ParticleInCell method [60] consists in solving the Vlasov equation using a particle method, i.e. advancing numerically the particle trajectories which are the characteristics of the Vlasov equation, using the equations of motion which are the ordinary differential equations defining the characteristics. The selffields are computed using a standard method on a structured or unstructured grid of physical space. The coupling between the field solve and the particle advance is done on the one hand by depositing the particle data on the grid to get the charge and current densities for Maxwell's equations and, on the other hand, by interpolating the fields at the particle positions. This coupling is one of the difficult issues and needs to be handled carefully.
Maxwell's equations in singular geometry
The solutions to Maxwell's equations are a priori defined in a function space such that the curl and the divergence are square integrable and that satisfy the electric and magnetic boundary conditions. Those solutions are in fact smoother (all the derivatives are square integrable) when the boundary of the domain is smooth or convex. This is no longer true when the domain exhibits nonconvex geometrical singularities (corners, vertices or edges).
Physically, the electromagnetic field tends to infinity in the neighborhood of the reentrant singularities, which is a challenge to the usual finite element methods. Nodal elements cannot converge towards the physical solution. Edge elements demand considerable mesh refinement in order to represent those infinities, which is not only time and memoryconsuming, but potentially catastrophic when solving timedependent equations: the CFL condition then imposes a very small time step. Moreover, the fields computed by edge elements are discontinuous, which can create considerable numerical noise when the Maxwell solver is embedded in a plasma (e.g. PIC) code.
In order to overcome this dilemma, a method consists in splitting the solution as the sum of a regular part, computed by nodal elements, and a singular part which we relate to singular solutions of the Laplace operator, thus allowing to calculate a local analytic representation. This makes it possible to compute the solution precisely without having to refine the mesh.
This Singular Complement Method (SCM) had been developed [56] and implemented [55] in plane geometry.
An especially interesting case is axisymmetric geometry. This is still a 2D geometry, but more realistic than the plane case; despite its practical interest, it had been subject to much fewer theoretical studies [59] . The nondensity result for regular fields was proven [62] , the singularities of the electromagnetic field were related to that of modified Laplacians [52] , and expressions of the singular fields were calculated [53] . Thus the SCM was extended to this geometry. It was then implemented by F. Assous (now at BarIlan University, Israel) and S. Labrunie in a PIC–finite element Vlasov–Maxwell code [54] .
As a byproduct, spacetime regularity results were obtained for the solution to timedependent Maxwell's equation in presence of geometrical singularities in the plane and axisymmetric cases [70] , [53] .