Verification of 2D × 2D and two-species Vlasov-Poisson solvers

In [2], 1D ×1D two-species Vlasov-Poisson simulations are performed by the semi-Lagrangian method. Thanks to a classical first order dispersion analysis, we are able to check the validity of their simulations; the extension to second order is performed and shown to be relevant for explaining further details. In order to validate multi-dimensional effects, we propose a 2D × 2D single species test problem that has true 2D effects coming from the sole second order dispersion analysis.Finally, we perform, in the same code, full 2D × 2D non linear two-species simulations with mass ratio 0:01, and consider the mixing of semi-Lagrangian and Particle-in-Cell methods.


Introduction
We consider the two-species Vlasov-Poisson system in 2D × 2D. We look for ion and electron distribution functions f s = f s (t, x, v), with s ∈ {e, i} and electric field E = E(t, x), satisfying (1) and subject to initial distributions f i (t = 0, x, v) and f e (t = 0, x, v). Here q is the charge, m s is the mass of the species s and ε 0 is the dielectric constant, t ∈ R + is the time, x = (x 1 , x 2 ) ∈ Ω x = R 2 /(L 1 Z × L 2 Z) the position and v ∈ R 2 the velocity. Φ = Φ(t, x) is the electric potential. The original aim of the PICSL Cemracs project is to develop a code that works both for Particle-in-Cell (PIC) and semi-Lagrangian (SL) method and that is able to solve equation (1). We focus here on some of the difficulties of kinetic simulations that are the multi-dimensionality (here 2D × 2D instead of 1D × 1D), multi-species (ions and electrons) and multi-methods (both PIC and SL) aspects. Extensions to higher dimensions, Vlasov-Maxwell (see [26] for such a recent work, that discusses also the pros/cons between PIC and semi-Lagrangian methods) and gyrokinetics (that includes the issue of using more complex geometries) is out of the scope of this paper and will be the subject of further research.
In the literature, works on single species 1D × 1D Vlasov-Poisson solvers is abundant. We refer here to [23,21,14,12,13,10,3,2,1,15,19,24] for works on 2D × 2D Vlasov-Poisson simulations and to [18,17] on multi-species simulations. This list is far from being exhaustive; there is a huge number of papers in plasma physics on the subject. Validation with respect to the dispersion relation is often performed for 1D × 1D simulations (see for example [22], where cross code comparison is also performed). The dispersion relation analysis, even if less used, permits also to study multi-dimensional and multi-species simulations. Using such an analysis, our first aim is to justify the two-species simulations of [18], following [25], and consisting in the linearization of the equations around the Maxwellian equilibrium. Note that the dispersion analysis dates back to Landau [16]. With respect to the usual single species case, two main modes play here a role and their relative weights have an importance, in order to catch the right behavior. A finer study permits to exhibit relevant nonlinear effects, thanks to a second order expansion. We then had in mind to extend the results to the multi-dimensional case. We focus there first on the single species case and are able to show a true multi-dimensional effect solely visible from a second order expansion. Note that in the literature, usual 2D ×2D test cases reduce to 1D ×1D (see previous references based on Landau and two-stream instabilities test cases; note that in [19], an effort has been put to get a two dimensional character). We then could study the two-species case, in the 2D × 2D setting, but such analysis is basically the superposition of the two previous analysis and we prefer here to focus on performing nonlinear simulations, where the analysis coming from the dispersion relation is anyway no more valid.
If a subsequent part of the work relies on the verification of the codes through the dispersion relation, another part of the work is based on the comparison and the mixing of semi-Lagrangian and PIC methods. As the two methodologies are quite different, getting the same answers for both codes is a further validation. On the other hand, defining a common code that works for both of them is a step to allow for a wider application range, as each method has his own benefits and drawbacks. We could for example consider in the future, the PIC method when a species is well localized (which permits to prevent from the discretization of the whole phase space) and the semi-Lagrangian method for other species. Note that the coupling of the two methods seems not to have been considered; at least, such approach is rare, as generally one method is privileged for a given simulation. The difficulties rely here on the definition of a common framework and on the need of expertise in both methods. In the numerical results, we will see that the present approach does not generate extra problems due to the coupling of the methods, which is encouraging for further developments.
In the sequel, we will first work on getting analytical results for our system, through the study of the dispersion relation. The analysis will be upgraded to second order expansion as done in [11], for example. This will enable us to propose the true 2D × 2D one species test case and also to study the two-species (1D × 1D) test case proposed by Badsi-Herda [18]. We will then present the numerical method and implementation details that permit to deal with both methods, before giving the numerical results, which are in accordance with the analytical results. We then present full non linear two-species results in 2D × 2D (an extension to two-species of a test case presented in [19]). Finally we give some details about the efficiency of the code in the context of parallel programming on supercomputers.

