SPARSE GRID RECONSTRUCTIONS FOR PARTICLE-IN-CELL METHODS

. In this article, we propose and analyse Particle-In-Cell (PIC) methods embedding sparse grid reconstructions such as those introduced in Ricketson and Cerfon [ Plasma Phys. Control. Fusion 59 (2017) 024002] and Muralikrishnan et al . [ J. Comput. Phys. X 11 (2021) 100094]. The sparse grid reconstructions offer a significant improvement on the statistical error of PIC schemes as well as a reduction in the complexity of the problem providing the electric field. Main results on the convergence of the electric field interpolant and conservation properties are provided in this paper. Besides, tailored sparse grid reconstructions, in the frame of the offset combination technique, are proposed to introduce PIC methods with improved efficiency. The methods are assessed numerically and compared to existing PIC schemes thanks to classical benchmarks with remarkable prospects for three dimensional computations.


Introduction
Particle-In-Cell (PIC) discretizations have been among the most used numerical methods in the simulation of kinetic plasmas for years [3,9,22,31] and are still topical [10,12,14,27].The method consists in a coupling between a Lagrangian method for the Vlasov equation, based on the integration of numerical particle trajectories, and a mesh-based discretization of Poisson's equation (or Maxwell's system) for the computation of the self-consistent field.Despite their simplicity, ease of parallelization and robustness, Particle-In-Cell schemes still contain a significant drawback: the statistical error originating from the sampling of the distribution function by numerical particles.This numerical noise decreases slowly with the increase of the average number of particles per cell.Therefore, a large number of particles may be required, necessitating tremendous computational resources, specifically for three dimensional simulations for which the desired precision may impose a number of cells as large as 10 9 , the number of particles exceeding 10 12 .Noise reduction strategies aim at maintaining the accuracy of computations with a reduced set of particles.They have therefore received a lot of attention with, for instance, variance reduction methods such as the  method [11] or the quiet start initialization procedure [31] as well as filtering methods in either Fourier domain [4], wavelet domain [17] or variants [33].
Sparse grid methods [6,13], originally developed for the interpolation of high dimensional functions, then extended to the approximation of partial differential equations [5,18,19,21], have recently been applied, in the framework of the so-called saprse grid combination technique [7,20], to Particle-In-Cell schemes [8,15,16,24,29].The aim here is to improve the properties of PIC methods with respect to the statistical error resulting from the particle sampling.In the sparse grid reconstructions, the numerical approximations are recomposed from partial representations carried out on a hierarchy of sparse grids with coarse resolutions.Compared to a regular Cartesian grid, the mean number of particles per cell is larger for any of the sparse grids.This crucial feature offers either a mitigation of the statistical noise or a decrease of the total number of numerical particles for a precision comparable to standard PIC discretizations.Besides, considering the thorough studies conducted during recent years to apply the combination technique to the resolution of PDEs, promising improvements in the computational efficiency are expected for the resolution of Poisson's equation (see [7,26,28]) providing the electric field in the Particle-In-Cell framework.Early implementations of sparse grids Particle-In-Cell schemes [15,16,29] show computational gains which are forecasted substantial for three dimensional applications.
The objective of the present paper is to provide an overview of the existing Particle-In-Cell discretizations implementing the sparse grid combination technique, to conduct a formal analysis to explain their merits and weaknesses and support the development of new methods with improved efficiency.The proofs of the results provided herein are limited to two dimensional geometries.Nonetheless, they are readily extendable to an arbitrary number of dimensions, using the same tools, however with an additional complexity of notation avoided within the present document.
The analyses of standard as well as sparse grid PIC discretizations reveal that, for both methods, the approximation error may be decomposed into three contributions.The precision of the methods is characterized by the accuracy of the most probable value of the statistics associated to the particle sampling, this component being referred to as the bias.This is a grid-based error related to both the mesh size (ℎ) and the smoothness of the solution with a component depending on the mixed derivatives of the solution and another contribution depending on non-mixed derivatives.The last error component is the so-called numerical noise or particle sampling error, providing the magnitude of the dispersion of the values attached to a sample of particles.The introduction of sparse grid reconstructions within PIC discretizations entails an increase of the grid error together with a significant mitigation of the statistical noise.This outlines the potential of these approaches: sparse grid reconstructions may be tailored to define different trades-off between the components of the error and finally mitigate the most detrimental one for the precision of PIC numerical approximations (the statistical noise).This leads to the derivation of the new sparse-grid methods introduced herein, with an improved numerical efficiency.
A specific attention is payed to the approximation of the electric field which can be computed thanks to two different approaches.The first one consists in computing the electric field on a refined Cartesian mesh thanks to the sparse interpolant of the charge density obtained by the combination technique.The second relies on a computation on each subgrid using the projected density.The electric field interpolant is then obtained by recombining the local approximants of the sub-grids.The analyses conducted in this paper are aimed at demonstrating the convergence of the electric field sparse grid interpolant and highlight the differences between these two approaches.The main results regarding the electric field, given by the Propositions 3.7 and 3.8, point the strong dependance of the error on the mixed derivatives of the solution, especially for the second approach, which has proven to be more dependant of the smoothness of the solution.This is a major contribution of the paper since no convergence properties have already been proposed so far for the electric field, which is the critical quantity when determining the overall accuracy of PIC discretizations [32].
The paper is organized as follows.In Section 2, the Particle-In-Cell scheme is outlined in its conventional framework with the definition of the grid-based and the particle sampling errors.In Section 3.1, sparse grids approximations are introduced in the specific framework of the saprse grid combination technique before being merged with Particle-In-Cell methods in Section 3.2.In Section 3.3, a generalization of the combination technique, referred to as offset combination technique, is proposed in order to tune the trade-off between the different components of the sparse grid approximation error.Section 3.4 is devoted to the introduction of two enhanced sparse PIC methods to improve the efficiency of the electric field computation.In Section 4, the merits of the different methods are investigated thanks to two dimensional computations relating classical plasma physics test cases: the linear and non linear Landau damping as well as the diocotron instability.The conclusions are drawn in Section 5, with an emphasize of the remarkable prospects of sparse grid Particle-In-Cell methods for three dimensional computations.

