Advances in the Implementation of the Exactly Energy Conserving Semi-Implicit (ECsim) Particle in Cell Method

The Energy Conserving semi-implicit method (ECsim), presented by Lapenta in 2017, is a Particle in Cell (PIC) algorithm for the simulation of plasmas. Energy conservation is achieved within a semi-implicit formulation that does not require any non-linear solver. A mass matrix is introduced to express linearly the particle-field coupling. With the mass matrix the algorithm preserves energy conservation to machine precision. The construction of the mass matrix is the central nature of the method and also the main cost of the computational cycle. We analyze here three methods that modify the construction of the mass matrix. First, we consider how the sub-cycling of the particle motion modifies the mass matrix. Second, we introduce a form of smoothing that reduces the noise while retaining exact energy conservation. Finally, we discuss an approximation of the mass matrix that transform the ECsim scheme to the implicit moment method.


I. INTRODUCTION
The Energy Conserving semi implicit method (ECSim) is an algorithm for plasma simulation based on the particle in cell (PIC) approach [19]. PIC methods can be explicit, semi-implicit or fully implicit. Plasmas are governed by two sets of equations: the equations for the motion of the particles and the equations for the evolution of the fields. The two sets of equations are coupled because the field equations need the sources (current and charge) from the particles and the particles need the fields to compute the force. As a particle moves, the fields are modified and as the fields change the forces on the particles are modified. This link is central to the physics of plasmas: plasmas are collective sets of particles interacting via the fields. The coupling between particles and fields is non-linear and to represent it in discretized equations in its fullness one needs fully implicit methods. In fully implicit methods, the particle equations of motion and the field equations are solved together within a non-linear solver, such as the Newton-Krylov approach [9,27]. In explicit methods, conversely, the coupling between particles and fields is suspended for a small time step [1]. In that small time interval one assumes that the known fields can be used unchanged for moving the particles and the particle information can be used unchanged to evolve the fields. This has three major consequences.
First, the explicit method is very simple, no iteration is needed and explicit PIC can be implemented as some of the most efficient algorithm known in computer science, consistently being a top performance achiever on any new computer architecture introduced. For example, PIC was one of the first applications to reach petascale performance [5]. Implicit PIC is much more complex in its implementation. Especially on massively parallel computers, reaching high efficiency is a challenge.
Second, in explicit PIC, the time step becomes limited by numerical stability considerations, requiring to use high resolution. The peculiarity of PIC is that the resolution needs not just be refined in time, to resolve the electron plasma frequency ω pe ∆t < 2 [14], but also in space, to avoid the so-called finite-grid instability [1]. This limitation is removed by the implicit approach that allows one to select grid spacing and time step based on the accuracy needed and not on the stability of the numerical algorithm [18].
Third, in explicit PIC, energy is not conserved. Using a good resolution, energy is acceptably maintained. There is a secular trend of energy increase always [1], but as the resolution is relaxed closer to the stability limit of the finite grid instability, the energy increase becomes more severe, until at the instability limit it starts to grow exponentially. This effect cannot be avoided but it can be improved by using smoothing and higher order interpolation techniques. Recent structure-preserving geometric particle-in-cell methods use symplectic integrators to ensure local energy conservation at small time steps [16]. The implicit PIC method, instead, conserves energy exactly, whatever resolution is used [9,27]. This feature is physically important and practically impactful. Physically, of course, confidence comes form using an algorithm that preserves one of the most established properties in physics: conservation of energy. If energy starts to spontaneously increase confidence in the results is shaken. Practically, lack of energy conservation requires a tedious and careful tuning of the parameters to make sure the simulation does not increase its energy excessively, leading in some situations to excessively resolved models that need to use much more resolution than the processes of interest require.
The semi-implicit PIC method tries to make a compromise and retain some of the advantages of both approaches. In semiimplicit methods the particles and the fields are still advanced together and an iteration is needed but the coupling is linearized and the iteration uses linear solvers. Different methods are used to linearize the coupling. The implicit moment method formulates the particle response to changes in the fields using the moment closure method [6]. The direct implicit method uses a formal linear expansion of the coupling operator [17]. In all these approaches the stability properties of the semi-implicit method are superior to explicit methods and allow a good compromise in the resolution needed [18]. However, energy is not conserved. Unlike explicit PIC, energy can either increase or decrease depending on the implementation because dissipation terms are included in the algorithm to suppress energy growth.
The ECsim approach is the first semi-implicit PIC to retain exact energy conservation as in the fully implicit PIC. ECsim uses a mathematical construct called mass matrix to express the coupling between particles and fields. With the mass matrix, the coupling is linear but energy is conserved exactly. Below we review how this was achieved [19].
Since its recent introduction, ECsim has found application in a number of applications in space [21,22,24,35,36] and in fusion [13,30]. An important improvement has removed the lack of charge conservation in the original scheme [10,31]. We report here some extensions of the method that can widen its practical applicability.
First, we describe a method to introduce smoothing to reduce noise, while retaining exact energy conservation. Smoothing is an affordable way to introduce effectively a higher order interpolation scheme. It removes the high frequency part of the spectrum. Without special attention, smoothing will tend to break energy conservation by just removing high frequency fluctuations from the system. We present, instead, a method that conserves energy.
Second, sub-cycling is a convenient approach in plasma simulation to address the faster scales seen by the particles. In some applications, the particle response (or that of a subset of the particle population) is faster than the evolution of the fields and it is beneficial to move the particles several times without advancing the fields. In explicit PIC, the operation is very simple since fields and particles are not advanced together. In implicit and semi-implicit PIC, instead, sub-cycling needs to be done also in the coupling requiring to modify the algorithm to compute the current used in the field solution.
Finally, we illustrate how the mass matrix formulation opens up the opportunity to approximate the mass matrix in certain limit cases to reduce the cost of the simulation.
The paper is organizes as follows. Section II recaps the key properties of the ECsim needed for the discussion. Section III introduces how smoothing can be implemented to reduce noise while retaining the property of energy conservation. Section 30 is dedicated to the algorithm for sub-cycling the particle motion within the ECsim scheme. Section V derives a limit case of the mass matrix formulation that transforms the ECsim algorithm into the standard implicit moment method [6]. Results are presented in Sect. VI and final conclusions are drawn in Sect. VII.

