Valentin Pagès
email: pagesv@ljll.math.upmc.fr

Martin Campos Pinto
email: martin.campos-pinto@ipp.mpg.de

Martin Campos

A semi-implicit electromagnetic FEM-PIC scheme with exact energy and charge conservation

Keywords: Semi-implicit particle-in-cell, Energy conservation, Charge conservation, Vlasov-Maxwell, Kinetic plasma simulation, Structure-preserving finite elements

, which essentially consists of improving the latter with appropriate charge-conserving current deposition terms and an explicit correction step. As such, the resulting scheme is free of nonlinear iterations, and it preserves both energy and charge exactly.

Introduction

The long-time accuracy of numerical simulations usually relies on the fact that fundamental invariants of the physical system, such as its total energy, are well preserved at the discrete level. In the case of electromagnetic particle solvers, it has been long recognized that Gauss's laws also represent key invariants [START_REF] Boris | Relativistic plasma simulation-optimization of a hybrid code[END_REF][START_REF] Marder | A method for incorporating Gauss' law into electromagnetic PIC codes[END_REF][START_REF] Villasenor | Rigorous charge conservation for local electromagnetic field solvers[END_REF], with spurious virtual charges building up when they are not preserved by the numerical method, see e.g. [START_REF] Barthelmé | Numerical charge conservation in particle-incell codes[END_REF][START_REF] Stock | Investigation of the purely hyperbolic maxwell system for divergence cleaning in discontinuous galerkin based particle-in-cell methods[END_REF][START_REF] Pinto | Charge-conserving FEM-PIC schemes on general grids[END_REF]. At the continuous level, the preservation of Gauss's laws by the time-dependent Maxwell-Vlasov equations hinges on two properties: (i) the charge and current densities satisfy a continuity equation, and (ii) the divergence of a curl is always zero. On a larger perspective, the key role played by Gauss's laws in the long-time stability of Maxwell's equations is now well understood in the scope of their de Rham structure [START_REF] Bossavit | Computational electromagnetism: variational formulations, complementarity, edge elements[END_REF][START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF][START_REF] Gross | Electromagnetic theory and computation: a topological approach[END_REF][START_REF] Campos | Compatible Maxwell solvers with particles I: conforming and non-conforming 2d schemes with a strong Ampere law[END_REF].

In Particle-in-Cell (PIC) methods [START_REF] Hockney | Computer simulation using particles[END_REF][START_REF] Birdsall | Plasma physics via computer simulation[END_REF] where Maxwell solvers are coupled to particle approximate sources, a variety of charge-(that is, Gauss-) preserving solvers have been proposed to adress this issue. They may be divided in two categories: field correction methods, where the electromagnetic field is modified to eliminate its deficiency in Gauss's laws [START_REF] Boris | Relativistic plasma simulation-optimization of a hybrid code[END_REF][START_REF] Munz | Divergence Correction Techniques for Maxwell Solvers Based on a Hyperbolic Model[END_REF][START_REF] Langdon | On enforcing Gauss' law in electromagnetic particle-incell codes[END_REF][START_REF] Marder | A method for incorporating Gauss' law into electromagnetic PIC codes[END_REF], and conservative deposition methods, where the current density is computed in a way that preserves a proper continuity equation at the discrete level [START_REF] Villasenor | Rigorous charge conservation for local electromagnetic field solvers[END_REF][START_REF] Zh | Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor[END_REF][START_REF] Eastwood | The virtual particle electromagnetic particle-mesh method[END_REF][START_REF] Umeda | A new charge conservation method in electromagnetic particle-in-cell simulations[END_REF][START_REF] Pinto | Charge-conserving FEM-PIC schemes on general grids[END_REF]. Similarly as on the continuous level, these methods also require that the Maxwell solver preserves some essential features of the underlying de Rham structure. For this reason they may be connected to the larger class of structure-preserving methods, which were extensively studied in the last decades : see e.g. [START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF][START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF][START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF][START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF] for the general properties of these methods, and [START_REF] Buffa | Isogeometric analysis in electromagnetics: B-splines approximation[END_REF][START_REF] Campos | Compatible Maxwell solvers with particles II: conforming and non-conforming 2d schemes with a strong Faraday law[END_REF][START_REF] Kraus | GEMPIC: geometric electromagnetic particle-in-cell methods[END_REF][START_REF] Kormann | Energy-conserving time propagation for a geometric particle-in-cell Vlasov-Maxwell solver[END_REF][START_REF] D Y Na | Axisymmetric charge-conservative electromagnetic particle simulation algorithm on unstructured grids: Application to microwave vacuum electronic devices[END_REF][START_REF] Moon | Exact chargeconserving scatter-gather algorithm for particle-in-cell simulations on unstructured grids: a geometric perspective[END_REF][START_REF] Mark | Discretizing Transient Current Densities in the Maxwell Equations[END_REF] for their applications in particular problems from plasma physics or electrical engineering.

Recently, several PIC codes have been proposed that preserve exactly the total energy of the discrete system. In [START_REF] Markidis | The energy conserving particlein-cell method[END_REF], Markidis and Lapenta have proposed an energy-conserving method based on an implicit time discretization. The resulting fully nonlinear time steps are solved by an iterative Newton Krylov scheme, and a particle enslavement technique is proposed to reformulate the nonlinear problems on the fields only. In this method the Gauss laws are satisfied approximately, by adding field correction terms as pseudo-currents in the nonlinear problems. In a related method [START_REF] Chacón | A charge-and energy-conserving implicit, electrostatic particle-in-cell algorithm on mapped computational meshes[END_REF], Chen, Chacón and Barnes have taken this approach a step further by taking advantage of the particle enslavement principle to deposit the current with a charge-conserving algorithm, which allows to preserve exactly the Gauss laws even when the iterative Newton scheme is not applied until convergence. At this point we may remind that the classical, explicit leap-frog time stepping [START_REF] Hockney | Computer simulation using particles[END_REF] coupled with a charge-conserving current deposition [START_REF] Eastwood | The virtual particle electromagnetic particle-mesh method[END_REF][START_REF] Pinto | Charge-conserving FEM-PIC schemes on general grids[END_REF] allows to preserve the Gauss laws, but it only preserves a modified discrete energy. As a consequence of its explicit nature, stable runs require an explicit CFL constraint which limits the time steps to approximately the cell propagation time of a light wave. In contrast, energy conserving schemes enjoy the possibility of using much larger time steps. However in both methods above the implicit steps are fully nonlinear.

This limitation is removed in the recent Energy-Conserving Semi-Implicit Method (ECSIM) proposed by Lapenta [START_REF] Lapenta | Exactly energy conserving semi-implicit particle in cell formulation[END_REF]. In this novel scheme, the particle enslavement principle is cleverly applied at the very level of the time stepping in such a way that the implicit problems are not only reduced to the field unknowns, but they also become linear. A drawback of this semi-implicit time stepping is that it is a priori not compatible with a charge-conserving current deposition algorithm. Indeed the latter involves time-averages on the particle trajectories, which cannot be expressed as a linear function of the electric field.

To address this issue, a few methods have been devised: in [START_REF] Chen | Gauss's Law satisfying Energy-Conserving Semi-Implicit Particle-in-Cell method[END_REF] the ECSIM scheme is completed by a correction on the particle positions based on a local linearization of the particle shapes, which can effectively reduce the Gauss' law error when applied in a recursive fashion. A different approach is followed in [START_REF] Kormann | Energy-conserving time propagation for a geometric particle-in-cell Vlasov-Maxwell solver[END_REF], where the authors consider the general framework of geometric (hamiltonian) and structure-preserving FEM-PIC methods [START_REF] Kraus | GEMPIC: geometric electromagnetic particle-in-cell methods[END_REF]. There a time stepping is proposed which preserves exactly the energy and the Gauss laws, albeit at the cost of a new nonlinear problem to be solved on the fields.

In this article we propose a new method where the semi-implicit time stepping is free of any nonlinear iterations, and where both the energy and Gauss laws are preserved exactly. Our strategy essentially consists of first endowing Lapenta's ECSIM scheme [START_REF] Lapenta | Exactly energy conserving semi-implicit particle in cell formulation[END_REF] with a semi-conservative current deposition, and of adding an explicit correction to complete the time step with a fully charge-conserving deposition. A velocity rescaling is eventually performed in order to restore the energy conservation of the original scheme. As a result, the computational cost of ChECSIM is close to that of ECSIM with an additional conservative current deposition. similar to what is done in the explicit charge-conserving leap-frog scheme [START_REF] Pinto | Charge-conserving FEM-PIC schemes on general grids[END_REF]. Moreover it is naturally formulated in the general framework of structure-preserving finite elements, as in [START_REF] Campos | Compatible Maxwell solvers with particles II: conforming and non-conforming 2d schemes with a strong Faraday law[END_REF][START_REF] Perse | Energy-conserving Implicit Time Discretisation for the GEMPIC Framework[END_REF].

This framework is important as it also allows unstructured meshes on complex geometries, with possible local grid adaptations. Indeed, one important field of applications is the industrial context of vacuum electronic devices, where designers are interested in performing numerous long-time simulations within a rather short time-frame. Fast solvers usually rely on simplified models, such as separating the representation of the fields in the beam area in resonant Klystron cavities, see e.g. [START_REF] Cooke | Validation of the large-signal klystron Simulation code TESLA[END_REF][START_REF] Gao | Development of a 2.5-Dimensional Particle-In-Cell Code for Efficient High-Power Klystron Design[END_REF]. Nevertheless, in the ever more challenging improvement of the radio-frequency power sources, simulations based on reduced models have shown difficulties to capture troublemaking phenomena that impede the development of new components, and consequently the demand for solvers of the complete time-dependent Maxwell-Vlasov model has arisen. To this end, stable schemes emancipating from CFL constraints and having a complexity that remains affordable for short runtime are attractive, provided that they are robust enough to show long-time accuracy.

The article is organized as follows. After recalling the model in Section 2 we summarize in Section 3 the principles of structure-preserving FEM-PIC spatial discretizations. Once this is settled, we provide in Section 4 a comparative description of the two algorithms that we use as reference FEM-PIC schemes, namely the explicit, charge-conserving leap-frog (CCLF) scheme [START_REF] Pinto | Charge-conserving FEM-PIC schemes on general grids[END_REF] and the ECSIM scheme of Lapenta [START_REF] Lapenta | Exactly energy conserving semi-implicit particle in cell formulation[END_REF], reformulated in the framework of finite elements. With these building blocks at hand we can introduce the ChECSIM scheme in Section 5. The guess step is described in Section 5.1 as a semi-conservative modification of the original ECSIM, and the corrective step is presented in Section 5.2. This section is completed by an a priori estimate which shows that the eventual velocity correction is possible under a mild condition on the time step that is linked to the steepness of the electric field, see Proposition 5.7. If necessary, a dynamic adaptation of the time step suffices to fulfill this criterion. In Section 6 we finally perform some numerical experiments on a couple of standard 2D test cases from plasma physics, which confirm the improved stability properties of the proposed method.