Notations
Let us introduce some notations for multi-variate functions and functional spaces.Let  be the dimension of the problem considered.For a multi-index  = ( 1 , . . .,   ) ∈ N  , a function  defined on the -dimensional unit interval Ω = [0, 1]  , we denote by     the partial derivative of  with respect to   and order   for  ∈ {1, . . ., }.We introduce the following functional spaces: where (Ω) denotes the space of continuous functions on Ω,   0 (Ω) and   0 (Ω) the spaces of functions vanishing on the boundary.Let us define order relations on multi-indexes by: and consider the supremum norm for a function  belonging to one of the previous spaces, the l 1 norm, l ∞ norm for a multi-index  ∈ N  be: For a vector function u : Ω → R  ,  ∈ N, we introduce the supremum norm notation: Let us also consider the family of -dimensional anisotropic grids on Ω called sub-grids, denoted Ω ℎ l and parametrized by the multi-index l ∈ : and: where ℎ l := (ℎ 1 , . . ., ℎ   ) := 2 −l is called the grid discretization.We consider also a regular isotropic grid, named Cartesian grid, denoted Ω (∞) ℎ , corresponding to a sub-grid of level n =  • 1 with discretization ℎ  in all directions, which is typically the underlying grid in PIC methods or standard interpolation: (2.9)

Standard Particle-In-Cell (PIC-Std)
In this section, the non-relativistic Vlasov-Poisson system with external magnetic field B is considered: In this problem,   (x, v, ) is the phase-space distribution function attached to the species ;   ,   are the corresponding charge and mass, E is the electric field and  is the charge density obtained from the phase-space distribution of each species: ℎ with finite differences according to: Considering a leap frog scheme and a second order finite difference scheme, the time discretization error is  (︀ ∆ 2 )︀ and the field solver error is  (︀ ℎ 2

𝑛
)︀ (see Prop. 3.7), where ∆ is the time step, ℎ  is the grid discretization corresponding to the cell size of the Cartesian grid.

Projection of the density onto the grid
The charge density of any species is approximated on a grid from a collection of numerical particles with a certain error that we shall make explicit here.For ease of presentation in the following, we consider one type of particle, ommit time dependance and let  =   / be the probability density function associated to the phase-space distribution   , where  = / is the total number of physical particles and  the total charge of the particles.Starting from the definition stated by equation (2.11), the density is recast into the following integral: (2.14) The rewriting in equation (2.14) allow us to define a numerical approximation of the density in the following.Let us consider a sub-grid of level l ∈ N  with a discretization ℎ l (corresponding to the cell width) and, as an ersatz of the convolution kernel, a -dimensional shape function, denoted  ℎ l and represented on Figure 1, which is constructed by tensor products of one dimensional hat functions: Substituting the convolution kernel with the shape function, we define the approximation of the density with: (2.17) In order to approximate this multidimensional integral, we consider a Monte Carlo framework using statistical sampling techniques.Let X be a random variable with the following probability density function, expected value and variance: The quantity  ℎ l (x − X) being a random variable and the function defined by  ↦ →  ℎ l (x − ) being a mesurable function of R  , owing to Transfer theorem [30], the integral in (2.17) can be recast into an expected value: Using an independent random sample (X 1 , . . ., X  ) from the random variable X, one can define an estimator of the integral by: where  is the number of random variables in the sample.To be able to make use of this result in particle simulations, we note that we can associate the values of the random sample (X 1 , . . ., X  ) with the particle positions [1] (x 1 , . . ., x  ), yielding the following statistical estimator: The expected value and variance of this estimator are defined from the independant random sample (X 1 , . . ., X  ) by: Definition 2.1.The local error between the density and its statistical estimator is recast into two components: . (2.23) The first, refered to as the particle sampling error and denoted  ,ℎ l , is a centered random variable corresponding to the error stemming from the variance of the sampling with a finite number of particles.The second is the bias of the estimator, denoted Bias(ρ ℎ l ) also reered to as the grid based error.It depends on both the mesh size and the solution smoothness and measures how close to the density the most probable value of the estimator is.
Proposition 2.2.Let  (•, v) ∈  2 (Ω), then the square root variance of the particle sampling error and the grid-based error are given by: Proof of Proposition 2.2.This result is a generalization to anisotropic grids of the estimates provided for instance in [30] for Cartesian meshes.
Remark 2.3.The estimate for a regular isotropic grid, used for regular PIC methods is recovered setting l =  • 1.In this section, the sparse grid notations will be introduced for the specific framework of the so-called sparse grid combination technique [7,20].Let (l, ) ∈ R  × N be the level and degree of basis functions defined by tensor products of unidmensional splines as follows:

Merging Particle-In-Cell with sparse grids
where j ∈  l is the index of the basis function.The unidimensional spline functions are defined recursively starting from the spline of degree 0 that we shall denoted  0 , and higher order splines of degree , for  ∈ N * defined by: In this paper, we consider only linear splines and ommit the superscript  in the following.It yields the following definition for the unidimensional linear splines: Let us introduce the spaces of -dimensional spline functions with respect to Ω ℎ l and Ω (∞) ℎ , denoted  ℎ l and  (∞) ℎ , and defined by: where { ℎ l ,j | j ∈  l } is called the nodal basis of the space  ℎ l and  l the nodal basis index set.Eventually, we define the interpolation of a function  onto the spaces  ℎ l ,  (∞) ℎ by: where the coefficients  l,j are determined by the resolution of a linear system and falls down to the nodal values of the function for linear splines ( l,j =  (jℎ l )).
Remark 3.1.The basis functions verify a partition of unit property: (3.6)

Sparse grid combination technique
The sparse grid combination technique [7,20] is a method of interpolation using evaluations of the function on the nodes of sub-grids.The interpolant is obtained by a linear combination of partial representations of the function on the sub-grids.Let  be a function and  ℎ l an approximation of this function in the space  ℎ l (e.g.  ℎ l  ), then a sparse grid reconstruction, denoted   ℎ is defined by linear combination of the contributions  ℎ l of each sub-grid (see Fig. 2): and  then: where for  = 2,  = 5 4 ‖ 1,2 (•; ℎ 1 , ℎ 2 )‖ ∞ .Proof of Proposition 3.2.See [20] for a proof in the two-dimensional and three-dimensional cases.
The combination technique is remarkable for the reduction of the number of interpolation nodes, explicited by the relations: while achieving nearly the same precision than the standard interpolation (with a negligible multiplicative term | log ℎ  | −1 ).
Remark 3.3.The sparse grid reconstruction of a nonnegative function is not nonnegative.

Application to Particle-In-Cell discretizations
In this section, two applications of the sparse grid combination technique to Particle-In-Cell methods are presented.In the regular PIC approximation, the main drawback is the statistical error decreasing with the number of particles per cell.Indeed, the particle sampling and the grid based errors scale respectively as

𝑛
)︀ , ℎ  being the mesh size of the Cartesian grid denoted Ω (∞) ℎ and  the total number of particles.These two error estimates together lead to the following onerous conditions for convergence of the scheme ℎ 2  ≪ 1,  ℎ   ≫ 1 which can require an extremely large number of particles, specifically for three dimensional simulations.The combination technique achieves a representation of a function using a sequence of sparse grids, coarser than the standard full mesh.This ends up in a reduced number of interpolation nodes and, accordingly, an increased mean number of particles per cell.This feature motivates the application of sparse grid techniques to PIC methods.
The bias of the estimator of the density on a sub-grid verifies the assumption of equation (3.8).Therefore, the combination of these projections onto the sub-grids according to equation (3.7) shall lead to cancellations for the grid-based error.This feature, resulting from the tensor product form of the shape function, is the motivation for the application of sparse grid combination to PIC methods.

Discretization on hybrid grids (PIC-Hg)
In this section, a first sparse grid application to PIC methods, introduced in [24], is presented.In order to take advantage of the cancellations exposed in the precedent section, a sparse grid interpolant of the density is constructed.The density is accumulated onto each sub-grid, achieving a reduction of the statistical error thanks to the large cells of the sub-grids, then evaluated onto the Cartesian grid with the combination technique.The electric field is obtained by resolving the Poisson equation on the Cartesian grid.Let us introduce the sparse grid reconstruction of the density, denoted   ℎ , and the electric field, denoted E ℎ , as: where ∇ ℎ and ∆ ℎ are the discrete gradient and discrete laplacian operators defined by finite differences.Let us introduce the corresponding scheme in Algorithm 2, that we shall name hybrid grid PIC scheme (owing to the Cartesian grid and sub-grids considered within the scheme).
Remark 3.5.The positivity of the charge density is not preserved by the sparse grid interpolation (see Rem. 3.3).
Proposition 3.6.Assuming  ∈  2 (Ω), then the following estimation holds: where  and  are constants related to the derivatives of  depending on the dimension of the problem.For two dimensional computations, it holds: ) The benefit of the sparse grid combination technique as a noise reduction strategy is pointed out by the equation (3.13).Indeed, owing to the projection of the density onto the sub-grids, the particle sampling error is significantly reduced: , where  ,ℎ (respectively  ,ℎ l ) is the particle sampling error of the density projection onto the Cartesian grid (respectively a sub-grid).The profit is significant for refined grids as well as three dimensional problems.
Conversely, an increase of the grid-based error results from the combination; the grid-based error is increased from , the loss may appear negligible, however the dominant component of the density error depends on 2th order derivatives.As a comparison, the grid-based error of regular PIC methods is dominated by component with dependencies on 2nd order derivatives.This feature indicates a drawback of the combination technique and may limit the efficiency of the method when the solution develops strong gradients not aligned with the Cartesian grid.A similar estimation can be stated for the electric field.Proposition 3.7.Assuming enough smoothness on the solutions, i.e.Φ ∈  5 0 (Ω),  ∈  2 (Ω) ∩  3 (Ω), then the reconstruction of the electric field verifies: where  and  are constants related to the derivatives of E and depends on the dimension of the problem.For two dimensional computations, the following estimation holds true: ) Proof of Proposition 3.4.Recasting the sums and integral, it holds true: where the two following relations have been used: Proof of Propositions 3.6 and 3.7.See Appendix A.