II. SUMMARY OF THE ENERGY CONSERVING SEMI-IMPLICIT (ECSIM) METHOD
The energy conserving semi-implicit method [19] is based on a formulation of the mover similar to the classic leap-frog scheme: where v p = (v n+1 p + v n p )/2 and E n+θ = θ E n+1 + (1 − θ )E n . This mover differs from both the explicit leap-frog mover [4] and the implicit θ -mover [7]: it combines the first equation from the explicit mover with the second equation from the implicit mover but with an important difference: the electric and magnetic fields are computed at the known position x n+1/2 p rather than at the unknown position x p . The important consequence is that the particle equations can be solved directly without any iteration needed among themselves. Instead in the standard θ -mover, a predictor-corrector iteration is required. Yet the mover is still implicit because the new fields are not known until the field equations are solved. The ECsim method retains the coupling between advanced fields and advanced particles, requiring the solution of a linear coupled system. However, the mover itself does not require any iteration, a substantial simplification.
In the ECsim scheme, the electric and magnetic fields are computed at the known position x n+1/2 p . In the θ -scheme, instead, they are computed at the unknown position x p . These two positions are conceptually the same, they express the particle position at the mid-time between the old and new evaluations of the velocity. But one is computed explicitly, in the leap-frog sense, while the other is computed as part of a predictor-corrector iteration [20,33]. Both methods are second order accurate but the ECsim scheme is simpler to compute. This simplicity is not just a virtue in itself but leads to an important consequence: the simplicity allows us to formulate the coupling with the fields in a way that insures exact energy conservation without requiring non-linear iterations. The θ -scheme can be made energy conserving but at the cost of a fully non-linear iteration requiring a non-linear solver [9,27].The ECsim scheme allows exact energy conservation without requiring any non-linear iteration.
The properties of stability are determined by the field-particle coupling and in this sense the method is still implicit. For this reason, not requiring any non-linear iteration but still requiring a liner solver to deal with the field-particle coupling, the method is semi-implicit. This nomenclature is to distinguish it from the fully implicit method that requires the non-linear iteration.
Note that the force term is written using the magnetic field at the initial time level B n (x n+1/2 p ) but the electric field is written at The reason for this choice is simplicity and the fact that the magnetic field does no work and using the old time level does not introduce any loss of energy conservation.
The coupling of particles and fields require to interpolate the fields to the particle positions: We used a generic index g for the grid. In the specific implementation within the iPic3D code [12,28] the electric field and magnetic field are not colocated and g label either centers (for B) or vertices (for E). Here we simplify the notation as: ). In our implementation, the interpolation function W is a b-splines of order = 1 [3]: This expression reduces trivially in 1D for the examples reported below. For the Maxwell's equation we use the standard θ -scheme [28]: The spatial operators in eq. (5) are discretized on the grid labelled by g introduced above.
The coupling of the field equations with the particles is expressed by the current for each species: where the summation is over the particles of the same species, labeled by s.
As with the θ -mover, the velocity equation can be rewritten in the equivalent form [34]: with: and the rotation matrix α n p given by: where I is the dyadic tensor (matrix with diagonal of 1) and β s = q p ∆t/2m p (independent of the particle weight and unique to a given species). The elements of the rotation matrix are indicated as α i j,n p with label i and j referring to the 3 components of the vector space (x, y, z).
Substituting then eq. (7) into eq. (6), we obtain without any approximation or linearization: where we shortened the notation W pg = W (x n+1/2 p − x g ) and the summation is intended over all particles of species s. Using eq. (8), the expression for the current becomes: where we defined: Computing then the electric field on the particles by interpolation form the grid as in eq. (3), it follows that: Exchanging the order of summation we obtain: where we have defined the actor in the leading role of the ECsim scheme: the mass matrix [8]: There are 3v (where v is the number of velocity directions) mass matrices and in matrix notation they can be written as M gg , that is without the indices i, j for the vector directions. The mass matrices M s,gg that are the most important aspect of the ECsim method and are also the most expensive part of the computation [12]. A number of symmetries can be used to reduce the cost. Speed up of the construction can be achieved using offloading to accelerator processors (e.g. graphical processing units, GPU) [2]. The mass matrices, eq. (14), provide an explicit linear link between the advanced current at the mid-point of the time step and the electric field at the advanced time. This linear relationship can be substituted into the discretized Maxwell's equations (5) to form a linear set of equations: where the total current is J g = ∑ s J sg and the species summed mass matrices, that written by elements are: The direct link provided by the mass matrix is analytically exact for the original set of discretized equations. Unlike the implicit moment method where the equations have to be approximated by Taylor series expansion [34], here the link is still exactly the same as in the original set of discretized equations. Having eliminated the need for any approximation or Taylor series expansion is the reason why ECsim conserves energy exactly.
We consider now how energy conservation can be shown in the case θ = 1/2. It is important for the derivations below to consider what key steps enable energy conservation. The inner product of the velocity equation (1) with the average speed, v p gives by summing over all particles: where the electric field is computed as average consistent with the choice θ = 1/2 and the magnetic field drops out as obvious from the properties of the cross product. Exchanging the summation over particles and cells leads to: where it is recognized that J g = ∑ p q p v p W pg .
Multiplying the first equation (5) by B g and the second by E g and summing them leads to: Assuming a mimetic grid discretization that preserves the continuum properties of the operators and summing over all grid points gives: This conservation law states that the variation of the magnetic and electric energy, as measured on the grid, equals the amount exchanged with the particles and carried by the grid-discretized divergence of the Poynting flux. For energy to be conserved in the system, the energy exchange term on the particle equations (eq. (19)) needs to be identical to that on the field equations, (eq. (21)). This term is indeed identically equal to ∑ g J g · E g in both equations. Energy conservation is enforced exactly, to round off.
Besides guaranteeing physical conservation of energy, a cornerstone in any physical model, the existence of this conservation constraint also guarantees a form of non-linear stability of the discretized equations [19] expanding the stability of the semiimplicit method compared with the moment implicit scheme [13,22].