General Framework

The continuous model

Vlasov-Maxwell's equations are among the most accurate models in plasma physics or in the design of radio-frequency sources. They consist of Vlasov's equation for the charged particles phase-space density (for the sake of simplicity, we restrict ourselves to the case of a single species of particles),

∂ t f + v • ∇ x f + q m (E + v × B) • ∇ v f = 0, (1)
which essentially expresses the transport of particles with elementary charge q and mass m along the electromagnetic Lorentz force, and of Maxwell's equations for the electromagnetic fields. These may also be decomposed in two sets of equations: an evolution system consisting of Ampère and Faraday's equations

∂ t E -c 2 curl B = - 1 0 J , (2)
∂ t B + curl E = 0, (3)
and Gauss's laws

div E = 1 0 ρ, (4)
div

B = 0. (5)
Here, the source terms are the charge and current densities carried by the particles' distribution, defined as the first and second velocity moments of f ,

ρ = ρ[f ] = q ˆR3 f dv, (6)
J = J [f ] = q ˆR3 vf dv. (7)
It should be noted that the evolution system composed of the Vlasov-Ampère-Faraday equations, namely (1)-( 3), defines a unique solution in time starting from an initial one. The Gauss laws (4)-( 5) may then be seen as divergence constraints on the fields, and it is an elementary exercice to verify that (i) they will be preserved by any solution to (1)-( 3) if the sources satisfy the so-called continuity equation,

∂ t ρ + div J = 0, (8)
and (ii) this continuity equation is satisfied by any source of the form ( 6)-( 7), as long as f solves a transport equation of the form (1).

As emphasized in the introduction, another key invariant is the total energy

E = 0 E 2 2 + B 2 2µ 0 + m 2 ˆR3 |v| 2 f dv (9)
which is preserved by exact solutions in the absence of external contributions.

Discrete FEM-PIC models

In finite element Particle-In-Cell approximations of the Vlasov-Maxwell system, the electromagnetic field (E, B) is represented by a finite element field (E h , B h ), typically a piecewise polynomial on a given mesh of the computational domain Ω ⊂ R 3 , that is determined as the solution to some discrete version of the evolution equations ( 2)-( 3), namely Ampère's and Faraday's equations. This will be specified in Section 3 and below. As for the plasma phase-space density f , it is represented by a weighted collection of N numerical particles (here with Dirac masses), of the form

f N (t, x, v) = N k=1 w k δ(x -x k (t))δ(v -v k (t)), (10)
with positions x k and velocities v k that follow discrete Lorentz force trajectories. The coupling of these discrete solutions then consists of specifying (i) how the finite element fields define the Lorentz force at the particle positions, see Section 3.3, and (ii) how the current and charge density sources in the discrete FEM Maxwell system are defined from their particle approximations,

ρ N := ˆR3 f N dv = N p=1 q k δ(x -x k (t)), (11)
J N := ˆR3 vf N dv = N k=1 q k v k (t)δ(x -x k (t)). ( 12
)
Here q k = qw k is the charge carried by the numerical particle of index k, and likewise m k = mw k will denote its mass (for simplicity we consider here a single particle species of charge q and mass m).

Since FEM-PIC solutions rely only on the time-domain equations ( 1)-( 3), preserving a discrete version of the Gauss laws (4)- [START_REF] Barthelmé | Numerical charge conservation in particle-incell codes[END_REF], as well as a discrete analog to the total energy [START_REF] Bossavit | Computational electromagnetism: variational formulations, complementarity, edge elements[END_REF], is clearly a non-straightforward feature to enforce. In particular, we already observe that preserving the Gauss laws will require a particular care in the discretization of the FEM current density source, indeed the latter plays a pivotal role in the charge-conserving coupling between Maxwell and Vlasov equations (see Section 3.3 below).

Spatial discretization with structure-preserving finite elements

In this section we recall the main features of a structure-preserving finite elements discretization of the Maxwell equations, and its coupling with point particles. Similarly as in [START_REF] Monk | An analysis of Nédélec's method for the spatial discretization of Maxwell's equations[END_REF][START_REF] Mark | Discretizing Transient Current Densities in the Maxwell Equations[END_REF][START_REF] Moon | Exact chargeconserving scatter-gather algorithm for particle-in-cell simulations on unstructured grids: a geometric perspective[END_REF][START_REF] Campos | Compatible Maxwell solvers with particles II: conforming and non-conforming 2d schemes with a strong Faraday law[END_REF][START_REF] Kraus | GEMPIC: geometric electromagnetic particle-in-cell methods[END_REF][START_REF] D Y Na | Axisymmetric charge-conservative electromagnetic particle simulation algorithm on unstructured grids: Application to microwave vacuum electronic devices[END_REF] we follow the general framework of Finite Element Exterior Calculus (FEEC), detailed in [START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF][START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF][START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF]. It essentially relies on finite elements spaces that preserve the de Rham structure of the continuous problem [START_REF] Bossavit | Computational electromagnetism: variational formulations, complementarity, edge elements[END_REF]. Remark 3.1. In this article we consider finite element discretizations for their high order approximation properties on general meshes. A similar derivation is possible in the context of finite difference schemes, by using the classical methods of Villasenor-Buneman [START_REF] Villasenor | Rigorous charge conservation for local electromagnetic field solvers[END_REF] or Esirkepov [START_REF] Zh | Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor[END_REF] instead of the finite element conservative deposition that will be recalled below.

Mixed weak-form of Maxwell equations

For simplicity, we consider here that the Maxwell equations are posed on a bounded open domain Ω of R 3 with Lipschitz-continuous boundary supplemented with boundary conditions of metallic type, see e.g. [START_REF] Monk | Finite Element Methods for Maxwell's Equations[END_REF],

n × E = 0 on ∂Ω ( 13
)
where n is the outward normal unit vector (for boundary conditions of absorbing type, we refer to Appendix C). This leads us to considering the following de Rham sequence,

V 0 V 1 V 2 V 3 grad curl div (14)
with Hilbert spaces

V 0 = H 1 0 (Ω), V 1 = H 0 (curl, Ω), V 2 = H 0 (div, Ω) and V 3 = L 2 (Ω)
, endowed with their standard norms [START_REF] Boffi | Mixed finite element methods and applications[END_REF][START_REF] Assous | Mathematical foundations of computational electromagnetism[END_REF].

The first step is to reformulate the time-dependent Maxwell equations in weak form, using only the spaces involved in the sequence [START_REF] Campos | Compatible Maxwell solvers with particles I: conforming and non-conforming 2d schemes with a strong Ampere law[END_REF]. Following e.g. [START_REF] Monk | An analysis of Nédélec's method for the spatial discretization of Maxwell's equations[END_REF], a common choice is to integrate by parts in the Ampère equation. Given the metallic boundary condition on the electric field [START_REF] Campos | Gauss-compatible Galerkin schemes for time-dependent Maxwell equations[END_REF], this gives

∂ t E(t), φ 1 -c 2 B(t), curl φ 1 = - 1 0 J (t), φ 1 , ∀φ 1 ∈ V 1 , (15)
∂ t B(t), φ 2 + curl E(t), φ 2 = 0, ∀φ 2 ∈ V 2 (16)
for all t ∈ [0, T ]. Here the bracket •, • denotes the standard scalar product in L 2 (Ω), and the solution (E, B) is sought for in the spaces

C 1 ([0, T ]; V 1 ) × C 1 ([0, T ]; V 2
). We observe that curl E ∈ C 0 ([0, T ]; V 2 ), so that (16) amounts to writing the Faraday equation in a strong sense, namely

∂ t B(t) + curl E(t) = 0 ∀t ∈ [0, T ]. ( 17
)

Structure-preserving finite elements

In order to preserve the de Rham structure at the discrete level, we consider structure-preserving finite element spaces which form a discrete sequence,