Discretization on sub-grids (PIC-Sg)
In this section, a second application of the sparse grid combination technique to PIC methods [29] is presented.This implementation does not use a projection of the density onto the Cartesian grid.The idea is to deposit the charge density, solve the Poisson equation, differentiate the electric potential on each sub-grid and eventually interpolate the electric field from the sub-grids at particle positions by means of the combination technique.Let us introduce the sparse grid reconstruction of the electric field, denoted E  ℎ , defined by the relation: where E ℎ l is the approximation of the electrid field on a sub-grid with finite differences.Solving the Poisson problem on each sub-grid rather than on the Cartesian grid, is likely to speed-up the computations.Indeed, the gain may be coarsely estimated by the reduced number of cells in all the sub-grids compared to the Cartesian mesh (see Eq. (3.10)), the result being a linear system with reduced size to solve for this scheme.Let us introduce the scheme in Algorithm 3, that we shall name the sub-grid PIC scheme in the following.
d d where The estimation of the error requires stronger assumptions on the smoothness of the electric potential (Φ ∈  4 0 (Ω)∩ 5 0 (Ω)) and density ( ∈  4 (Ω)∩ 5 (Ω)) than for the PIC-Hg scheme (Φ ∈  5 0 (Ω),  ∈  2 (Ω)∩ 3 (Ω)), especially on mixed derivatives.The significant differences with the PIC-Hg scheme are the additional dominant terms depending on higher mixed derivatives of the solution: Φ‖ ∞ which are (4 − 1)th order derivatives of the density and (4 + 1)th order derivatives of the potential.In comparison, the dominant terms in the PIC-Hg scheme estimation in equation (3.15) are (2 + 1)th order for the density, and the negligible terms are fifth order for the potential.
In order to prove the Proposition 3.8, we need to introduce some notations and lemma.Let us introduce ∇ ℎ l and ∆ ℎ l the discrete second order finite difference operators defined on the component grid Ω ℎ l by: are left-sided, right-sided differences and  0 ,ℎ   is centered difference: and i  ∈ N 2 is the index whose value is 1 along the th coordinate and 0 elsewhere.Let us explicit the notation { 1 , . . .,   } ∈ I := {∅, {1}, {2}, {1, 2}}, for 0 ≤  ≤ 2, standing for (3.25) The quantities (functions, domains, operators, etc.) associated with this notation correspond either to continuous ones for  = 0, semi-discrete ones for  = 1 or discrete ones for  = 2.Then, we introduce the semi-discrete operators defined by: ∆ where  is defined on the following hyperplane: Ω (1,...,) for  ≥ 1 and Ω ∅ ℎ l := Ω.In order to prove the Proposition 3.8, we need the following lemma: of the interpolation onto  ℎ l is:  with: , where { 1 , . . .,   } ∈ I and 0 ≤  ≤ 2, verifying the semi-discrete problem: then the following majoration holds: Proof of Lemmas 3.10-3.12.See Appendix A.
Proof of Proposition 3.8.The proof is an outcome of Proposition 3.2, we therefore need that the local error (E  ℎ − E) verify the assumption stated by equation (3.8).Unlike the problem discretized on the Cartesian grid (see proof of Proposition 3.7), the resolution of the Poisson equation on the sub-grids requires more work because of the anisotropy of the grids, the operator of the problem depending on the discretizations in each direction, which are different.Thus, invertibility of the operator cannot provide directly an error expansion of the form of equations (3.8).This problem has been thoroughly studied in [7,28].We therefore use the framework introduced in [28] for this proof.
First the local error can be recast into: (3.34) In the following of this proof, we ommit the dependance upon the grid discretization in the coefficients of the form  1,..., (︀ •; ℎ  1 , . . ., ℎ   )︀ for simplicity of notation.Applying Lemma 3.11 for  = 0 and introducing the estimator of the density, one gets: Let us introduce the semi-discrete problems: The equation (3.35) can be recast into: where the Lemma 3.11 has been used on the right hand side and owing to the Lemmas 3.11 and 3.12, it holds for  = 1, 2:

.37)
Eventually applying the operators ∆ −1 ℎ l , ∇ ℎ l which commute, Lemma 3.12 and a Taylor expansion: where ) ) On the other side, a Taylor expansion gives us the truncation error of the discrete gradient: (3.43) Lemma 3.10, applied first for  = 1 and  = 0, gives us: where where Finally applying Proposition 3.2 and a argument similar to the proof of Proposition 3.6 for the particle sampling error, we get the result.
Proof of Proposition 3.9.The total momentum of the system is defined as where E  ℎ is the sparse grid reconstruction of the electric field evaluated at particle positions, given by: . (3.54) Exchanging the sums, we get:

𝑙2
)︃ • A change of index together with the periodicity of the problem yields: Thus, owing to equivalent relations in  2 , we get the result.
Remark 3.13.The proof on the conservation of the total momentum cannot be applied to the PIC-Hg scheme.Indeed, the source term of the Poisson equation does not appear in the field expression in equation (3.55).