III. SMOOTHING WITH THE MASS MATRIX FORMULATION
Smoothing can be designed to be compatible with the energy conserving properties of the mass matrix. We choose to smooth only the electric field, since the magnetic field tends to be much smoother in PIC simulations and smoothing is not needed.
From the proof of energy conservation recapped above, it is clear that for energy conservation the smoothing of the current must be done in the same way as that of the electric field in the mover. Starting from the mover and calling S gg the smoothing operator, we define a smoothed electric field on the grid as: From eq, (22), the smoothed electric field acting on a particle can be computed as: The second equation of motion, eq. (1), then uses the smoothed electric field as: where B n p is still computed as above. From eq. (24), we can compute again the current as The energy exchange term for the particles then becomes Applying now smoothing to the source term of the second of the Maxwell eq. (16), we have: where the last term is expressed from eq. (25). For the fields, the energy integral then becomes: For the two energy integrals to be the same the exchange term seen by the particles must be equal to that seen by the fields. The right-hand sides of eq.(26) must then equal that of eq. (28): Switching g with g (just names) in the right-hand side, the equivalence above holds when the smoothing operator is symmetric (i.e. the matrix repenting it is symmetric), a common property shared by many smoothing operators [9]. Note that we smooth the electric field but not the magnetic field that tends to be less noisy by its nature.