V 0 h V 1 h V 2 h V 3 h . grad curl div (18)
As studied in [START_REF] Campos | Compatible Maxwell solvers with particles I: conforming and non-conforming 2d schemes with a strong Ampere law[END_REF][START_REF] Campos | Compatible Maxwell solvers with particles II: conforming and non-conforming 2d schemes with a strong Faraday law[END_REF], a key feature of such discretization spaces is to provide discrete Gauss laws that allow for a rigorous long-time stability analysis of chargeconserving Maxwell solvers. Moreover, it is known that structure-preserving FEM schemes associated with stable commuting diagrams are free of spurious eigenvalues for the Maxwell evolution operator, see e.g. [START_REF] Hiptmair | Finite elements in computational electromagnetism[END_REF][START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF][START_REF] Campos | Gauss-compatible Galerkin schemes for time-dependent Maxwell equations[END_REF]. Several sequences of finite elements may be chosen here. A classical option on unstructured meshes consists in using the Lagrange-Nédélec-Raviart-Thomas sequence [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms[END_REF][START_REF] Boffi | Mixed finite element methods and applications[END_REF] studied and generalized in [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF], see also [START_REF] Campos | Gauss-compatible Galerkin schemes for time-dependent Maxwell equations[END_REF]. Another option is to follow [START_REF] Kraus | GEMPIC: geometric electromagnetic particle-in-cell methods[END_REF] and use the compatible spline spaces defined on tensor-product domains [START_REF] Buffa | Isogeometric analysis in electromagnetics: B-splines approximation[END_REF][START_REF] Buffa | Isogeometric discrete differential forms in three dimensions[END_REF][START_REF] Ratnani | An Arbitrary High-Order Spline Finite Element Solver for the Time Domain Maxwell Equations[END_REF], or mimetic spectral elements [START_REF] Gerritsma | An Introduction to a Compatible Spectral Discretization Method[END_REF][START_REF] Palha | Physics-compatible discretization techniques on single and dual grids, with application to the Poisson equation of volume forms[END_REF]. Using structure-preserving finite elements we approximate the equations above by a semi-discrete system

∂ t E h (t), φ 1 h -c 2 B h (t), curl φ 1 h = - 1 0 J h (t), φ 1 h , ∀φ 1 h ∈ V 1 h , ( 19
)
∂ t B h (t) + curl E h (t) = 0, (20)
again for all t ∈ [0, T ]. Here the solution (E h , B h ) is sought for in the spaces

C 1 ([0, T ]; V 1 h )×C 1 ([0, T ]; V 2 h ), and J h is a source current density in C 0 ([0, T ]; V 1 h ).
For the subsequent study of the charge-conservation properties, it will be convenient to rewrite the weak Ampère equation as

∂ t E h -c 2 curl h B h = - 1 0 J h , (21)
where the discrete operator

curl h : V 2 h → V 1 h ( 22
)
is defined as the discrete adjoint of the operator curl :

V 1 h → V 2 h , i.e.,

by the relations

curl h φ 2 h , φ 1 h = φ 2 h , curl φ 1 h , ∀φ 1 h ∈ V 1 h and φ 2 h ∈ V 2 h .
Following [START_REF] Campos | Compatible Maxwell solvers with particles I: conforming and non-conforming 2d schemes with a strong Ampere law[END_REF][START_REF] Campos | Compatible Maxwell solvers with particles II: conforming and non-conforming 2d schemes with a strong Faraday law[END_REF], we know that this discrete functional setting allows to derive long-time stability estimates when the associated discrete Gauss laws, namely

div h E h = 1 0 ρ h (23)
div B h = 0 [START_REF] Gao | Development of a 2.5-Dimensional Particle-In-Cell Code for Efficient High-Power Klystron Design[END_REF] are preserved in time. Here the magnetic Gauss law is written in strong form and is obviously preserved since Faraday's equation is also solved in a strong form [START_REF] Eastwood | The virtual particle electromagnetic particle-mesh method[END_REF], whereas [START_REF] Frenkel | Understanding Molecular Simulation[END_REF] involves the discrete weak divergence

div h : V 1 h → V 0 h ( 25
)
defined as the adjoint of the operator -grad : V 0 h → V 1 h , and ρ h ∈ V 0 h is defined as the L 2 projection of the particle charge density [START_REF] Buffa | Isogeometric discrete differential forms in three dimensions[END_REF], that is,

ρ h (t), φ 0 h = ρ N (t), φ 0 h = N k=1 q k φ 0 h (x k (t)), ∀φ 0 h ∈ V 0 h . (26)
In particular, the electric Gauss law [START_REF] Frenkel | Understanding Molecular Simulation[END_REF] is equivalently rewritten in a weak finite element sense,

-E h , grad φ 0 h = 1 0 ρ h , φ 0 h , ∀φ 0 h ∈ V 0 h . (27)
In order to write System ( 19)- [START_REF] Eastwood | The virtual particle electromagnetic particle-mesh method[END_REF] in matrix form, we choose a basis Λ i , i = 1, . . . , N for each space V h , = 0, . . . , 3 and we denote by σ i the associated degrees of freedom. Again, examples of degrees of freedom can be found in [START_REF] Boffi | Mixed finite element methods and applications[END_REF]. Note that the particular choice of degrees of freedom will affect the mass and operator matrices, but not the resulting solutions of system ( 19)- [START_REF] Eastwood | The virtual particle electromagnetic particle-mesh method[END_REF]. The corresponding mass matrices are then

M = Λ i , Λ j 1≤i,j≤N , (28)
and we also denote by

C = σ 2 i (curl Λ 1 j ) 1≤i≤N2,1≤j≤N1 , the N 2 × N 1 matrix of the operator curl : V 1 h → V 2 h
, and by

G = σ 1 i (grad Λ 0 j ) 1≤i≤N1,1≤j≤N0
, the N 1 × N 0 matrix of the operator grad : V 0 h → V 0 h . Finally we write the time-dependent field coefficients as column vectors,

e(t) = σ 1 i (E h (t)) 1≤i≤N1 and b(t) = σ 2 i (B h (t)) 1≤i≤N2 .
System ( 19)-( 20) reads then

M 1 d dt e -c 2 (M 2 C) T b = - 1 0 j (29)
d dt b + Ce = 0 ( 30
)
where j contains the moments of the finite element current source, namely

j(t) = J h (t), Λ 1 i 1≤i≤N1 ,
and discrete Gauss' law [START_REF] Gross | Electromagnetic theory and computation: a topological approach[END_REF] reads

(M 1 G) T e = 1 0 M 0 r (31)
where r contains the moments of the finite element charge density, namely

r(t) = ρ h (t), Λ 0 i 1≤i≤N0 ,

Energy-and charge-conserving coupling with particles

At the time-continuous level, we consider the same particle-field coupling as done in, e.g., [START_REF] Pinto | Charge-conserving FEM-PIC schemes on general grids[END_REF]. Using discrete point particles, the trajectories read

   ẋk (t) = v k (t), vk (t) = q k m k E h (t, x k (t)) + v k (t) × B h (t, x k (t)) , (32)
and the FEM current source

J h ∈ V 1
h is defined as the generalized L 2 projection of the particle current density [START_REF] Pinto | Charge-conserving FEM-PIC schemes on general grids[END_REF], namely

J h (t), φ 1 h = J N (t), φ 1 h = N k=1 q k v k (t) • φ 1 h (x k (t)), ∀φ 1 h ∈ V 1 h . (33)
In order to reduce the numerical noise due to Dirac particle shapes, one option is to consider smoothed charge and current densities of the form ρ S N := S * ρ N and J S N := S * J N , where S is a local smoothing kernel, such as a B-spline. The extension of our method to such cases can be performed along the lines of schemes studied in [START_REF] Ameres | Variational framework for structure-preserving electromagnetic particle-in-cell methods[END_REF] In the following it will be convenient to use a compact notation for the particles positions and velocities,

X = x k 1≤k≤N and V = v k 1≤k≤N .
Proposition 3.2. The time-continuous FEM-PIC equations (29)-( 33) preserve the discrete energy

E h (V , e, b) = N k=1 m k 2 |v k | 2 + 0 E h 2 2 + B h 2 2µ 0 ( 34
)
and the discrete Gauss laws (23)- [START_REF] Gao | Development of a 2.5-Dimensional Particle-In-Cell Code for Efficient High-Power Klystron Design[END_REF].

Proof. This result is well known. The conservation of energy is verified by taking the time derivative of (34) and using [START_REF] Langdon | On enforcing Gauss' law in electromagnetic particle-incell codes[END_REF] together with ( 19)- [START_REF] Eastwood | The virtual particle electromagnetic particle-mesh method[END_REF],

d dt E = k m k v k • vk + 0 E h , ∂ t E h + 1 µ 0 B h , ∂ t B h (35) = k q k v k • E h (x k ) + c 2 0 curl E h , B h -E h , J h - 1 µ 0 B h , curl E h = 0 (36)
where we have used 0 µ 0 = 1/c 2 , and the definition of J h in the last equality, see [START_REF] Lapenta | Exactly energy conserving semi-implicit particle in cell formulation[END_REF]. Turning to the Gauss laws, we first observe that the strong magnetic Gauss law readily follows from the Faraday equation being satisfied in a strong sense. To verify next that the electric Gauss law is preserved, we differentiate in time [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms[END_REF],

∂ t ρ h , φ 0 h = N k=1 q k v k • grad φ 0 h (x k ) = J h , grad φ 0 h ( 37
)
which shows that the FE sources satisfy the semi-discrete continuity equation

∂ t ρ h + div h J h = 0 ( 38
)
involving the weak divergence operator div h , see [START_REF] Gerritsma | An Introduction to a Compatible Spectral Discretization Method[END_REF]. Applying this operator to the discrete (weak) Ampere equation ( 21) yields then

∂ t div h E h = div h c 2 curl h B h - 1 0 J h = ∂ t 1 0 ρ h . ( 39
)
This shows that the discrete Gauss law ( 23) is preserved.

Two reference FEM-PIC schemes: CCLF and ECSIM

In this section we describe the two fully discrete FEM-PIC schemes that will serve us as reference methods in this work, corresponding to two different time integrations of the equations ( 29)- [START_REF] Lapenta | Exactly energy conserving semi-implicit particle in cell formulation[END_REF]. The first one is the charge-conserving leap-frog scheme (CCLF) studied in [START_REF] Pinto | Charge-conserving FEM-PIC schemes on general grids[END_REF], and the second one is a finite element version of the Energy conserving semi-implicit (ECSIM) scheme of Lapenta [START_REF] Lapenta | Exactly energy conserving semi-implicit particle in cell formulation[END_REF], similar to the version derived in [START_REF] Perse | Energy-conserving Implicit Time Discretisation for the GEMPIC Framework[END_REF]. Although they are naturally expressed on staggered times, we will present them on the same integer time-step [t n , t n+1 ] to facilitate comparisons, i.e.

(X n , V n , e n , b n ) → (X n+1 , V n+1 , e n+1 , b n+1 ).

The charge-conserving leap-frog scheme (CCLF)

The explicit leap-frog scheme [START_REF] Hockney | Computer simulation using particles[END_REF] consists of integrating both the particle trajectories [START_REF] Langdon | On enforcing Gauss' law in electromagnetic particle-incell codes[END_REF] and the field equations ( 19)-( 20) in a staggered fashion. On an integer time-step, this reads

             v n+ 1 2 k -v n k = ∆t 2 q k m k E n h (x n k ) + v n k × B n h (x n k ) , x n+1 k -x n k = ∆tv n+ 1 2 k , v n+1 k -v n+ 1 2 k = ∆t 2 q k m k E n+1 h (x n+1 k ) + v n+1 k × B n+1 h (x n+1 k ) , (40)
for the particle trajectories (note that

v n k = 1 2 (v n+ 1 2 k + v n-1 2 k
)), and

             b n+ 1 2 -b n = - ∆t 2 Ce n , M 1 (e n+1 -e n ) = ∆t c 2 (M 2 C) T b n+ 1 2 - 1 0 j n+ 1 2 , b n+1 -b n+ 1 2 = - ∆t 2 Ce n+1 , (41)
for the Maxwell equations. Following [START_REF] Eastwood | The virtual particle electromagnetic particle-mesh method[END_REF][START_REF] Pinto | Charge-conserving FEM-PIC schemes on general grids[END_REF] we consider here a time averaged current source, namely

j n+ 1 2 = k q k v n+ 1 2 k • {Λ 1 i } n k 1≤i≤N1
where the curly brackets denote time averages over interpolated trajectories,

{g} n k := 1 ∆t ˆtn+1 t n g x k (t) dt with x k (t) = x n k + (t -t n )v n+ 1 2 k . (42)
Here, a key point is that the v

n+ 1 2 k
are constant on the time-step, so that they can be placed either inside or in front of the the curly brackets {•} n k . We observe that (41) may be expressed in terms of the finite element fields and operators defined in Section 3.2, as

             B n+ 1 2 h -B n h = - ∆t 2 curl E n h , E n+1 h -E n h = c 2 ∆t curl h B n+ 1 2 h - ∆t 0 J n+ 1 2 h , B n+1 h -B n+ 1 2 h = - ∆t 2 curl E n+1 h , (43)
with a FEM current source

J h ∈ V 1 h defined as J n+ 1 2 h , φ 1 h = k q k v n+ 1 2 k • {φ 1 h } n k , ∀φ 1 h ∈ V 1 h . (44)
An important property of this method is that it preserves exactly the proper discrete Gauss laws, see [START_REF] Pinto | Charge-conserving FEM-PIC schemes on general grids[END_REF][START_REF] Campos | Compatible Maxwell solvers with particles II: conforming and non-conforming 2d schemes with a strong Faraday law[END_REF].

Proposition 4.1. The CCLF scheme (40)-( 41) preserves the discrete Gauss laws (23)-( 24) evaluated at the discrete times, namely

div h E n h = 1 0 ρ n h (45) div B n h = 0 ( 46
)
with a FEM charge density defined as the instant value (26

), i.e., ρ n h := ρ h (t n ).
Proof. The argument is the same as for the time continuous result (Prop. 3.2), the key ingredient being that the time-averaged current (44) naturally leads to a proper discrete continuity equation. Indeed the interpolated trajectories (42) have a constant velocity over the time-step [t n , t n+1 ], so that the above current density corresponds to

J n+ 1 2 h = 1 ∆t ˆtn+1 t n J h (t) dt,
where the time-continuous current is defined as in [START_REF] Lapenta | Exactly energy conserving semi-implicit particle in cell formulation[END_REF]. In particular, a time integration of (38) readily yields

1 ∆t ρ n+1 h -ρ n h + div h J n+ 1 2 h = 0, (47)
and applying the discrete div h operator to the FEM Ampere equation [START_REF] Perse | Energy-conserving Implicit Time Discretisation for the GEMPIC Framework[END_REF] shows that the discrete Gauss law ( 45) is indeed preserved. As for [START_REF] Mark | Discretizing Transient Current Densities in the Maxwell Equations[END_REF], it readily follows from the strong form of the discrete Faraday equation in [START_REF] Perse | Energy-conserving Implicit Time Discretisation for the GEMPIC Framework[END_REF].

Energy-Conserving Semi-Implicit Method (ECSIM)

As is well known, the above leap-frog scheme does not preserve the energy of the discrete solution,

E h (V n , e n , b n ) = k m k 2 |v n k | 2 + 0 E n h 2 2 + B n h 2 2µ 0 . ( 48
)
In the source-free case for example, the leap-frog Maxwell scheme only preserves a pseudo-energy, which controls the energy (48) when a CFL time-stepping condition is satisfied, see e.g. [START_REF] Taflove | Computational electrodynamics : the finite-difference timedomain method[END_REF]. As recalled in the introduction, a time scheme that exactly preserves energy at the cost of a linear implicitation in the Maxwell solver has been recently proposed by Lapenta [START_REF] Lapenta | Exactly energy conserving semi-implicit particle in cell formulation[END_REF]. Here we describe this energy-conserving semi-implicit method (ECSIM) in our FEM setting. The particle pusher reads

             x n+ 1 2 k -x n k = ∆t 2 v n k , v n+1 k -v n k = ∆t q k m k E n+ 1 2 h (x n+ 1 2 k ) + v n+ 1 2 k × B n h (x n+ 1 2 k
) ,

x n+1 k -x n+ 1 2 k = ∆t 2 v n+1 k , (49)
with short-hand notations

v n+ 1 2 k := 1 2 (v n k + v n+1 k ) and E n+ 1 2 h := 1 2 (E n h + E n+1 h ),
and the discrete Maxwell solver takes the following implicit form

   M 1 (e n+1 -e n ) = ∆t c 2 (M 2 C) T b n+ 1 2 - 1 0 j n+ 1 2 , b n+1 -b n = -∆tCe n+ 1 2 , ( 50
)
with e n+ 1 2 := 1 2 (e n + e n+1 ), b n+ 1 2 := 1 2 (b n + b n+1
), and a discrete current now defined through point values at the particle positions,

j n+ 1 2 = k q k v n+ 1 2 k • Λ 1 i (x n+ 1 2 k ) 1≤i≤N1 . (51)
As will be verified below, this scheme preserves exactly the total energy [START_REF] Umeda | A new charge conservation method in electromagnetic particle-in-cell simulations[END_REF]. One apparent difficulty is that, as the current source (51) involves the unknown velocities v n+1 k , it is a priori fully implicit. However, just as its original FD version [START_REF] Lapenta | Exactly energy conserving semi-implicit particle in cell formulation[END_REF], it has the remarkable property that it can be decomposed into a fully explicit particle pusher and a linearly-implicit field solver. In order to see it we rewrite the velocity kick in [START_REF] Villasenor | Rigorous charge conservation for local electromagnetic field solvers[END_REF] as

v n+ 1 2 k = 1 2 (v n k + v n+1 k ) = v n k + q k ∆t 2m k E n+ 1 2 h (x n+ 1 2 k ) + B n k v n+ 1 2 k
, where the matrix

B n k :=   0 B n k,3 -B n k,2 -B n k,3 0 B n k,1 B n k,2 -B n k,1 0   with B n k := B n h (x n+ 1 2 k ) is such that B n k v = v × B n h (x n+ 1 2 k ) for all v ∈ R 3 . Denoting η k = q k ∆t 2m k and A n k := I -η k B n k -1 = 1 1 + β 2   1 + β 2 1 -β 3 + β 1 β 2 β 2 + β 1 β 3 β 3 + β 1 β 2 1 + β 2 2 -β 1 + β 2 β 3 -β 2 + β 1 β 3 β 1 + β 2 β 3 1 + β 2 3   (52)
here with β := η k B n h (x

n+ 1 2 k
) as a short-hand notation, this gives

v n+ 1 2 k = A n k v n k + η k E n+ 1 2 h (x n+ 1 2 k
) , which allows us to decompose the fully implicit current (51) in two parts,

j n+ 1 2 = j n exp + j n imp [E n+ 1 2 h ], (53)
with explicit and implicit parts defined by the moment arrays

                 j n exp = k q k A n k v n k • Λ 1 i (x n+ 1 2 k ) 1≤i≤N1 , j n imp [E n+ 1 2 h ] = ∆t 2 k q 2 k m k A n k E n+ 1 2 h (x n+ 1 2 k ) • Λ 1 i (x n+ 1 2 k ) 1≤i≤N1 = ∆t 2 L n e n + e n+1 2 .
(54)

Here we have introduced a finite element "Lapenta" matrix [START_REF] Lapenta | Exactly energy conserving semi-implicit particle in cell formulation[END_REF],

L n := k q 2 k m k A n k Λ 1 j (x n+ 1 2 k ) • Λ 1 i (x n+ 1 2 k ) 1≤i,j≤N1 . (55)
Using ( 53) and (54) we then rewrite the discrete Maxwell system (50) as

M 1 + ∆t 2 4 0 L n -c 2 ∆t 2 (M 2 C) T ∆t 2 C I e n+1 b n+1 = M 1 -∆t 2 4 0 L n c 2 ∆t 2 (M 2 C) T -∆t 2 C I e n b n - ∆t 0 j n exp 0 , (56)
and summarize the above steps as follows.

Proposition 4.2. The ECSIM scheme (49)-( 51) is equivalently performed by the following steps:

(i) push the particles between t n and t n+ 1 2 as in (49),

(ii) compute the explicit current term (54) and the Lapenta matrix (55),

(iii) solve the linearly-implicit Maxwell system (56),

(iv) update the velocities as in (49), (v) push the particles between t n+ 1 2 and t n+1 as in [START_REF] Villasenor | Rigorous charge conservation for local electromagnetic field solvers[END_REF].

As claimed above, this scheme enjoys an exact energy conservation property.

Proposition 4.3. The ECSIM scheme (49)-( 51) preserves the discrete energy [START_REF] Umeda | A new charge conservation method in electromagnetic particle-in-cell simulations[END_REF].

Proof. In terms of finite element fields, the Maxwell solver (50) reads :

     E n+1 h -E n h , φ 1 h = c 2 ∆t B n+ 1 2 h , curl φ 1 h - ∆t 0 J n+ 1 2 h , φ 1 h , ∀φ 1 h ∈ V 1 h , B n+1 h -B n h = -∆t curl E n+ 1 2 h , (57)
with a current source defined by its moments :

J n+ 1 2 h , φ 1 h = k q k v n+ 1 2 k • φ 1 h (x n+ 1 2 k ) for all φ 1 h ∈ V 1 h . (58)
The energy variation from time-step n to time-step n + 1 is then the sum of three contributions:

B n+1 h 2 -B n h 2 2µ 0 = 1 µ 0 B n+1 h -B n h , B n+ 1 2 h = - 1 µ 0 ∆t B n+ 1 2 h , curl E n+ 1 2 h , 0 E n+1 h 2 -E n h 2 2 = 0 E n+1 h -E n h , E n+ 1 2 h = 1 µ 0 ∆t B n+ 1 2 h , curl E n+ 1 2 h -∆t J n+ 1 2 h , E n+ 1 2 h , k m k 2 (|v n+1 k | 2 -|v n k | 2 ) = k m k (v n+1 k -v n k ) • v n+ 1 2 k = k q k ∆tv n+ 1 2 k • E n+ 1 2 h (x n+ 1 2 k ) = ∆t J n+ 1 2 h , E n+ 1 2 h .
and those three terms cancel out.

Before we close this section, let us state the following result on the solvability of the implicit Maxwell linear system involved in the finite element ECSIM scheme.

Proposition 4.4. The linear system (56) is always invertible.

We refer the reader to Appendix A for a proof.

ChECSIM: a charge-and energy-conserving semi-implicit scheme

In the above sections, we have presented two FEM-PIC schemes which preserve either a proper set of discrete Gauss laws, or the discrete energy of the system. However, none of them preserves both quantities. On the other hand, it is possible to take a time-step free from CFL constraint with ECSIM's Maxwell solver but not with CCLF's. We wish to combine all the advantageous properties in a new scheme that keeps ECSIM's low computational cost, that is, a linear-implicit field iteration and an explicit pusher.

In ECSIM, when performing the discrete Ampère's equation [START_REF] Weibel | Spontaneously Growing Transverse Waves in a Plasma Due to an Anisotropic Velocity Distribution[END_REF], it is not possible to straightforwardly apply charge-conserving current deposition of the form of (44) without resorting to a non-linear resolution. Indeed, the trajectory of particles is unknown over the second half of the time-step, since only the first operation in [START_REF] Villasenor | Rigorous charge conservation for local electromagnetic field solvers[END_REF] has been performed at this point. To circumvent this obstacle, we propose a prediction/correction two-stepped scheme : in the first step we use a modified ECSIM scheme to compute a temporary solution :

(X n , V n , e n , b n ) → (X n+ 1 2 , Ṽ n+1 , ẽn+1 , b n+1 ).
Then we perform a corrective step,

(X n+ 1 2 , Ṽ n+1 , ẽn+1 , b n+1 ) → (X n+1 , V n+1 , e n+1 , b n+1 ),
in such a way that the resulting ChECSIM scheme is both Gauss-and energypreserving.

ChECSIM guess step

The first step mimics an ECSIM iteration, with a few modifications: the trajectories are now defined by :

x n+ 1 2 k -x n k = ∆t 2 v n k , (59)
ṽn+1 k -v n k = ∆t q k m k { Ẽn+ 1 2 h } n,- k + ṽn+ 1 2 k × B n h (x n+ 1 2 k ) , (60)
with short-hand notations :

ṽn+ 1 2 k := 1 2 (v n k + ṽn+1 k ) and Ẽn+ 1 2 h := 1 2 (E n h + Ẽn+1 h ),
and the curly brackets with a minus sign denoting a time average over the first half time-step :

{g} n,- k := 2 ∆t ˆtn+ 1 2 t n g x k (t) dt with x k (t) = x n k + (t -t n )v n k . (61)
The Maxwell equations become

   M 1 (ẽ n+1 -e n ) = ∆t c 2 (M 2 C) T b n+ 1 2 - 1 0 n+ 1 2 , b n+1 -b n = -∆tCẽ n+ 1 2 , (62)
with ẽn+ 1 2 := 1 2 (e n + ẽn+1 ) and b n+ 1 2 := 1 2 (b n + b n+1 ). Here the current term is n+

1 2 = k q k ṽn+ 1 2 k • {Λ 1 i } n,- k 1≤i≤N1 . ( 63
)
We observe that the main difference with ECSIM is that in several places we have used time-averages corresponding to the first half time-step, i.e.

{•} n,- k , instead of point values at x n+ 1 2 k
. Namely, the ChECSIM counter part of ( 49) is ( 59)-(60) and that of (51) is (63). This leaves the energy conservation untouched.

Proposition 5.1. The above guess step preserves the discrete energy (48), i.e.,

E h ( Ṽ n+1 , ẽn+1 , b n+1 ) = E h (V n , e n , b n ). ( 64
)
Proof. Similarly as in Proposition 4.3 we rewrite the Maxwell solver (62) in terms of finite element fields,

     Ẽn+1 h -E n h , φ 1 h = c 2 ∆t B n+ 1 2 h , curl φ 1 h -∆t 1 0 J n+ 1 2 h , φ 1 h , B n+1 h -B n h = -∆t curl Ẽn+ 1 2 h , (65)
with a current term now defined as

J n+ 1 2 h , φ 1 h = k q k ṽn+ 1 2 k • {φ 1 h } n,- k , for all φ 1 h ∈ V 1 h . ( 66
)
The same computations then apply, the only difference lying in the fact that point values of the electric field are replaced by time-averages over trajectories.

For the kinetic energy this gives

k m k 2 (|ṽ n+1 k | 2 -|v n k | 2 ) = k m k ṽn+ 1 2 k • (ṽ n+1 k -v n k ) (67) = k q k ∆tṽ n+ 1 2 k • { Ẽn+ 1 2 h } n,- k (68) = ∆t J n+ 1 2 h , Ẽn+ 1 2 h (69)
and this cancels out the electric energy update, similarly as in Prop. 4.3.

As in the ECSIM scheme, the guess step above may be decomposed as an explicit particle pusher and a linearly implicit field solver. To picture that, we first rewrite the discrete Maxwell system (62) as a block-matrix equation,

M 1 -c 2 ∆t 2 (M 2 C) T ∆t 2 C I ẽn+1 b n+1 = M 1 c 2 ∆t 2 (M 2 C) T -∆t 2 C I e n b n - ∆t 0 n+ 1 2 0 , (70)
and similarly as for ECSIM, we rewrite the mid-point velocity as

ṽn+ 1 2 k = 1 2 (v n k + ṽn+1 k ) = A n k v n k + η k { Ẽn+ 1 2 h } n,- k , with η k = ∆tq k 2m k and A n k = A n k (b n , X n+ 1
2 ) the matrix defined in (52). This allows us to decompose the current in two parts,

n+ 1 2 = n exp + n imp [ Ẽn+ 1 2 h ], (71)
with explicit and implicit parts defined by the moment arrays

                 n exp = k q k A n k v n k • {Λ 1 i } n,- k 1≤i≤N1 , n imp [ Ẽn+ 1 2 h ] = ∆t 2 k q 2 k m k A n k { Ẽn+ 1 2 h } n,- k • {Λ 1 i } n,- k 1≤i≤N1 = ∆t 2 Ln e n + ẽn+1 2 . ( 72
)
and a time-averaged analog to our finite element Lapenta matrix (55),

Ln = k q 2 k m k A n k {Λ 1 j } n,- k • {Λ 1 i } n,- k 1≤i,j≤N1 . (73)
Plugging this decomposition into (70) then gives

M 1 + ∆t 2 4 0 Ln -c 2 ∆t 2 (M 2 C) T ∆t 2 C I ẽn+1 b n+1 = M 1 -∆t 2 4 0 Ln c 2 ∆t 2 (M 2 C) T -∆t 2 C I e n b n - ∆t 0 n exp 0 , (74)
which leads us to the following computational steps.

Proposition 5.2. The ChECSIM guess step (59)-( 60)-( 66) is equivalently performed by the following steps:

(i) push the particles between t n and t n+ 1 2 as in (59),

(ii) compute the explicit current term (72) and the new Lapenta matrix (73),

(iii) solve the linearly-implicit Maxwell system (74),

(iv) update the velocities as in (60).

Finally we verify that the arguments of Proposition 4.4 readily apply to the modified Maxwell system (74). Again, we refer to Appendix Appendix A for a detailed proof.

Proposition 5.3. The linear system (74) is always invertible.

Remark 5.4. During this guess step, we choose to modify ECSIM with time averages in the Lorentz force (60), and current term (63). This is not out of necessity : one could as well perform the predictive step with the regular ECSIM scheme from Section 4.2. Nonetheless, using time averged terms has two advantages :

(i) Since the error in Gauss law appears when using a point-valued current deposition instead of trajectory integrals, our choice results in a better guess field, in the sense that it has a smaller charge error. (ii) In order to perform the full correction, the time averages have to be computed on both half time-steps. Since the first one is available after pushing the particles in the predictive step (59), it is simpler to use it at this point.

ChECSIM correction step

As previously remarked, the guess step above does not preserve the Gauss laws, which is due to the fact that the linear decoupling procedure makes it impossible to average over the trajectories over the full time-step. However, since we now have a guess velocity for second half step, we may correct the deposited current so that it becomes charge-conserving. Specifically, decompose the guess current source (66) in two parts,

J n+ 1 2 h = J n,- h + J n,- h 2 , ( 75
)
where J n,- h is a charge-conserving term defined as an average current density,

J n,- h , φ 1 h = k q k v n k • {φ 1 h } n,- k = 2 ∆t ˆtn+ 1 2 t n J h (t), φ 1 h dt, (76)
see ( 33), (61), and a remainder term

J n,- h , φ 1 h = k q k ṽn+1 k • {φ 1 h } n,- k (77)
that is not charge-conserving. The aim of the correction step is then to replace this term by a conserving current for the second half time-step, defined as

J n,+ h , φ 1 h = k q k ṽn+1 k • {φ 1 h } n,+ k , (78)
where the curly brackets with a plus sign denote time averages over the second half time-step,

{g} n,+ k := 2 ∆t ˆtn+1 t n+ 1 2 g x k (t) dt with x k (t) = x n+ 1 2 k +(t-t n+ 1 2 )ṽ n+1 k . (79)
Once again, the velocities of the particles are constants on both halves of the time step, so in all current expressions above, the term v n k can be placed either inside or in front of curly brackets {•} n,- k , and same goes for ṽn+1 k and {•} n,+ k . Once this current correction is performed, the energy conservation of Proposition 5.1 is lost, since we change the electric field from ẽn+1 to e n+1 , and

E h ( Ṽ n+1 , e n+1 , b n+1 ) = E h ( Ṽ n+1 , ẽn+1 , b n+1 ).
To restore it, we rescale the velocities. The equations are then as follows.

           x n+1 k -x n+ 1 2 k = ∆t 2 ṽn+1 k , M 1 (e n+1 -ẽn+1 ) = - ∆t 2 0 j n,+ -n,-, v n+1 k = λṽ n+1 k . ( 80
)
Here the current sources correspond to the moments of ( 77) and (78), namely

j n,+ -n,-= k q k ṽn+1 k • {Λ 1 i } n,+ k -{Λ 1 i } n,- k 1≤i≤N1
, and λ > 0 is a global correction parameter given by the requirement that

0 2 E n+1 h 2 - Ẽn+1 h 2 = k m k 2 |ṽ n+1 k | 2 -|v n+1 k | 2 = (1-λ 2 ) k m k 2 |ṽ n+1 k | 2 ,
(81) which yields

λ := 1 - 0 ( E n+1 h 2 - Ẽn+1 h 2 ) k m k |ṽ n+1 k | 2 1 2 . ( 82
)
Below we will give a sufficient condition for this parameter to be well-defined. Assuming this for now, we can state the main properties of our scheme.

Theorem 5.5. The ChECSIM scheme defined by the guess step (59)-( 60)-( 61) and the correction step (80)-(82) preserves the discrete Gauss laws (45)-( 46) as well as the discrete energy [START_REF] Umeda | A new charge conservation method in electromagnetic particle-in-cell simulations[END_REF].

Proof. Both properties are easily verified. By construction, the final electric field E n+1 h satisfies a discrete Ampère equation

E n+1 h -E n h = c 2 ∆t curl h B n+ 1 2 h -∆t 1 0 J n+ 1 2 h
with a current density corresponding to

J n+ 1 2 h = J n+ 1 2 h + 1 2 J n,+ h - J n,- h = 1 2 J n,- h + J n,+ h = 1 ∆t ˆtn+1 t n J h (t) dt
where again, J h (t) is the time-continuous current density defined by [START_REF] Lapenta | Exactly energy conserving semi-implicit particle in cell formulation[END_REF]. Similarly as for the charge-conserving leap-frog scheme (Prop. 4.1), this shows that the discrete sources satisfy a continuity equation [START_REF] Taflove | Computational electrodynamics : the finite-difference timedomain method[END_REF] and since the discrete Maxwell equations are solved by structure-preserving finite elements, this is enough to show that the Gauss equations are preserved. The energy conservation follows from Prop. 5.1 and the requirement (81).

Remark 5.6. In contrast to Marder's pseudo-current correction method [START_REF] Marder | A method for incorporating Gauss' law into electromagnetic PIC codes[END_REF], the current correction in (80) is both local and exactly charge-conserving.

We now provide a sufficient condition for the existence of a correction parameter λ. It will be convenient to denote by

v n m := k m k |v k | 2 1 2
the weighted 2 norm involved in the kinetic energy.

2∆t γ n m ≤ v n m , (83)
then the field correction satisfies

0 E n+1 h 2 - Ẽn+1 h 2 ≤ C∆t 2 ṽn+1 2 m ∇ Ē 1,∞ (84)
with

C = q m (1 + √ 5) and ∇ Ē 1,∞ := i,j ∂ j Ēi 2 L ∞ (Ω) 1 2 with Ē := E n+1 h + Ẽn+1 h 2 .
In particular, the correction step (80)-( 82) is well defined for ∆t < ∆t * , where

∆t * := min v n m 2 γ n m , q m (1 + √ 5) ∇ Ē 1,∞ -1 2 . ( 85
)
We refer the reader to Appendix B for a proof.

Remark 5.8. Velocity rescaling techniques are commonly used in molecular dynamics simulations to preserve the structural stability of the discrete systems, and fine methods based on thermostat or variational principles have been devised that allow a localized treatment of the particles, see e.g. [START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Maggs | Adding an energy-like conservation law to the leapfrog integrator[END_REF]. In (80) we have considered a global rescaling for simplicity, as it gives satisfactory results for our numerical test cases. For real-life problems however, one would need to apply a finer strategy, either by adapting thermostating methods, or by using a localized velocity rescaling v n+1 k = λ k ṽn+1 k with a rescaling factor λ k that is close to 1 in the regions where the current correction is small. In practice a convenient approach is to consider a factor of the form

λ 2 k := 1 -γ J n,+ h - J n,- h (x n+1 k ) 2 , ( 86
)
where γ is determined so that energy conservation is restored, yielding

γ = 0 ( E n+1 h 2 - Ẽn+1 h 2 ) k m k |ṽ n+1 k | 2 J n,+ h - J n,- h (x n+1 k ) 2 . ( 87
)
In the numerical section of the paper, we investigate this local velocity correction in the diode test-case (see Section 6.2).

Numerical experiments

We now investigate the validity of the new ChECSIM algorithm, and compare its properties with the charge-conserving leap-frog (CCLF) and the EC-SIM schemes, recalled in Sections 4.1 and 4.2. These numerical experiments are performed with a 2D research code called Coffee which implements Lagrange-Nédélec finite elements on triangular meshes. We note that in a 2D setting the de Rham sequence is reduced. Corresponding to the choice of a weak Ampère equation, it reads

V 0 h V 1 h V 2 h grad curl (88)
see, e.g. [START_REF] Campos | Compatible Maxwell solvers with particles II: conforming and non-conforming 2d schemes with a strong Faraday law[END_REF]. Here,

V 2 h = P p-1 (T h ) := {u ∈ L 2 (Ω) : u| T ∈ P p-1 (T ) for T ∈ T h } (89)
corresponds to broken (discontinuous) piecewise polynomials of degree p -1,

V 1 h = u ∈ H 0 (curl, Ω) : u| T ∈ P p-1 (T ) 2 + -y x P p-1 (T ) for T ∈ T h (90)
is the (first-kind) Nédélec finite element space of maximal degree p, and

V 0 h = L p (T h ) := P p (T h ) ∩ C 0 (Ω) (91)
is the Lagrange finite element space of degree p ≥ 1. We refer e.g. to [START_REF] Boffi | Mixed finite element methods and applications[END_REF][START_REF] Girault | Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms[END_REF][START_REF] Ern | Theory and Practice of Finite Elements[END_REF] for a complete description of these spaces. Before describing our test cases, we may point out an interesting result from our simulations. As we have seen above, in the ChECSIM algorithm the velocity correction step (80) may not be feasible if the kinetic energy is lower than the field energy correction, and when that is not the case Proposition 5.7 guarantees that an adaptation of the time-step allows to have a well-posed correction, such as the one proposed in Appendix E. A first numerical observation is that for all the runs presented in this section, the correction step was always possible and we never had to adapt the time-step.

Weibel instability

In order to certify the validity of our alogrithms, we first simulate the electromagnetic Weibel instability [START_REF] Weibel | Spontaneously Growing Transverse Waves in a Plasma Due to an Anisotropic Velocity Distribution[END_REF]. This phenomenon consists in the exponential growth of the self-induced magnetic field B = B z in a plasma where the velocity distribution of electrons is anisotropic, with higher kinetic energy in one of the physical dimensions. In our 2D2V simulation code, we reproduce the experiment documented in [START_REF] Crouseilles | Hamiltonian splitting for the vlasov-maxwell equations[END_REF] and [START_REF] Kraus | GEMPIC: geometric electromagnetic particle-in-cell methods[END_REF].

Here, we consider Equations ( 1)-( 2) in a non-dimensional setting where the speed of light is set to c = 1, as well as the electric permittivity of vacuum are considered to be fixed over time. At initial time, the electrons are distributed uniformly in space and have an anisotropic Maxwellian velocity distribution

f 0 (x, v) = 1 2πσ x σ y exp 1 2 v 2 x σ 2 x + v 2 y σ 2 y ( 92
)
with σ x = 0.02 √ 2 and σ y = 0.02 √ 6, hence a higher thermal velocity in y. The initial magnetic field B 0 h ∈ V 2 h is computed as an approximation of

B 0 (x) = β cos(kx) ( 93
)
with β = -10 -6 , and the initial electric field E 0 h ∈ V 1 h is obtained from the initial particle distribution by solving the corresponding finite element Poisson problem and taking the gradient of the resulting discrete potential φ 0 h ∈ V 0 h , in order to satisfy the discrete Gauss law [START_REF] Stock | Investigation of the purely hyperbolic maxwell system for divergence cleaning in discontinuous galerkin based particle-in-cell methods[END_REF] for n = 0.

We recall that the Weibel instability corresponds to a magnetic mode of the form B k (t, x) = B k e i(kx-ωt) , with ω the most unstable solution of the dispersion relation [START_REF] Weibel | Spontaneously Growing Transverse Waves in a Plasma Due to an Anisotropic Velocity Distribution[END_REF]

D(ω, k) := ω 2 -k 2 + σ y σ x 2 -1 - σ y σ x 2 Φ ω σ x k ω σ x k = 0, (94)
i.e., the solution with largest positive imaginary part γ := (ω) > 0. Here Φ(z) = exp(-1/2z 2 ) ´z -i∞ exp(1/2ξ 2 )dξ, and the analytic growth rate of the instability is γ ≈ 0.02784, see eg [START_REF] Kraus | GEMPIC: geometric electromagnetic particle-in-cell methods[END_REF]. We shall verify whether this rate is obtained in the numerical simulations. We simulate the plasma during a time T = 450. The domain is discretized by a conforming triangular mesh T h composed of 74 elements of typical diameter h = 0.5585. The order of the finite element spaces (88) is p = 2, and the number of numerical macro-particles used in the PIC approximation is N = 2.96 • 10 5 . The time-step used is ∆t = 0.05, corresponding to approximately 80% of the maximum stable time-step for the explicit leap-frog scheme for the resulting Finite Element Maxwell solver.

On Figure 1, we present the time evolution of the magnetic energy 1 2 B n h 2 , computed with each of the three schemes presented in this article. Then, we present on Figure 2 the time evolution of the relative error on the conservation of the total discrete energy [START_REF] Umeda | A new charge conservation method in electromagnetic particle-in-cell simulations[END_REF],

|E h (V n , e n , b n ) -E h (V 0 , e 0 , b 0 )| |E h (V 0 , e 0 , b 0 )| , ( 95
)
and the absolute error on the electric Gauss law [START_REF] Kraus | GEMPIC: geometric electromagnetic particle-in-cell methods[END_REF],

(M 1 G) T e n - 1 0 r n 2 . ( 96
)
This latter error is divided by the typical physical scale of the charge density, yielding a non-dimensional value which expresses the error on Gauss' law relatively to its characteristic value. Finally, on the left of Figure 3, we present the evolution of ChECSIM's velocity rescaling factor λ, while on the right of the same figure is the time evolution of the momentum conservation error

|P h (V n , e n , b n ) -P h (V 0 , e 0 , b 0 )|, ( 97
)
where the discrete total momentum is defined by In Table 1, we show the execution times of the schemes relative to the CCLF one.

P h (V n , e n , b n ) = N k=1 m k v n k + ˆΩ 0 E n h × B n h dx. ( 98
)
On Figure 1, the exponential growth of the magnetic energy is clearly visible, and it is comparable for the three schemes. Since this energy depends on the square of the norm of B, its growth rate is twice that of the field. A linear fit of the ChECSIM curve between t = 250 and 310 yields a growth rate of γ ≈ 0.024, in good agreement with the analytical value.

On the left of figure 2, we observe that the relative energy conservation (95) is close to machine accuracy for ECSIM and ChECSIM (a slight increase is visible for ChECSIM, but still on the order of roundoff errors). With CCLF the energy is also rather stable, but the conservation error is higher by about 6 orders of magnitude. On the right of the same figure, we observe that both ChECSIM and CCLF preserve Gauss's law with very high accuracy. For ECSIM, we observe that the error on the electric Gauss law is significantly larger, again by about 5 orders of magnitude.

On Figure 3, left, we observe that throughout the simulation, the rescaling factor λ stays very close to one, deviating by an order of 10 -7 maximum. On the right, one sees that all methods behave similarly regarding momentum conservation errors.

In Table 1, we notice that ECSIM is more than two times faster than ChEC-SIM. Our new scheme is also 140% slower than CCLF. This is not unexpected, since each two-stepped iteration in ChECSIM loops twice on all the particles, and computes time averages on the trajectories, which is by far the most timeconsuming process. Nevertheless, we shall see in the next numerical experiment that ChECSIM has other advantages that might temper this drawback, allowing quicker runtime.

This provides us with a first numerical validation of the basic properties of the ChECSIM scheme, namely the conservation of energy and charge for an isolated system.

Electron beam in a simple diode

To assess the properties of the different schemes in a more realistic setting where the system is no longer isolated, we now consider a second test-case consisting of an electron beam in a simple, academic diode.

Following [START_REF] Stock | Investigation of the purely hyperbolic maxwell system for divergence cleaning in discontinuous galerkin based particle-in-cell methods[END_REF], we consider a diode modeled by a square domain Ω of dimensions L x = L y = 0.1 m. In this domain, electrons are emitted at a cathode lying in the face of abscissa x = 0, and accelerated towards the anode, lying in the face of abscissa x = 0.1, by a constant exterior electric field E ext = -(U ext /L x )x with an operating voltage U ext = 10 5 V (see Appendix D for the modeling of exterior fields in the ECSIM and ChECSIM schemes). More precisely, the electrons are injected in the domain through a portion of the cathode : a strip centered at y = 0.05 and of fixed width L inj = 0.03786 m. The current density at the emission surface is I = -5536 A.m -2 , corresponding to 75% of the space-charge limiting current density.

The electron beam is approximated by numerical macro-particles emitted at rate r = 9000 ns -1 . In this 2D model, each numerical particle carries a linear density w of electrons in the z dimension. Thus, the the total charge injected during a time ∆t reads Q inj (∆t) := IL z L inj ∆t = r∆twq e L z where L z is the height of the 3D diode, q e is the charge of an electron and w is the uniform weight of the 2D macro-particles, which has the physical dimension of a linear density. We then obtain w ≈ 1.45

• 10 8 m -1 .
The domain is discretized by an unstructured conforming triangular mesh T h , composed of 2158 elements, the typical diameter of a cell being h = 0.005 m. The electromagnetic fields are approached by the structure-preserving finite element complex (88), with a precision order p = 2.

We simulate the flow of electrons within the diode during a time of T = 50 ns, using the three schemes presented in this article : the charge-conserving leapfrog scheme (CCLF), ECSIM and ChECSIM. We note that for this test-case the energy is not expected to remain constant, as the system is not isolated. Up to T stat = 2 ns approximately, the system is in a transitory regime : the electronic current sets up in the initially empty diode. Then, from T stat to T , the phase-space density of particles and the electromagnetic fields are approximately constant, and the system is in a stationary regime. These two phases are clearly visible on the left panel of Figure 4, where we plot the time evolution of the total energy (48) in the diode.

Specifically, Figure 4 displays two sets of energy curves: on the left panel the runs are obtained using a constant time-step ∆t CFL = 1.25 • 10 -3 ns, corresponding to 80% of the CFL time-step of the leap-frog Maxwell solver. This constraint is needed for the stability of the CCLF scheme. On the right panel, CCLF ECSIM ChECSIM / CCLF 1 0.6637 1.3977 the curves correspond to ECSIM and ChECSIM runs using a larger time step. Indeed these methods are not constrained by an explicit CFL condition, so it is interesting to compare their results when the time step is increased. In order to have a similar accuracy for the transitory regime, we have kept the time-step ∆t CFL up to time T stat , and then switched to a longer ∆t stat for the stationary phase.

Then, on Figures 5 and6 we plot the profiles of the numerical electron beam and of the E x field at an itermediate time t = 16 ns. Again, several values are used for the time-step ∆t stat , as described above.

On the left of Figure 7 we represent the time evolution of λ, the global velocity rescaling factor of ChECSIM. Once more, different values are used for the time-step ∆t stat starting from time T stat , so we plot the three corresponding values for λ. On the right of the same figure, we plot the momentum conservation error of all three schemes with the smallest time-step.

On Figure 8, we investigate the local correction factor λ k defined by ( 86)-(87). At the position of each particle k, the value of |1 -λ k | is represented with a color scale. On the left is the state of the beam in the diode as we described it above, with a uniform exterior electric field. On the right is a similar diode in which we introduce a slight modification : the accelerating electric field E ext is not uniform in space, but rather reads as follows,

E ext (x, y) =      -(U ext /L x )x, y < 0.050 m, -[1 + 20(y -0.050)](U ext /L x )x, 0.050 m y 0.055 m, -2(U ext /L x )x, y > 0.055 m.
(99) That is, the accelerating field is continuous, has constant magnitude in the lower part of the diode, then an affine growth in magnitude from the center towards the upper part, where it is constant again, at twice the previous amplitude. The aim of such localized variation of the field is to artificially amplify the prediction errors on the current, to verify that the local correction is indeed preeminent in that region.

Lastly, in Table 2 we compare the execution time for the different schemes. For ECSIM and ChECSIM, when longer time-steps are used, we also present the corresponding runtime.

On the left of Figure 4, we observe a stationary total energy from t = T stat to T with all three schemes. On the right, one notices that for the longest choosen time-step ∆t stat = 10∆t CFL , a significant fraction of the energy is lost by ECSIM. Although this may seem surprising, we emphasize that it does not contradict the theoretical properties of the energy conserving scheme, since they are only guaranteed for an isolated system. Here we are not in this situation : particles enter and leave the domain, they gain energy from the external field, and absorbing boundary conditions on the horizontal boundaries allow some of the field energy to be dissipated. For the intermediate time-step ∆t stat = 2∆t CFL (not shown here) the energy profile is similar for ECSIM and ChECSIM, but we have observed that the total energy also decreases significantly with ECSIM in longer simulations with T = 50 ns. For ChECSIM we observe that Here t is expressed in nanoseconds while the energy is normalized. On the right panel the colored zone starting at t = 2 ns emphasizes the change of time-step during the second part of the simulation, as described in the text.

the energy is well preserved during the stationary phase, even when longer time-steps are used. We should point out that this behaviour is not a direct consequence of the sole energy conservation properties of the scheme since the system is not isolated, but rather an indication of the improved stability of the ChECSIM algorithm.

We now turn to the profiles of the particle beam and the electric field on Figures 5 and6. For the smallest time-step ∆t stat = ∆t CFL , the CCLF profiles on panels 5a and 6a show a good qualitative behavior, very similar to those obtained with ChECSIM (not shown here). Here we note that the moiré patterns on the beams are due to the pseudo-random sampling of particles in the injection zone. The profiles for the ECSIM run in panels 5b and 6b are close, despite some subtle filamentation developping in the upper half of the particle beam (a typical result of small space-charge errors due to the lack of charge-conservation, see e.g. [START_REF] Barthelmé | Numerical charge conservation in particle-incell codes[END_REF]), and a slightly stronger electric field at the cathode.

For the time-step ∆t stat = 2∆t CFL , the ChECSIM profiles in panels 5c and 6c are virtually identical, but in the ECSIM run the numerical artifacts grow stronger: the filamentation of the beam is more pronounced in panel 5d, and the corresponding field in panel 6d is now higher at the cathode by around 15%. With the longest time-step ∆t stat = 10∆t CFL , the profiles obtained with ChECSIM on panels 5e and 6e remain very close to the previous ones, despite a higher level of noise visible in the field. This is by far not the case with the ECSIM run, where the accumulation of charge errors leads to significant errors: on panel 5f we see that the particle beam is badly scattered, and on panel 6f the field takes clearly unphysical values at the cathode, and is also disrupted by strong spurious oscillations close to the anode.

Turning to Figure 7, we notice that the global rescaling parameter λ remains close to one, even when using longer time-steps. In these latter cases, the rescaling factor is bigger, with a deviation from unity attaining an order 10 -5 while it remains below 10 -7 for short time-step. This is expected since a longer time-step comes with a bigger discrepancy between pointwise current and timeaveraged current, hence a bigger correction is needed to catch this up. On the left of the same figure, one verifies that all three schemes behave similarly on the momentum conservation : once the stationnary regime is attained, the total momentum stays roughly constant. Now, when the local correction factor λ k defined by (86)-( 87) is used, we see on the left pannel of Figure 8 that the values of the local factors are quite uniformly distributed in the beam, around an average value which seems to be close to that of the global factor. On the right panel where the case of the nonuniform accelerating electric field is presented, we observe that the less regular the electric field is, the more important the local correction is, which legitimates this method as a finer strategy than the global correction. Finally, although the sufficient criterion of Proposition 5.7 no longer applies to the local factor, it was never necessary in the expermient to dynamically adapt the time-step for all factors λ k to be well-defined.

To finish with our observations, we compare the runtimes. For similar timesteps, ChECSIM is slower, as expected. Nevertheless, its good stability properties allow to take longer time-steps, where it becomes the quickest method.

As a result, we find that this test-case completes the numerical evidence of Section 6.1, and confirms that the combined energy and charge-conserving properties of the ChECSIM algorithm lead to an improved numerical stability for long-time simulations.

CCLF

Conclusion

In this article, we have reformulated the ECSIM algorithm in a structurepreserving finite element setting and we have used the properties of this general discrete framework to complement the original energy-conserving method with an explicit, local charge-conserving correction step based on FEM-PIC conservative current deposition techniques. The resulting ChECSIM scheme is still linearly implicit in the fields and explicit in the particles, and it now preserves exactly the total energy and the discrete charge in a sense that is compatible with the structure-preserving finite element framework. A condition for our correction step to be well-defined is that the kinetic energy is higher than the field correction energy, and we have provided a rigorous criterion on the time step, under which this condition is always fulfilled. In practice one may need to dynamically adapt the time-step to meet this criterion, but this was not the case in the numerical experiments performed in this article. Furthermore, we proposed to adapt locally the particle correction, to avoid situations where velocities are rescaled in parts of the domain where there is little error on the current in the guess step. The numerical experiments conducted in this paper have allowed us to exhibit a clear improvement of the stability properties of the new scheme, in particular for configurations involving large time steps and a model system that is not isolated from exterior energy sources.

Appendix B. Proof of the well-posedness of the rescaling factor λ

Here we give a proof of Proposition 5.7 stating that there exist a time-step ∆t small enough for the definition (82) of λ to be well-posed.

From (60) we compute Assuming that the finite element field Ē is Lipschitz, the first term is bounded by and similarly for the second one,

v n 2 m -ṽn+1 2 m = -∆t k m k γ n k • (ṽ n+1 k + v n k ) ≤ ∆t γ n m ( ṽn+1 m + v n m ) ≤
ˆtn+1 t n+ 1 2 Ē(x k (t)) -Ē(x
ˆtn+ 1 2 t n Ē(x n+ 1 2 k ) -Ē(x k (t)) dt ≤ ∆t 2 2 |v n k | ∇ Ē 1,∞ .
The net electric energy variation induced by the correction (80) reads then :

E n+1 h 2 - Ẽn+1 h 2 ≤ ∆t 2 0 k q k (|ṽ n+1 k | 2 + |ṽ n+1 k ||v n k |) ∇ Ē 1,∞ ≤ ∆t 2 0 q m ṽn+1 m ( ṽn+1 m + v n m ) ∇ Ē 1,∞ ≤ ∆t 2 0 q m ṽn+1 2 m (1 + √ 5) ∇ Ē 1,∞
where we have used (B.1). This shows (84) and ends the proof.

Appendix C. Modelling of absorbing boundary conditions

Absorbing boundaries may be modelled by decomposing the boundary Γ := ∂Ω into a part Γ M with metallic boundary conditions and an part Γ A with absorbing Silver-Müller conditions [START_REF] Assous | Mathematical foundations of computational electromagnetism[END_REF], so that ( 13) is replaced by

n × E = 0 on Γ M -n × (n × cB) on Γ A . (C.1)
Using this boundary condition in the integration by parts, the weak Ampère equation [START_REF] Campos | Compatible Maxwell solvers with particles II: conforming and non-conforming 2d schemes with a strong Faraday law[END_REF] becomes

∂ t E(t), φ 1 -B(t), curl φ 1 + c n × E, n × φ 1 Γ A = -J (t), φ 1 (C.2)
for all φ 1 ∈ V 1 , with an electric field now sought for in the space C 1 ([0, T ]; V 1 ) with

V 1 = H 0,Γ M (curl, Ω) = {u ∈ H(curl, Ω) : n × u = 0 on Γ M }.
Accordingly, the finite element space V 1 h should have non-homogeneous boundary conditions on Γ A , and the time-continous matrix Ampère equation ( 29) becomes

M 1 d dt e -(M 2 C) T b + cT Γ A e = -j (C.3)
where T A is a boundary mass-matrix involving the tangential traces on the absorbing boundary Γ A ,

T Γ A = n × Λ 1 i , n × Λ 1 j Γ A 1≤i,j≤N1 . k = 1 2 (v n k + v n+1 k ) = v n k + q∆t 2m (E n+ 1 2 h + E ext )(x n+ 1 2 k ) + B n k v n+ 1 2 k with a matrix B n k such that B n k v = v × (B n h + B ext )(x n+ 1 2 k
). Therefore we need to modify the matrix A n k accordingly,

A n k := I -η k B n k -1
and to decompose the fully implicit current (51) into

j n+ 1 2 = j n exp + j n imp [E n+ 1 2 h + E ext ]
with an explicit part j n exp defined similarly as in (54), and an implicit part which now involves the exterior field. Namely, we have

                 j n exp = k q k A n k v n k • Λ 1 i (x n+ 1 2 ) 1≤i≤N1 , j n imp [E n+ 1 2 h + E ext ] = ∆t 2 k q 2 k m k A n k (E n+ 1 2 h + E ext )(x n+ 1 2 k ) • Λ 1 i (x n+ 1 2 ) 1≤i≤N1 = ∆t 2
L n e n + e n+1 2 + e ext .

(D.1) where the matrix L n is defined as in (55), with the proper correction on the magnetic matrices. In particular, the linearly-implicit Maxwell system (56) becomes 2) where we see that the exterior electric field also appears as an additional source term. In the ChECSIM guess step, the modification is exactly the same.

M 1 + ∆t 2 4 0 L n -c 2 ∆t 2 (M 2 C) T

Appendix E. Dynamic adaptation of the time-step

In the ChECSIM algorithm, the velocity correction step (80) is only welldefined when the kinetic energy is larger than the variation of the field energy due to the charge-conserving current correction. Although this was always the case in the numerical experiments conducted for this article, it may no longer be true when the particles have a low velocity or when the field correction is important. Since Proposition 5.7 guarantees that the correction step is wellposed for a time-step small enough, in practice one may use a dynamic recursive splitting of the time-step, such as the following one.

Proposition 5 . 7 .

 57 Let γ n k terms in the guess step (60). If the time-step is such that

Figure 1 :

 1 Figure 1: Weibel test-case : time evolution of the field energy using ChECSIM (left) and magnetic field B at t = 349.22ω -1 pe (right).

Figure 2 :

 2 Figure 2: Weibel test-case : time evolution of (left) the energy conservation relative error (95) and (right) the electric Gauss's law absolute error (96).

Figure 3 :

 3 Figure 3: Weibel test-case : (left) time evolution of |λ -1|, where λ is the global rescaling factor (82) of ChECSIM, (right) errors on the momentum conservation (97) in arbitrary units.

Figure 4 :

 4 Figure 4: Diode test-case : time evolution of the total energy (48) for two values of ∆tstat.Here t is expressed in nanoseconds while the energy is normalized. On the right panel the colored zone starting at t = 2 ns emphasizes the change of time-step during the second part of the simulation, as described in the text.

Figure 5 :

 5 Figure 5: Diode test-case: positions of the numerical particles showing the beam profile at t = 16 ns for various schemes and time-steps. The values on the x and y axes are in meters.

Figure 6 :

 6 Figure 6: Diode test-case: profile of the Ex field at t = 16 ns for various schemes and timesteps. The values on the x and y axes are in meters, those on the z axis are in volt/meter.

Figure 7 :

 7 Figure 7: Diode test-case : (left) time evolution of |λ -1| where λ is the global rescaling factor (82) of ChECSIM, and (right) momentum conservation errors (97), for ∆tstat = ∆t CFL , expressed in kg.m.s -1 .

Figure 8 :

 8 Figure 8: Diode test-case with local velocity rescaling: particle beam at t = 5 ns in a diode where (left) the exterior electric field is uniform and (right) where it is non-uniform, see (99). The colors represent the values of |1 -λ k | at the particle positions, where λ k is the local correction factor (86)-(87) of the particles. In both cases, the timestep is ∆tstat = 10∆t CFL .

2 m 1 ))

 21 where we have used (83) and Young's inequality in the last equation. This gives a bound for the guess velocities ṽn+1 Integrating next the field correction in (80) against Ē and using the definition of the respective current terms yields then set δ n k ( Ē) := { Ē} n,+ k -{ Ē} n,- k . By definition of the time averages (61) and (79) we next compute that ∆t -Ē(x k (t)) dt.

Algorithm 1

 1 ChECSIM with dynamic time-step splittingfor n = 0, • • • , #(iterations) -1 do τ = ∆t s = 0 while s < ∆t do τ = min(τ, ∆t -s)Perform guess step (59)-(60)-(63) with time-step τ Push the particles and correct the electric field (

Table 1 :

 1 Weibel test-case : comparison of the execution times, relative to CCLF's.

  ECSIM ChECSIM ∆t stat = ∆t CFL ∆t CFL 2∆t CFL 10∆t CFL ∆t CFL 2∆t CFL 10∆t CFL

	/ CCLF	1	0.7397 0.4673	0.1868	1.086	0.6679	0.3193

Table 2 :

 2 Diode test-case : comparison of the execution times, relative to CCLF's.

  80) with time-step τ if λ is well-defined by (82) then Perform velocity correction step in (80) s = s + τ if s + 2τ ≤ ∆t then

	τ = 2τ
	end if
	else
	Rollback particles and fields to n∆t + s
	τ = τ /2
	end if
	end while
	end for