Offset combination technique
From the analysis conducted in the precedent sections, the following conclusions may be stated.Sparse grid reconstructions define a different trade off between the two components of the errors (grid-based and particle sampling errors) as compared to standard PIC approximations.The particle sampling error is significantly mitigated thanks to sparse grid approximations.This is an important feature since this component of the error is generally the most detrimental for the computations and ultimately limits the precision of the approximation.Contrariwise, the grid-based error is less favorable for the sparse grid approximations: the dominant term depends on mixed derivatives of the solutions (see Tabs. 1a and 1b for a comparison of the methods).In this section, a general framework, that we shall referred to as offset combination technique, is introduced in order to reduce the grid based error of sparse grid reconstructions.The offset combination technique is motivated by the property of the dominant term error to be an increasing function of the number of sub-grids involved in the combination ( ) as well as the combined sum of the errors of the partial estimators.The offset combination consists therefore in both reducing the number of sub-grids considered within the combination and using sub-grids with increased minimum levels.In this respect, the offset combination borrows some ideas to the so-called truncated combination [2,24].However, a more subtle strategy is implemented within the offset combination in order to select efficiently the subset of sub-grids.
The main drawback of the truncated method is the additional statistical error introduced in the simulation, because of the smaller cells of the grids considered in the combination, the mean number of particles per cell is reduced.The offset combination is aimed at mitigating the increase of the statistical noise as well as the dominant component of the grid based error in return of a deterioration of the error component depending on the non-mixed derivatives of the solution.The tuning of the balance between the different components of the error is implemented thanks to the two parameters  0 ,  1 ∈ N. The index  0 , which is the truncation parameter, is used to parametrize the minimum discretization level for the sub-grids, with the aim of discarding the most anisotropic sub-grids from the combination.The parameter  1 , which is the offset parameter, sets the loss of discretization in the directions aligned with an axis (contributing to the negligible terms of the grid-based error with non-mixed derivatives) as illustrated on Figure 3.For  1 = ( − 1) 0 the offset method boils down to the truncated combination technique introduced in [2,24].The set of sub-grids for the offset combination is defined by: where the constants are the same constants as in Proposition 3.6 and: Compared to the estimation of the PIC-Hg scheme (Eq.(3.13)) both the grid-based component and the particle sampling component of the error are mitigated in the offset method.This is obtained thanks to the reduced number of sub-grids in the combination which scales with Besides, the use of sub-grids with higher levels, i.e. parametrized by indices l with larger |l| 1 , entails an offset for both the particle sampling error and the grid-based error, if  2 ≥ .Conversely, the negligible components in the grid-based error are increased by the elimination of the more anisotropic grids, however this is of no consequence on the accuracy of the numerical method since these components are negligible.
Remark 3.15.The offset combination technique can be applied to either the PIC-Hg scheme or the PIC-Sg scheme without the loss of their conservativity properties (total charge, total momentum, etc.).

Objectives
The analyses conducted within Section 3.2 highlight a deterioration of the electric field approximations carried out by sparse grid PIC methods (see Tab. 1b) compared to standard PIC methods.Similarly to the density interpolant, this altered precision, more important for the PIC-SG scheme, is due to an increase of the grid based error component depending on high order cross derivatives of the solution.The corrections proposed herein address this specific issue.

Oversampled hybrid grid (PIC-OHg)
The PIC-Hg scheme does not take advantage of sparse-grid techniques for the resolution of the electric field, the potential being carrying out on a Cartesian mesh, unlike the PIC-Sg scheme, which makes the resolution more expansive than the latter.This computation may be expansive for refined discretizations and particularly for three dimensional problems.The benefit of this approach is a reduced dependency of the electric field approximation to the solution cross derivatives (compared to the PIC-Sg scheme, see Tab. 1b).The strategy introduced to alleviate the numerical cost of the electric field computation, consists in using grids with different resolutions for the charge density deposition and the electric field computation.Precisely, the electric field is carried out on a full mesh Ω (∞) ℎ while the density is projected onto a sequence of sub-grids associated to a more refined (oversampled) full mesh Ω (∞) ℎñ , where ℎ ñ = ℎ  • ℎ Δ , ∆ ∈ N. The corrected scheme, that we shall name the oversampled hybrid grid PIC scheme, is similar to the PIC-Hg scheme, except that the sub-grids are considered in the following index set: where   ,  ñ correspond to the volume of a cell of the grid considered: This correction consists in enhancing the sub-grids for the resolution of the electric field by introducing subgrids more refined than those used for the projection of the density.The sub-grids carrying the electric field are considered with discretization ℎ l = ℎ l • ℎ Δ , ∆ ∈ N being a parameter denoting the additional depth of these enhanced sub-grids.The corrected scheme, that we shall name enhanced sub-grid PIC scheme is similar to the PIC-Sg scheme except that after the projection, on each sub-grid, the partial representation of the density is interpolated to the enhanced sub-grid Ω ℎl with standard interpolation: ρℎ˜l (iℎ l) := ∑︁ j∈ l  l,j  ℎ l ,j (iℎ l), for i ∈  l, (3.64)where the coefficients  l,j are determined by interpolation conditions.Eventually, the electric field is computated on the enhanced sub-grids and reconstructed on the non-enhanced sub-grids in a way similar to equation (3.61).