IV. SUB-CYCLING WITH THE MASS MATRIX FORMULATION
A mass matrix can be defined also in presence of sub-cycling or orbit averaging movers. The velocity update of ECsim can be reformulated for sub-cycling as: We assume that the time step ∆t between field updates is subdivided into N ν not necessarily equal sub-steps ∆t ν . The positions for the field evaluations, x ν p , during the sub-cycle can be computed in different ways. The simplest is to assume a straight orbit within ∆t, similar to the leap-frog approach: that can all be computed at once since the same velocity is used for all points along the trajectory. The first step starts from the old position: . The fields are assumed to be those computed at the time level θ within the field update time step ∆t. Another promising approach is to recall that most often in plasma physics particles are not moving in straight lines but rather they are frozen into the field lines, moving in cyclotron orbits with drifts due to the in-homogeneity of the fields. In the spirit of gyro-averaging, certain applications of implicit method might need to step over the gyration time scale and the positions of the particles used in eq.(30) would then be chosen to achieve accurate gyro-averaging, for example taking N ν positions along the gyro-orbit of a particle [26] to compute an average force on the particle's center of gyration.
The example of the two strategies above for computing the intermediate positions can be made in a single explicit step that generates all positions at once: in this case each substep contribution to the mass matrix and the moments can be computed in parallel, greatly improving the parallel performance. In practice, the N ν operations required by the substepping algorithm can all be done in parallel in an embarrassingly parallel approach that scales ideally on supercomputers: no communication between the particles and between the substeps is needed.
The equation (30) can be inverted with the same vector manipulations used for eq. (1), to obtain : where hatted quantities have been rotated by the magnetic field computed at the location x ν p : via a rotation matrix α n p defined as in the case of a singe step but the magnetic field computed at the last substep position: From eq. (32) we can obtain directly from its definition (6) the mean current over all sub-steps, without any further approximation or linearization: where W ν pg = W (x ν p − x g ). Using now the definitions of the hatted quantities, eq. (33), we can cast the mean current in the same form as in the single step formulation: but with the new definition for and the mass matrix for a sub-cycled trajectory defined by: Note that the definition of the sub-cycled mass matrix (eq. (38)) and the sub-cycled hatted current (eq. (37)) treats each subinterval for each particle as if they were independent: the equations see each sub-interval as a particle. It is as if the system has N p N ν particles made by each particle for each sub-interval ∆ ν . This feature lends itself to a simpler computing implementation where each particle is spawned into N ν treated as independent particles in the interpolation, mass matrix computation and current gathering step. A valuable approach in parallel and vectorized computer architectures such as GPUs.
Regardless of how the positions for the particles during the sybcycling are chosen, if the current defined in eq. (36) is computed with the mass matrices defined in eq. (38) energy is conserved. In fact, the term ∑ g J g · E g is again identical when computed from the equations for the particles and for the fields.

V. SIMPLIFICATION OF THE MASS MATRIX: THE LIMIT OF THE IMPLICIT MOMENT METHOD
There is a specific case where the mass matrix formulation takes a much simplified form: in case of the nearest grid point (NGP) interpolation. In that case, the interpolation function W pg is simple: 1 for the nearest grid point, that we label as g p , and 0 everywhere else: In this cases the mass matrix becomes: M i j s,gg = ∑ p q p α i j,n p δ g p g δ g p g When substituted in the expression for the current this leads to: Given the properties of Kronecker's delta, the summation over g can be done first: and substituting: where we have used the fact that the α's are computed using the magnetic fields of the nearest grid point, consistent with the NGP interpolation. Using again the properties of Kronecker's delta, the summation over the particles can be done directly: to obtain: This is the same expression that links the electric field and the current in the implicit moment method [32] used in Venus2D [6], in Celeste3D [20] and iPic3D [28]. It is then worth considering this limit expression for the mass matrix formulation. Most often, the NGP cannot be used in practice for its well known excessive noise, but we can still use the expression 45 even in presence of other interpolation schemes. Eq. (45) is exact only for the NGP scheme where it still leads to exact energy conservation. When, instead, other orders of interpolation are used, energy conservation is lost but it still provides a meaningful link between current and electric field: it is in fact the same used for decades by the implicit moment method. If eq. (45) is used, ECsim becomes more similar to the implicit moment method but it still differs in one key aspect: the interpolations are computed using the position provided by the leap-frog algorithm for the particle position. This is known explicitly and does not require any iteration. In the implicit moment method, the mover requires to iterate between velocity update and position update using a predictor-corrector scheme. This is not required in ECsim where the particle position is known explicitly from the previous time step.

