2.1 Non-canonical Hamiltonian structure

The Vlasov–Maxwell system possesses a noncanonical Hamiltonian structure. The system of equations (2.1)–(2.3) can be obtained from the following Poisson bracket, a bilinear, anti-symmetric bracket that satisfies Leibniz’ rule and the Jacobi identity:

where $[f,g]=\unicode[STIX]{x1D735}_{\boldsymbol{x}}f\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{v}}g-\unicode[STIX]{x1D735}_{\boldsymbol{x}}g\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{v}}f$ . This bracket was introduced in Morrison (Reference Morrison1980), with a term corrected in Marsden & Weinstein (Reference Marsden and Weinstein1982) (see also Weinstein & Morrison Reference Weinstein and Morrison1981; Morrison Reference Morrison1982), and its limitation to divergence-free magnetic fields first pointed out in Morrison (Reference Morrison1982). See also Chandre et al. (Reference Chandre, Guillebon, Back, Tassi and Morrison2013) and Morrison (Reference Morrison2013), where the latter contains the details of the direct proof of the Jacobi identity

The time evolution of any functional ${\mathcal{F}}[\,f_{s},\boldsymbol{E},\boldsymbol{B}]$ is given by

with the Hamiltonian ${\mathcal{H}}$ given as the sum of the kinetic energy of the particles and the electric and magnetic field energies,

In order to obtain the Vlasov equations, we consider the functional

for which the equations of motion (2.10) are computed as

For the electric field, we consider

so that from (2.10) we obtain Ampère’s equation,

where the current density $\boldsymbol{J}$ is given by

And for the magnetic field, we consider

and obtain the Faraday equation,

Our aim is to preserve this noncanonical Hamiltonian structure and its features at the discrete level. This can be done by taking only a finite number of initial positions for the particles instead of a continuum and by taking the electromagnetic fields in finite-dimensional subspaces of the original function spaces. A good candidate for such a discretization is the finite element particle-in-cell framework. In order to satisfy the continuity equation as well as the identities from vector calculus and thereby preserve Gauss’ law and the divergence of the magnetic field, the finite element spaces for the different fields cannot be chosen independently. The right framework is given by FEEC.

Before describing this framework in more detail, we shortly want to discuss some conservation laws of the Vlasov–Maxwell system. In Hamiltonian systems, there are two kinds of conserved quantities, Casimir invariants and momentum maps.

2.2 Invariants

A family of conserved quantities are Casimir invariants (Casimirs), which originate from the degeneracy of the Poisson bracket. Casimirs are functionals ${\mathcal{C}}(f_{s},\boldsymbol{E},\boldsymbol{B})$ which Poisson commute with every other functional ${\mathcal{G}}(f_{s},\boldsymbol{E},\boldsymbol{B})$ , i.e. $\{{\mathcal{C}},{\mathcal{G}}\}=0$ . For the Vlasov–Maxwell bracket, there are several such Casimirs (Morrison Reference Morrison1987; Morrison & Pfirsch Reference Morrison and Pfirsch1989; Chandre et al. Reference Chandre, Guillebon, Back, Tassi and Morrison2013). First, the integral of any real function $h_{s}$ of each distribution function $f_{s}$ is preserved, i.e.

This family of Casimirs is a manifestation of Liouville’s theorem and corresponds to conservation of phase space volume. Further we have two Casimirs related to Gauss’ law (2.4) and the divergence-free property of the magnetic field (2.5),

where $h_{E}$ and $h_{B}$ are arbitrary real functions of $\boldsymbol{x}$ . The latter functional, ${\mathcal{C}}_{B}$ , is not a true Casimir but should rather be referred to as pseudo-Casimir. It acts like a Casimir in that it Poisson commutes with any other functional, but the Jacobi identity is only satisfied when $\text{div}\,B=0$ (see Morrison Reference Morrison1982, Reference Morrison2013).

A second family of conserved quantities are momentum maps $\unicode[STIX]{x1D6F7}$ , which arise from symmetries that preserve the particular Hamiltonian ${\mathcal{H}}$ , and therefore also the equations of motion. This means that the Hamiltonian is constant along the flow of $\unicode[STIX]{x1D6F7}$ , i.e.

From Noether’s theorem it follows that the generators $\unicode[STIX]{x1D6F7}$ of the symmetry are preserved by the time evolution, i.e.

If the symmetry condition (2.22) holds, this is obvious by the anti-symmetry of the Poisson bracket as

Therefore $\unicode[STIX]{x1D6F7}$ is a constant of motion if and only if $\{\unicode[STIX]{x1D6F7},{\mathcal{H}}\}=0$ .

The complete set of constants of motion, the algebra of invariants, will be discussed elsewhere. However, as an example of a momentum map we shall consider here the total momentum

By direct computations, assuming periodic boundary conditions, it can be shown that

defining $Q(\boldsymbol{x}):=\unicode[STIX]{x1D70C}-\text{div}\,\boldsymbol{E}$ , which is a local version of the Casimir ${\mathcal{C}}_{E}$ . Therefore, if at $t=0$ the Casimir $Q\equiv 0$ , then momentum is conserved. If at $t=0$ the Casimir $Q\not \equiv 0$ , then momentum is not conserved and it changes in accordance with (2.26). For a multi-species plasma $Q\equiv 0$ is equivalent to the physical requirement that Poisson’s equation be satisfied. If for some reason it is not exactly satisfied, then we have violation of momentum conservation.

For a single species plasma, say electrons, with a neutralizing positive background charge $\unicode[STIX]{x1D70C}_{B}(\boldsymbol{x})$ , say ions, Poisson’s equation is

The Poisson bracket for this case has the local Casimir

and it does not recognize the background charge. Because the background is stationary, the total momentum is

and it satisfies

We will verify this relation in the numerical experiments of § 7.5.

3.1 Maxwell’s equations

When Maxwell’s equations are used in some material medium, they are best understood by introducing two additional fields. The electromagnetic properties are then defined by the electric and magnetic fields, usually denoted by $\boldsymbol{E}$ and $\boldsymbol{B}$ , the displacement field $\boldsymbol{D}$ and the magnetic intensity $\boldsymbol{H}$ . For simple materials, the electric field is related to the displacement field and the magnetic field to the magnetic intensity by

where $\unicode[STIX]{x1D73A}$ and $\unicode[STIX]{x1D741}$ are the permittivity and permeability tensors reflecting the material properties. In vacuum they become the scalars $\unicode[STIX]{x1D700}_{0}$ and $\unicode[STIX]{x1D707}_{0}$ , which are unity in our scaled variables, while for more complicated media such as plasmas they can be nonlinear operators (Morrison Reference Morrison2013). The Maxwell equations with the four fields read

The mathematical interpretation of these fields become clearer when interpreting them as differential forms: $\boldsymbol{E}$ and $\boldsymbol{H}$ are 1-forms, $\boldsymbol{D}$ and $\boldsymbol{B}$ are 2-forms. The charge density $\unicode[STIX]{x1D70C}$ is a 3-form and the current density $\boldsymbol{J}$ a 2-form. Moreover, the electrostatic potential $\unicode[STIX]{x1D719}$ is a 0-form and the vector potential $\boldsymbol{A}$ a 1-form. The $\text{grad}$ , $\text{curl}$ , $\text{div}$ operators represent the exterior derivative applied respectively to 0-forms, 1-forms and 2-forms. To be more precise, there are two kinds of differential forms, depending on the orientation. Straight differential forms have an intrinsic (or inner) orientation, whereas twisted differential forms have an outer orientation, defined by the ambient space. Faraday’s equation and $\text{div}\,\boldsymbol{B}=0$ are naturally inner oriented, whereas Ampère’s equation and Gauss’ law are outer oriented. This knowledge can be used to define a natural discretization for Maxwell’s equations. For finite difference approximations a dual mesh is needed for the discretization of twisted forms. This can already be found in Yee’s scheme (Yee Reference Yee1966). In the finite element context, only one mesh is used, but dual operators are used for the twisted forms. As an implication, the charge density $\unicode[STIX]{x1D70C}$ will be treated as a 0-form and the current density $J$ as a 1-form, instead of a (twisted) 3-form and a (twisted) 2-form, respectively. Another consequence is that Ampère’s equation and Gauss’ law are being treated weakly while Faraday’s equation and $\text{div}\,\boldsymbol{B}=0$ are treated strongly. A detailed description of this formalism can be found, e.g. in Bossavit’s lecture notes (Bossavit Reference Bossavit2006).

3.2 Finite element spaces of differential forms

The full mathematical theory for the finite element discretization of differential forms is due to Arnold et al. (Reference Arnold, Falk and Winther2006, Reference Arnold, Falk and Winther2010) and is called finite element exterior calculus (see also Monk Reference Monk2003, Christiansen et al. Reference Christiansen, Munthe-Kaas and Owren2011). Most finite element spaces appearing in this theory were known before, but their connection in the context of differential forms was not made clear. The first building block of FEEC is the following commuting diagram:

where $\unicode[STIX]{x1D6FA}\subset \mathbb{R}^{3}$ , $\unicode[STIX]{x1D6EC}^{k}(\unicode[STIX]{x1D6FA})$ is the space of $k$ -forms on $\unicode[STIX]{x1D6FA}$ that we endow with the inner product $\langle \unicode[STIX]{x1D6FC},\unicode[STIX]{x1D6FD}\rangle =\int \unicode[STIX]{x1D6FC}\wedge \star \unicode[STIX]{x1D6FD}$ , $\star$ is the Hodge operator and $\unicode[STIX]{x1D625}$ is the exterior derivative that generalizes the gradient, curl and divergence. Then we define

and the Sobolev spaces of differential forms

Obviously in a three-dimensional manifold the exterior derivative of a 3-form vanishes so that $H\unicode[STIX]{x1D6EC}^{3}(\unicode[STIX]{x1D6FA})=L^{2}(\unicode[STIX]{x1D6FA})$ . This diagram can also be expressed using the standard vector calculus formalism:

The first row of (3.9) represents the sequence of function spaces involved in Maxwell’s equations. Such a sequence is called a complex if at each node, the image of the previous operator is in the kernel of the next operator, i.e. $\text{Im}(\text{grad})\subseteq \text{Ker}(\text{curl})$ and $\text{Im}(\text{curl})\subseteq \text{Ker}(\text{div})$ . The power of the conforming finite element framework is that this complex can be reproduced at the discrete level by choosing the appropriate finite-dimensional subspaces $V_{0}$ , $V_{1}$ , $V_{2}$ , $V_{3}$ . The order of the approximation is dictated by the choice made for $V_{0}$ and the requirement of having a complex at the discrete level. The projection operators $\unicode[STIX]{x1D6F1}_{i}$ are the finite element interpolants, which have the property that the diagram is commuting. This means for example, that the $\text{grad}$ of the projection on $V_{0}$ is identical to the projection of the $\text{grad}$ on $V_{1}$ . As proven by Arnold, Falk and Winther, their choice of finite elements naturally leads to stable discretizations.

There are many known sequences of finite element spaces that fit this diagram. The sequences proposed by Arnold, Falk and Winther are based on well-known finite element spaces. On tetrahedra these are $H^{1}$ conforming $\mathbb{P}_{k}$ Lagrange finite elements for $V_{0}$ , the $H(\text{curl})$ conforming Nédélec elements for $V_{1}$ , the $H(\text{div})$ conforming Raviart–Thomas elements for $V_{2}$ and discontinuous Galerkin elements for $V_{3}$ . A similar sequence can be defined on hexahedra based on the $H^{1}$ conforming $\mathbb{Q}_{k}$ Lagrange finite elements for $V_{0}$ .

Other sequences that satisfy the complex property are available. Let us in particular cite the mimetic spectral elements (Kreeft, Palha & Gerritsma Reference Kreeft, Palha and Gerritsma2011; Gerritsma Reference Gerritsma2012; Palha et al. Reference Palha, Rebelo, Hiemstra, Kreeft and Gerritsma2014) and the spline finite elements (Buffa, Sangalli & Vázquez Reference Buffa, Sangalli and Vázquez2010; Buffa et al. Reference Buffa, Rivas, Sangalli and Vázquez2011; Ratnani and Sonnendrücker Reference Ratnani and Sonnendrücker2012) that we shall use in this work, as splines are generally favoured in PIC codes due to their smoothness properties that enable noise reduction.

3.3 Finite element discretization of Maxwell’s equations

This framework is enough to express discrete relations between all the straight (or primal forms), i.e. $\boldsymbol{E}$ , $\boldsymbol{B}$ , $\boldsymbol{A}$ and $\unicode[STIX]{x1D719}$ . The commuting diagram yields a direct expression of the discrete Faraday equation. Indeed projecting all the components of the equation onto $V_{2}$ yields

which is equivalent, due to the commuting diagram property, to

Denoting with an $h$ index the discrete fields, $\boldsymbol{B}_{h}=\unicode[STIX]{x1D6F1}_{2}\boldsymbol{B}$ , $\boldsymbol{E}_{h}=\unicode[STIX]{x1D6F1}_{1}\boldsymbol{E}$ , this yields the discrete Faraday equation,

In the same way, the discrete electric and magnetic fields are defined exactly as in the continuous case from the discrete potentials, thanks to the compatible finite element spaces,

so that automatically we get

On the other hand, Ampère’s equation and Gauss’ law relate expressions involving twisted differential forms. In the finite element framework, these should be expressed on the dual complex to (3.9). But due to the property that the dual of an operator in $L^{2}(\unicode[STIX]{x1D6FA})$ can be identified with its $L^{2}$ adjoint via an inner product, the discrete dual spaces are such that $V_{0}^{\ast }=V_{3}$ , $V_{1}^{\ast }=V_{2}$ , $V_{2}^{\ast }=V_{1}$ and $V_{3}^{\ast }=V_{0}$ , so that the dual operators and spaces are not explicitly needed. They are most naturally used seamlessly by keeping the weak formulation of the corresponding equations. The weak form of Ampère’s equation is found by taking the dot product of (2.2) with a test function $\bar{\boldsymbol{E}}\in H(\text{curl},\unicode[STIX]{x1D6FA})$ and applying a Green identity. Assuming periodic boundary conditions, the weak solution of Ampère’s equation $(\boldsymbol{E},\boldsymbol{B})\in H(\text{curl},\unicode[STIX]{x1D6FA})\times H(\text{div},\unicode[STIX]{x1D6FA})$ is characterized by

The discrete version is obtained by replacing the continuous spaces by their finite-dimensional subspaces. The approximate solution $(\boldsymbol{E}_{h},\boldsymbol{B}_{h})\in V_{1}\times V_{2}$ is characterized by

In the same way the weak solution of Gauss’ law with $\boldsymbol{E}\in H(\text{curl},\unicode[STIX]{x1D6FA})$ is characterized by

its discrete version for $\boldsymbol{E}_{h}\in V_{1}$ being characterized by

The last step for the finite element discretization is to define a basis for each of the finite-dimensional spaces $V_{0},V_{1},V_{2},V_{3}$ , with $\dim V_{k}=N_{k}$ and to find equations relating the coefficients on these bases. Let us denote by $\{\unicode[STIX]{x1D6EC}_{i}^{0}\}_{i=1\ldots N_{0}}$ and $\{\unicode[STIX]{x1D6EC}_{i}^{3}\}_{i=1\ldots N_{3}}$ a basis of $V_{0}$ and $V_{3}$ , respectively, which are spaces of scalar functions, and $\{\unicode[STIX]{x1D726}_{i,\unicode[STIX]{x1D707}}^{1}\}_{i=1\ldots N_{1},\unicode[STIX]{x1D707}=1\ldots 3}$ a basis of $V_{1}\subset H(\text{curl},\unicode[STIX]{x1D6FA})$ and $\{\unicode[STIX]{x1D726}_{i,\unicode[STIX]{x1D707}}^{2}\}_{i=1\ldots N_{2},\unicode[STIX]{x1D707}=1\ldots 3}$ a basis of $V_{2}\subset H(\text{div},\unicode[STIX]{x1D6FA})$ , which are vector valued functions,

Let us note that the restriction to a basis of this form is not strictly necessary and the generalization to more general bases is straightforward. However, for didactical reasons we stick to this form of the basis as it simplifies some of the computations and thus helps to clarify the following derivations. In order to keep a concise notation, and by slight abuse of the same, we introduce vectors of basis functions $\unicode[STIX]{x1D726}^{k}=(\unicode[STIX]{x1D726}_{1,1}^{k},\unicode[STIX]{x1D726}_{1,2}^{k},\ldots ,\unicode[STIX]{x1D726}_{N_{k},3}^{k})^{\text{T}}$ for $k=1,2$ , which are indexed by $I=3(i-1)+\unicode[STIX]{x1D707}=1\ldots 3N_{k}$ with $i=1\ldots N_{k}$ and $\unicode[STIX]{x1D707}=1\ldots 3$ , and $\unicode[STIX]{x1D6EC}^{k}=(\unicode[STIX]{x1D6EC}_{1}^{k},\unicode[STIX]{x1D6EC}_{2}^{k},\ldots ,\unicode[STIX]{x1D6EC}_{N_{k}}^{k})^{\text{T}}$ for $k=0,3$ , which are indexed by $i=1\ldots N_{k}$ .

We shall also need for each basis the dual basis, which in finite element terminology corresponds to the degrees of freedom. For each basis $\unicode[STIX]{x1D6EC}_{i}^{k}$ for $k=0,3$ and $\unicode[STIX]{x1D726}_{I}^{k}$ for $k=1,2$ , the dual basis is denoted by $\unicode[STIX]{x1D6F4}_{i}^{k}$ and $\unicode[STIX]{x1D72E}_{I}^{k}$ , respectively, and defined by