General settings
In the following, we consider the electrons immersed in a uniform, immobile, background of ions: A series of numerical tests is carried out in two dimensional geometries: a Landau damping in both the linear and the non-linear regimes as well as a diocotron instability.These tests are performed with the standard Particle-In-Cell scheme (PIC-Std), the sub-grid Particle-In-cell scheme (PIC-Sg), the hybrid grid Particle-In-cell scheme (PIC-Hg), the enhanced sub-grid Particle-In-Cell (PIC-ESg) scheme and the oversampled hybrid grid Particle-In-Cell scheme (PIC-OHg), respectively presented in Section 2.2, 3.2.1,3.2.2,3.4.3,3.4.2.The schemes are implemented either with the classical, when not specified, or the offset combination technique.The results are compared without any filtering methods for both sparse and standard schemes.Throughout this section we will refer to the mean number of particle per cell, denoted   , relating the amount of statistical noise in the simulation and introduced in [29].This quantity depends on the underlying Cartesian grid or the sub-grids used in the combination.
where the grid discretization now depends on the domain size (ℎ  = 2 − ) and  1 ,  2 are defined in equations (3.59).The mean number of particles per cell is therefore provided by equations (4.1)-(4.3)for the PIC-Std scheme, the PIC-Sg, PIC-Hg, PIC-ESg, PIC-OHg schemes with the classical combination technique and the offset combination technique.
In any of the following test cases, the total charge of the density is exactly conserved (to machine precision ≈ 10 −16 ) for the PIC-Std, PIC-Hg scheme and PIC-OHg schemes with any combination technique which confirms the Propositions 3.4 and 3.16.Recalling that within the PIC-Sg or the PIC-ESg scheme, because the density is never fully projected onto a grid since we only proceed a partial projection onto each sub-grid and combine the resulting partial electric fields, the conservation of density is therefore not assessed numerically for these schemes.

Landau damping
As a first test case, we consider the evolution in time of a perturbation known as the Landau damping, in the linear and non-linear regimes.When a plasma is slightly perturbed from an equilibrum state, it returns to its equilibrium with an exponential damping.A perturbation in the electron distribution of an equilibrium state is considered: where ‖v‖ is the magnitude and  1,2 is the period of the perturbation in each dimension.

Linear regime
The perturbation considered in the distribution of electrons has to be small enough so that a linear approximation is valid.Under this assumption, the electric field decreases exponentially fast in time according to a damping rate [23].The motivation here is to recover this damping rate with the different schemes.Let us parametrize the perturbation with and the final time  = 25 −1  .The grid discretization is chosen so that ℎ  ≃ 0.34  , the configuration of the grid and particles is indicated in Table 2. Interpolations in the combination technique are done with linear splines for all the methods.The numerical results are represented in Figure 4a.Both the PIC-Hg, PIC-Sg, PIC-Std schemes agrees well with the damping rate.Besides, there are four times less total particles in the simulation for the sparse grid schemes.