A. Effects of sub-cycling
To check how the innovations described above perform in practice, we consider first the effect of sub-cycling. The goal is twofold. First, we want to verify that the energy is indeed conserved and sub-cycling does not break the energy conservation. Second, we test on a specific case how far one can push the sub-cycling before the physics deteriorates. The first task above is just a confirmation of the rigorously exact calculations above and it is merely a verification test for the code implemented here. The second task is only an illustration because the usefulness of sub-cycling is highly problem dependent.
A two-stream instability is initiated by considering a domain of L/d e = 2π divided into 64 cells with 10,000 particles. A mini-app is useful to study readily implementations in different computer architectures: it is intended merely as a test not as a real world problem. Space is normalized in electron skin depth, d e = c/ω pe , and time in electron plasma frequency, ω pe t. The ions are kept immobile while the electrons have a thermal speed of v the /c = 0.02 and the two beams have a net speed of V 0 /c = ±0.1. The reference case without sub-cycling (N ν = 1) has a CFL = V 0 ∆t/∆x = 0.1. We then proceed to do 10 runs with N ν = 1 − 10 and ∆t is increased by the same factor. That is we keep the particle ∆t fixed and increase the ∆t of the field. As the number of sub-cycles is increased the ∆t increases but the total time remains the same: ω pe t = 50. A perturbation of the particle velocity is added with mode number m = 5: with an amplitude V 1 = v 0 /10. We use a MATLAB implementation that we test on a MAC OSX 10.14.6 with processor Intel Core i7 2,6GHz complemented with 16GB DDR3 1600 MHz memory.
The results of the study are summarized in Table I and on Fig. 1 -2. Table I reports the reduction of cost of the simulation as the number of sub-cycles is increased. The MATLAB implementation is well vectorized and efficient in handling the particle projections, so much so that as the number of sub-cycles is increased the cost per ∆t hardly increases at all. If one compares the cost of a time step without sub-cyclig (0.08395s) with that with 10 sub-cycles (0.0958s), the increase is minimal even though particles have been moved and projected to the moments and the construction of the mass matrix 10 times more. While this result is specific to the vectorization done via MATLAB, other forms of vectorizations are also possible on GPUs using OpenACC and CUDA, opening up a similar opportunity for this type of optminization also in other implementations.
As required by the theoretical derivations above, energy should be conserved in all cases. And this is indeed the case in Table  I, a confirmation that our MATLAB implementation is bug free.
But if sub-cycling is beneficial in reducing the cost of the simulation, it introduces a degradation in the fidelity of the physics. Particles are moved with the same time step in each sub-cycle regardless of N ν but as N ν increases the overall time step ∆t for the field recalculation is increased. In a number of problem still it is beneficial to accept this compromise. Figure 1 and Fig. 2 show the degradation of the results with N ν . Unfortunately the instability modeled is a complex non-linear process and the linear phase is not important: an easy quantitative metric of accuracy is not available. We start already from a significant perturbation. If we do not add that, several modes develop at the same time and still the linear phase is muddied by the interaction of many modes. There is no single metric that can easily summarize the quality of the evolution. We have to rely on qualitative visual comparison. For this reason we look carefully at the velocity space distribution (Fig. 1) and phase space (Fig. 2). As N ν is increased the correct physics is progressively lost.
The most characteristic feature of the two stream instability is the formation of a flat top distribution, a distribution where for a range of velocities the distribution remains flat. This feature is prominent without sub-cycling, and it is still reasonably well represented at higher sub-cycling numbers. However, the tail of higher energy particles is reduced when sub-cycling is increased, a reflection of the fact that the electric field responsible for particle acceleration [23] is less accurately computed. These are non physical artifacts of excessive sub-cycling.
The same conclusion is reached analysing phase space. In this case the most important feature is the formation of electron holes, regions of phase space depleted of electrons and in fact completely void of them. These features, observed also in experimental measurements, are distorted and expelled to the edge of the velocity range as N ν is increased.
sub-cycling is compatible with the mass matrix formulation, retaining the property of exact energy conservation. It can be implemented very effectively, increasing only minimally the cost of the computational cycle, despite the increase in the number of particle operations, while reducing the number of time steps for the same total time. However, it must be used with care because not updating the fields after the particles are moved in a sub-cycle introduces physical errors even though energy is still exactly conserved.