for the scalar valued bases $\unicode[STIX]{x1D6EC}_{i}^{k}$ , and

for the vector valued bases $\unicode[STIX]{x1D726}_{I}^{k}$ , where $\left<\boldsymbol{\cdot },\boldsymbol{\cdot }\right>$ denotes the $L^{2}$ inner product in the appropriate space and $\unicode[STIX]{x1D6FF}_{IJ}$ is the Kronecker symbol, whose value is unity for $I=J$ and zero otherwise. We introduce the linear functionals $L^{2}\unicode[STIX]{x1D6EC}^{k}(\unicode[STIX]{x1D6FA})\rightarrow \mathbb{R}$ , which are denoted by $\unicode[STIX]{x1D70E}_{i}^{k}$ for $k=0,3$ and by $\unicode[STIX]{x1D748}_{I}^{k}$ for $k=1,2$ , respectively. On the finite element space they are represented by the dual basis functions $\unicode[STIX]{x1D6F4}_{i}^{k}$ and $\unicode[STIX]{x1D72E}_{I}^{k}$ and defined by

and

so that $\unicode[STIX]{x1D70E}_{i}^{k}(\unicode[STIX]{x1D6EC}_{j}^{k})=\unicode[STIX]{x1D6FF}_{ij}$ and $\unicode[STIX]{x1D748}_{I}^{k}(\unicode[STIX]{x1D726}_{J}^{k})=\unicode[STIX]{x1D6FF}_{IJ}$ for the appropriate $k$ . Elements of the finite-dimensional spaces can be expanded on their respective bases, e.g. elements of $V_{1}$ and $V_{2}$ , respectively, as

denoting by $\boldsymbol{e}=(e_{1,1},\,e_{1,2},\ldots ,e_{N_{1},3})^{\text{T}}\in \mathbb{R}^{3N_{1}}$ and $\boldsymbol{b}=(b_{1,1},\,b_{1,2},\ldots ,b_{N_{2},3})^{\text{T}}\in \mathbb{R}^{3N_{2}}$ the corresponding degrees of freedom with $e_{i,\unicode[STIX]{x1D707}}=\unicode[STIX]{x1D748}_{i,\unicode[STIX]{x1D707}}^{1}(\boldsymbol{E}_{h})$ and $b_{i,\unicode[STIX]{x1D707}}=\unicode[STIX]{x1D748}_{i,\unicode[STIX]{x1D707}}^{2}(\boldsymbol{B}_{h})$ , respectively.

Denoting the elements of $\boldsymbol{e}$ by $\boldsymbol{e}_{I}$ and the elements of $\boldsymbol{b}$ by $\boldsymbol{b}_{I}$ , we have that $\boldsymbol{e}_{I}=\unicode[STIX]{x1D748}_{I}^{1}(\boldsymbol{E}_{h})$ and $\boldsymbol{b}_{I}=\unicode[STIX]{x1D748}_{I}^{2}(\boldsymbol{B}_{h})$ , respectively, and can re-express (3.25) as

Henceforth we will use both notations in parallel, choosing whichever is more practical at any given time.

Due to the complex property we have that $\text{curl}\,\boldsymbol{E}_{h}\in V_{2}$ for all $\boldsymbol{E}_{h}\in V_{1}$ , so that $\text{curl}\,\boldsymbol{E}_{h}$ can be expressed in the basis of $V_{2}$ by

Let us also denote by $\boldsymbol{c}=(c_{1,1},\,c_{1,2},\ldots ,c_{N_{2},3})^{\text{T}}$ , so that $\text{curl}\,\boldsymbol{E}_{h}$ can also be written as

On the other hand

Denoting by $\mathbb{C}$ the discrete curl matrix,

the degrees of freedom of $\text{curl}\,\boldsymbol{E}_{h}$ in $V_{2}$ are related to the degrees of freedom of $\boldsymbol{E}_{h}$ in $V_{1}$ by $\boldsymbol{c}=\mathbb{C}\boldsymbol{e}$ . In the same way we can define the discrete gradient matrix $\mathbb{G}$ and the discrete divergence matrix $\mathbb{D}$ , given by

respectively. Denoting by $\unicode[STIX]{x1D753}=(\unicode[STIX]{x1D711}_{1},\ldots ,\unicode[STIX]{x1D711}_{N_{0}})^{\text{T}}$ and $\boldsymbol{a}=(\mathit{a}_{1,1},\,\mathit{a}_{1,2},\,\ldots ,\,\mathit{a}_{N_{1},3})^{\text{T}}$ the degrees of freedom of the potentials $\unicode[STIX]{x1D719}_{h}$ and $\boldsymbol{A}_{h}$ , with $\unicode[STIX]{x1D711}_{i}=\unicode[STIX]{x1D70E}_{i}^{0}(\unicode[STIX]{x1D719}_{h})$ for $1\leqslant i\leqslant N_{0}$ and $\boldsymbol{a}_{I}=\unicode[STIX]{x1D748}_{I}^{1}(\boldsymbol{A}_{h})$ for $1\leqslant I\leqslant 3N_{1}$ , the relation (3.13) between the discrete fields (3.25) and the potentials can be written using only the degrees of freedom as

Finally, we need to define the so-called mass matrices in each of the discrete spaces $V_{i}$ , which define the discrete Hodge operator linking the primal complex with the dual complex. We denote by $(\mathbb{M}_{0})_{ij}=\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D6EC}_{i}^{0}(\boldsymbol{x})\,\unicode[STIX]{x1D6EC}_{j}^{0}(\boldsymbol{x})\,\text{d}\boldsymbol{x}$ with $1\leqslant i,j\leqslant N_{0}$ and $(\mathbb{M}_{1})_{IJ}=\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D726}_{I}^{1}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D726}_{J}^{1}(\boldsymbol{x})\,\text{d}\boldsymbol{x}$ with $1\leqslant I,J\leqslant 3N_{1}$ the mass matrices in $V_{0}$ and $V_{1}$ , respectively, and similarly $\mathbb{M}_{2}$ and $\mathbb{M}_{3}$ the mass matrices in $V_{2}$ and $V_{3}$ . Using these definitions as well as $\unicode[STIX]{x1D754}=(\unicode[STIX]{x1D71A}_{1},\ldots ,\unicode[STIX]{x1D71A}_{N_{0}})^{\text{T}}$ and $\boldsymbol{j}=(j_{1,1},\,j_{1,2},\ldots ,j_{N_{1},3})^{\text{T}}$ with $\unicode[STIX]{x1D71A}_{i}=\unicode[STIX]{x1D70E}_{i}^{0}(\unicode[STIX]{x1D70C}_{h})$ for $1\leqslant i\leqslant N_{0}$ and $\boldsymbol{j}_{I}=\unicode[STIX]{x1D748}_{I}^{1}(\boldsymbol{J}_{h})$ for $1\leqslant I\leqslant 3N_{1}$ (recall that the charge density $\unicode[STIX]{x1D70C}$ is treated as a 0-form and the current density $\boldsymbol{J}$ as a 1-form), we obtain a system of ordinary differential equations for each of the continuous equations, namely

It is worth emphasizing that $\text{div}\,\boldsymbol{B}=0$ is satisfied in strong form, which is important for the Jacobi identity of the discretized Poisson bracket (cf. § 4.4). The complex properties can also be expressed at the matrix level. The primal sequence being

with $\text{Im}\mathbb{G}\subseteq \text{Ker}\mathbb{C}$ , $\text{Im}\mathbb{C}\subseteq \text{Ker}\mathbb{D}$ , and the dual sequence being

with $\text{Im}\mathbb{D}^{\text{T}}\subseteq \text{Ker}\mathbb{C}^{\text{T}}$ , $\text{Im}\mathbb{C}^{\text{T}}\subseteq \text{Ker}\mathbb{G}^{\text{T}}$ .

3.4 Example: B-spline finite elements

In the following, we will use so-called basic splines, or B-splines, as bases for the finite element function spaces. B-splines are piecewise polynomials. The points where two polynomials connect are called knots. The $j$ th basic spline (B-spline) of degree $p$ can be defined recursively by

where

and

with the knot vector $\unicode[STIX]{x1D6EF}=\{x_{i}\}_{1\leqslant i\leqslant N+k}$ being a non-decreasing sequence of points. The knot vector can also contain repeated knots. If a knot $x_{i}$ has multiplicity $m$ , then the B-spline is $C^{p-m}$ continuous at $x_{i}$ . The derivative of a B-spline of degree $p$ can easily be computed as the difference of two B-splines of degree $p-1$ ,

For convenience, we introduce the following shorthand notation for differentials,

In the case of an equidistant grid with grid step size $\unicode[STIX]{x0394}x=x_{j+1}-x_{j}$ , this simplifies to

Using $D_{j}^{p}$ the recursion formula (3.39) becomes

In more than one dimension, we can define tensor-product B-spline basis functions, e.g. for three dimensions as

The bases of the differential form spaces will be tensor products of the basis functions $N_{i}^{p}$ and the differentials $D_{j}^{p}$ . The discrete function spaces are given by