Non-linear regime
When the perturbation of the equilibrium state is considered large enough to invalidate the linear approximation, the precedent analytic damping rate is not available anymore.In order to assess the efficiency of the methods, the results will be compared to the PIC-Std scheme in the same configuration of grid discretization.A larger perturbation than for the linear case is considered with  1 = 0.2,  2 = 0.15,  1 = 4,  2 = 3 in equation (4.4) and let  = 60  , ∆ = 1 20  −1  .The grid discretization is chosen so that ℎ  ≃ 0.47  , the configuration of the grid, particles and schemes is indicated in Table 3. Interpolations in the combination technique are done with linear splines for each of the methods.
A first series of simulations is performed with different grid resolutions, numbers of particles per cell and total number of particles in order to get a comparison of the projection error between the methods at initial time.It is therefore possible to assess precisely the precision of the density projected onto the grid for the different methods by comparison with the analytic expression of the initial electron density ( ex ).The error of the density in  2 or supremum norm is defined as: where  ex is the analytic density at initial time.The integrals in the  2 norm expressions are approximated with a Riemann sum on the Cartesian grid.The numerical results are represented in Figure 5.The results for the PIC-Sg scheme are not represented because at initial time the scheme is equivalent to the PIC-Hg scheme.
On the left figure, for any of the methods, the precision is limited by the particle sampling error for quite coarse grid resolutions, as soon as the number of cells is larger than 2 5 in each direction.Though the total number of particles is increased with the mesh refinement, the mean number of particles per cell remains unchanged, which explains the non decreasing error observed on the plots of Figure 5 (proportional to 1 √  for the standard method).This highlights that, though the number of particles per cell is consequent (five hundred or one thousand), it should be increased further to obtain an optimal precision and reduce the particle sampling error to a value comparable to the grid based error.This outlines an important characteristic of PIC methods: the grid based error is marginal compared to the particle error.This proves that gains may be expected from numerical methods with a better control of the statistical noise, hence the interest for sparse grid reconstructions.On the right figure the grid discretization is frozen so that the grid based error is constant, proportional to ℎ 2  for all the methods which testify that the multiplicative term | log ℎ  | −1 in the error expression of the sparse grid methods is negligible.When the number of particle increases, the particle sampling error converges (in the  2 norm) to zero at the rate  − 1 2 .The global error converges to the grid-based error.For refined grids with more than 2 6 cells in each direction, the amount of statistical noise of the PIC-Hg scheme with one hundred particles per cell is equal to that of the PIC-Std scheme run with five hundred particles per cell.This amounts to a total number of particles 32 times less for the sparse method compared to the standard method.
The second series of simulation for the non-linear Landau damping is dedicated to the time evolution of the perturbation.First, the conservation of the total momentum is investigated for the different methods.To this end, let us introduce the following error for the momentum:    where v  () is the velocity of the th particle at time  and  ℎ := √︀ 2    /  is the thermal velocity of the electrons.The default of momentum conservation is represented as a function of time in Figures 6a and 6b with   = 500.As predicted by the analysis conducted before, the total momentum is exactly conserved (to machine precision) for the PIC-Sg with any of the classical combination technique or the offset combination technique, as well as the PIC-Std scheme.The conservation default for the PIC-Hg, PIC-OHg and PIC-ESg schemes is observed to remain marginal (≈10 −6 ) and bounded independently of time.
The projection onto the Cartesian grid and a section in the -direction of the electron density are proposed in Figures 7 and 8 after two periods of oscillation of the electric field (at time  = 6.3 −1  ) for the different configurations of Table 3.Since the density is never projected onto the Cartesian grid within the PIC-Sg and PIC-ESg schemes, we perform an interpolation on the Cartesian grid according to the combination technique for the diagnostics.
First, it appears that the reduction of the numerical noise is manifest for all the methods using sparse grid reconstructions.It is an essential property since, this error, due to the undersampling of the distribution function, is the most detrimental in the precision of numerical methods.This better control of the numerical noise is obvious on the density plots displayed on Figures 7 and 8. Though the estimated mean number of particles per cell is equal to 1000 for the PIC-Hg, PIC-Sg, PIC-OHg (offset combination) and PIC-ESg (offset), the magnitude of the dispersion is observed to be comparable, or even less, to that of the PIC-Std scheme with 4000 particles per cell which amounts to a total number of particles 25, or 36 times greater (see Tab. 3).It is also noticeable that sparse grid approximations do not introduce too much numerical diffusion, the extrema of the numerical approximations being comparable whatever the scheme.Second, the grid-based error can be observed, particularly on the plots of the PIC-Sg and PIC-Hg methods Figure 8.This error reproduces the patterns of the coarsest grid levels.This error is specific to sparse grid approximations of the density and analysed to be dominated by the component scaling with  (︀ ℎ 2

𝑛
)︀ .This latter term results from the accumulation of error on the different levels of sub-grids and is reduced thanks to the offset combination technique.The smearing of the grid structure on the error plots related to the computations performed with the offset combination technique can be observed on Figure 8. Indeed, in these computations the number of sub-grids considered for the reconstruction is reduced considering  1 <  (from 2 − 1 = 13 sub-grids to 2 1 − 1 = 7 sub-grids).The control of the numerical noise (with approximately the same number of total particles), though expected deteriorated, due to the use of more refined grids compared to the original methods, remains very effective improving significantly the quality of the numerical results.For these methods, still a good control of the numerical diffusion shall be pointed out.Similarly, the increase of the negligible grid-based component (see Prop. 3.14) reveals to remain marginal on the precision of the numerical approximation.

Diocotron instability
In this test case, we consider a hollow profile in the electron distribution, confined by a uniform magnetic field B [25], with the following Maxwellian distribution of electrons: ,  s.t.

∫︁ ∫︁
The external magnetic field is considered uniform along the -axis B = (0, 0,   ) (  = 2.5 × 10 −5 T) and strong enough so that the electron dynamics is dominated by advection in the self-consistent field E × B. The instability caused by the magnetic field deforms the initially axisymmetric electron density distribution, leading, in the nonlinear phase, to the formation of a discrete number of vortices.As we expect fine scale structure to form in the process, a high discretization in space is required to reproduce these structures.Again, this test case defines a demanding benchmark for sparse grid approximations very likely to outline the grid error introduced with these reconstructions.
Let the parameters be  = 22  , ∆ = 0.1 −1  , the system is observed at time  = 54 −1  .The numericals results will be compared to the PIC-Std scheme with a grid composed of 256 × 256 cells and   = 200 particles per cells as indicated in Table 4.The interpolations are implemented using linear splines while splines of degree two are used for the density visualization.Following steps of proof of Proposition 3.6, one can see that the use of spline of degree two, even restricted to the visualization of the density on the grid, provides a better representation of the density.
The projection of the electron density onto the Cartesian grid and a section in the -direction (at  = 2 3 ) at time  for the different configurations of Table 4 are represented in Figures 9 and 10.The numerical results are presented only for the PIC-OHg scheme and not for the PIC-Hg scheme because the resolution of the Poisson problem for the latter is exceedingly costly in comparison to the others methods.Besides, the same accuracy is achieved with either the PIC-OHg scheme or the HG scheme, hilighting the benefit of the oversampled correction of the scheme.The discretizations of the sparse grid schemes with the classical combination technique (ℎ  = 2 −10 , ℎ ñ = 2 −10 ) are chosen higher than the standard ones (ℎ  = 2 −8 ) because the sparse grid schemes fail to reproduce the fine-scale structures depending on cross derivative terms.Indeed, we have shown in Section 3