B. Effects of smoothing
To test the effects of smoothing we consider another type of streaming instability: the transverse electromagnetic instability driven by counter streaming beams [11]. We consider again a 1D plasma (x is the only spatial variable) but we consider now all three components of the particle velocity. The two counter streaming beams are directed along y. Using the speed of light c for normalization, we consider a case with the beams having v th /c = 0.01 counter streaming with speed v 0y /c = 0.2. We add no initial perturbation and let the natural noise initiate the instability. The transverse electromagnetic streaming instability is well known [11] and it is often used in the context of understanding how magnetic fields can be generated in the Universe [25], it is a form of dynamo. We use a setup similar to that reported in Innocenti et al. [15]. A characteristic of this instability is that it segregates in phase space particles with opposite signs of v y , initially residing in the two different beams. Figure 3 shows the evolution of the phase space cross section (x, v x ): the red particles have v y < 0 and the blue particles v y > 0. The evolution initially retains their separation but in time phase space mixing takes over. More details of the physics of this simulation can be found in Innocenti et al. [15]. We are focusing here on comparing the normal ECsim case without smoothing with one where the smoothing kernel S = [1/4, 1/2, 1/4] is applied by convolution 3 times. Figure 3 compares the smoothed and non smoothed case. Smoothing is not altering in any profound ways the evolution but it affects it. Since the simulation is started from its natural noise, no two simulations will be identical and of course the smoothed simulation will have less noise. The location of the islands formed in phase space is not the same but the overall features are similar.
The energy exchange is also similar in the smoothed and non smoothed run, as shown in Fig. 4. In particular the kinetic energy is lost primarily to produce magnetic energy. The transverse counter streaming instability is a form of magnetic dynamo that spontaneously creates a magnetic field by using the kinetic energy of the counter streaming beams. The total energy is conserved to machine precision in both runs. Smoothing, as shown in Sect. III theoretically, indeed conserves energy exactly. The main effect of smoothing is to eliminate the high frequency part of the spectrum. Figure 5 and 6 compares the k − ω spectrum of selected fields. The noise at high values of k is reduced. The spectrum identifies the low k part of the spectrum as dominant. The characteristic arch of the electromagnetic waves (light waves) is also prominent.

VII. DISCUSSION
The ECsim method allows within a semi-implicit approach to make PIC simulation that conserve energy exactly. This feature besides having its own intrinsic value, leads also to improved stability allowing, for example, to consider realistic plasma conditions in the heliosphere, including the colder solar wind where the Debye length is very small compared with the scales of interest.
We report here three new developments. First, we introduced a method for smoothing the electric field in an algorithm that preserves energy conservation. When it can be avoided, smoothing should be avoided, but when there is a need for it, the algorithm presented here achieves smoothing without breaking energy conservation. This can be for example the case when excessive noise alters the correct transport properties of a plasma, an issue of great importance for example in fusion energy studies [29].
Second, we introduced a method for computing the mass matrix in presence of particle sub-cycling. Sub-cycling can be useful in a number of situations. Most often if it can be avoided, it should be avoided because particles and fields move together, the frozen-in condition being a cardinal property of plasmas. However, frozen in is valid at the large scales typical of MHD. There are many situations where the fields evolve slowly while particles move quickly. For example in gyro motion. Sub-cycling can be used together with gyroaveraging. The method presented here allows to use sub-cycling while constructing a mass matrix that continues to preserve the exact energy conservation. Finally, we discussed a limit case when the mass matrix calculation becomes especially simple and show that in this limit the plasma particle response in ECsim method becomes identical to that of the implicit moment method (IMM). This has two implications. First, it allows to use the same code either in the full energy conserving ECsim mode or in the less expensive IMM mode. Comparisons can then be made more readily. Second, the derivation clarifies the theoretical links between ECsim and IMM, revealing in what limit the two become identical.
These three new steps have both a theoretical significance and practical value, bringing understanding of the properties of ECsim and broadening its range of applications.

ACKNOWLEDGMENTS
The work reported received funding from the KULeuven Bijzonder Onderzoeksfonds (BOF) under the C1 project TRACESpace and from the European Union project DEEP-SEA (grant agreement 955606).