=

1, and the electron charge and mass q e = -1 and m e = 1. The computational domain Ω is a periodic rectangle of length L x = 2π/k, with k = 1.25, and width L y = L x /4. The domain is occupied by a plasma composed of electrons and a uniform neutralizing background of ions. The ions

Acknowledgements

The authors would like to thank Benedikt Perse, Katharina Kormann and Eric Sonnendrücker for fruitful discussions.

Appendix A. Well-posedness of the Maxwell implicit linear system Here we detail the invertibility of Maxwell linear system 56. The matrix on the left-hand side has a block structure

where we have removed scalar parameters such as c or ∆t for the sake of simplicity, since they do not play any role in the argument below. This matrix is invertible if and only if there exists matrices P, Q, R, S such that

which is equivalent to

Hence, the invertiblity of M is equivalent to the invertibility of the matrices

Here the mass matrix M 1 is obviously positive definite, since for any function w h = 0 of V 1 h with coeficient vector denoted by w, one has

Finally, using (55) and (52) we write

) and

). This shows that the Lapenta matrix L n is also positive semi-definite, so that both U and V are positive definite. They are therefore invertible, which proves that the corresponding implicit Maxwell system (56) is well-posed.

Note that all the arguments in this proof apply to the slightly modified system (74) used in ChECSIM.

Note that if Γ A = ∅, then T Γ A = 0 and we recover Equation [START_REF] Hockney | Computer simulation using particles[END_REF]. Using an semi-implicit time-stepping for the absorbing boundary terms, the Maxwell solver in the leap-frog scheme then becomes

.

(C.4) Similarly, the fully implicit Maxwell solver [START_REF] Weibel | Spontaneously Growing Transverse Waves in a Plasma Due to an Anisotropic Velocity Distribution[END_REF] in the ECSIM scheme takes the form

and the corresponding linearly implicit formulation (56) becomes

(C.6) The same applies for the ChECSIM guess step, where the linearly implicit formulation (74) becomes

(C.7) Finally in order to obtain an electric field e n+1 satisfying the Ampère matrix equation

with the charge-conserving current term j n+ 1 2 = 1 2 (j n,-+ j n,+ ), the field correction step (80) rewrites as

Appendix D. Influence of applied exterior fields

In the case where an exterior field (E ext , B ext ) is applied, we point out that the linearly-implicit formulation (56) of the ECSIM scheme needs to be modified in a non-trivial way. Indeed, the velocity kick in [START_REF] Villasenor | Rigorous charge conservation for local electromagnetic field solvers[END_REF] becomes