𝑛
)︀ for the standard scheme).Despite the use of a more refined grid discretization, we observe on the plots of the section in the -direction (Fig. 10) that the PIC-Ohg and PIC-Sg schemes with the classical combination technique still fail to reproduce correctly the fine-scale structure of the density.Indeed, where the solution has steep gradients ( ≈ 5.5,  ≈ 9 in Fig. 10) the sparse grid schemes show a bit of numerical diffusion.The correction of the schemes with the offset combination technique achieves a fair representation of the density and even reproduce some fine-scale structures that are faded by the statistical noise of the PIC-Std scheme with   = 40 (see the zooms on Fig. 10).The improvement is manifest on the plot of the section in the -direction (Fig. 10) where the numerical diffusion introduced by the sparse grid schemes has been mitigated.Here again the sparse grid techniques achieve an improvement on the statistical noise with the same mean number of particles per cell (and thus less total particles in the simulation) than the regular PIC approximation (see the plot of the section in the -direction for the density in Fig. 10).Though this test case is a very demanding benchmark for sparse grid approximations a reduction of the total number of particles is achieved (see Tab. 4) and all the fine-scale structures appearing in the reference solution are well reproduced by the corrected schemes.

Conclusions
In this paper, we have presented, analysed and proposed Particle-In-Cell numerical methods embedding sparse-grid reconstructions by means of the combination technique.These methods have been numerically experienced and compared against regular PIC methods.Sparse-grid PIC approaches offer a reduction of the memory cost of the method (see Fig. 10) thanks to a better control of the statistical noise which entails a decrease of the particle number.A reduction of cost for the computation of the electric field is also accessible, owing to the reduced number of cells composing sparse-grids (

𝑛
)︀ for a Cartesian mesh).The analyses conducted within this document show that the approximation error may be decomposed into a particle error and a grid based error.Two components define the grid based error, one depending on mixed derivatives and the other depending on non-mixed derivatives of the solution, with a balance between these two contributions depending on the numerical methods.The particle error is related to the sampling of the distribution function by particles and characterizes the dispersion of the sampling.One can conclude that the Sub-grid (PIC-Sg) and Hybrid-grid (PIC-Hg) schemes achieve both a fair representation of the density with an improved control of the statistical noise compared to regular PIC schemes.Nonetheless, the grid based error is deteriorated, due to the increase of the component depending on the solution high order cross derivatives.This increase is more substantial for three dimensional computations.Similar conclusions may be drawn from the formal analysis stating the first convergence (and rate of convergence) results for the electric field sparse grid approximation established within this document (see Props.3.7 and 3.8 for the PIC-Sg and PIC-Hg schemes).Furthermore, the PIC-Sg scheme is proved to be compliant with the conservation properties of standard PIC methods (total charge and momentum), except positivity.The offset technique is introduced to decrease the dominant component of the grid based error in sparse grid PIC methods.This framework permits to tune the balance between the different components of the error and improve the quality of PIC sparse grid approximations.
-Particle sampling error: the following estimate bound holds: Let us look at the covariance of a pair of sub-grids.Since any of the random variables   ℎ l  ,ℎ l are centered, equations (2.23), (3.5) and the Cauchy-Schwarz inequality gives: where, owing to the linearity of the expected value: From equations (2.22) and (2.23), the linearity of the expected value, then recasting the sum upon (, ) ∈ 1,  2 into two sums, whether if  is equal to  or not, which contain  and  ( − 1) terms, and by independancy of the random variables of the sample in the latter case, we have: Finally, owing to Lemma A.1 we get:  Note that the -terms with odd exponent no longer vanish because Ψ is not symmetrical.We have: Eventually, using the partition of unit property of the basis function of equation (3.6)Then adding the two relations, dividing by ℎ 2   and summing upon the dimensions , one gets the result.

Figure 2 .
Figure 2. Illustration of two dimensional sparse grids used in the combination approximation.

.62) Remark 3 . 16 . 3 .
The total charge of the density is conserved by the projection onto Ω Enhanced sub-grids (PIC-ESg) The domain is a square Ω := [0, ]  , whose dimensions depend on the Debye length  ∝   ,   = √︀  0   /   0 , with the following charge, mass and temperature for the electrons   = 1.602 × 10 −19 C,   = 9.109 × 10 −31 kg,   = 1 eV.Periodic boundary conditions are considered for the particles, the electric potential and the electric field.The time discretization depends on the plasma frequency: , ∆ ∝  −1  ,   = √︀    0 /   0 .The particles are pushed with a leap-frog temporal scheme following the equations of motion.

Figure 4 .
Figure 4. Linear Landau damping: evolution of the electric field.

Figure 7 .
Figure 7. Non linear Landau damping: density profile of the electron density.

Figure 10 .
Figure 10.Number of particles and memory storage.
.11)The particle distribution   is represented by a collection of  numerical particles whose positions and velocities are denoted (x , v  ), for  ∈ 1,  .The standard PIC scheme consists of four steps repeated at each iteration in time (see Algorithm 1).Algorithm 1. PIC-Std scheme.Require: Particle positions and velocities (x, v), time step Δ, external magnetic field B. for each Δ do Accumulate the charge density  onto Ω (∞) ℎ (see details in the following).Compute E from  on Ω (∞) the notation for the binomial coefficient.

Table 1 .
Grid-based error: dominant and marginal terms with dependances on  or Φ derivatives.

Table 4 .
Diocotron instability: Configuration of the numerical methods.