These choices appear quite natural when one considers the action of the gradient on 0-forms, the curl on 1-forms and the divergence on 2-forms. In the following, we will exemplify this using the potentials and fields of Maxwell’s equations. The semi-discrete potentials are written in the respective spline basis (3.47) as

and

Computing the gradient of the semi-discrete scalar potential $\unicode[STIX]{x1D719}_{h}$ , we find

Assuming periodic boundary conditions, the sums can be re-arranged to give

Similarly the curl of the semi-discrete vector potential $\boldsymbol{A}_{h}$ is computed as

Again, assuming periodic boundary conditions, the sums can be re-arranged to give

Given the above, we determine the semi-discrete electromagnetic fields $\boldsymbol{E}_{h}$ and $\boldsymbol{B}_{h}$ . The electric field $\boldsymbol{E}_{h}=-\text{grad}\,\unicode[STIX]{x1D719}_{h}-\unicode[STIX]{x2202}\boldsymbol{A}_{h}/\unicode[STIX]{x2202}t$ is computed as

and the magnetic field $\boldsymbol{B}_{h}=\text{curl}\boldsymbol{A}_{h}$ as

Now it becomes clear why definitions (3.47) are the natural choice for the spline bases and it is straightforward to verify that $\text{div}\,\boldsymbol{B}_{h}=0$ .

4.1 Discretization of the functional field derivatives

Upon inserting (3.25), any functional ${\mathcal{F}}[\boldsymbol{E}_{h}]$ can be considered as a function $F(\boldsymbol{e})$ of the finite element coefficients,

Therefore, we can write the first variation of ${\mathcal{F}}[\boldsymbol{E}]$ ,

as

with

Let $\unicode[STIX]{x1D72E}_{I}^{1}(\boldsymbol{x})$ denote the dual basis of $\unicode[STIX]{x1D726}_{I}^{1}(\boldsymbol{x})$ with respect to the $L^{2}$ inner product (3.24) and let us express the functional derivative on this dual basis, i.e.

Using (4.5) for $\bar{\boldsymbol{e}}=(0,\ldots ,0,1,0,,0)^{\text{T}}$ with $1$ at the $I$ th position and $0$ everywhere else, so that $\bar{\boldsymbol{E}}_{h}=\unicode[STIX]{x1D726}_{I}^{1}$ , we find that

and can therefore write

On the other hand, expanding the dual basis in the original basis,

and taking the $L^{2}$ inner product with $\unicode[STIX]{x1D726}_{J}^{1}(\boldsymbol{x})$ , we find that the matrix $\mathbb{A}=(\boldsymbol{a}_{IJ})$ verifies $\mathbb{A}\mathbb{M}_{1}=\mathbb{I}_{3N_{1}}$ , where $\mathbb{I}_{3N_{1}}$ denotes the $3N_{1}\times 3N_{1}$ identity matrix, so that $\mathbb{A}$ is the inverse of the mass matrix $\mathbb{M}_{1}$ and

In full analogy we find

Next, using (3.30), we find

Finally, we can re-express the following term in the Poisson bracket

The other terms in the bracket involving functional derivatives with respect to the fields are handled similarly. In the next step we need to discretize the distribution function and the corresponding functional derivatives.

4.2 Discretization of the functional particle derivatives

We proceed by assuming a particle-like distribution function for $N_{p}$ particles labelled by  $a$ ,

with mass $m_{a}$ , charge $q_{a}$ , weights $w_{a}$ , particle positions $\boldsymbol{x}_{a}$ and particle velocities $\boldsymbol{v}_{a}$ . Here, $N_{p}$ denotes the total number of particles of all species, with each particle carrying its own mass and signed charge. Functionals of the distribution function, ${\mathcal{F}}[f]$ , can be considered as functions of the particle phase space trajectories, $F(\boldsymbol{X},\boldsymbol{V})$ , upon inserting (4.15),

Variation of the left-hand side and integration by parts gives,

Upon equating this expression with the variation of the right-hand side of (4.16),

we obtain

Considering the kinetic part of the Poisson bracket (2.8), the discretization proceeds in two steps. First, replace $f_{s}$ with $f_{h}$ to get

Then insert (4.19) in order to obtain the discrete kinetic bracket.

4.3 Discrete Poisson bracket

Replacing all functional derivatives in (2.8) as outlined in the previous two sections, we obtain the semi-discrete Poisson bracket

with the curl matrix $\mathbb{C}$ as given in (3.30). In order to express the semi-discrete Poisson bracket (4.21) in matrix form, we denote by $\unicode[STIX]{x1D726}^{1}(\boldsymbol{X})$ the $3N_{p}\times 3N_{1}$ matrix with generic term $\unicode[STIX]{x1D726}_{I}^{1}(\boldsymbol{x}_{a})$ , where $1\leqslant a\leqslant N_{p}$ and $1\leqslant I\leqslant 3N_{1}$ , and by $\mathbb{B}(\boldsymbol{X},\boldsymbol{b})$ the $3N_{p}\times 3N_{p}$ block diagonal matrix with generic block

Further, let us introduce a mass matrix $\unicode[STIX]{x1D648}_{p}$ and a charge matrix $\unicode[STIX]{x1D648}_{q}$ for the particles. Both are diagonal $N_{p}\times N_{p}$ matrices with elements $(\unicode[STIX]{x1D648}_{p})_{aa}=m_{a}w_{a}$ and $(\unicode[STIX]{x1D648}_{q})_{aa}=q_{a}w_{a}$ , respectively. Additionally, we will need the $3N_{p}\times 3N_{p}$ matrices

where $\mathbb{I}_{3}$ denotes the $3\times 3$ identity matrix. This allows us to rewrite

Here, the derivatives are represented by the $3N_{p}$ vector

and correspondingly for $\unicode[STIX]{x2202}G/\unicode[STIX]{x2202}\boldsymbol{V}$ , $\unicode[STIX]{x2202}F/\unicode[STIX]{x2202}\boldsymbol{e}$ , $\unicode[STIX]{x2202}F/\unicode[STIX]{x2202}\boldsymbol{b}$ , etc. Thus, the discrete Poisson bracket (4.21) becomes

The action of this bracket on two functions $F$ and $G$ can also be expressed as

denoting by $\mathit{D}$ the derivative with respect to the dynamical variables

and by $\mathbb{J}$ the Poisson matrix, given by

We immediately see that $\mathbb{J}(\boldsymbol{u})$ is anti-symmetric, but it remains to be shown that it satisfies the Jacobi identity.

4.4 Jacobi identity

The discrete Poisson bracket (4.26) satisfies the Jacobi identity if and only if the following condition holds (see e.g. (Morrison Reference Morrison1998, § IV) or (Hairer, Lubich & Wanner Reference Hairer, Lubich and Wanner2006, § VII.2, Lemma 2.3)):

Here, all indices $i,j,k,l$ run from $1$ to $6N_{p}+3N_{1}+3N_{2}$ . To simplify the verification of (4.30), we start by identifying blocks of the Poisson matrix $\mathbb{J}$ whose elements contribute to the above condition. Therefore, we write

The Poisson matrix $\mathbb{J}$ only depends on $\boldsymbol{X}$ and $\boldsymbol{b}$ , so in (4.30) we only need to sum $l$ over the corresponding indices, $1\leqslant l\leqslant 3N_{p}$ and $6N_{p}+3N_{1}<l\leqslant 6N_{p}+3N_{1}+3N_{2}$ , respectively. Considering the terms $\mathbb{J}_{li}(\boldsymbol{u})$ , $\mathbb{J}_{lj}(\boldsymbol{u})$ and $\mathbb{J}_{lk}(\boldsymbol{u})$ , we see that in the aforementioned index ranges for $l$ , only $\unicode[STIX]{x1D645}_{12}=\mathbb{M}_{p}^{-1}$ and $\unicode[STIX]{x1D645}_{43}=-\mathbb{M}_{2}^{-1}\mathbb{C}\mathbb{M}_{1}^{-1}$ are non-vanishing, so that we have to account only for those two blocks, i.e. for $1\leqslant l\leqslant 3N_{p}$ we only need to consider $3N_{p}<i,j,k\leqslant 6N_{p}$ and for $6N_{p}+3N_{1}<l\leqslant 6N_{p}+3N_{1}+3N_{2}$ we only need to consider $6N_{p}<i,j,k\leqslant 6N_{p}+3N_{1}$ . Note that $\unicode[STIX]{x1D645}_{12}$ is a diagonal matrix, therefore $(\unicode[STIX]{x1D645}_{12})_{ab}=(\unicode[STIX]{x1D645}_{12})_{aa}\unicode[STIX]{x1D6FF}_{ab}$ with $1\leqslant a,b\leqslant N_{p}$ . Further, only $\unicode[STIX]{x1D645}_{22}$ , $\unicode[STIX]{x1D645}_{23}$ and $\unicode[STIX]{x1D645}_{32}$ depend on $\boldsymbol{b}$ and/or $\boldsymbol{X}$ , so only those blocks have to be considered when computing derivatives with respect to $\boldsymbol{u}$ . In summary, we obtain two conditions. The contributions involving $\unicode[STIX]{x1D645}_{22}$ and $\unicode[STIX]{x1D645}_{12}$ are

