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

Giovanni Lapenta  Center for mathematical Plasma Astrophysics, Department of
Mathematics, University of Leuven, KULeuven, Belgium

###### Abstract

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 [27,
9]. 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 $\omega_{pe}\Delta 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 [27, 9]. 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 semi-implicit 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 [22, 35, 36, 21, 24] 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:

$\begin{array}[]{c}\displaystyle\mathbf{x}_{p}^{n+1/2}=\mathbf{x}_{p}^{n-1/2}+\Delta
t\mathbf{v}_{p}^{n}\\\ \\\
\displaystyle\mathbf{v}_{p}^{n+1}=\mathbf{v}_{p}^{n}+\frac{q_{p}\Delta
t}{m_{p}}\left(\mathbf{E}^{n+\theta}(\mathbf{x}_{p}^{n+1/2})+\overline{\mathbf{v}}_{p}\times\mathbf{B}^{n}(\mathbf{x}_{p}^{n+1/2})\right)\end{array}$
(1)

where $\overline{\mathbf{v}}_{p}=(\mathbf{v}_{p}^{n+1}+\mathbf{v}_{p}^{n})/2$
and $\mathbf{E}^{n+\theta}=\theta\mathbf{E}^{n+1}+(1-\theta)\mathbf{E}^{n}$.

