1. Introduction
The energy-conserving semi-implicit (ECSim) method is an algorithm for plasma simulation based on the particle-in-cell (PIC) approach [
1]. 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 [
2,
3]. In explicit methods, conversely, the coupling between particles and fields is suspended for a small time step [
4]. 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 straightforward with no iteration is needed, and explicit PIC can be implemented as one of the most efficient algorithms 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 the use of high resolution. The peculiarity of PIC is that the resolution needs to not just be refined in time to resolve the electron plasma frequency:
[
6], where
denotes the electron plasma frequency and
is the time step, but also in space to avoid the so-called finite-grid instability [
4]. This limitation is removed by the implicit approach that allows one to select grid spacing and a time step based on the accuracy needed and not on the stability of the numerical algorithm [
7].
Third, in explicit PIC, energy is not conserved. Using a good resolution, energy is acceptably maintained. There is always a secular trend of energy increase [
4], 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 [
8]. The implicit PIC method, instead, conserves energy exactly, whatever resolution is used [
2,
3]. This feature is physically important and practically impactful. Physically, of course, confidence comes from 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 [
9]. The direct implicit method uses a formal linear expansion of the coupling operator [
10]. 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 [
7]. 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 a mass matrix to express the coupling between particles and fields. With the mass matrix, the coupling is linear, but energy is conserved exactly. Below, it is reviewed how this was achieved [
1].
Since its recent introduction, ECsim has found application in a number of applications in space [
11,
12,
13,
14,
15] and in fusion [
16,
17]. An important improvement has removed the lack of charge conservation in the original scheme [
18,
19]. The current paper reports some extensions of the method that can widen its practical applicability.
First, a method to introduce smoothing to reduce noise while retaining exact energy conservation is described. Smoothing is an affordable way to effectively introduce 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. Instead, a method that conserves energy is presented.
Second, sub-cycling is a convenient approach in plasma simulations 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 straightforward since fields and particles are not advanced together. In implicit and semi-implicit PIC, instead, sub-cycling also needs to be done in the coupling, requiring the modification of 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 organized as follows.
Section 2 recaps the key properties of the ECsim needed for the discussion.
Section 3 introduces how smoothing can be implemented to reduce noise while retaining the property of energy conservation.
Section 4 is dedicated to the algorithm for sub-cycling the particle motion within the ECsim scheme.
Section 5 derives a limit case of the mass matrix formulation that transforms the ECsim algorithm into the standard implicit moment method [
9]. Results are presented in
Section 6, and final conclusions are drawn in
Section 7.
2. Summary of the ECsim Method
The ECsim method [
1] is based on a formulation of the mover similar to the classic leap-frog scheme:
where
denotes the particle (labeled with
p index) space position,
the particle velocity,
is particle charge,
is the particle mass,
n denotes the number of a particle, and
and
are electrical and magnetic fields, respectively,
is the averaged speed, and
.
This mover differs from both the explicit leap-frog mover [
20] and the implicit
-mover [
21]: 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
rather than at the unknown position
. 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
. In the
-scheme, instead, they are computed at the unknown position
. 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. However, one is computed explicitly, in the leap-frog sense, while the other is computed as part of a predictor-corrector iteration [
22,
23]. 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 [
2,
3].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 expressed using the magnetic field, , at the initial time level, but the electric field is considered at the advanced intermediate level: . 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 requires us to interpolate the fields to the particle positions:
A generic index
g is used for the grid. In the specific implementation within the iPic3D code [
24,
25], the electric field and magnetic field are not colocated, and
g labels either centers (for
) or vertices (for
). Here, the notation is simplified:
and
. In the implementation here, the interpolation function,
W, is a b-spline of order
[
26]:
This expression reduces trivially in 1D (one-dimensional) for the examples reported below.
For Maxwell’s equation, the standard
-scheme [
24] is used:
The spatial operators in Equation (
5) are discretized on the grid labeled 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, and
is the cell volume.
As with the
-mover, the velocity equation can be rewritten in the equivalent form [
27],
with
and the rotation matrix
given by
where
is the dyadic tensor (a matrix with a diagonal of 1) and
(independent of the particle weight and unique to a given species). The elements of the rotation matrix are indicated as
, with labels
i and
j referring to the 3 components of the vector space (
x,
y,
z).
Substituting Equation (
7) into Equation (
6), one obtains without any approximation or linearization,
where the notation is shortened,
, and the summation is intended over all particles of species
s.
Using Equation (
8), the expression for the current becomes
where
Computing then the electric field on the particles by interpolation form the grid as in Equation (
3), it follows that
Exchanging the order of summation, one obtains:
where the actor in the leading role of the ECsim scheme—the mass matrix [
28] is defined:
There are (where v is the number of velocity directions) mass matrices, and in matrix notation, they can be written as , that is, without the indices for the vector directions.
The mass matrices
are the most important aspect of the ECsim method and are also the most expensive part of the computation [
25]. A number of symmetries can be used to reduce the cost. Speeding up the construction can be achieved using offloading to accelerator processors (e.g., graphical processing units, GPU) [
29]. The mass matrices, as shown in Equation (
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 Equation (
5) to form a linear set of equations:
where the total current is
, and the species summed mass matrices that are 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 [
27], 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.
Let us consider now how energy conservation can be shown in the case . This is important for the derivations below to consider what key steps enable energy conservation.
The inner product of the velocity from Equation (
1) with the average speed,
, gives, by summing over all particles,
where the electric field is computed as an average consistent with the choice
, 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
.
Multiplying the first Equation (
5) by
and the second by
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 (Equation (
19)) needs to be identical to that on the field equations, (Equation (
21)). This term is indeed identically equal to
in both equations. Energy conservation is enforced exactly to round off.
In addition to 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 [
1], expanding the stability of the semi-implicit method compared with the moment implicit scheme [
11,
16].
3. Smoothing with the Mass Matrix Formulation
Smoothing can be designed to be compatible with the energy conserving properties of the mass matrix. Only the electric field is chosen to be smoothed, 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
the smoothing operator, a smoothed electric field on the grid is defined as:
From Equation (
22), the smoothed electric field acting on a particle can be computed as:
The second equation of motion, Equation (
1), then uses the smoothed electric field as:
where
is still computed as above.
From Equation (
24), one can again compute the current as
The energy exchange term for the particles then becomes
Applying smoothing to the source term of the second of the Maxwell Equation (
16), one obtains:
where the last term is expressed from Equation (
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 Equation (
26) must then equal that of Equation (
28):
Switching
g with
(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 [
3].
Note that the electric field is smoothed but not the magnetic field, which tends to be less noisy by its nature.
4. Sub-Cycling with the Mass Matrix Formulation
A mass matrix can also be defined in the presence of sub-cycling or orbit-averaging movers. The velocity update of ECsim can be reformulated for sub-cycling as
It is assumed assume that the time step between field updates is subdivided into , not necessarily equal sub-steps .
The positions for the field evaluations,
, during the sub-cycle can be computed in different ways. The simplest is to assume a straight orbit within
, similar to the leap-frog approach:
which 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
and old velocity
, and the last step leads to the final position
and final velocity
. The fields are assumed to be those computed at the time level
within the field update time step
.
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 inhomogeneity of the fields. In the spirit of gyro-averaging, certain applications of the implicit method might need to step over the gyration time scale, and the positions of the particles used in Equation (
30) would then be chosen to achieve accurate gyro-averaging, for example taking
positions along the gyro-orbit of a particle [
30] 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 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.
Equation (
30) can be inverted with the same vector manipulations used for Equation (
1) to obtain:
where hatted quantities have been rotated by the magnetic field computed at the location
:
via a rotation matrix,
, defined as in the case of a singe step, but with the magnetic field computed at the last substep position:
From Equation (
32), one can obtain directly from its definition (
6) the mean current over all sub-steps without any further approximation or linearization:
where
. Using now the definitions of the hatted quantities, Equation (
33), one 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 (Equation (
38)) and the sub-cycled hatted current (Equation (
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
particles made by each particle for each sub-interval
. This feature lends itself to a simpler computing implementation where each particle is spawned into
treated as independent particles in the interpolation, mass matrix computation and current gathering step. This is 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 Equation (
36) is computed with the mass matrices defined in Equation (
38), energy is conserved. In fact, the term
is again identical when computed from the equations for the particles and for the fields.
5. 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 the case of the interpolation of the nearest grid point (NGP). In that case, the interpolation function
is simple: 1 for the nearest grid point, which is labeled here as
, and 0 for everywhere else:
where
is the Kronecker delta.
In this case, the mass matrix becomes
When substituted in the expression for the current, this leads to
Given the properties of Kronecker delta, the summation over
can be done first:
and substituting:
where the point that the
s are computed using the magnetic fields of the nearest grid point is used, consistent with the NGP interpolation. Using again the properties of Kronecker 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 [
31] used in Venus2D [
9], in Celeste3D [
23] and iPic3D [
24]. 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 the expression (
45) still can be used even in the presence of other interpolation schemes.
Equation (
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 the current and electric field: it is, in fact, the same used for decades by the implicit moment method. If Equation (
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 the 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.
7. Discussion
The ECsim method allows to make a PIC simulation within a semi-implicit approach that conserves energy exactly. This feature, besides having its own intrinsic value, also leads 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.
This paper reports on three new developments.
First, a method for smoothing the electric field in an algorithm that preserves energy conservation is introduced. 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 [
38].
Second, a method for computing the mass matrix in the presence of particle sub-cycling is introduced. 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, with the frozen-in condition being a cardinal property of plasmas. However, frozen-in properties are 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 us to use sub-cycling while constructing a mass matrix that continues to preserve the exact energy conservation.
Finally, a limit case when the mass matrix calculation becomes especially simple is discussed and it is shown that, in this limit, the plasma particle response in the ECsim method becomes identical to that of the implicit moment method (IMM). This has two implications. First, it allows us 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.