for $1\leqslant b,c,d\leqslant 3N_{p}$ , which corresponds to (4.30) for $3N_{p}<i,j,k\leqslant 6N_{p}$ . Inserting the actual values for $\unicode[STIX]{x1D645}_{12}$ and $\unicode[STIX]{x1D645}_{22}$ and using that $\mathbb{M}_{p}$ is diagonal, equation (4.32) becomes

All outer indices of this expression belong to the inverse matrix $\mathbb{M}_{p}^{-1}$ . As this matrix is constant, symmetric and positive definite, we can contract the above expression with $\mathbb{M}_{p}$ on all indices, to obtain

If in the first term of (4.34) one picks a particular index $k$ , then this selects the $\unicode[STIX]{x1D70E}$ component of the position $\boldsymbol{x}_{a}$ of some particle. At the same time, in the second and third terms, this selects a block of (4.22), which is evaluated at the same particle position $\boldsymbol{x}_{a}$ . This means that the only non-vanishing contributions in (4.34) will be those indexed by $b$ and $c$ with $\boldsymbol{X}_{b}$ and $\boldsymbol{X}_{c}$ corresponding to components $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D708}$ of the same particle position $\boldsymbol{x}_{a}$ . Therefore, the condition (4.22) reduces further to

where $\widehat{B}_{\unicode[STIX]{x1D707}\unicode[STIX]{x1D708}}$ denotes the components of the matrix in (4.22). When all three indices are equal, this corresponds to diagonal terms of the matrix $\widehat{\unicode[STIX]{x1D63D}}_{h}(\boldsymbol{x}_{a},t)$ which vanish. When two of the three are equal, it cancels because of the skew symmetry of the same matrix and for all three indices distinct, this condition corresponds to $\text{div}\,\boldsymbol{B}_{h}=0$ . Choosing initial conditions such that $\text{div}\,\boldsymbol{B}_{h}(\boldsymbol{x},0)=0$ and using a discrete deRham complex guarantees $\text{div}\,\boldsymbol{B}_{h}(\boldsymbol{x},t)=0$ for all times $t$ . Note that this was to be expected, because it is the discrete version of the $\text{div}\,\boldsymbol{B}=0$ condition for the continuous Poisson bracket (2.8) (Morrison Reference Morrison1982). Further note that in the discrete case, just as in the continuous case, the Jacobi identity requires the magnetic field to be divergence free, but it is not required to be the curl of some vector potential. In the language of differential forms this is to say that the magnetic field as a 2-form needs to be closed but does not need to be exact.

From the contributions involving $\unicode[STIX]{x1D645}_{22}$ , $\unicode[STIX]{x1D645}_{23}$ , $\unicode[STIX]{x1D645}_{32}$ and $\unicode[STIX]{x1D645}_{43}$ we have

for $1\leqslant b,c\leqslant 3N_{p}$ and $1\leqslant I\leqslant 3N_{1}$ , which corresponds to (4.30) for $3N_{p}<i,j\leqslant 6N_{p}<k\leqslant 6N_{p}+3N_{1}$ . Writing out (4.36) and using that $\mathbb{M}_{p}$ is diagonal gives

Again, we can contract this with the matrix $\mathbb{M}_{p}$ on the indices $b$ and $c$ , in order to remove $\mathbb{M}_{p}^{-1}$ , and with $\mathbb{M}_{1}$ on the index $I$ , in order to remove $\mathbb{M}_{1}^{-1}$ . Similarly, $\mathbb{M}_{q}$ can be removed by contracting with $\mathbb{M}_{q}^{-1}$ , noting that this matrix is constant and therefore commutes with the curl. This results in the simplified condition

The $\boldsymbol{b}_{J}$ derivative of $\mathbb{B}$ results in the $3N_{p}\times 3N_{p}$ block diagonal matrix $\unicode[STIX]{x1D726}_{J}^{2}(\boldsymbol{X})$ with generic block

so that (4.38) becomes

This condition can be compactly written as

That the charge matrices $\mathbb{M}_{q}$ can be eliminated can be seen as follows. Similarly as before, picking a particular index $c$ in the first term selects the $\unicode[STIX]{x1D708}$ component of the position $\boldsymbol{x}_{a}$ of some particle. At the same time, in the second term, this selects the $\unicode[STIX]{x1D708}$ component of $\unicode[STIX]{x1D726}^{1}$ , evaluated at the same particle position $\boldsymbol{x}_{a}$ . The only non-vanishing derivative of this term is therefore with respect to components of the same particle position, so that $\boldsymbol{X}_{b}$ denotes the $\unicode[STIX]{x1D707}$ component of $\boldsymbol{x}_{a}$ . Hence, condition (4.40) simplifies to

for $1\leqslant a\leqslant N_{p}$ , $1\leqslant \unicode[STIX]{x1D707},\unicode[STIX]{x1D708}\leqslant 3$ and $1\leqslant I\leqslant 3N_{1}$ . The charge $q_{a}$ is the same on both sides and can therefore be removed. This conditions states how the curl of the 1-form basis, evaluated at some particle’s position, is expressed in the 2-form basis, evaluated at the same particle’s position, using the curl matrix $\mathbb{C}$ . For spaces which build a deRham complex, this is always satisfied. This concludes the verification of condition (4.30) for the Jacobi identity to hold for the discrete bracket (4.26).

4.5 Discrete Hamiltonian and equations of motion

The Hamiltonian is discretized by inserting (4.15) and (3.25) into (2.11),

which in matrix notation becomes

The semi-discrete equations of motion are given by

which are equivalent to

With $\mathit{D}H(\boldsymbol{u})=(0,\,\mathbb{M}_{p}\boldsymbol{V},\,\mathbb{M}_{1}\boldsymbol{e},\,\mathbb{M}_{2}\boldsymbol{b})^{\text{T}}$ , we obtain

where the first two equations describe the particle dynamics and the last two equations describe the evolution of the electromagnetic fields. Note that $\mathbb{M}_{p}$ and $\mathbb{M}_{q}$ are diagonal matrices and that $\mathbb{M}_{p}^{-1}\mathbb{M}_{q}$ is nothing but the factor $q_{a}/m_{a}$ for the particle labelled by $a$ . The purpose of introducing these matrices is solely to obtain a compact notation and to display the Poisson structure of the semi-discrete system. However, these matrices are neither explicitly constructed nor is there a need to invert them.

4.6 Discrete Gauss’ law

Multiplying (4.47c ) by $\mathbb{G}^{\text{T}}\mathbb{M}_{1}$ on the left, we get

As $\mathbb{C}\mathbb{G}=0$ from (3.37), the first term on the right-hand side vanishes. Observe that

and using $\text{d}\boldsymbol{x}_{a}/\text{d}t=\boldsymbol{v}_{a}$ we find that

for any $\unicode[STIX]{x1D6F9}_{h}\in V_{0}$ with $\text{grad}\,\unicode[STIX]{x1D6F9}_{h}\in V_{1}$ , so that we obtain

where $\unicode[STIX]{x1D7D9}_{N_{p}}$ denotes the column vector with $N_{p}$ terms all being unity, needed for the sum over the particles when there is no velocity vector. This shows that the discrete Gauss’ law is conserved,

Moreover, note that (4.51) also contains a discrete version of the continuity equation (2.7),

with the discrete charge and current density given by

The conservation of this relation in the fully discrete system depends on the choice of the time stepping scheme.

4.7 Discrete Casimir invariants

Let us now find the Casimir invariants of the semi-discrete Poisson structure. In order to obtain the discrete Casimir invariants, we assume that the discrete spaces in (3.9) not only form a complex, so that $\text{Im}(\text{grad})\subseteq \text{Ker}(\text{curl})$ and $\text{Im}(\text{curl})\subseteq \text{Ker}(\text{div})$ , but that they form an exact sequence, so that $\text{Im}(\text{grad})=\text{Ker}(\text{curl})$ and $\text{Im}(\text{curl})=\text{Ker}(\text{div})$ , i.e. at each node in (3.9), the image of the previous operator is not only a subset of the kernel of the next operator but is exactly equal to the kernel of the next operator. We will then see that this requirement is not necessary for the identified functionals to be valid Casimir invariants. However, it is a useful assumption in their identification.

The Casimir invariants of the semi-discrete system are functions $C(\boldsymbol{X},\boldsymbol{V},\boldsymbol{e},\boldsymbol{b})$ such that $\{C,F\}=0$ for any function $F$ . In terms of the Poisson matrix $\mathbb{J}$ , this can be expressed as $\mathbb{J}(\boldsymbol{u})\,\mathit{D}C(\boldsymbol{u})=0$ . Upon writing this for each of the lines of $\mathbb{J}$ of (4.29), we find for the first line