Description of the equations and test cases
In this section, we describe the test cases and motivate our work.

A 1D × 1D two-species test case
The first test case has been studied by [18]. We look for f i , f e satisfying with ε = me mi , the root of the mass ratio between ions and electrons. ( The initial functions are given by with k = 2π L , and A the amplitude of the perturbation. The phase-space domain is We will take here σ = 1 2 and L = 21, as in [18]. This is a first example of two-species simulation. Our goal is to reproduce these results [18] from the literature with our commonly used methods (as in [20] for example) and also to provide a dispersion analysis, which permits to further validate the code. Note that this test is 1D × 1D, but it will be simulated in the 2D × 2D code; this enables to have a first check of the code.

A 2D × 2D one-species test case
We focus then on 2D × 2D phase space. We look for f satisfying with initial function We take L 1 = L 2 = 4π, the phase-space domain is [0, 4π] 2 × [−v max , v max ] 2 and v max = 10. This test permits to capture the interaction of different modes and reveals 2D space features, which would not be visible by 1D × 1D codes. The authors are not aware of such a test case in the literature; it seems not to be standard.

A 2D × 2D two-species test case
We look for f i , f e satisfying with ε = me mi , the root of the mass ratio between ions and electrons.
The initial functions are given by with k 1 = 2π L1 , k 2 = 2π L2 , the perturbation amplitudes A 1 , A 2 , the velocity drift v d and the thermal velocities Our aim is to develop 2D × 2D two-species simulations; so, here is such an example. It is a generalization of the first test to the 2D × 2D framework; we use a 2D × 2D initial function for the ions that was developed in [19]. We will take here v d = 2.4, A 1 = 0.005, A 2 = 0.25, σ 1 = 0.5, σ 2 = 1, k 1 = k 2 = 0.2 together with v max = 10.

Dispersion analysis
In this section, we perform the dispersion analysis for the first and second test cases (Subsection 2.1 and 2.2). The dispersion analysis consists in studying a perturbation of an equilibrium solution, which leads formally to a solution of the linearized problem. For the first test case, we obtain in particular the behavior in the mass ratio limit ε → 0. A finer analysis using second order expansion permits to explain an instability that cannot be captured through first order expansion. Such an expansion will also be a key in the development of the second test case (Subsection 2.2).

First order expansion
We first recall the first order expansion with respect to A, in the simplest case and then look how to adapt it to multi-dimensional and multi-species cases. The procedure is quite standard. We give here the general form and explicit the computations for Maxwellian type initial functions.

The one dimensional case with one species.
We first recall the 1D × 1D dispersion analysis (see [25] for details). Taking an initial condition of the form of the Vlasov-Poisson equation neglecting O(A 2 ) terms. We then get an equation for f 1 which reads Considering f 1 (t = 0, x, v) = f (t = 0, v)e ikx with k ∈ 2π L Z * , we get a solution of (8) which is with E 1 = E(t)e ikx . In order to express E, we introduce D k and N k defined by and get thanks to a Laplace transform and the theorem of residuals Res which leads to an expression of E whose L 2 norm will be computed as a diagnostic in the numerical results.

3.1.2
The one dimensional case with m species.
with general numbers a j , b j , z j . We get similarly the solution of the linearized problem which writes j dv, and are given by the expressions Res ω e −iωt , with , We then can further explicit these computations, for Maxwellian type functions, thanks to the plasma dispersion function of Fried and Conte (see [25]) In (11), we typically can have terms like We have from Appendix B of [6] Z with the derivatives The first terms are Z 0 = Z and In fact, we can check that we have the following formula: In particular, for the first two-species test case (described in Subsection 2.1), we define We get This leads to In the sequel, we will write D k (ω, ε), instead of D k (ω) to specify the dependence of this quantity to ε, through ξ e = ωε √ 2k .

Asymptotic behavior for the two-species case
In the previous example, we remark that, in practice, it may be difficult to determinate the zeros of D k (·, ε) when ε is small. As a consequence, we focus here on the asymptotic behavior of the zeros of D k (·, ε), when ε goes to 0.
In this study, we consider D k (·, ε) as a continuous family of entire functions. By the classical holomorphic function theory, there exists a family of continuous functions (ω n ) of the variable ε, indexed by a set of natural integers S , such that for all ε, {ω n (ε) | n ∈ S } is the set of the zeros of D k (·, ε). We can determine the asymptotic behavior of these functions as ε → 0.
considering here D k as a function of ξ e and ξ i . The latter limit is obtained thanks to the limits coming from the asymptotic expansion We have more precisely the following result, which gives the asymptotic location of the zeros of D k when ε → 0.

The multi-dimensional case.
When we consider the d-dimensional case (d ≥ 1) with one species case 1 , we take an initial condition of the form and look for a solution We introduce the component of the velocity along the mode k and its orthogonal: and {e k , e ⊥,j k , j = 1, . . . , d − 1} is an orthogonal basis of R d . Note that this analysis is already present in the original work of Landau [16]. There he considered, e k is along the x-direction. We detail here the computations, without assuming this simplification, which permits later (see Subsection 3.2) to consider several modes, and their interactions. The analog to (9) is Note that there is no evolution in v ⊥ . We now integrate the equation in v ⊥ . Introducing the notation < . > ⊥ defined by 1 a similar expression is straightforward for the multi-species case We finally can reuse the one-dimensional analysis and the first order dispersion relation rewrites Res ω e −iωt e k , with

Introduction
We focus here on second order expansion with respect to A; such expansion has been developed in [11] for example for the Landau damping. Through the previous first order analysis, we can get the superposition of 1D modes. Indeed we can consider an initial function of the form for a given set of modes K and get a linearized solution of the form Res ω e −iωt e ik·x e k .
Note that in this expression, we do not see the multi-dimensional (non-linear) interaction of these modes. Thanks to a second order expansion, we are able to study the interaction of different modes in space and design test cases more relevant to the dimension 2 (second test case, Subsection 2.2). As a by product, it will also be useful to better qualitatively describe the first 1D two-species test case of Subsection 2.1.

Results
We look for solutions of the one species Vlasov-Poisson system set for x ∈ [0, We assume that f 2 (0, x, v) = 0, and take initially Let g = S(t)u be the solution of the homogeneous system with the definition Proceeding as usual by linearization (see the previous analysis) we find the classical equation for f 1 and E 1 . Then, we obtain f 2 and E 2 as solutions of First, we deduce from (18) that we can look for the modes k = 0 for E 2 . Indeed, since f 2 dvdx = 0, which is a consequence of the hypothesis f 2 (0, x, v) = 0, the Poisson equation in (18) admits, modulo a constant, a unique solution for the electric potential, which leads to E 2 dx = 0.
Then, thanks to the dispersion analysis we can compute the source term E 1 · ∇ v f 1 above. Therefore, by the Duhamel formula we can write As f 2 (0, ·, ·) = 0, we get Inserting E 1 and f 1 which are of the form Using (17), we get where k = −k is due to the condition of k = 0 above. We look for the dependence in time of this expression.
We write The expansion ofμ k +k makes appear terms like e −iωt and . A precise result related to this expansion is rather technical and will be detailed elsewhere.

Numerical method
In order to solve the Vlasov equation, we develop a code that is able to use both Particle-in-Cell (PIC) and semi-Lagrangian methods, in the framework of the Selalib 2 library. For the Poisson solver, we classically use the FFT; time and space (semi-Lagrangian or PIC) will be further detailed thereafter. The framework is such that we can use PIC for the two-species, semi-Lagrangian for the two-species, or PIC for one species and semi-Lagrangian for the other species.

Time discretization
We consider two types of time discretizations. The first one is based on a splitting by direction, as in [18] and the second one is a splitting by species.

Splitting first by direction
The algorithm can be sketched as in Figures 1 and 2 for the classical Strang splitting. This scheme can be generalized to higher order splitting; we will here use the classical 6 th order splitting of Blanes and Moan [5], as in [7].
Parameters : ∆t, the time step. ncx × ncy, the size of the spatial grid. Advection of f electrons in x over ∆t/2 ∂tf electrons + v · ∇xf electrons = 0 4 Compute ρ from f electrons (integration in v for SL, deposit for PIC) 5 Compute E from ρ Poisson solver 6 Advection of f electrons in v over ∆t ∂tf electrons − E · ∇vf electrons = 0 7 Advection of f electrons in x over ∆t/2 ∂tf electrons + v · ∇xf electrons = 0 8 End Foreach  Advection of f electrons in x over ∆t/2 Advection of f ions in x over ∆t/2 ∂tf ions + v · ∇xf ions = 0 5 Compute ρ electrons from f electrons and ρ ions from f ions (integration in v for SL, deposit for PIC) 6 Compute E from ρ = ρ ions − ρ electrons Poisson solver 7 Advection of f electrons in v over ∆t ∂tf electrons − (1/ε)E · ∇vf electrons = 0 8 Advection of f ions in v over ∆t ∂tf ions + E · ∇vf ions = 0 9 Advection of f electrons in x over ∆t/2 ∂tf electrons + (1/ε)v · ∇xf electrons = 0 10 Advection of f ions in x over ∆t/2 ∂tf ions + v · ∇xf ions = 0 11 End Foreach Figure 2: Two-species pseudo-code, splitting by direction.

Splitting first by species
In order to have the possibility of dealing with the species differently, we developed another scheme based on a splitting by species. It is depicted in Figure 3.
Parameters : ∆t, the time step. ncx × ncy, the size of the spatial grid. ε = m electrons m ions , the square root of the mass ratio. Compute ρ species from f 9 Compute E from ρ = ρ ions − ρ electrons 10 Advection of f species in v over time step 11 Advection of f species in x over time step/2 12 Compute ρ species from f species Once more, this scheme can be generalized to higher order splitting. Nevertheless, we cannot use the splitting coefficients from the classical 6 th order splitting, because they are not suitable for a splitting by species.

Semi-Lagrangian discretization
We use a classical backward semi-Lagrangian (BSL) method, consisting here in solving successive constant advection equations on a uniform 1D periodic mesh [9]. Centered Lagrange interpolation of degree 9 is used for the interpolation; see for example [20]. Concerning the splitting by species, we use Strang splitting on each species for the corresponding solving of the Vlasov-Poisson equation.

PIC discretization
A Particle-in-Cell (PIC) method consists in discretizing (sampling) the distribution function by a collection of N macro-particles that move in the phase space following the characteristics of the Vlasov equation. We use the classical PIC method, following the lines of [4], with linear or cubic splines for the deposition of the charge and for the interpolation of the electric field. The macro-particles are initialized randomly, which ensures a stochastic convergence in 1 √ N . All the results of this paper are shown with the random initialization.
Time schemes presented in Figures 1, 2 and 3 are still valid for the PIC method. These schemes are used when running simulations using PIC for one species and BSL for the other. However, to ensure efficiency when running simulations only with the PIC method, a leap-frog scheme is used (second order in time).

Numerical results
We now give the numerical results for the test cases described in Section 2.

First test case: 1D × 1D two-species
This test case is described in Subsection 2.1. We take here ε = 1, as we first want here to validate the twospecies feature; this permits to have a first example taken from the literature [18] that is here justified with the dispersion relation analysis and that can be cheaply reproduced in this one dimensional context. On Figure 4 (left, logarithmic scale; right, standard scale), we represent, for the perturbation A = 0.0001, the electric energy defined by 1 2 L 0 |E| 2 dx vs time t and also the absolute value of the first and second Fourier modes multiplied by 1 2 , in order to be comparable to the electric energy. We represent also theoretical results, coming from the study of the dispersion analysis developed in Section 3. The theoretical first mode is here the expression It comes from the first order dispersion relation whose more precise expression, using the two first relevant zeros, is We remark that this analytical expression permits to describe precisely, up to time t = 80, the behavior of the first mode that is simulated and also the whole electric energy, as this first mode is dominant. For the simulation, we have used the BSL method on a 1024 × 2048 grid, with ∆t = 0.1. The study of the linear analysis at order 2 developed in Subsection 3.2 permits to explain the behavior of the electric energy up to time t = 90, and the behavior of the second Fourier mode from initial time to time t = 90. We have used here the following analytical expression for the second Fourier mode Here the coefficients 0.3 and 0.7 are chosen to fit the numerical results, 0.089 is an approximation of γ 1 and 0.145 is the rounding of 0.144982725814, coming from the dispersion analysis of 3.2. The theoretical electric energy is then given in the figures by E 2 1 + E 2 2 . Note that after time 100, we are in the non linear phase and the dispersion relation analysis is no more valid. On Figure 5, we take A = 0.01, as in [18]. We take here as parameters, the BSL method on a grid 128 × 256 with ∆t = 0.02. The behavior is similar. As the perturbation is bigger, the non linear phase appears sooner. We can note also that the first mode does not have time to develop and that the instability is essentially explained by the second order expansion. Then, we study the influence of the numerical parameters, on Figure 6. We see that for A = 0.0001, the grid 64 × 256 is quite good, as the difference with the refined run on a grid 1024 × 2048 (similar to 2048 × 4096) is only visible at the end of the simulation, around T = 150. For A = 0.01, we get converged results until T around 80 − 100; then for longer times, we see that the results start to differ, and the grid 64 × 256 seems not fine enough. For the time step, it seems that ∆t = 0.1 is a good choice, as the results are very similar between ∆t = 0.02 or ∆t = 0.1. On Figures 13, 14, 15, 16, we can appreciate the convergence on the diagnostics of conservation of L 1 , L 2 norms; the mass is conserved up to machine precision.
x − vx cut permits to measure the structures and the filaments, here for A = 0.01. It is confirmed that at time T = 80, the mesh 64 × 256 correctly describes the ions (see Figures 11 and 12). At time T = 150, however, as already seen on the electric energy (Figure 6 right), we see the differences between the fine run (512 × 2048 grid) and the coarse one (64 × 256) for the ions (see Figures 9 and 10); for the electrons the differences are smaller.

Second test case: 2D × 2D one-species
This test case is described in Subsection 2.2.
On Figure 18 (left, logarithmic scale; right, standard scale), we represent, for the perturbation A = 0.1, the electric energy defined by 1 2 L 0 |E| 2 dx versus time t and also the absolute value of the first and second Fourier modes multiplied by 1 2 , in order to be comparable to the electric energy. We represent also theoretical results, coming from the study of the dispersion analysis developed in Section 3. The theoretical first mode is here the expression We have here only used the theoretical values 1.416, −0.1533 and fitted the two other coefficients. We remark that this analytical expression permits to describe precisely, up to time t = 12, the behavior of the first mode that is simulated and also the whole electric energy, as this first mode is dominant. For the simulation, we have used the BSL method on a 32 × 32 × 256 × 256 grid, with ∆t = 0.1. The study of the linear analysis at order 2 developed in Subsection 3.2 permits to explain the behavior of the electric energy up to time t = 25, and the behavior of the second Fourier mode from initial time to time t = 25. We have used here the following analytical expression for the second Fourier mode E 2 = 0.0028 exp(0.259t) Here the coefficients 0.0028 is chosen to fit the numerical results, 0.259 is coming from the dispersion analysis of 3.2. The theoretical electric energy is then given in the figures by E 2 1 + E 2 2 . Note that after time 25 − 30, we are in the non linear phase and the dispersion relation analysis is no more valid. We then study the convergence on the diagnostic of the electric energy on Figures 19 and 20 (left). We notice that both PIC and BSL methods converge to the same state in the non linear phase, which permits to validate the results, from this cross comparison. We see on Figure 20 (right) the time evolution of the L 2 norm; we notice that the conservation is clearly improved by refining the grid in space. On Figures 21 and 22, we see x − vx and y − vy cuts; the first looks similar to two-stream instability and the second to Landau damping simulations. The filaments seem to be well resolved thanks to a relatively high number of points in the velocity directions. On Figure 23, we see the contour plots of ρ at different times; we clearly see the behavior of the modes: first the mode (0, 1) dominates and then it is the mode (1, 0).

Third test case: 2D × 2D two-species
This test case is described in Subsection 2.3. We take here ε = 0.01. This leads to a more oscillatory behavior. We focus here on the electric energy. On Figure 24, we give the electric energy for BSL using the 6-th order scheme, for ∆t = 0.02 and ∆t = 0.01. We remark that there is a lot of oscillations. We see that the results are very similar, which is a mark of the fact that the scheme is converged in time. We then do the comparison with other methods and numerical parameters. The same quantity is plotted for other numerical parameters on Figures 25,26,27,28,29. On Figure 25, we see that the result is equivalent with using the Strang scheme with ∆t = 0.0025. On Figure 26, we see that the convergence is not complete when passing from a grid 32 × 256 × 32 × 512 to a grid 32 × 512 × 32 × 1024, which means that high resolution in y − vy is needed. On Figure 27, we see that on the contrary, 32 points in x seem sufficient, as the curve for the 32 × 512 × 32 × 1024 and 64 × 512 × 32 × 1024 well match, and we see that going to ∆t = 0.005 in the Strang splitting case changes more the solution; so that it seems to be a little better to stick to ∆t = 0.0025. On Figure 28, we see that more clearly that high resolution in y − vy is needed: the grid 128 in y and 512 in v y is clearly not sufficient. On Figure 29, we see simulations using a splitting first by species. The time step for the ions is ∆t i = 0.1; for the electrons the time step is ∆t e = 0.01; BSL (resp. PIC) is used for the electrons on the left (resp. right) figure.
The results are converged (they are compared to a "reference" solution: BSL with 6-th order time scheme and ∆t = 0.01 on grid 32 × 512 × 32 × 2048). Thus, we validate the splitting by species using BSL for ions and BSL or PIC for electrons, with sub-steps for the electrons. This opens the door to use specific PIC (or BSL) schemes that are designed for capturing high oscillations (see [8]). On Figure 30, we compare the total energy conservation between Strang and the 6-th order splitting; we remark that the conservation is really improved with the 6-th order splitting, which is coherent with [7], where such a splitting is also used for a single species. Then, on Figures 31, 32, 33, 34, 35, we give some 2D plots. On Figure 31, we see the x − vx cut for the electrons (left) and the ions (right). We note that this picture does not change much with time; in particular, luckily, a two-stream instability is here not developed, which permit to keep a resolution small in these directions.  Figure 32, we see on the contrary, that for the y − vy cut for the electrons, very fine structures appear; this confirms the fact that high resolution is here needed. On Figure 33, we see a Landau damping behavior in for the y − vy cut for the ions. On Figure 34, we see the time evolution of ρ e = f e (x, y, v x , v y )dv x dv y , and on Figure 35, the time evolution of ρ i = f i (x, y, v x , v y )dv x dv y . We see the rapid change of ρ e with respect to time; we remark also some structures in x and the amplitude of ρ − 1 is small.

Efficiency results
We present below the time needed to achieve the simulations shown in this paper, along with strong scaling results of our code on the second test case, for both the PIC method and the BSL method. All the experiments were run on the supercomputer Curie 3 : 5,040 nodes, each node has 2 sockets, each socket being a Intel Xeon E5-2680 @2.7 GHz (SandyBridge), with 64 GB of RAM, 4 memory channels and 8 cores. The code is compiled with ifort 16.0.3.210 from Intel.
For the BSL method, we show efficiency as the number of cells updated by second. It can be computed by the following formula : ncx × ncy × ncvx × ncvy × num iteration × num split step × 2 num processes × execution time where ncx×ncy ×ncvx×ncvy are the grid sizes, num iteration is the number of time iterations, num split step is the number of steps of the time splitting (3 for the Strang splitting, 23 for the classical 6 th order splitting), and 2 because each 2D advection is split as two 1D advections.  The first test case can be simulated very fast, as it is only 1D × 1D. However, our code being designed for 2D × 2D, we took 4 points in the y and the vy directions, with uniform distribution across those directions. Execution time and efficiency (in millions of cells updated per second per processor) for each discretization are shown in Table 1.

Test case 2
For the second test case, execution time and efficiency (in millions of cells updated per second per processor) for each discretization are shown in Table 2. The scalability of our code was tested using this test case. Results for the BSL code are shown in Figure 36, for the 64 × 64 × 512 × 512 grid size (minimum grid size needed for convergence).
What can be seen in this figure, as well as in the previous tables, could be predicted. As long as there is enough computations to do for each processor, the code is efficient. As soon as the number of cells per processor is too low, the communications become a bottleneck, and the efficiency is lowered.  Results for the PIC code are shown in Figure 37, for 800 million particles with a 128 × 128 grid size (minimum particle number needed for convergence). The PIC code is made parallel using hybrid MPI + OpenMP parallelism, to use as much as possible the shared memory from the available architecture. On Curie, there are 8 cores per socket, which means that the best efficiency can be reached using 1 MPI process per socket and 8 threads per MPI process.  The left part shows the performance of the code on 1 socket, using only OpenMP parallelism. The scaling is ideal up to 4 threads and deteriorates on 8 threads. The reason bases on the reduced number of memory channels per socket, 4 precisely, and on the well-known fact that PIC codes are demanding in memory bandwidth. The scalability on 8 threads is thus limited by the number of memory channels. In all, the code processes 39 million particles per second per thread and 238 million particles per second when using 8 threads.
The right part shows that, as long as enough computations are deployed per process, the code has a very good scalability. However, when we reach 128 MPI processes, only 6 million particles are distributed per MPI process. The number of the corresponding computations is clearly not high and therefore the execution time of the simulation is dominated by the communication time between processes.

Conclusion
We have performed two-species 2D × 2D simulations with both PIC and semi-Lagrangian methods. Validation of the code is done through dispersion analysis and/or cross comparisons between the results and the numerical parameters. We notice that a high order splitting method in time is relevant, when so many time steps have to be used in this time oscillatory problem that comes from the mass anisotropy. The sixth order splitting scheme leads to better energy conservation with respect to the classical Strang splitting. The time splitting by direction does however not permit to treat differently the ions and the electrons. So, another splitting scheme by species has been introduced; the latter one permits to use different time steps for the ions and the electrons. In the future, we plan to treat the electrons differently using specific methods (several methods, already proposed in the literature in particular in the PIC case could be tried) in order to speed its computation.