This mover differs from both the explicit leap-frog mover [4] and the implicit
$\theta$-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 $\mathbf{x}_{p}^{n+1/2}$ rather than at the unknown position
$\overline{\mathbf{x}}_{p}$. The important consequence is that the particle
equations can be solved directly without any iteration needed among
themselves. Instead in the standard $\theta$-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 $\mathbf{x}_{p}^{n+1/2}$. In the $\theta$-scheme, instead, they
are computed at the unknown position $\overline{\mathbf{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 [33, 20]. 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
$\theta$-scheme can be made energy conserving but at the cost of a fully non-
linear iteration requiring a non-linear solver [27, 9].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 $\mathbf{B}^{n}(\mathbf{x}_{p}^{n+1/2})$ but the electric field is
written at the advanced intermediate level
$\mathbf{E}^{n+\theta}(\mathbf{x}_{p}^{n+1/2})$. 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:

$\displaystyle\mathbf{E}_{p}^{n+\theta}$
$\displaystyle=\mathbf{E}^{n+\theta}(\mathbf{x}_{p}^{n+1/2})=\sum_{g}\mathbf{E}_{g}^{n+\theta}W(\mathbf{x}_{p}^{n+1/2}-\mathbf{x}_{g})$
(2) $\displaystyle\mathbf{B}_{p}^{n}$
$\displaystyle=\mathbf{B}^{n+\theta}(\mathbf{x}_{p}^{n+1/2})=\sum_{g}\mathbf{B}_{g}^{n}W(\mathbf{x}_{p}^{n+1/2}-\mathbf{x}_{g})$
(3)

We used a generic index $g$ for the grid. In the specific implementation
within the iPic3D code [28, 12] the electric field and magnetic field are not
colocated and $g$ label either centers (for $\mathbf{B}$) or vertices (for
$\mathbf{E}$). Here we simplify the notation as:
$\mathbf{B}_{p}^{n}=\mathbf{B}^{n}(\mathbf{x}_{p}^{n+1/2})$ and
$\mathbf{E}_{p}^{n+\theta}=\mathbf{E}^{n+\theta}(\mathbf{x}_{p}^{n+1/2})$. In
our implementation, the interpolation function $W$ is a b-splines of order
$\ell=1$ [3]:

$W(\mathbf{x}_{p}-\mathbf{x}_{g})=b_{\ell}(x_{p}-x_{g})b_{\ell}(y_{p}-y_{g})b_{\ell}(z_{p}-z_{g})$
(4)

This expression reduces trivially in 1D for the examples reported below.

For the Maxwell’s equation we use the standard $\theta$-scheme [28]:

$\begin{array}[]{ccc}\displaystyle\nabla_{g}\times\mathbf{E}^{n+\theta}+\frac{1}{c}\frac{\mathbf{B}^{n+1}_{g}-\mathbf{B}^{n}_{g}}{\Delta
t}=0\\\ \\\
\displaystyle\nabla_{g}\times\mathbf{B}^{n+\theta}-\frac{1}{c}\frac{\mathbf{E}^{n+1}_{g}-\mathbf{E}^{n}_{g}}{\Delta
t}=\frac{4\pi}{c}\overline{\mathbf{J}}_{g}\end{array}$ (5)

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:

$\overline{\mathbf{J}}_{sg}=\frac{1}{V_{g}}\sum_{p\in
s}q_{p}\overline{\mathbf{v}}_{p}W(\mathbf{x}_{p}^{n+1/2}-\mathbf{x}_{g})$ (6)

where the summation is over the particles of the same species, labeled by $s$.

As with the $\theta$-mover, the velocity equation can be rewritten in the
equivalent form [34]:

$\overline{\mathbf{v}}_{p}=\widehat{\mathbf{v}}_{p}+\beta_{s}\widehat{\mathbf{E}}_{p}$
(7)

with:

$\begin{array}[]{c}\widehat{\mathbf{v}}_{p}={\alpha}^{n}_{p}\mathbf{v}^{n}_{p}\\\
\\\
\widehat{\mathbf{E}}_{p}={\alpha}^{n}_{p}\mathbf{E}_{p}^{n+\theta}\end{array}$
(8)

and the rotation matrix ${\alpha}_{p}^{n}$ given by:

${\alpha}_{p}^{n}=\frac{1}{1+(\beta_{s}B_{p}^{n})^{2}}\left(\mathbb{I}-\beta_{s}\mathbb{I}\times\mathbf{B}_{p}^{n}+\beta_{s}^{2}\mathbf{B}_{p}^{n}\mathbf{B}_{p}^{n}\right)$
(9)

where $\mathbb{I}$ is the dyadic tensor (matrix with diagonal of 1) and
$\beta_{s}=q_{p}\Delta t/2m_{p}$ (independent of the particle weight and
unique to a given species). The elements of the rotation matrix are indicated
as ${\alpha}^{ij,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:

$\overline{\mathbf{J}}_{sg}=\frac{1}{V_{g}}\sum_{p}q_{p}\widehat{\mathbf{v}}_{p}W_{pg}+\frac{\beta_{s}}{V_{g}}\sum_{p}q_{p}\widehat{\mathbf{E}}_{p}^{n+\theta}W_{pg}$
(10)

where we shortened the notation
$W_{pg}=W(\mathbf{x}_{p}^{n+1/2}-\mathbf{x}_{g})$ and the summation is
intended over all particles of species $s$.

Using eq. (8), the expression for the current becomes:

$\overline{\mathbf{J}}_{sg}=\widehat{\mathbf{J}}_{sg}+\frac{\beta_{s}}{V_{g}}\sum_{p}q_{p}{\alpha}^{n}_{p}\mathbf{E}_{p}^{n+\theta}W_{pg}$
(11)

where we defined:

$\widehat{\mathbf{J}}_{sg}=\frac{1}{V_{g}}\sum_{p}q_{p}\widehat{\mathbf{v}}_{p}W_{pg}$
(12)

Computing then the electric field on the particles by interpolation form the
grid as in eq. (3), it follows that:

$\overline{\mathbf{J}}_{sg}=\widehat{\mathbf{J}}_{sg}+\frac{\beta_{s}}{V_{g}}\sum_{p}\sum_{g^{\prime}}q_{p}{\alpha}^{n}_{p}\mathbf{E}_{g^{\prime}}^{n+\theta}W_{pg^{\prime}}W_{pg}$
(13)

Exchanging the order of summation we obtain:

$\overline{\mathbf{J}}_{sg}=\widehat{\mathbf{J}}_{sg}+\frac{\beta_{s}}{V_{g}}\sum_{g^{\prime}}M_{s,gg^{\prime}}\mathbf{E}_{g^{\prime}}^{n+\theta}$
(14)

where we have defined the actor in the leading role of the ECsim scheme: the
mass matrix [8]:

$M_{s,gg^{\prime}}^{ij}=\sum_{p}q_{p}{\alpha}^{ij,n}_{p}W_{pg^{\prime}}W_{pg}$
(15)

There are $3v$ (where $v$ is the number of velocity directions) mass matrices
and in matrix notation they can be written as $M_{gg^{\prime}}$, that is
without the indices $i,j$ for the vector directions.

The mass matrices $M_{s,gg^{\prime}}$ 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:

$\left\\{\begin{array}[]{l}\displaystyle\nabla_{g}\times\mathbf{E}^{n+\theta}+\frac{1}{c}\frac{\mathbf{B}^{n+1}-\mathbf{B}^{n}}{\Delta
t}=0\\\ \\\
\displaystyle\nabla_{g}\times\mathbf{B}^{n+\theta}-\frac{1}{c}\frac{\mathbf{E}^{n+1}-\mathbf{E}^{n}}{\Delta
t}=\frac{4\pi}{c}\left(\widehat{\mathbf{J}}_{g}+\sum_{g^{\prime}}M_{gg^{\prime}}\mathbf{E}_{g^{\prime}}^{n+\theta}\right)\end{array}\right.$
(16)

where the total current is
$\widehat{\mathbf{J}}_{g}=\sum_{s}\widehat{\mathbf{J}}_{sg}$ and the species
summed mass matrices, that written by elements are:

$M_{gg^{\prime}}^{ij}=\sum_{s}\frac{\beta_{s}}{V_{g}}M_{s,gg^{\prime}}^{ij}$
(17)

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 $\theta=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, $\overline{\mathbf{v}}_{p}$ gives by summing over all
particles:

$\frac{1}{2}\sum_{p}\left(m_{p}(\mathbf{v}_{p}^{n+1})^{2}-(\mathbf{v}_{p}^{n})^{2}\right)=\Delta
t\sum_{p}\left(q_{p}\sum_{g}\overline{\mathbf{v}}_{p}\cdot\overline{\mathbf{E}}_{g}W_{pg}\right)$
(18)

where the electric field is computed as average consistent with the choice
$\theta=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:

$\frac{1}{2}\sum_{p}\left(m_{p}(\mathbf{v}_{p}^{n+1})^{2}-(\mathbf{v}_{p}^{n})^{2}\right)=\Delta
t\sum_{g}\overline{\mathbf{J}}_{g}\cdot\overline{\mathbf{E}}_{g}$ (19)

where it is recognized that
$\overline{\mathbf{J}}_{g}=\sum_{p}q_{p}\overline{\mathbf{v}}_{p}W_{pg}$.

Multiplying the first equation (5) by $\overline{\mathbf{B}}_{g}$ and the
second by $\overline{\mathbf{E}}_{g}$ and summing them leads to:

$\begin{array}[]{l}\displaystyle\frac{(\mathbf{B}_{g}^{n+1})^{2}-(\mathbf{B}_{g}^{n})^{2}}{2c}+\frac{(\mathbf{E}_{g}^{n+1})^{2}-(\mathbf{E}_{g}^{n})^{2}}{2c}=\\\
\\\ \displaystyle\Delta
t\left(\frac{4\pi}{c}\overline{\mathbf{J}}_{g}\cdot\overline{\mathbf{E}}_{g}+\mathbf{E}_{g}\cdot\nabla_{g}\times\mathbf{B}-\mathbf{B}_{g}\cdot\nabla_{g}\times\mathbf{E}\right)\end{array}$
(20)

Assuming a mimetic grid discretization that preserves the continuum properties
of the operators and summing over all grid points gives:

$\begin{array}[]{l}\displaystyle\sum_{g}\frac{(\mathbf{B}_{g}^{n+1})^{2}-(\mathbf{B}_{g}^{n})^{2}}{4\pi}+\sum_{g}\frac{(\mathbf{E}_{g}^{n+1})^{2}-(\mathbf{E}_{g}^{n})^{2}}{4\pi}=\\\
\\\ \displaystyle\Delta
t\sum_{g}\overline{\mathbf{J}}_{g}\cdot\overline{\mathbf{E}}_{g}+\frac{c\Delta
t}{4\pi}\sum_{g}\nabla_{g}\cdot(\mathbf{E}_{g}\times\mathbf{B}_{g})\end{array}$
(21)

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
$\sum_{g}\overline{\mathbf{J}}_{g}\cdot\overline{\mathbf{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 semi-implicit method compared with the moment implicit scheme
[22, 13].

##  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 $\mathbb{S}_{gg^{\prime}}$ the smoothing operator, we define a
smoothed electric field on the grid as:

$\mathbf{E}^{SM}_{g}=\sum_{g^{\prime}}\mathbb{S}_{gg^{\prime}}\mathbf{E}^{n+\theta}_{g^{\prime}}$
(22)

From eq, (22), the smoothed electric field acting on a particle can be
computed as:

$\mathbf{E}^{SM}_{p}=\sum_{g}\mathbf{E}^{SM}_{g}W(\mathbf{x}_{p}^{n+1/2}-\mathbf{x}_{g})$
(23)

The second equation of motion, eq. (1), then uses the smoothed electric field
as:

$\mathbf{v}_{p}^{n+1}=\mathbf{v}_{p}^{n}+\frac{q_{p}\Delta
t}{m_{p}}\left(\mathbf{E}^{SM}_{p}+\overline{\mathbf{v}}_{p}\times\mathbf{B}^{n}_{p}\right)$
(24)

where $\mathbf{B}^{n}_{p}$ is still computed as above.

From eq. (24), we can compute again the current as

$\overline{\mathbf{J}}_{sg}=\widehat{\mathbf{J}}_{sg}+\frac{\beta_{s}}{V_{g}}\sum_{g^{\prime}}M_{s,gg^{\prime}}\mathbf{E}^{SM}_{g^{\prime}}$
(25)

The energy exchange term for the particles then becomes

$\frac{1}{2}\sum_{p}\left(m_{p}(\mathbf{v}_{p}^{n+1})^{2}-(\mathbf{v}_{p}^{n})^{2}\right)=\Delta
t\sum_{s}\sum_{g}\left(\widehat{\mathbf{J}}_{sg}+\frac{\beta_{s}}{V_{g}}\sum_{g^{\prime}}M_{s,gg^{\prime}}\mathbf{E}^{SM}_{g^{\prime}}\right)\cdot\mathbf{E}^{SM}_{g}$
(26)

Applying now smoothing to the source term of the second of the Maxwell eq.
(16), we have:

$\nabla_{g}\times\mathbf{B}^{n+\theta}-\frac{1}{c}\frac{\mathbf{E}^{n+1}_{g}-\mathbf{E}^{n}_{g}}{\Delta
t}=\frac{4\pi}{c}\sum_{s}\sum_{g^{\prime}}\mathbb{S}_{gg^{\prime}}\left(\widehat{\mathbf{J}}_{sg^{\prime}}+\frac{\beta_{s}}{V_{g}}\sum_{g^{\prime\prime}}M_{s,g^{\prime}g^{\prime\prime}}\mathbf{E}^{SM}_{g^{\prime\prime}}\right)$
(27)

where the last term is expressed from eq. (25).

For the fields, the energy integral then becomes:

$\begin{array}[]{l}\displaystyle\sum_{g}\frac{(\mathbf{B}_{g}^{n+1})^{2}-(\mathbf{B}_{g}^{n})^{2}}{4\pi}+\sum_{g}\frac{(\mathbf{E}_{g}^{n+1})^{2}-(\mathbf{E}_{g}^{n})^{2}}{4\pi}-\frac{c\Delta
t}{4\pi}\sum_{g}\nabla_{g}\cdot(\mathbf{E}_{g}\times\mathbf{B}_{g})=\\\ \\\
\displaystyle\Delta
t\sum_{s}\sum_{g}\sum_{g^{\prime}}\mathbb{S}_{gg^{\prime}}\left(\widehat{\mathbf{J}}_{sg^{\prime}}+\frac{\beta_{s}}{V_{g^{\prime}}}\sum_{g^{\prime\prime}}M_{s,g^{\prime}g^{\prime\prime}}\mathbf{E}^{SM}_{g^{\prime\prime}}\right)\cdot\overline{\mathbf{E}}_{g}\end{array}$
(28)

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):

$\begin{split}\sum_{g}\sum_{g^{\prime}}\mathbb{S}_{gg^{\prime}}\left(\widehat{\mathbf{J}}_{sg^{\prime}}+\frac{\beta_{s}}{V_{g^{\prime}}}\sum_{g^{\prime\prime}}M_{s,g^{\prime}g^{\prime\prime}}\mathbf{E}^{SM}_{g^{\prime\prime}}\right)\cdot\overline{\mathbf{E}}_{g}=\\\
\sum_{g}\left(\widehat{\mathbf{J}}_{sg}+\frac{\beta_{s}}{V_{g}}\sum_{g^{\prime\prime}}M_{s,gg^{\prime\prime}}\mathbf{E}^{SM}_{g^{\prime\prime}}\right)\cdot\sum_{g^{\prime}}\mathbb{S}_{gg^{\prime}}\overline{\mathbf{E}}_{g^{\prime}}\end{split}$
(29)

Switching $g$ with $g^{\prime}$ (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:

$\displaystyle\mathbf{v}_{p}^{\nu+1}=\mathbf{v}_{p}^{\nu}+\frac{q_{p}\Delta
t_{\nu}}{m_{p}}\left(\mathbf{E}^{n+\theta}(\mathbf{x}_{p}^{\nu})+\frac{\mathbf{v}_{p}^{\nu+1}+\mathbf{v}_{p}^{\nu}}{2}\times\mathbf{B}^{n}(\mathbf{x}_{p}^{\nu})\right)$
(30)

We assume that the time step $\Delta t$ between field updates is subdivided
into $N_{\nu}$ not necessarily equal sub-steps $\Delta t_{\nu}$.

The positions for the field evaluations, $\mathbf{x}_{p}^{\nu}$, during the
sub-cycle can be computed in different ways. The simplest is to assume a
straight orbit within $\Delta t$, similar to the leap-frog approach:

$\mathbf{x}_{p}^{\nu}=\mathbf{x}_{p}^{n-1/2}+\mathbf{v}_{p}^{n}\sum_{\nu^{\prime}=0}^{\nu^{\prime}=\nu}\Delta
t_{\nu^{\prime}}$ (31)

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:
$\mathbf{x}_{p}^{\nu=0}=\mathbf{x}_{x}^{n}$ and old velocity
$\mathbf{v}_{p}^{\nu=0}=\mathbf{v}_{x}^{n}$ and the last step leads to the
final position $\mathbf{x}_{p}^{\nu=N_{\nu}}=\mathbf{x}_{p}^{n+1}$ and final
velocity $\mathbf{v}_{p}^{\nu=N_{\nu}}=\mathbf{v}_{p}^{n+1}$ . The fields are
assumed to be those computed at the time level $\theta$ within the field
update time step $\Delta 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_{\nu}$ 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_{\nu}$ 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 :

$\frac{\mathbf{v}_{p}^{\nu+1}+\mathbf{v}_{p}^{\nu}}{2}=\widehat{\mathbf{v}}_{p}^{\nu}+\beta_{s}\widehat{\mathbf{E}}_{p}^{\nu}$
(32)

where hatted quantities have been rotated by the magnetic field computed at
the location $\mathbf{x}_{p}^{\nu}$:

$\begin{array}[]{c}\widehat{\mathbf{v}}_{p}^{\nu}={\alpha}^{\nu}_{p}\mathbf{v}^{\nu}_{p}\\\
\\\
\widehat{\mathbf{E}}_{p}^{\nu}={\alpha}^{\nu}_{p}\mathbf{E}^{n+\theta}(\mathbf{x}_{p}^{\nu})\end{array}$
(33)

via a rotation matrix ${\alpha}_{p}^{n}$ defined as in the case of a singe
step but the magnetic field computed at the last substep position:

${\alpha}_{p}^{\nu}=\frac{1}{1+(\beta_{s}B^{n}(\mathbf{x}_{p}^{\nu}))^{2}}\left(\mathbb{I}-\beta_{s}\mathbb{I}\times\mathbf{B}^{n}(\mathbf{x}_{p}^{\nu})+\beta_{s}^{2}\mathbf{B}^{n}(\mathbf{x}_{p}^{\nu})\mathbf{B}^{n}(\mathbf{x}_{p}^{\nu})\right)$
(34)

From eq. (32) we can obtain directly from its definition (6) the mean current
over all sub-steps, without any further approximation or linearization:

$\overline{\mathbf{J}}_{sg}=\frac{1}{V_{g}}\sum_{p}q_{p}\sum_{\nu}\frac{\Delta
t_{\nu}}{\Delta
t}\left(\widehat{\mathbf{v}}_{p}^{\nu}+\beta_{s}\widehat{\mathbf{E}}_{p}^{\nu}\right)W_{pg}^{\nu}$
(35)

where $W_{pg}^{\nu}=W(\mathbf{x}_{p}^{\nu}-\mathbf{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:

$\overline{\mathbf{J}}_{sg}=\widehat{\mathbf{J}}_{sg}+\frac{\beta_{s}}{V_{g}}\sum_{g^{\prime}}M_{s,gg^{\prime}}\mathbf{E}_{g^{\prime}}^{n+\theta}$
(36)

but with the new definition for

$\widehat{\mathbf{J}}_{sg}=\frac{1}{V_{g}}\sum_{p}q_{p}\sum_{\nu}\frac{\Delta
t_{\nu}}{\Delta t}\widehat{\mathbf{v}}_{p}^{\nu}W_{pg}^{\nu}$ (37)

and the mass matrix for a sub-cycled trajectory defined by:

$M_{s,gg^{\prime}}^{ij}=\sum_{p}q_{p}\sum_{\nu}\frac{\Delta t_{\nu}}{\Delta
t}{\alpha}^{ij,\nu}_{p}W_{pg^{\prime}}^{\nu}W_{pg}^{\nu}$ (38)

Note that the definition of the sub-cycled mass matrix (eq. (38)) and the sub-
cycled hatted current (eq. (37)) treats each sub-interval 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_{\nu}$ particles made by each particle for
each sub-interval $\Delta_{\nu}$. This feature lends itself to a simpler
computing implementation where each particle is spawned into $N_{\nu}$ 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
$\sum_{g}\overline{\mathbf{J}}_{g}\cdot\overline{\mathbf{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:

$W_{pg}=\delta_{g_{p}g}$ (39)

In this cases the mass matrix becomes:

$M_{s,gg^{\prime}}^{ij}=\sum_{p}q_{p}{\alpha}^{ij,n}_{p}\delta_{g_{p}g^{\prime}}\delta_{g_{p}g}$
(40)

When substituted in the expression for the current this leads to:

$\overline{\mathbf{J}}_{sg}=\widehat{\mathbf{J}}_{sg}+\frac{\beta_{s}}{V_{g}}\sum_{p}\sum_{g^{\prime}}q_{p}{\alpha}^{n}_{p}\mathbf{E}_{g^{\prime}}^{n+\theta}\delta_{g_{p}g^{\prime}}\delta_{g_{p}g}$
(41)

Given the properties of Kronecker’s delta, the summation over $g^{\prime}$ can
be done first:

$\sum_{g^{\prime}}\mathbf{E}_{g^{\prime}}^{n+\theta}\delta_{g_{p}g^{\prime}}=\mathbf{E}_{g_{p}}^{n+\theta}$
(42)

and substituting:

$\overline{\mathbf{J}}_{sg}=\widehat{\mathbf{J}}_{sg}+\frac{\beta_{s}}{V_{g}}\sum_{p}q_{p}{\alpha}^{n}_{g_{p}}\mathbf{E}_{g_{p}}^{n+\theta}\delta_{g_{p}g}$
(43)

where we have used the fact that the $\alpha$’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:

$\frac{1}{V_{g}}\sum_{p}q_{p}=\rho_{sg}$ (44)

to obtain:

$\overline{\mathbf{J}}_{sg}=\widehat{\mathbf{J}}_{sg}+\beta_{s}\rho_{sg}{\alpha}^{n}_{g}\mathbf{E}_{g}^{n+\theta}$
(45)

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.

##  VI Results

###  VI.1 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\pi$ 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/\omega_{pe}$, and time in electron
plasma frequency, $\omega_{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=\pm 0.1$. The reference case without sub-cycling
($N_{\nu}=1$) has a CFL = $V_{0}\Delta t/\Delta x=0.1$. We then proceed to do
10 runs with $N_{\nu}=1-10$ and $\Delta t$ is increased by the same factor.
That is we keep the particle $\Delta t$ fixed and increase the $\Delta t$ of
the field. As the number of sub-cycles is increased the $\Delta t$ increases
but the total time remains the same: $\omega_{pe}t=50$. A perturbation of the
particle velocity is added with mode number $m=5$:

$u_{p}=u_{p}+V_{1}\sin\left(\frac{2\pi mx_{p}}{L}\right);$ (46)

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 1 and on Fig. 1 -2.

$N_{\nu}$ | Time (s) | Seconds/ $\Delta t$ | $\Delta E$
---|---|---|---
1.0 | 42.73 | 0.08395 | $1.75\,{10}^{-16}$
2.0 | 20.94 | 0.08244 | $1.1\,{10}^{-16}$
3.0 | 14.63 | 0.08657 | $1.015\,{10}^{-16}$
4.0 | 10.88 | 0.08567 | $7.388\,{10}^{-17}$
5.0 | 8.97 | 0.08881 | $6.764\,{10}^{-17}$
6.0 | 7.36 | 0.08762 | $8.342\,{10}^{-17}$
7.0 | 6.46 | 0.08972 | $7.988\,{10}^{-17}$
8.0 | 5.74 | 0.09111 | $8.492\,{10}^{-17}$
9.0 | 5.27 | 0.09411 | $8.14\,{10}^{-17}$
10.0 | 4.79 | 0.0958 | $8.61\,{10}^{-17}$
Table 1: Effect of sub-cycling on a two-stream instability run. Varying the
number of sub-cycles from 1 to 10, while correspondingly increasing by the
same factor $\Delta t$, the duration of the simulation in seconds, the seconds
spent per $\Delta t$ and the energy error are reported.  Figure 1: Particle
velocity distribution at the end of 10 different simulations of same total
duration varying the number of sub-cycles $N_{\nu}=1-10$ (from top,
$N_{\nu}=1$, to bottom, $N_{\nu}=10$) and correspondingly increasing $\Delta
t$ by the same factor. Figure 2: Scatter plot on the phase space for the
particles at the end of 10 different simulations of same total duration
varying the number of sub-cycles $N_{\nu}=1-10$ (from top, $N_{\nu}=1$, to
bottom, $N_{\nu}=10$) and correspondingly increasing $\Delta t$ by the same
factor.

Table 1 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 $\Delta 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 1, 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_{\nu}$ but as
$N_{\nu}$ increases the overall time step $\Delta 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_{\nu}$. 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_{\nu}$ 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_{\nu}$ 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.

###  VI.2 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 $\mathbf{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.

No smoothing | Smoothing
---|---
a) $\omega_{pe}t=6.35$ | b) $\omega_{pe}t=6.35$
|
c) $\omega_{pe}t=25$ | d) $\omega_{pe}t=25$
|
e) $\omega_{pe}t=43.75$ | f) $\omega_{pe}t=43.75$
|
g) $\omega_{pe}t=62.5$ | h) $\omega_{pe}t=62.5$
|
Figure 3: Evolution of the phase space section $(x,v_{x})$ in the transverse
counter streaming instability. The red particles have $v_{y}<0$ and the blue
particles $v_{y}>0$. Three times are shown and the left (right) panels report
the non smoothed (smoothed) run.

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 $\mathbb{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.

No smoothing | Smoothing
---|---
|
Figure 4: Evolution of the energy in the transverse counter streaming
instability. The left panels show the case without smoothing and the right
panels with smoothing. The tope panel show the exchange between kinetic and
magnetic energy with a minority contribution going to the electric field
energy. The bottom panel shows the total energy conservation.

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-\omega$ 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.

|  |
---|---|---
Figure 5: Transverse counter streaming instability in absence of smoothing. Three fields are shown from left to right: $E_{x}$, $E_{z}$ and $B_{y}$. The top row shows the spatiotemporal plane $(x,t)$ and the bottom row the spectral plane $(k,\omega)$.  |  |
---|---|---
Figure 6: Transverse counter streaming instability in presence of smoothing.
Three fields are shown from left to right: $E_{x}$, $E_{z}$ and $B_{y}$. The
top row shows the spatiotemporal plane $(x,t)$ and the bottom row the spectral
plane $(k,\omega)$.

##  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.

###### Acknowledgements.

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).

## References

  * Birdsall & Langdon [2004] Birdsall, C., & Langdon, A. 2004, Plasma Physics Via Computer Simulation (London: Taylor & Francis)
  * Boella et al. [2022] Boella, E., Innocenti, M. E., Bettencourt, M., et al. 2022, in OpenACC, Location: Virtual
  * Boor [1978] Boor, C. D. 1978, A practical guide to splines (Springer)
  * Boris et al. [1970] Boris, J. P., et al. 1970, in Proc. Fourth Conf. Num. Sim. Plasmas, 3–67
  * Bowers et al. [2009] Bowers, K. J., Albright, B. J., Yin, L., et al. 2009, Journal of Physics Conference Series, 180, 012055,
  * Brackbill & Forslund [1982] Brackbill, J., & Forslund, D. 1982, J. Computat. Phys., 46, 271
  * Brackbill & Cohen [1985] Brackbill, J. U., & Cohen, B. I., eds. 1985, Multiple time scales.
  * Burgess et al. [1992] Burgess, D., Sulsky, D., & Brackbill, J. 1992, Journal of Computational Physics, 103, 1
  * Chen et al. [2011] Chen, G., Chacón, L., & Barnes, D. C. 2011, Journal of Computational Physics, 230, 7018
  * Chen & Tóth [2019] Chen, Y., & Tóth, G. 2019, Journal of Computational Physics, 386, 632
  * Fried [1959] Fried, B. D. 1959, The Physics of Fluids, 2, 337.
  * Gonzalez-Herrero et al. [2018] Gonzalez-Herrero, D., Boella, E., & Lapenta, G. 2018, Computer Physics Communications, 229, 162
  * Gonzalez-Herrero et al. [2019] Gonzalez-Herrero, D., Micera, A., Boella, E., Park, J., & Lapenta, G. 2019, Computer Physics Communications, 236, 153
  * Hockney & Eastwood [1988] Hockney, R., & Eastwood, J. 1988, Computer simulation using particles (Taylor & Francis)
  * Innocenti et al. [2011] Innocenti, M., Lazar, M., Markidis, S., Lapenta, G., & Poedts, S. 2011, Physics of Plasmas, 18, 052104
  * Jianyuan et al. [2018] Jianyuan, X., Hong, Q., & Jian, L. 2018, Plasma Science and Technology, 20, 110501
  * Langdon et al. [1983] Langdon, A., Cohen, B., & Friedman, A. 1983, J. Computat. Phys., 51, 107
  * Lapenta [2012] Lapenta, G. 2012, Journal of Computational Physics, 231, 795
  * Lapenta [2017] —. 2017, Journal of Computational Physics, 334, 349
  * Lapenta et al. [2006] Lapenta, G., Brackbill, J. U., & Ricci, P. 2006, Phys. Plasmas, 13, 055904
  * Lapenta et al. [2020] Lapenta, G., El Alaoui, M., Berchem, J., & Walker, R. 2020, Journal of Geophysical Research: Space Physics, 125, no
  * Lapenta et al. [2017] Lapenta, G., Gonzalez-Herrero, D., & Boella, E. 2017, Journal of Plasma Physics, 83
  * Lapenta et al. [2007] Lapenta, G., Markidis, S., Marocchino, A., & Kaniadakis, G. 2007, The Astrophysical Journal, 666, 949
  * Lapenta et al. [2022] Lapenta, G., Schriver, D., Walker, R. J., et al. 2022, Journal of Geophysical Research: Space Physics, 127, e2021JA030241
  * Lazar et al. [2009] Lazar, M., Schlickeiser, R., Wielebinski, R., & Poedts, S. 2009, The Astrophysical Journal, 693, 1133
  * Lee [1987] Lee, W. 1987, Journal of Computational Physics, 72, 243
  * Markidis & Lapenta [2011] Markidis, S., & Lapenta, G. 2011, Journal of Computational Physics, 230, 7037
  * Markidis et al. [2010] Markidis, S., Lapenta, G., & Rizwan-uddin. 2010, Mathematics and Computers and Simulation, 80, 1509
  * Nevins et al. [2005] Nevins, W., Hammett, G., Dimits, A. M., Dorland, W., & Shumaker, D. 2005, Physics of plasmas, 12, 122305
  * Park et al. [2019] Park, J., Lapenta, G., Gonzalez-Herrero, D., & Krall, N. A. 2019, Frontiers in Astronomy and Space Sciences, 6, 74
  * Pinto & Pagès [2022] Pinto, M. C., & Pagès, V. 2022, Journal of Computational Physics, 453, 110912
  * Ricci et al. [2002] Ricci, P., Lapenta, G., & Brackbill, J. 2002, J. Computat. Phys., 183, 117
  * Vu & Brackbill [1995] Vu, H., & Brackbill, J. 1995, Journal of Computational Physics, 116, 384
  * Vu & Brackbill [1992] Vu, H. X., & Brackbill, J. U. 1992, Comput. Phys. Comm., 69, 253
  * Walker et al. [2019] Walker, R. J., Lapenta, G., Berchem, J., El-Alaoui, M., & Schriver, D. 2019, Journal of Plasma Physics, 85
  * Zhou et al. [2019] Zhou, H., Tóth, G., Jia, X., Chen, Y., & Markidis, S. 2019, Journal of Geophysical Research: Space Physics, 124, 5441