This implies that $C$ does not depend on $\boldsymbol{V}$ , which we shall use in the sequel. Then the third line simply becomes

Then, because of the exact sequence property, there exists $\bar{\boldsymbol{b}}\in \mathbb{R}^{N_{3}}$ such that $\mathit{D}_{\boldsymbol{b}}C=\mathbb{D}^{\text{T}}\bar{\boldsymbol{b}}$ . Hence all functions of the form

are Casimirs, which means that $\mathbb{D}\boldsymbol{b}$ , the matrix form of $\text{div}\,\boldsymbol{B}_{h}$ , is conserved.

The fourth line, using that $C$ does not depend on $\boldsymbol{V}$ , becomes

Because of the exact sequence property there exists $\bar{\boldsymbol{e}}\in \mathbb{R}^{N_{1}}$ , such that $\mathit{D}_{\boldsymbol{e}}C=\mathbb{M}_{1}\mathbb{G}\bar{\boldsymbol{e}}$ . Finally, the second line couples $\boldsymbol{e}$ and $\boldsymbol{X}$ and reads, upon multiplying by $M_{p}$ ,

using the expression for $\mathit{D}_{\boldsymbol{e}}C$ derived previously and (4.49). It follows that all functions of the form

are Casimirs, so that $\mathbb{G}^{\text{T}}\mathbb{M}_{1}\boldsymbol{e}+\unicode[STIX]{x1D726}^{0}(\boldsymbol{X})^{\text{T}}\mathbb{M}_{q}\unicode[STIX]{x1D7D9}_{N}$ is conserved. This is the matrix form of Gauss’ law (4.52).

Having identified the discrete Casimir invariants (4.57) and (4.60), it is easy to see by plugging them into the discrete Poisson bracket (4.26) that they are valid Casimir invariants, even if the deRham complex is not exact, because all that is needed for $\mathbb{J}(\boldsymbol{u})\,\mathit{D}C(\boldsymbol{u})=0$ is the complex property $\text{Im}\mathbb{D}^{\text{T}}\subseteq \text{Ker}\mathbb{C}^{\text{T}}$ and $\text{grad}\unicode[STIX]{x1D726}^{0}(\boldsymbol{X})=\unicode[STIX]{x1D726}^{1}(\boldsymbol{X})\mathbb{G}$ .

5.1 Solution of the subsystems

The discrete equations of motion for $H_{E}$ are

For initial conditions $(\boldsymbol{X}(0),\boldsymbol{V}(0),\boldsymbol{e}(0),\boldsymbol{b}(0))$ the exact solutions at time $\unicode[STIX]{x0394}t$ are given by the map $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,E}$ defined as The discrete equations of motion for $H_{B}$ are For initial conditions $(\boldsymbol{X}(0),\boldsymbol{V}(0),\boldsymbol{e}(0),\boldsymbol{b}(0))$ the exact solutions at time $\unicode[STIX]{x0394}t$ are given by the map $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,B}$ defined as The discrete equations of motion for $H_{p}$ are For general magnetic field coefficients $\boldsymbol{b}$ , this system cannot be exactly integrated (He et al. Reference He, Qin, Sun, Xiao, Zhang and Liu2015). Note that each component $\dot{\boldsymbol{V}}_{\unicode[STIX]{x1D707}}$ of the equation for $\dot{\boldsymbol{V}}$ does not depend on $\boldsymbol{V}_{\unicode[STIX]{x1D707}}$ , where $\boldsymbol{V}_{\unicode[STIX]{x1D707}}=(\mathit{v}_{1,\unicode[STIX]{x1D707}},\mathit{v}_{2,\unicode[STIX]{x1D707}},\ldots ,\mathit{v}_{N_{p},\unicode[STIX]{x1D707}})^{\text{T}}$ , etc., with $1\leqslant \unicode[STIX]{x1D707}\leqslant 3$ . Therefore we can split this system once more into

with

For concise notation we introduce the $N_{p}\times N_{1}$ matrix $\unicode[STIX]{x1D726}_{\unicode[STIX]{x1D707}}^{1}(\boldsymbol{X})$ with generic term $\unicode[STIX]{x1D726}_{i,\unicode[STIX]{x1D707}}^{1}(\boldsymbol{x}_{a})$ , and the $N_{p}\times N_{p}$ diagonal matrix $\unicode[STIX]{x1D726}_{\unicode[STIX]{x1D707}}^{2}(\boldsymbol{b},\boldsymbol{X})$ with entries $\sum _{i=1}^{N_{2}}b_{i}(t)\unicode[STIX]{x1D726}_{i,\unicode[STIX]{x1D707}}^{2}(\boldsymbol{x}_{a})$ , where $1\leqslant \unicode[STIX]{x1D707}\leqslant 3$ , $1\leqslant a\leqslant N_{p}$ , $1\leqslant i\leqslant N_{1}$ , $1\leqslant j\leqslant N_{2}$ . Then, for $H_{p_{1}}$ we have

for $H_{p_{2}}$ we have and for $H_{p_{3}}$ we have For $H_{p_{1}}$ and initial conditions $(\boldsymbol{X}(0),\boldsymbol{V}(0),\boldsymbol{e}(0),\boldsymbol{b}(0))$ the exact solutions at time $\unicode[STIX]{x0394}t$ are given by the map $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,p_{1}}$ defined as for $H_{p_{2}}$ by the map $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,p_{2}}$ defined as and for $H_{p_{3}}$ by the map $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,p_{3}}$ defined as respectively, where all components not specified are constant. The only challenge in solving these equations is the exact computation of line integrals along the trajectories (Campos Pinto et al. Reference Campos Pinto, Jund, Salmon and Sonnendrücker2014; Squire et al. Reference Squire, Qin and Tang2012; Moon et al. Reference Moon, Teixeira and Omelchenko2015). However, because only one component of the particle positions $\boldsymbol{x}_{a}$ is changing in each step of the splitting and, moreover, the trajectory during one sub step of the splitting is approximated by a straight line, this is not very complicated. Compared to standard particle-in-cell methods for the Vlasov–Maxwell system, the exact integration causes slightly increased computational costs. These are, however, comparable to existing charge-conserving algorithms like that of Villasenor & Buneman (Reference Villasenor and Buneman1992) or the Boris correction method (Boris Reference Boris1970).

5.2 Splitting methods

Given initial conditions $\boldsymbol{u}(0)=(\boldsymbol{X}(0),\boldsymbol{V}(0),\boldsymbol{e}(0),\boldsymbol{b}(0))^{\text{T}}$ , a numerical solution of the discrete Vlasov–Maxwell equations (4.47a )–(4.47d ) at time $\unicode[STIX]{x0394}t$ can be obtained by composition of the exact solutions of all the subsystems. A first-order integrator can be obtained by the Lie–Trotter composition (Trotter Reference Trotter1959),

A second-order integrator can be obtained by the symmetric Strang composition (Strang Reference Strang1968),

where $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,L}^{\ast }$ denotes the adjoint of $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,L}$ , defined as

Explicitly, the Strang splitting can be written as

Let us note that the Lie splitting $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,L}$ and the Strang splitting $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,S2}$ are conjugate methods by the adjoint of $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,L}$ (McLachlan & Quispel Reference McLachlan and Quispel2002), i.e.

The last equality holds by the group property of the flow, but is only valid when the exact solution of each subsystem is used in the composition (and not some general symplectic or Poisson integrator). This implies that the Lie splitting shares many properties with the Strang splitting which are not found in general first-order methods.

Another second-order integrator with a smaller error constant than $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,S2}$ can be obtained by the following composition,

Here, $\unicode[STIX]{x1D6FC}$ is a free parameter which can be used to reduce the error constant. A particularly small error is obtained for $\unicode[STIX]{x1D6FC}=0.1932$ (McLachlan Reference McLachlan1995). Fourth-order time integrators can easily be obtained from a second-order integrator like $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,S}$ by the following composition (Suzuki Reference Suzuki1990; Yoshida Reference Yoshida1990),

with

Alternatively, we can compose the first-order integrator $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,L}$ together with its adjoint $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t,L}^{\ast }$ as follows (McLachlan Reference McLachlan1995),

with

For higher-order composition methods see e.g. Hairer et al. (Reference Hairer, Lubich and Wanner2006) and McLachlan & Quispel (Reference McLachlan and Quispel2002) and references therein.

5.3 Backward error analysis

In the following, we want to compute the modified Hamiltonian for the Lie–Trotter composition (5.17). For a splitting in two terms, $H=H_{A}+H_{B}$ , the Lie–Trotter composition can be written as

where $\boldsymbol{X}_{A}$ and $\boldsymbol{X}_{B}$ are the Hamiltonian vector fields corresponding to $H_{A}$ and $H_{B}$ , respectively, and $\tilde{\boldsymbol{X}}$ is the modified vector field corresponding to the modified Hamiltonian $\tilde{H}$ . Following Hairer et al. (Reference Hairer, Lubich and Wanner2006, Chapter IX), the modified Hamiltonian $\tilde{H}$ is given by

with

In order to compute the modified Hamiltonian for splittings with more than two terms, we have to apply the Baker–Campbell–Hausdorff formula recursively. Denoting by $\boldsymbol{X}_{E}$ the Hamiltonian vector field corresponding to $H_{E}$ , that is $\boldsymbol{X}_{E}=\mathbb{J}(\boldsymbol{\cdot },\mathit{D}H_{E})$ , and similar for $\boldsymbol{X}_{B}$ and $\boldsymbol{X}_{p_{i}}$ , the Lie–Trotter splitting (5.17) can be expressed in terms of these Hamiltonian vector fields as

Let us split $\unicode[STIX]{x1D711}_{\unicode[STIX]{x0394}t}=\exp (\unicode[STIX]{x0394}t\,\tilde{\boldsymbol{X}})$ into

where the corresponding Hamiltonians are given by

The Hamiltonian $\tilde{H}$ corresponding to $\tilde{X}$ is given by

with the first-order correction $\tilde{H}_{1}$ obtained as

The various Poisson brackets are computed as follows,

where $\mathbb{B}_{\unicode[STIX]{x1D707}}(\boldsymbol{X},\boldsymbol{b})$ denotes $N_{p}\times N_{p}$ diagonal matrix with elements $B_{h,\unicode[STIX]{x1D707}}(\boldsymbol{x}_{a})$ . The Lie–Trotter integrator (5.17) preserves the modified energy $\tilde{H}=H+\unicode[STIX]{x0394}t\,\tilde{H}_{1}$ to ${\mathcal{O}}(\unicode[STIX]{x0394}t^{2})$ , while the original energy $H$ is preserved only to ${\mathcal{O}}(\unicode[STIX]{x0394}t)$ .

6.1 Non-canonical Hamiltonian structure

Under the assumptions of (6.1) and (6.2), the bracket of (2.8) reduces to

where now

The Hamiltonian (2.11) becomes

The bracket of (6.9) with Hamiltonian (6.11), generates (6.4), (6.5) and (6.6), with the $\boldsymbol{J}$ of (6.8).

6.2 Reduced Jacobi identity

The bracket of (6.9) can be shown to satisfy the Jacobi identity by direct calculations. However since (2.8) satisfies the Jacobi identity for all functionals, it must also satisfy it for all reduced functionals. Here we have closure, i.e. if ${\mathcal{F}}$ and ${\mathcal{G}}$ are reduced functionals, then $\{{\mathcal{F}},{\mathcal{G}}\}$ is a reduced functional, where the bracket is the full bracket of (2.8).

More specifically, to understand this closure, observe that the reduced functionals have the form

where in (6.12) we assumed (6.1) and (6.2). In the second equality of (6.12) the integrations over $x_{2}$ , $x_{3}$ and $v_{3}$ are easily performed because the integrand is independent of these variables or it has been performed with an explicit dependence on $x_{2}$ , $x_{3}$ and $v_{3}$ that makes the integrals converge. Any constant factors resulting from the integrations are absorbed into the definition of $\bar{{\mathcal{F}}}$ . The closure condition on $\{{\mathcal{F}},{\mathcal{G}}\}$ amounts to the statement that given any two functionals ${\mathcal{F}}$ , ${\mathcal{G}}$ of the form of (6.12), their bracket is again a reduced functional of this form. This follows from the fact that the bracket of two such functionals reduces (2.8) to (6.9), which of course is reduced.

That not all reductions of functionals have closure can be seen by considering ones of the form ${\mathcal{F}}[\boldsymbol{E},f]$ , i.e. ones for which dependence on $\boldsymbol{B}$ is absent. The bracket of two functions of this form gives the bracket of (2.8) with the absence of the last term. Clearly this bracket depends on $\boldsymbol{B}$ and thus there is no closure. A consequence of this is that the bracket (2.8) with the absence of the last term does not satisfy the Jacobi identity. We note, however, that adding a projector can remedy this, as was shown in Chandre et al. (Reference Chandre, Guillebon, Back, Tassi and Morrison2013).

6.3 Discrete deRham complex in one dimension

Here, we consider the components of the electromagnetic fields separately and we have that $E_{1}$ is a 1-form, $E_{2}$ is a 0-form and $B_{3}$ is again a 1-form. We denote the semi-discrete fields by $D_{h}$ , $E_{h}$ and $B_{h}$ respectively, and write

Next we introduce an equidistant grid in $x$ and denote the spline of degree $p$ with support starting at $x_{i}$ by $N_{i}^{p}$ . We can express the derivative of $N_{i}^{p}$ as follows

In the finite element field solver, we represent $E_{h}$ by an expansion in splines of degree $p$ and $D_{h}$ , $B_{h}$ by an expansion in splines of degree $p-1$ , that is

6.4 Discrete Poisson bracket

The discrete Poisson bracket for this reduced system reads

Here, we denote by $\mathbb{M}_{p}=\unicode[STIX]{x1D648}_{p}$ and $\mathbb{M}_{q}=\unicode[STIX]{x1D648}_{q}$ the $N_{p}\times N_{p}$ diagonal matrices holding the particle masses and charges, respectively. We denote by $\unicode[STIX]{x1D726}^{0}(\boldsymbol{X})$ the $N_{p}\times N_{0}$ matrix with generic term $\unicode[STIX]{x1D6EC}_{i}^{0}(\mathit{x}_{a})$ , where $1\leqslant a\leqslant N_{p}$ and $1\leqslant i\leqslant N_{0}$ , and by $\unicode[STIX]{x1D726}^{1}(\boldsymbol{X})$ the $N_{p}\times N_{1}$ matrix with generic term $\unicode[STIX]{x1D6EC}_{i}^{1}(\mathit{x}_{a})$ , where $1\leqslant a\leqslant N_{p}$ and $1\leqslant i\leqslant N_{1}$ . Further, $\mathbb{B}(\boldsymbol{X},\boldsymbol{b})$ denotes the $N_{p}\times N_{p}$ diagonal matrix with entries

The reduced bracket can be shown to satisfy the Jacobi identity by direct proof in full analogy to the proof for the full bracket shown in § 4.4. However, one can also follow along the lines of § 6.2 in order to arrive at the same result.

6.5 Discrete Hamiltonian and equations of motion

The discrete Hamiltonian is given in terms of the reduced set of degrees of freedom $\boldsymbol{u}=(\boldsymbol{X},\boldsymbol{V}_{1},\boldsymbol{V}_{2},\boldsymbol{d},\boldsymbol{e},\boldsymbol{b})$ by

and the equations of motion are obtained as

which is seen to be in direct correspondence with (6.3)–(6.7).

6.6 Hamiltonian splitting

The solution of the discrete equations of motion for $H_{D}+H_{E}$ at time $\unicode[STIX]{x0394}t$ is

The solution of the discrete equations of motion for $H_{B}$ is

The solution of the discrete equations of motion for $H_{p_{1}}$ is

and for $H_{p_{2}}$ it is

respectively.

7.1 Weibel instability

We consider the Weibel instability studied in Weibel (Reference Weibel1959) in the form simulated in Crouseilles et al. (Reference Crouseilles, Einkemmer and Faou2015). We study a reduced 1d2v model with a perturbation along $x_{1}$ , a magnetic field along $x_{3}$ and electric fields along the $x_{1}$ and $x_{2}$ directions. Moreover, we assume that the distribution function is independent of $v_{3}$ . The initial distribution and fields are of the form

and $E_{1}(x,t=0)$ is computed from Poisson’s equation. In our simulations, we use the following choice of parameters, $\unicode[STIX]{x1D70E}_{1}=0.02/\sqrt{2}$ , $\unicode[STIX]{x1D70E}_{2}=\sqrt{12}\unicode[STIX]{x1D70E}_{1}$ , $k=1.25$ , $\unicode[STIX]{x1D6FC}=0$ and $\unicode[STIX]{x1D6FD}=-10^{-4}$ . Note that these are the same parameters as in Crouseilles et al. (Reference Crouseilles, Einkemmer and Faou2015) except for the fact that we sample from the Maxwellian without perturbation in $x_{1}$ .

The dispersion relation from Weibel (Reference Weibel1959) applied to our model reads

where $\unicode[STIX]{x1D719}(z)=\exp (-1/2z^{2})\int _{-\text{i}\infty }^{z}\exp (1/2\unicode[STIX]{x1D709}^{2})\,\text{d}\unicode[STIX]{x1D709}$ . For our parameter choice, this gives a growth rate of 0.02784. In figure 1, we show the electric and magnetic energies together with the analytic growth rate. We see that the growth rate is verified in the numerical solution. This simulation was performed with 100 000 particles, 32 grid points, splines of degree 3 and 2 and $\unicode[STIX]{x0394}t=0.05$ . Note that we have chosen a very large number of particles in order to obtain a solution of very high quality. In practice, the Weibel instability can also be simulated with much fewer particles (cf. § 7.5).

In table 1, we show the conservation properties of our splitting with various orders of the splitting (cf. § 5.2) and compare them also to the Boris–Yee scheme. The other numerical parameters are kept as before.

We can see that Gauss’ law is satisfied in each time step for the Hamiltonian splitting. This is a Casimir (cf. § 2.2) and therefore naturally conserved by the Hamiltonian splitting. On the other hand, this is not the case for the Boris–Yee scheme.

We can also see that the energy error improves with the order of the splitting; however, the Hamiltonian splitting method as well as the Boris–Yee scheme are not energy conserving. The time evolution of the total energy error is depicted in figure 2 for the various methods.

Figure 1.

Weibel instability: the two electric and the magnetic energies together with the analytic growth rate.

Figure 2.

Weibel instability: difference of total energy and its initial value as a function of time for various integrators.

Table 1.

Weibel instability: maximum error in the total energy and Gauss’ law until time 500 for simulation with various integrators: Lie–Trotter splitting from (5.17) (Lie), Strang splitting from (5.18) (Strang), second-order splitting with 4 Lie parts defined in (5.22) (2nd, 4 Lie), fourth-order splitting with 3 Strang parts defined in (5.23) (4th, 3 Strang) and 10 Lie parts defined in (5.25) (4th, 10 Lie).

7.2 Streaming Weibel instability

As a second test case, we consider the streaming Weibel instability. We study the same reduced model as in the previous section, but following Califano, Pegoraro & Bulanov (Reference Califano, Pegoraro and Bulanov1997), Cheng et al. (Reference Cheng, Gamba, Li and Morrison2014c ) the initial distribution and fields are prescribed as

and $E_{1}(x,t=0)$ is computed from Poisson’s equation.

We set the parameters to the following values $\unicode[STIX]{x1D70E}=0.1/\sqrt{2}$ , $k=0.2$ , $\unicode[STIX]{x1D6FD}=-10^{-3}$ , $v_{0,1}=0.5$ , $v_{0,2}=-0.1$ and $\unicode[STIX]{x1D6FF}=1/6$ . The parameters are chosen as in the case 2 of Cheng et al. (Reference Cheng, Gamba, Li and Morrison2014c ). The growth rate of energy of the second component of the electric field was determined to be 0.03 in Califano et al. (Reference Califano, Pegoraro and Bulanov1997). In figure 3, we show the electric and magnetic energies together with the analytic growth rate. We see that the growth rate is verified in the numerical solution. This simulation was performed on the domain $x\in [0,2\unicode[STIX]{x03C0}/k)$ with 20 000 000 particles, 128 grid points, splines of degree 3 and 2 and $\unicode[STIX]{x0394}t=0.01$ . Observe that the energy of the  $E_{1}$ component of the electric field starts to increase at times earlier than in Califano et al. (Reference Califano, Pegoraro and Bulanov1997), which is caused by particle noise.

As for the Weibel instability, we compare the conservation properties in table 2 for various integrators. Again we see that the Hamiltonian splitting conserves Gauss’ law as opposed to the Boris–Yee scheme. The energy conservation properties of the various schemes show approximately the same behaviour as in the previous test case (see also figure 4 for the time evolution of the energy error).

Table 2.

Streaming Weibel instability: maximum error in the total energy and Gauss’ law until time 200 for simulation with various integrators.

Figure 3.

Streaming Weibel instability: the two electric and the magnetic energies together with the analytic growth rate.

Figure 4.

Streaming Weibel instability: difference of total energy and its initial value as a function of time for various integrators.

7.3 Strong Landau damping

Table 3.

Damping and growth rates for strong Landau damping.

Finally, we also study the electrostatic example of strong Landau damping with initial distribution

The physical parameters are chosen as $\unicode[STIX]{x1D70E}=1$ , $k=0.5$ , $\unicode[STIX]{x1D6FC}=0.5$ and the numerical parameters as $\unicode[STIX]{x1D6E5}t=0.05$ , $n_{x}=32$ and 100 000 particles. The fields $B_{3}$ and $E_{2}$ are initialized to zero and remain zero over time. In this example, we essentially solve the Vlasov–Ampère equation with results equivalent to the Vlasov–Poisson equations. In figure 5 we show the time evolution of the electric energy associated with $E_{1}$ . We have also fitted a damping and growth rate (using the marked local maxima in the plot). These are in good agreement with other codes (see table 3). Again the energy conservation for the various method is visualized as a function of time in figure 6. And again we see that the fourth-order methods give excellent energy conservation.

Figure 5.

Landau damping: electric energy with fitted damping and growth rates.

Figure 6.

Landau damping: total energy error.

7.4 Backward error analysis

For the Lie–Trotter splitting, the error in the Hamiltonian $H$ is of order $\unicode[STIX]{x0394}t$ . However, using backward error analysis (cf. § 5.3), modified Hamiltonians can be computed, which are preserved to higher order. Accounting for first-order corrections $\tilde{H}_{1}$ , the error in the modified Hamiltonian,

is of order $(\unicode[STIX]{x0394}t)^{2}$ . For the 1d2v example, this correction is obtained as

with the various Poisson brackets computed as

In figure 7, we show the maximum and $\ell _{2}$ error of the energy and the corrected energy for the Weibel instability test case with the parameters in § 7.1. The simulations were performed with 100 000 particles, 32 grid points, splines of degree 3 and 2 and $\unicode[STIX]{x0394}t=0.01,0.02,0.05$ . We can see that the theoretical convergence rates are verified in the numerical experiments. Figure 8 shows the energy as well as the modified energy for the Weibel instability test case simulated with a time step of $0.05$ .

Figure 7.

Weibel instability: error (maximum norm and $\ell _{2}$ norm) in the total energy for simulations with $\unicode[STIX]{x0394}t=0.01,0.02,0.05$ .

Figure 8.

Weibel instability: energy and first-order corrected energy for simulation with $\unicode[STIX]{x0394}t=0.05$ .

7.5 Momentum conservation

Finally, let us discuss conservation of momentum for our one species 1d2v equations. In this case, equation (2.30) becomes

For $P_{e,1}(t)$ from (6.4) we get

Since in general $\int j_{1}(x,t)\,\text{d}x\neq 0$ , momentum is not conserved. As a means of testing momentum conservation, Crouseilles et al. (Reference Crouseilles, Einkemmer and Faou2015) replaced (6.4) by

Thus, for this artificial dynamics, momentum will be conserved if $\int E_{1}(x,0)\,\text{d}x=0$ . Instead, we do not modify the equations but check the validity of (7.12). We define a discrete version of (7.12), integrated over time, in the following way:

Our numerical scheme does not conserve momentum exactly. However, the error in momentum can be kept rather small during the linear phase of the simulation. Note that in all our examples, $\int j_{1,2}(x,t)\,\text{d}x=0$ . For a Gaussian initial distribution, the antithetic sampling ensures that $\int j_{1,2}(x,t)\,\text{d}x=0$ holds in the discrete sense. Figure 9 shows the momentum error in a simulation of the Weibel instability as considered in § 7.1, but up to time 2000 and with 25 600 particles sampled from pseudo-random numbers and Sobol numbers, for both plain and antithetical sampling. For the plain sampling, momentum error is a bit smaller for Sobol numbers compared to pseudo-random numbers during the linear phase. For the antithetic sampling, we can see that the momentum error is very small until time 200 (linear phase). However, when nonlinear effects start to dominate, the momentum error slowly increases until it has reached the same level as the momentum error for plain Sobol number sampling. The level depends on the number of particles (cf. table 4). Note that the sampling does not seem to have an influence on the energy conservation as can be seen in figure 10 that compares the energy error for the various sampling techniques. The curves show that the energy error is related to the increase in potential energy during the linear phase but does not further grow during the nonlinear phase.

Figure 9.

Weibel instability: error in the first component of the momentum for plain and antithetic Sobol sampling.

Figure 10.

Weibel instability: total energy error for plain and antithetic Sobol sampling.

Table 4.

Weibel instability: maximum error in both components of the momentum for simulations until time 2000 with various numbers of particles and 32 grid cells.

For the streaming Weibel instability on the other hand, we have a sum of two Gaussians in the second component of the velocity. Since in our sampling method we draw the particles based on Sobol quasi-random numbers, the fractions drawn from each of the Gaussians are not exactly given by $\unicode[STIX]{x1D6FF}$ and $1-\unicode[STIX]{x1D6FF}$ as in (7.5). Hence $\int j_{2}(x,t)\,\text{d}x$ is small but non-zero. Comparing the discrete momentum defined by (7.16) with the discretization of its definition (2.25),

we see a maximum deviation over time of $1.6\times 10^{-16}$ for a simulation with Strang splitting, $20\,000\,000$ particles, 128 grid points and $\unicode[STIX]{x0394}t=0.01$ . This shows that our discretization with a Strang splitting conserves (7.12) in the sense of (7.16) to very high accuracy.