Gauge invariance and conservation laws in the variational formulation of macro-particle plasma models

Higuet, Adam J.; Shadwick, B. A.

doi:10.1063/5.0245210

Recently, there has been significant interest in formulations of macro-particle models using variational methods. This is attractive because many of the inherent pathologies of traditional particle-in-cell methods are avoided. Broadly, there are two approaches to constructing these models that on the surface appear incompatible: a Lagrangian method based on the Low Lagrangian and a Hamiltonian method based on the noncanonical Vlasov–Maxwell Poisson bracket. In both cases, a reduction is performed on the distribution function, which is replaced by a finite sum of macro-particles with a fixed spatial structure and definite momentum. The distribution is, thus, replaced as dynamical variable by a collection of macro-particle positions and momenta. Here, we resolve the apparent discrepancy between the two approaches and show (1) that the noncanonical Hamiltonian formulation can be transformed into a canonical Hamiltonian system; (2) that this canonical Hamiltonian system is identical to that obtained from the Lagrangian by Legendre transformation; and (3) that introducing a grid in the Lagrangian setting requires careful treatment and an unusual structure in the variational principle to account for the breaking of gauge invariance.

I. INTRODUCTION

Recently, there has been renewed interest in the variational formulation of macro-particle kinetic plasma models. Early work^1,2 involved a reduction of the Low Lagrangian³ using point particles that suffer from significant computational noise. More modern efforts have been based on the phase space action,⁴ the Low Lagrangian,^5–7 or the noncanonical bracket.^8,9 Benefits of the variational approach include exact energy conservation, which precludes grid heating,^10–12 reduced noise,⁷ and a rigorous theoretical foundation based on a direct connection to the underlying Maxwell–Vlasov system. While representing the potentials (or fields) with a truncated basis set is compatible with any of these formulations, discretizing space using a grid, which is straightforward in the Lagrangian formulation, seems incompatible with the noncanonical Hamiltonian approach. The use of a grid is essential for computationally performant implementations on large multi-node computer systems. An outstanding question was what are phenomenological differences between macro-particle models obtained from the Lagrangian vs the noncanonical Hamiltonian formulation? Here, we show that once a macro-particle reduction is performed, the noncanonical system can be transformed into a canonical Hamiltonian system and then to a Lagrangian system that is identical to that obtained from the Low Lagrangian by the same macro-particle reduction. Thus, we conclude all formulations are identical. Furthermore, we demonstrate spatial discretization using a grid in both the Lagrangian and canonical Hamiltonian formulations. The latter is of significant interest as it allows for the use of symplectic integrators¹³ on the combined macro-particle and field phase space.^6,14 We show that the grid-based discretization breaks gauge invariance and, consequently, special care is needed in handling the variational principle. Finally, we give a general, minimalist basis formulation that preserves gauge invariance. As a special case, we demonstrate the use of a truncated Fourier representation for the potentials that preserves gauge invariance and conserves momentum.

II. LAGRANGIAN MACRO-PARTICLE REDUCTION

We consider a single mobile species with charge q, mass m, and phase space distribution f governed by the Vlasov–Maxwell equations and immobile ions with charge density

ρ^{ion}

(the extension to multiple mobile species is obvious). We replace the distribution function with a finite set of macro-particles^5–7 with fixed spacial structure, definite momentum, and weight

w_{α}

⁠,

f = \sum_{α} w_{α} S (r - ξ^{α}) δ (v - {\dot{ξ}}^{α}),

(1)

where

ξ^{α} (t)

and

{\dot{ξ}}^{α} (t)

are the macro-particle position and velocity, respectively, and

\int d^{3} r S (r) = 1,

(2)

so that

\sum_{α} w_{α}

is the total number of particles. While this reduction appears ad hoc, it can be viewed as a spatially “smoothed” version of the Klimontovich distribution.¹⁵ In this way, we replace f as the dynamical variable in favor of

ξ^{α} (t)

and

{\dot{ξ}}^{α} (t)

⁠. As we will see, the

ξ^{α}

are not simply characteristic curves of the Vlasov–Maxwell equation. Applying this reduction to the Low Lagrangian³ yields^5–7

\begin{matrix} L = - m c^{2} \sum_{α} w_{α} \sqrt{1 - {| \frac{{\dot{ξ}}^{α}}{c} |}^{2}} + q \sum_{α} w_{α} \int d^{3} r S (r - ξ^{α}) (\frac{{\dot{ξ}}^{α}}{c} \cdot A - φ) - \int d^{3} r ρ^{ion} φ \\ + \frac{1}{8 π} \int d^{3} r [\frac{1}{c^{2}} {| \frac{\partial A}{\partial t} |}^{2} + \frac{2}{c} \nabla φ \cdot \frac{\partial A}{\partial t} + {| \nabla φ |}^{2} - {| \nabla \times A |}^{2}], \end{matrix}

(3)

where

φ

and

A

are the scalar and vector potentials, respectively, with the electric and magnetic fields given in the usual way as follows:

E = - \frac{1}{c} \frac{\partial A}{\partial t} - \nabla φ,

(4a)

B = \nabla \times A .

(4b)

We can interpret

ρ = q \sum_{α} w_{α} S (r - ξ^{α})

(5a)

and

J = q \sum_{α} w_{α} {\dot{ξ}}^{α} S (r - ξ^{α})

(5b)

as the macro-particle charge and current densities, respectively, and can see that the continuity equation

\frac{\partial ρ}{\partial t} + \nabla \cdot J = 0

(6)

is satisfied by virtue of the form of the reduction . The reduced Lagrangian inherits energy and momentum conservation and gauge-invariance from the underlying Vlasov–Maxwell system.

Taking the action

I = \int d t L

to be stationary yields the equations of motion in the usual way,

\frac{d}{d t} (γ_{α} m {\dot{ξ}}^{α}) = q \int d^{3} r S (r - ξ^{α}) [- \frac{1}{c} \frac{\partial A}{\partial t} - \nabla φ + \frac{{\dot{ξ}}^{α}}{c} \times (\nabla \times A)],

(7a)

\frac{\partial^{2} A}{\partial t^{2}} + c^{2} \nabla \times (\nabla \times A) + c \nabla \frac{\partial φ}{\partial t} = 4 π c J,

(7b)

and

- \nabla^{2} φ - \frac{1}{c} \nabla \cdot \frac{\partial A}{\partial t} = 4 π (ρ + ρ^{ion}),

(7c)

where

γ_{α} = \frac{1}{\sqrt{1 - {| {\dot{ξ}}^{α} / c |}^{2}}} .

(8)

We can see that is simply the Lorentz force averaged over the spatial extent of the macro-particle while is Ampère's Law and is Gauss's Law. Had we taken the macro-particles to be points in space,¹ i.e.,

S (r - ξ^{α}) = δ^{3} (r - ξ^{α})

⁠, then the macro-particle equations would be the characteristic curves of the Vlasov–Maxwell system. Using point particles, however, has unusably high noise levels [similar to the nearest grid-point charge deposition in the standard particle-in-cell (PIC) algorithm^16,17].

It is straightforward to show that leads to conservation of energy

\begin{matrix} E = m c^{2} \sum_{α} w_{α} γ_{α} + \frac{1}{8 π} \int d^{3} r ({| \frac{1}{c} \frac{\partial A}{\partial t} + \nabla φ |}^{2} + {| \nabla \times A |}^{2}) \end{matrix}

(9)

and total momentum

\begin{matrix} P = m \sum_{α} w_{α} γ_{α} {\dot{ξ}}^{α} - \frac{1}{c} \int d^{3} r ρ^{ion} A \\ - \frac{1}{4 π c} \int d^{3} r (\frac{1}{c} \frac{\partial A}{\partial t} + \nabla φ) \times (\nabla \times A) . \end{matrix}

(10)

Finally, as will become important below, we see that

p^{α}

⁠, the momentum conjugate to

ξ^{α}

⁠, is given by

p^{α} = \frac{\partial L}{\partial {\dot{ξ}}^{α}} = w_{α} γ_{α} m {\dot{ξ}}^{α} + \frac{q}{c} w_{α} \int d^{3} r S (r - ξ^{α}) A,

(11)

while the momentum conjugate to

A

is given by

\begin{matrix} Π = \frac{δ L}{δ (\partial A / \partial t)} = \frac{1}{4 π c} (\frac{1}{c} \frac{\partial A}{\partial t} + \nabla φ) = - \frac{E}{4 π c} . \end{matrix}

(12)

Of course, there is no momentum conjugate to

φ

which complicates the transformation to canonical coordinates below.

III. REDUCTION OF THE NONCANONICAL SYSTEM

The Vlasov–Maxwell system has a well-known noncanonical structure^18–21 with Poisson bracket

\begin{matrix} {\tilde{F}, \tilde{G}} = \int d^{3} r d^{3} p f [\frac{δ \tilde{F}}{δ f}, \frac{δ \tilde{G}}{δ f}] \\ + 4 π q \int d^{3} r d^{3} p \nabla_{p} f \cdot (\frac{δ \tilde{F}}{δ E} \frac{δ \tilde{G}}{δ f} - \frac{δ \tilde{G}}{δ E} \frac{δ \tilde{F}}{δ f}) \\ + 4 π c \int d^{3} r (\frac{δ \tilde{F}}{δ E} \cdot \nabla \times \frac{δ \tilde{G}}{δ B} - \frac{δ \tilde{G}}{δ E} \cdot \nabla \times \frac{δ \tilde{F}}{δ B}) \\ + \frac{q}{c} \int d^{3} r d^{3} p f B \cdot (\nabla_{p} \frac{δ \tilde{F}}{δ f} \times \nabla_{p} \frac{δ \tilde{G}}{δ f}), \end{matrix}

(13)

where

[a, b] = [\frac{\partial a}{\partial r} \cdot \frac{\partial b}{\partial p} - \frac{\partial a}{\partial p} \cdot \frac{\partial b}{\partial r}]

(14)

and

\nabla_{p} a = \frac{\partial a}{\partial p}

(15)

and Hamiltonian

H = \int d^{3} r d^{3} p γ m c^{2} f + \frac{1}{8 π} \int d^{3} r ({| E |}^{2} + {| B |}^{2}) .

(16)

The evolution of the dynamical variables is given by

\begin{matrix} \frac{\partial f}{\partial t} = {f, H}, \\ \frac{\partial E}{\partial t} = {E, H}, \\ \frac{\partial B}{\partial t} = {B, H} . \end{matrix}

(17)

This system has the Casimir invariants [null vectors of ]

C_{f} = \int d^{3} r d^{3} p α_{f} (f),

(18)

C_{E} = \int d^{3} r d^{3} p α_{E} (r) [\nabla \cdot E - 4 π (ρ + ρ^{ion})],

(19)

and

C_{B} = \int d^{3} r d^{3} p α_{B} (r) \nabla \cdot B

(20)

for arbitrary functions

α_{f}

⁠,

α_{E}

⁠, and

α_{B}

⁠. For to satisfy the Jacobi identity,

\nabla \cdot B = 0

and thus the only acceptable value is

C_{B} = 0

⁠. Since

E

and

B

are constrained by Gauss's Law and

\nabla \cdot B = 0

⁠, respectively, these constraints should be added to the Hamiltonian with Lagrange multipliers.²² These terms would be identical to the Casimirs

C_{E}

and

C_{B}

and, thus, could not contribute to the equations of motion and so can be ignored.

A. Macro-particle reduction

We now again reduce the phase space distribution function by introducing macro-particles as follows:

f = \sum_{α} f_{α} = \sum_{α} w_{α} q S (r - ξ^{α}) δ (p - γ_{α} m {\dot{ξ}}^{α}) .

(21)

To proceed, it appears necessary to treat each macro-particle as a separate species. Here, we replace functionals of f by corresponding functions of the macro-particle variables. We can view the macro-particle variables as moments of the distribution function

w_{α} = \int d^{3} r d^{3} p f_{α},

(22a)

ξ^{α} = \frac{1}{w_{α}} \int d^{3} r d^{3} p r f_{α},

(22b)

and

π^{α} = \int d^{3} r d^{3} p p f_{α},

(22c)

from which we find

\frac{δ w_{α}}{δ f_{α}} = 1,

(23a)

\frac{δ ξ^{α}}{δ f_{α}} = \frac{r - ξ^{α}}{w_{α}},

(23b)

and

\frac{δ π^{α}}{δ f_{α}} = p .

(23c)

For any functional

\tilde{F}

of the distribution function and the fields, there exists an equivalent hybrid object

\hat{F}

that is a function of the macro-particle variables and a functional of the field variables

\tilde{F} [f, E, B] = \hat{F} (w_{α}, ξ^{α}, π^{α}; E, B) .

(24)

The functional derivative with respect to the phase space distribution function can be expressed as

\begin{matrix} \frac{δ \tilde{F}}{δ f_{α}} = \frac{δ w_{α}}{δ f_{α}} \frac{\partial \hat{F}}{\partial w_{α}} + \frac{δ ξ^{α}}{δ f_{α}} \cdot \frac{\partial \hat{F}}{\partial ξ^{α}} + \frac{δ π^{α}}{δ f_{α}} \cdot \frac{\partial \hat{F}}{\partial π^{α}} \\ = \frac{\partial \hat{F}}{\partial w_{α}} + \frac{r - ξ^{α}}{w_{α}} \cdot \frac{\partial \hat{F}}{\partial ξ^{α}} + p \cdot \frac{\partial \hat{F}}{\partial π^{α}} . \end{matrix}

(25)

Now,

\begin{matrix} [\frac{δ \tilde{F}}{δ f_{α}}, \frac{δ \tilde{G}}{δ f_{α}}] = [\frac{\partial \hat{F}}{\partial w_{α}} + \frac{r - ξ^{α}}{w_{α}} \cdot \frac{\partial \hat{F}}{\partial ξ^{α}} + p \cdot \frac{\partial \hat{F}}{\partial π^{α}}, \frac{\partial \hat{G}}{\partial w_{α}} + \frac{r - ξ^{α}}{w_{α}} \cdot \frac{\partial \hat{G}}{\partial ξ^{α}} + p \cdot \frac{\partial \hat{G}}{\partial π^{α}}] \\ = \frac{\partial \hat{F}}{\partial ξ^{α}} \cdot [\frac{r - ξ^{α}}{w_{α}}, p] \cdot \frac{\partial \hat{G}}{\partial π^{α}} + \frac{\partial \hat{F}}{\partial π^{α}} \cdot [p, \frac{r - ξ^{α}}{w_{α}}] \cdot \frac{\partial \hat{G}}{\partial ξ^{α}} \\ = \frac{1}{w_{α}} (\frac{\partial \hat{F}}{\partial ξ^{α}} \cdot \frac{\partial \hat{G}}{\partial π^{α}} - \frac{\partial \hat{G}}{\partial ξ^{α}} \cdot \frac{\partial \hat{F}}{\partial π^{α}}) \\ = \frac{1}{w_{α}} {[\hat{F}, \hat{G}]}_{(ξ^{α}, π^{α})}, \end{matrix}

(26)

giving

\begin{matrix} \int d^{3} r d^{3} p f [\frac{δ \tilde{F}}{δ f}, \frac{δ \tilde{G}}{δ f}] = \sum_{α} \int d^{3} r d^{3} p f_{α} [\frac{δ \tilde{F}}{δ f_{α}}, \frac{δ \tilde{G}}{δ f_{α}}] \\ = \frac{1}{w_{α}} \sum_{α} {[\hat{F}, \hat{G}]}_{(ξ^{α}, π^{α})} \int d^{3} r d^{3} p f_{α} = \sum_{α} {[\hat{F}, \hat{G}]}_{(ξ^{α}, π^{α})} . \end{matrix}

(27)

Furthermore, from , we see

\nabla_{p} \frac{δ \tilde{F}}{δ f} = \frac{\partial \hat{F}}{\partial π^{α}}

(28)

and, thus,

\begin{matrix} \int d^{3} r d^{3} p \nabla_{p} f \cdot \frac{δ \tilde{F}}{δ E} \frac{δ \tilde{G}}{δ f} = - \int d^{3} r d^{3} p f \frac{δ \tilde{F}}{δ E} \cdot \nabla_{p} \frac{δ \tilde{G}}{δ f} \\ = - \int d^{3} r S (r - ξ^{α}) \frac{δ \tilde{F}}{δ E} \cdot \frac{\partial \hat{G}}{\partial π^{α}} . \end{matrix}

(29)

Similarly,

\begin{matrix} \int d^{3} r d^{3} p f B \cdot (\nabla_{p} \frac{δ \tilde{F}}{δ f} \times \nabla_{p} \frac{δ \tilde{G}}{δ f}) \\ = \sum_{α} w_{α} (\frac{\partial \hat{F}}{\partial π^{α}} \times \frac{\partial \hat{G}}{\partial π^{α}}) \cdot \int d^{3} r S (r - ξ^{α}) B . \end{matrix}

(30)

Combining , , and , we find the reduced bracket to be

\begin{matrix} {\hat{F}, \hat{G}} = \sum_{α} {[\hat{F}, \hat{G}]}_{(ξ^{α}, π^{α})} \\ - 4 π q \sum_{α} w_{α} \int d^{3} r S (r - ξ^{α}) (\frac{δ \hat{F}}{δ E} \cdot \frac{\partial \hat{G}}{\partial π^{α}} - \frac{δ \hat{G}}{δ E} \cdot \frac{\partial \hat{F}}{\partial π^{α}}) \\ + \frac{q}{c} \sum_{α} w_{α} (\frac{\partial \hat{F}}{\partial π^{α}} \times \frac{\partial \hat{G}}{\partial π^{α}}) \cdot \int d^{3} r S (r - ξ^{α}) B \\ + 4 π c \int d^{3} r (\frac{δ \hat{F}}{δ E} \cdot \nabla \times \frac{δ \hat{G}}{δ B} - \frac{δ \hat{G}}{δ E} \cdot \nabla \times \frac{δ \hat{F}}{δ B}) . \end{matrix}

(31)

Under this reduction, the Casimirs all become proportional to

\sum_{α} w_{α}

[due to the structure of , the only admissible

α_{f} (f)

are linear in f].

From , we see

π^{α} = w_{α} γ_{α} m {\dot{ξ}}^{α}

(32)

and, thus, we can write the reduced Hamiltonian as

H = m c^{2} \sum_{α} w_{α} γ_{α} + \frac{1}{8 π} \int d^{3} r ({| E |}^{2} + {| B |}^{2}),

(33)

where

γ_{α} = \sqrt{1 + \frac{{| π^{α} |}^{2}}{w_{α}^{2} m^{2} c^{2}}} .

(34)

Using and , we find

{\dot{ξ}}^{α} = \frac{π^{α}}{w_{α} γ_{α} m},

(35a)

\begin{matrix} {\dot{π}}^{α} = w_{α} q \int d^{3} r S (r - ξ^{α}) (E + \frac{π^{α}}{m c w_{α} γ_{α}} \times B) \\ = w_{α} q \int d^{3} r S (r - ξ^{α}) (E + \frac{{\dot{ξ}}^{α}}{c} \times B), \end{matrix}

(35b)

\begin{matrix} \frac{\partial E}{\partial t} = c \nabla \times B - 4 π \frac{q}{m} \sum_{α} \frac{π^{α}}{γ_{α}} S (r - ξ^{α}) = c \nabla \times B - 4 π J, \end{matrix}

(35c)

and

\frac{\partial B}{\partial t} = - c \nabla \times E .

(35d)

These equations of motion are identical to the equations of motion obtained from the reduced Lagrangian where Gauss's Law follows from Ampère's Law and the continuity equation .

B. Transformation to canonical coordinates

Previously, it has been shown⁵ that in one-dimension, the macro-particle reduction of the nonrelativistic Vlasov–Ampère system yields dynamical equations equivalent to those obtained from the Lagrangian formulation. Here, we show the equivalence of the two formulations in the general case. The presence of the magnetic field greatly complicates the analysis.

We begin by replacing the electric field with

Π

using and magnetic fields with the vector potential. Thus, for any

\hat{F}

there is an equivalent

\bar{F}

that is a function of the macro-particle quantities and a functional of

A

and

Π,

\hat{F} (w_{α}, ξ^{α}, π^{α}; E, B) = \bar{F} (w_{α}, ξ^{α}, π^{α}; A, Π) .

(36)

Now,

\frac{δ \hat{F}}{δ E} = - \frac{1}{4 π c} \frac{δ \bar{F}}{δ Π}

(37)

and

\nabla \times \frac{δ \hat{F}}{δ B} = \frac{δ \bar{F}}{δ A} .

(38)

Writing in terms of our new functionals gives

\begin{matrix} {\bar{F}, \bar{G}} = \sum_{α} {[\bar{F}, \bar{G}]}_{(ξ^{α}, π^{α})} \\ + \frac{q}{c} \sum_{α} w_{α} \int d^{3} r S (r - ξ^{α}) (\frac{δ \bar{F}}{δ Π} \cdot \frac{\partial \bar{G}}{\partial π^{α}} - \frac{δ \bar{G}}{δ Π} \cdot \frac{\partial \bar{F}}{\partial π^{α}}) \\ + \frac{q}{c} \sum_{α} w_{α} (\frac{\partial \bar{F}}{\partial π^{α}} \times \frac{\partial \bar{G}}{\partial π^{α}}) \cdot \int d^{3} r S (r - ξ^{α}) \nabla \times A \\ + \int d^{3} r (\frac{δ \bar{F}}{δ A} \cdot \frac{δ \bar{G}}{δ Π} - \frac{δ \bar{G}}{δ A} \cdot \frac{δ \bar{F}}{δ Π}) . \end{matrix}

(39)

In these variables,

C_{E}

, is no longer a Casimir and, of course,

C_{B} = 0

⁠. We still have the Gauss's Law constraint on

Π

⁠; thus, we must add a Lagrange multiplier to the Hamiltonian which will be adjusted to enforce the constraint. In the new coordinates, the Hamiltonian is

\begin{matrix} H = m c^{2} \sum_{α} w_{α} γ_{α} + \frac{1}{8 π} \int d^{3} r ({| \nabla \times A |}^{2} + 16 π^{2} c^{2} {| Π |}^{2}) \\ + \int d^{3} r α_{E} (c \nabla \cdot Π + ρ + ρ^{ion}), \end{matrix}

(40)

where

α_{E}

is the Lagrange multiplier.

Using and , we find

{\dot{ξ}}^{α} = \frac{π^{α}}{w_{α} γ_{α} m},

(41a)

{\dot{π}}^{α} = w_{α} q \int d^{3} r S (r - ξ^{α}) [- 4 π c Π + \frac{π^{α}}{m c w_{α} γ_{α}} \times (\nabla \times A)],

(41b)

\frac{\partial Π}{\partial t} = \frac{1}{4 π} \nabla \times (\nabla \times A) - \frac{q}{c} \sum_{α} \frac{π^{α}}{m γ_{α}} S (r - ξ^{α}),

(41c)

and

\frac{\partial A}{\partial t} = 4 π c^{2} Π - c \nabla α_{E} .

(41d)

Using the definition of

Π

, we see that becomes

\frac{\partial E}{\partial t} = c \nabla \times B - 4 π J,

(42)

leads to

E = - \frac{1}{c} \frac{\partial A}{\partial t} - \nabla α_{E},

(43)

and becomes

\frac{d}{d t} (γ_{α} m {\dot{ξ}}^{α}) = \int d^{3} r S (r - ξ^{α}) (E + \frac{{\dot{ξ}}^{α}}{c} \times B) .

(44)

The Lagrange multiplier is determined from the Gauss's Law constraint and to be

\nabla^{2} α_{E} + \frac{1}{c} \nabla \cdot \frac{\partial A}{\partial t} = - 4 π (ρ + ρ^{ion}) .

(45)

Clearly, the Lagrange multiplier plays the role of the scalar potential. It is only through the Gauss's Law constraint that

φ

enters as there is no corresponding dynamical equation. This leads to the natural interpretation that the role of the scalar potential it to ensure that Gauss's Law is satisfied.

We now make one last change of variables, replacing

π^{α}

with

p^{α}

⁠. For any

\bar{F}

⁠, there is an equivalent F that is a function of the macro-particle canonical coordinates and momenta and a functional of

A

and

Π

⁠,

\bar{F} (w_{α}, ξ^{α}, π^{α}; A, Π) = F (w_{α}, ξ^{α}, p^{α}; A, Π) .

(46)

From , we have

p^{α} = π^{α} + \frac{w_{α} q}{c} \int d^{3} r S (r - ξ^{α}) A,

(47)

giving

\begin{matrix} δ p^{α} = δ π^{α} + \frac{w_{α} q}{c} \int d^{3} r [A δ S (r - ξ^{α}) + S (r - ξ^{α}) δ A] . \end{matrix}

(48)

Now,

\begin{matrix} \int d^{3} r A δ S (r - ξ^{α}) = \int d^{3} r A δ ξ_{i}^{α} \frac{\partial S (r - ξ^{α})}{\partial ξ_{i}^{α}} \\ = - \int d^{3} r A δ ξ_{i}^{α} \partial_{i} S (r - ξ^{α}) \\ = \int d^{3} r S (r - ξ^{α}) δ ξ_{i}^{α} \partial_{i} A, \end{matrix}

(49)

where there is an implied sum over repeated roman indices and

\partial_{j} F

are the components of

\nabla F

⁠. Thus, becomes

δ p^{α} = δ π^{α} + \frac{w_{α} q}{c} \int d^{3} r S (r - ξ^{α}) (\partial_{i} A δ ξ_{i}^{α} + δ A) .

(50)

Here, the chain rule is significantly more complicated so comparing

δ \bar{F}

to

δ F

provides a more straightforward means to relate the derivatives. To this end, we compute

\begin{matrix} δ \bar{F} = \sum_{α} (\frac{\partial \bar{F}}{\partial ξ^{α}} \cdot δ ξ^{α} + \frac{\partial \bar{F}}{\partial π^{α}} \cdot δ π^{α}) + \int d^{3} r (\frac{δ \bar{F}}{δ A} \cdot δ A + \frac{δ \bar{F}}{δ Π} \cdot δ Π) \end{matrix}

(51)

and

\begin{matrix} δ F = \sum_{α} (\frac{\partial F}{\partial ξ^{α}} \cdot δ ξ^{α} + \frac{\partial F}{\partial p^{α}} \cdot δ p^{α}) + \int d^{3} r (\frac{δ F}{δ A} \cdot δ A + \frac{δ F}{δ Π} \cdot δ Π) . \end{matrix}

(52)

We now substitute into to obtain

\begin{matrix} δ F = \sum_{α} \frac{\partial F}{\partial ξ_{i}^{α}} δ ξ_{i}^{α} + \sum_{α} \frac{\partial F}{\partial p_{i}^{α}} [δ π_{i}^{α} + \frac{w_{α} q}{c} \int d^{3} r S (r - ξ^{α}) (\partial_{j} A_{i} δ ξ_{j}^{α} + δ A_{i})] + \int d^{3} r (\frac{δ F}{δ A} \cdot δ A + \frac{δ F}{δ Π} \cdot δ Π) \\ = \sum_{α} [\frac{\partial F}{\partial ξ_{i}^{α}} + \frac{w_{α} q}{c} \frac{\partial F}{\partial p_{j}^{α}} \int d^{3} r S (r - ξ^{α}) \partial_{i} A_{j}] δ ξ_{i}^{α} + \sum_{α} \frac{\partial F}{\partial p^{α}} \cdot δ π^{α} + \int d^{3} r [\frac{δ F}{δ A} + \frac{q}{c} \sum_{α} w_{α} S (r - ξ^{α}) \frac{\partial F}{\partial p^{α}}] \cdot δ A + \int d^{3} r \frac{δ F}{δ Π} \cdot δ Π, \end{matrix}

(53)

which is compared to to determine

\frac{\partial \bar{F}}{\partial ξ_{i}^{α}} = \frac{\partial F}{\partial ξ_{i}^{α}} + \frac{w_{α} q}{c} \frac{\partial F}{\partial p_{j}^{α}} \int d^{3} r S (r - ξ^{α}) \partial_{i} A_{j},

(54a)

\frac{\partial \bar{F}}{\partial π^{α}} = \frac{\partial F}{\partial p^{α}},

(54b)

\frac{δ \bar{F}}{δ A} = \frac{δ F}{δ A} + \frac{q}{c} \sum_{α} w_{α} S (r - ξ^{α}) \frac{\partial F}{\partial p^{α}},

(54c)

\frac{δ \bar{F}}{δ Π} = \frac{δ F}{δ Π} .

(54d)

We use the above to transform . We find that

\begin{matrix} {[\bar{F}, \bar{G}]}_{(ξ^{α}, π^{α})} = [\frac{\partial F}{\partial ξ_{i}^{α}} + \frac{w_{α} q}{c} \frac{\partial F}{\partial p_{j}^{α}} \int d^{3} r S (r - ξ^{α}) \partial_{i} A_{j}] \frac{\partial G}{\partial p_{i}^{α}} \\ - [\frac{\partial G}{\partial ξ_{i}^{α}} + \frac{w_{α} q}{c} \frac{\partial G}{\partial p_{j}^{α}} \int d^{3} r S (r - ξ^{α}) \partial_{i} A_{j}] \frac{\partial F}{\partial p_{i}^{α}} \\ = {[F, G]}_{(ξ^{α}, p^{α})} + \frac{w_{α} q}{c} \frac{\partial F}{\partial p_{i}^{α}} \frac{\partial G}{\partial p_{j}^{α}} \int d^{3} r S (r - ξ^{α}) (\partial_{j} A_{i} - \partial_{i} A_{j}) . \end{matrix}

(55)

The second term in transforms trivially as

\begin{matrix} \int d^{3} r S (r - ξ^{α}) (\frac{δ \bar{F}}{δ Π} \cdot \frac{\partial \bar{G}}{\partial π^{α}} - \frac{δ \bar{G}}{δ Π} \cdot \frac{\partial \bar{F}}{\partial π^{α}}) \\ = \int d^{3} r S (r - ξ^{α}) (\frac{δ F}{δ Π} \cdot \frac{\partial G}{\partial p^{α}} - \frac{δ G}{δ Π} \cdot \frac{\partial F}{\partial p^{α}}) . \end{matrix}

(56)

The third term also transforms trivially

\begin{matrix} (\frac{\partial \bar{F}}{\partial π^{α}} \times \frac{\partial \bar{G}}{\partial π^{α}}) \cdot \int d^{3} r S (r - ξ^{α}) \nabla \times A \\ = (\frac{\partial F}{\partial p^{α}} \times \frac{\partial G}{\partial p^{α}}) \cdot \int d^{3} r S (r - ξ^{α}) \nabla \times A . \end{matrix}

(57)

The last term becomes

\begin{matrix} \int d^{3} r (\frac{δ \bar{F}}{δ A} \cdot \frac{δ \bar{G}}{δ Π} - \frac{δ \bar{G}}{δ A} \cdot \frac{δ \bar{F}}{δ Π}) \\ = \int d^{3} r {\frac{δ G}{δ Π} \cdot [\frac{δ F}{δ A} + \frac{q}{c} \sum_{α} w_{α} S (r - ξ^{α}) \frac{\partial F}{\partial p^{α}}] - \frac{δ F}{δ Π} \cdot [\frac{δ G}{δ A} + \frac{q}{c} \sum_{α} w_{α} S (r - ξ^{α}) \frac{\partial G}{\partial p^{α}}]} \\ = \int d^{3} r (\frac{δ F}{δ A} \cdot \frac{δ G}{δ Π} - \frac{δ G}{δ A} \cdot \frac{δ F}{δ Π}) + \frac{q}{c} \sum_{α} w_{α} \int d^{3} r S (r - ξ^{α}) (\frac{\partial F}{\partial p^{α}} \cdot \frac{δ G}{δ Π} - \frac{\partial G}{\partial p^{α}} \cdot \frac{δ F}{δ Π}) . \end{matrix}

(58)

Putting the pieces together, becomes

\begin{matrix} {F, G} = \sum_{α} {[F, G]}_{(ξ^{α}, p^{α})} + \frac{q}{c} \sum_{α} w_{α} \frac{\partial F}{\partial p_{i}^{α}} \frac{\partial G}{\partial p_{j}^{α}} \int d^{3} r S (r - ξ^{α}) (\partial_{j} A_{i} - \partial_{i} A_{j}) \\ + \frac{q}{c} \sum_{α} w_{α} \int d^{3} r S (r - ξ^{α}) (\frac{δ F}{δ Π} \cdot \frac{\partial G}{\partial p^{α}} - \frac{δ G}{δ Π} \cdot \frac{\partial F}{\partial p^{α}}) + \frac{q}{c} \sum_{α} w_{α} (\frac{\partial \bar{F}}{\partial p^{α}} \times \frac{\partial \bar{G}}{\partial p^{α}}) \cdot \int d^{3} r S (r - ξ^{α}) \nabla \times A \\ + \int d^{3} r (\frac{δ F}{δ A} \cdot \frac{δ G}{δ Π} - \frac{δ G}{δ A} \cdot \frac{δ F}{δ Π}) + \frac{q}{c} \sum_{α} w_{α} \int d^{3} r S (r - ξ^{α}) (\frac{\partial F}{\partial p^{α}} \cdot \frac{δ G}{δ Π} - \frac{\partial G}{\partial p^{α}} \cdot \frac{δ F}{δ Π}) \\ = \sum_{α} {[F, G]}_{(ξ^{α}, p^{α})} + \int d^{3} r (\frac{δ F}{δ A} \cdot \frac{δ G}{δ Π} - \frac{δ G}{δ A} \cdot \frac{δ F}{δ Π}) \\ + \frac{q}{c} \sum_{α} w_{α} [\frac{\partial F}{\partial p_{i}^{α}} \frac{\partial G}{\partial p_{j}^{α}} \int d^{3} r S (r - ξ^{α}) (\partial_{j} A_{i} - \partial_{i} A_{j}) + (\frac{\partial \bar{F}}{\partial p^{α}} \times \frac{\partial \bar{G}}{\partial p^{α}}) \cdot \int d^{3} r S (r - ξ^{α}) \nabla \times A] . \end{matrix}

(59)

The last term in has the following form:

\begin{matrix} a_{i} b_{j} (\partial_{j} F_{i} - \partial_{i} F_{j}) + ε_{ijk} a_{i} b_{k} ε_{ilm} \partial_{l} F_{m} \\ = a_{i} b_{j} (\partial_{j} F_{i} - \partial_{i} F_{j}) + (δ_{j l} δ_{k m} - δ_{j m} δ_{k l}) a_{i} b_{k} \partial_{l} F_{m} \\ = a_{i} b_{j} (\partial_{j} F_{i} - \partial_{i} F_{j}) + a_{j} b_{k} \partial_{j} F_{k} - a_{j} b_{k} \partial_{k} F_{j} = 0, \end{matrix}

(60)

where

ε_{ijk}

is the usual 3D Levi–Civita symbol and, thus, becomes

\begin{matrix} {F, G} = \sum_{α} {[F, G]}_{(ξ^{α}, p^{α})} + \int d^{3} r (\frac{δ F}{δ A} \cdot \frac{δ G}{δ Π} - \frac{δ G}{δ A} \cdot \frac{δ F}{δ Π}), \end{matrix}

(61)

which is clearly canonical. Since we have performed a series of exact coordinate transformations leading to a canonical bracket, it must be the case that and are “good” Poisson brackets and, therefore, must satisfy the Jacobi identity. The Hamiltonian is given by but with the understanding that

π^{α}

is a function of

ξ^{α}

and

p^{α}

and functional of

A

through

π^{α} = p^{α} - \frac{w_{α} q}{c} \int d^{3} r S (r - ξ^{α}) A .

(62)

The equations of motion for the fields are

\begin{matrix} \frac{\partial A}{\partial t} = \frac{δ H}{δ Π} = 4 π c^{2} Π - c \nabla α_{E} \end{matrix}

(63a)

and

\begin{matrix} \frac{\partial Π}{\partial t} = - \frac{δ H}{δ A} = - \frac{1}{4 π} \nabla \times (\nabla \times A) + \frac{q}{c} \sum_{α} \frac{π^{α}}{m γ_{α}} S (r - ξ^{α}) . \end{matrix}

(63b)

Using , , and , these become

E = - \frac{1}{c} \frac{\partial A}{\partial t} - \nabla α_{E}

(64a)

and

\frac{\partial E}{\partial t} = c \nabla \times (\nabla \times A) - 4 π J .

(64b)

The macro-particle equations of motion are more involved. As we expect,

{\dot{ξ}}^{α} = \frac{\partial H}{\partial p^{α}} = \frac{π^{α}}{γ_{α} m w_{α}} .

(65)

Now,

\begin{matrix} {\dot{p}}^{α} = - \frac{\partial H}{\partial ξ^{α}} = \frac{q}{c} \frac{π_{k}^{a}}{m γ_{α}} \int d^{3} r \frac{\partial S (r - ξ^{α})}{\partial ξ^{α}} A_{k} \\ - q w_{α} \int d^{3} r α_{E} \frac{\partial S (r - ξ^{α})}{\partial ξ^{α}} \\ = \frac{q}{c} \frac{π_{k}^{a}}{m γ_{α}} \int d^{3} r S (r - ξ^{α}) \nabla A_{k} \\ - q w_{α} \int d^{3} r S (r - ξ^{α}) \nabla α_{E} . \end{matrix}

(66)

From , we have

{\dot{p}}^{α} = {\dot{π}}^{α} + \frac{q}{c} w_{α} \int d^{3} r [\frac{\partial S (r - ξ^{α})}{\partial ξ_{i}^{α}} {\dot{ξ}}_{i}^{α} A + S (r - ξ^{α}) \frac{\partial A}{\partial t}] .

(67)

Combining these equations, we have

\begin{matrix} \frac{d}{d t} (m γ_{α} {\dot{ξ}}^{α}) = q \int d^{3} r S (r - ξ^{α}) (- \frac{1}{c} \frac{\partial A}{\partial t} - \nabla α_{E}) + \frac{q}{c} \int d^{3} r S (r - ξ^{α}) ({\dot{ξ}}_{i}^{α} \nabla A_{i} - {\dot{ξ}}_{i}^{α} \partial_{i} A) \\ = q \int d^{3} r S (r - ξ^{α}) (- \frac{1}{c} \frac{\partial A}{\partial t} - \nabla α_{E}) + \frac{q}{c} \int d^{3} r S (r - ξ^{α}) {\dot{ξ}}^{α} \times (\nabla \times A) . \end{matrix}

(68)

The Lagrange multiplier is determined from the constraint by taking the divergence of to give

\nabla \cdot (- \nabla α_{E} - \frac{1}{c} \frac{\partial A}{\partial t}) = 4 π (ρ + ρ^{ion}) .

(69)

It is easy to see that these equations of motion are identical to with the identification of

α_{E}

with

φ

⁠.

C. Connection to the Lagrangian formulation

We now take the Legendre transform

L = \sum_{α} {\dot{ξ}}^{α} \cdot p^{α} + \int d^{3} r \frac{\partial A}{\partial t} \cdot Π - H

(70)

to obtain a Lagrangian, which we can then compare to . From the Hamiltonian, we have

{\dot{ξ}}^{α} = \frac{p^{α}}{w_{α} γ_{α} m} - \frac{q}{γ_{α} m c} \int d^{3} r S (r - ξ^{α}) A

(71a)

and

\frac{\partial A}{\partial t} = 4 π c^{2} Π + 4 π c \nabla α_{E},

(71b)

which we use to eliminate the canonical momenta, viz.,

p^{α} = w_{α} γ_{α} m {\dot{ξ}}^{α} + \frac{w_{α} q}{c} \int d^{3} r S (r - ξ^{α}) A

(72)

and

Π = \frac{1}{4 π c} (\frac{1}{c} \frac{\partial A}{\partial t} + \nabla α_{E}),

(73)

where

γ_{α}

is given by . A straightforward calculation yields

\begin{matrix} L = - m c^{2} \sum_{α} w_{α} \sqrt{1 - {| \frac{{\dot{ξ}}^{α}}{c} |}^{2}} \\ + q \sum_{α} w_{α} \int d^{3} r S (r - ξ^{α}) (\frac{{\dot{ξ}}^{α}}{c} \cdot A - α_{E}) \\ + \frac{1}{8 π} \int d^{3} r (\frac{1}{c^{2}} {| \frac{\partial A}{\partial t} |}^{2} + \frac{2}{c} \nabla α_{E} \cdot \frac{\partial A}{\partial t} + {| \nabla α_{E} |}^{2} - {| \nabla \times A |}^{2}), \end{matrix}

(74)

which we see is identical to the reduced Lagrangian if we choose

α_{E} = φ .

(75)

This choice re-introduces the scalar potential to the Lagrangian system derived from the noncanonical formulation and completes the equivalence between the two systems. This reinforces the view that the scalar potential can be interpreted of as a Lagrange multiplier rather than a dynamical quantity.

D. Constructing the primary Hamiltonian

We have shown that the noncanonical system is identical to the system with Lagrangian . We now perform a Legendre transform to obtain a canonical Hamiltonian system from the Lagrangian to demonstrate a well-known²² but nonetheless interesting quirk of electrodynamics. Determining the canonical momenta for the Lagrangian in the standard way reveals immediately that the momentum conjugate to the scalar potential,

Π_{φ}

⁠, is identically zero. Thus, we cannot write

Π_{φ}

as a function of

\dot{φ}

and we have a constrained system, which requires the use of the Dirac–Bergmann algorithm^22–25 to transform the canonical system even though no gauge condition has been chosen. Here,

Π_{φ} = 0

acts as a primary constraint. In the usual way, we replace the time-derivatives of the coordinates and fields using . We construct the constrained Hamiltonian as

H_{c} = \sum_{α} p^{α} \cdot {\dot{ξ}}^{α} + \int d^{3} r Π \cdot \frac{\partial A}{\partial t} + \int d^{3} r Π_{φ} \dot{φ} - L,

(76)

which is evaluated according to our primary constraint and to give

\begin{matrix} H_{c} = m c^{2} \sum_{α} w_{α} γ_{α} + \frac{1}{8 π} \int d^{3} r (16 π^{2} c^{2} {| Π |}^{2} + {| \nabla \times A |}^{2}) \\ + \int d^{3} r φ (c \nabla \cdot Π + ρ + ρ^{ion}) . \end{matrix}

(77)

Now,

H_{c}

is restricted to the sub-manifold of the phase space defined by the primary constraint; notice that this Hamiltonian matches the Hamiltonian . Furthermore, we have the Gauss's Law (which plays the role of a constitutive relation) constraining the canonical field momentum

c \nabla \cdot Π + ρ + ρ^{ion} = 0.

(78)

To define a Hamiltonian on all phase space, we continue the Dirac–Bergmann algorithm by adding the constraints with Lagrange multipliers to

H_{c}

to obtain the primary Hamiltonian

\begin{matrix} H_{p} = H_{c} + \int d^{3} r μ_{1} (r) Π_{φ} \\ + \int d^{3} r μ_{2} (r) (c \nabla \cdot Π + ρ + ρ^{ion}) . \end{matrix}

(79)

Explicitly, we have

\begin{matrix} H_{p} = m c^{2} \sum_{α} w_{α} γ_{α} \\ + \frac{1}{8 π} \int d^{3} r (16 π^{2} c^{2} {| Π |}^{2} + {| \nabla \times A |}^{2}) \\ + \int d^{3} r μ_{1} Π_{φ} + \int d^{3} r (φ + μ_{2}) [c \nabla \cdot Π + ρ + ρ^{ion}] . \end{matrix}

(80)

We also have to extend the Poisson bracket to full phase space in accordance with Dirac–Bergmann,

\begin{matrix} {F, G} = \sum_{α} {[F, G]}_{(ξ^{α}, p^{α})} + \int d^{3} r (\frac{δ F}{δ A} \cdot \frac{δ G}{δ Π} - \frac{δ G}{δ A} \cdot \frac{δ F}{δ Π}) \\ + \int d^{3} r (\frac{δ F}{δ φ} \frac{δ G}{δ Π_{φ}} - \frac{δ F}{δ Π_{φ}} \frac{δ G}{δ φ}) . \end{matrix}

(81)

This system yields equations of motion identical to those obtained from and .

The primary Hamiltonian (80) and extended bracket (81) is the correct canonical formulation according to the Dirac–Bergmann algorithm, yet the noncanonical system only produced us the constrained Hamiltonian and its corresponding bracket. Since the noncanonical formulation begins with $E$ and $B$ ⁠, it has no need to represent the entire phase space for the potentials and, as can be seen by its Casimirs, is only valid on the constrained phase space. In the noncanonical formulation, we do not need $φ$ to be a dynamical variable but simply “some scalar function that makes Gauss's Law true.” When starting from the Lagrangian, where the system is derived in terms of potentials, the actual dynamics of $φ$ are significant and can be carried into the Hamiltonian picture via the primary Hamiltonian.

IV. GAUGE INVARIANCE DISCRETIZATION: TRUNCATED FOURIER BASIS

Previous work has shown that using a Fourier representation for the potential results in exact momentum conservation in the electrostatic case.⁵ It has also been shown,⁶ in one-dimension, that if the potentials are expanded in a truncated basis that exactly represents the derivative operator, then a discrete from of the continuity holds and the Lagrangian is gauge invariant and momentum is exactly conserved.

We now introduce a truncated Fourier series representation of the potentials^5,6,14,26 which, as expected, preserves both gauge and translation invariance. This is a special case of the general basis expansion shown in the Appendix. To simplify the exposition, we take the domain to be a cube of size L and take the same number of modes in each direction. The generalization to an arbitrary Cartesian domain with different truncations in each direction is straightforward. We define the Fourier representation over a periodic domain of length L as

F (r) = \sum_{I} {\tilde{F}}_{I} (k_{I}) e^{i k_{I} \cdot r},

(82)

where I represents a triple of integers from the set

N = [- N, N] \times [- N, N] \times [- N, N]

with the sum ranging over all possible triples and

k_{I} = (k_{x}, k_{y}, k_{z}) = \frac{2 π}{L} (l, m, n),

(83)

where

(l, m, n) \in N

⁠. The amplitudes are defined

{\tilde{F}}_{I} (k_{I}) = \frac{1}{L^{3}} {\int \int \int}_{- L / 2}^{L / 2} d x d y d z e^{- i k_{I} \cdot r} F (r) .

(84)

Since we are dealing only with real-valued functions,

{\tilde{F}}_{I} = {\tilde{F}}_{- I}^{*} .

(85)

Using this representation, the Lagrangian becomes

\begin{matrix} L = - m c^{2} \sum_{α} w_{α} \sqrt{1 - {| \frac{{\dot{ξ}}^{α}}{c} |}^{2}} - L^{3} \sum_{I} [({\tilde{ρ}}_{I}^{*} + {\tilde{ρ}}_{I}^{ion}^{*}) {\tilde{φ}}_{I} \\ - \frac{1}{c} {\tilde{J}}_{I}^{*} \cdot {\tilde{A}}_{I}] + \frac{L^{3}}{8 π} \sum_{I} [\frac{1}{c^{2}} {| {\dot{\tilde{A}}}_{I} |}^{2} + \frac{2}{c} i k_{I} {\tilde{φ}}_{I}^{*} \cdot {\dot{\tilde{A}}}_{I} \\ + k_{I}^{2} {| {\tilde{φ}}_{I} |}^{2} - {| k_{I} \times {\tilde{A}}_{I} |}^{2}], \end{matrix}

(86)

where

{\tilde{ρ}}_{I} = q \sum_{α} w_{α} {\tilde{S}}_{I} e^{- i k_{I} \cdot ξ^{α}}

(87a)

and

{\tilde{J}}_{I} = q \sum_{α} w_{α} {\dot{ξ}}^{α} {\tilde{S}}_{I} e^{- i k_{I} \cdot ξ^{α}} .

(87b)

The equations of motion are

\begin{matrix} \frac{d}{d t} (m γ_{α} {\dot{ξ}}^{α}) = q L^{3} \sum_{I} {\tilde{S}}_{I}^{*} e^{i k_{I} \cdot ξ^{α}} [- \frac{1}{c} {\dot{\tilde{A}}}_{I} - i k_{I} {\tilde{φ}}_{I} \\ + \frac{{\dot{ξ}}^{α}}{c} \times (i k_{I} \times {\tilde{A}}_{I})], \end{matrix}

(88a)

{\ddot{\tilde{A}}}_{I} + i c k_{I} {\dot{\tilde{φ}}}_{I} + c^{2} k_{I} \times (k_{I} \times {\tilde{A}}_{I}) = 4 π c {\tilde{J}}_{I},

(88b)

and

k_{I}^{2} {\tilde{φ}}_{I} - \frac{i}{c} k_{I} \cdot {\dot{\tilde{A}}}_{I} = 4 π ({\tilde{ρ}}_{I} + {\tilde{ρ}}_{I}^{ion}) .

(88c)

These equations are identical to with the truncated expansion substituted for the potentials. In addition, if we compare the time-derivative of with the dot product of

i k_{I} / c

and , we find

{\dot{\tilde{ρ}}}_{I} + i k_{I} \cdot {\tilde{J}}_{I} = 0

(89)

[which also follows immediately from ]. This is what we expect based on the analysis in Sec. V; this is a discretization that preserves gauge invariance and, thus, in the absence of an imposed gauge condition, the equations of motion should imply the appropriate discrete continuity equation.

The Lagrangian retains the translation symmetry of and, therefore, there is a conserved momentum. The symmetry transformation is

ξ^{α} \to ξ^{α} + Δ,

(90a)

{\tilde{A}}_{I} \to {\tilde{A}}_{I} {exp}^{- i k_{I} \cdot Δ} \approx (1 - i k_{I} \cdot Δ) {\tilde{A}}_{I},

(90b)

and

{\tilde{φ}}_{I} \to {\tilde{φ}}_{I} {exp}^{- i k_{I} \cdot Δ} \approx (1 - i k_{I} \cdot Δ) {\tilde{φ}}_{I} .

(90c)

The conserved momentum given by Noether's Theorem is

\begin{matrix} P = \sum_{α} \frac{\partial L}{\partial {\dot{ξ}}^{α}} + i \sum_{I} k_{I} [{\tilde{A}}_{I} \cdot \frac{\partial L}{\partial {\dot{\tilde{A}}}_{I}} + {\tilde{φ}}_{I} \frac{\partial L}{\partial {\dot{\tilde{φ}}}_{I}}] \\ = \sum_{α} w_{α} (γ_{α} m {\dot{ξ}}^{α} + \frac{q L^{3}}{c} \sum_{I} {\tilde{S}}_{I}^{*} {\tilde{A}}_{I} e^{i k_{I} \cdot ξ^{α}}) \\ + i \frac{L^{3}}{4 π c} \sum_{I} k_{I} {\tilde{A}}_{I}^{*} \cdot (\frac{1}{c} {\dot{\tilde{A}}}_{I} + i k_{I} {\tilde{φ}}_{I}) . \end{matrix}

(91)

Using the equations of motion, this can be simplified to give

\begin{matrix} P = \sum_{α} w_{α} γ_{α} m {\dot{ξ}}^{α} - \frac{L^{3}}{c} \sum_{I} {\tilde{A}}_{I}^{*} {\tilde{ρ}}_{I}^{ion} \\ + i \frac{L^{3}}{4 π c} \sum_{I} (\frac{1}{c} {\dot{\tilde{A}}}_{I} + i k_{I} {\tilde{φ}}_{I}) \times (k_{I} \times {\tilde{A}}_{I}^{*}), \end{matrix}

(92)

in agreement with .

V. CHARGE CONSERVATION AND GAUGE SELECTION IN NON-GAUGE INVARIANT SYSTEMS

For the Lagrangian (3), the continuity equation is structural in the reduction (1) and not dynamic as it is, say, in the general Vlasov system. With the noncanonical formulation seen in Sec. III B, Gauss's Law does not appear as an equation of motion but rather as a constraint on the system and instead follows from the equations of motion as a consequence of $ρ$ and $J$ satisfying the continuity equation. While discretizations such as that in Sec. IV can maintain this property, representing the fields on a grid, which is desirable for computational performance, generically results in the discrete version of the continuity equation not being exactly satisfied.⁶ As we will see below, the failure of the continuity equation to hold implies the system is not gauge-invariant.

A. Varying a non-gauge invariant Lagrangian

Variations of the action obtained from with respect to

A

and

φ

yield

\begin{matrix} δ I = - \frac{1}{4 π c^{2}} \int d t \int d^{3} r {c^{2} δ φ [\nabla^{2} φ + \frac{1}{c} \nabla \cdot \frac{\partial A}{\partial t} \\ + 4 π (ρ + ρ^{ion})] + δ A \cdot [\frac{\partial^{2} A}{\partial t^{2}} - c^{2} \nabla^{2} A \\ + c^{2} \nabla (\nabla \cdot A) + c \nabla \frac{\partial φ}{\partial t} - 4 π c J]} \end{matrix}

(93)

(the variation with respect to

ξ^{α}

has been ignored as it is irrelevant to gauge invariance). These variations yield Gauss's Law

\nabla^{2} φ + \frac{1}{c} \nabla \cdot \frac{\partial A}{\partial t} = - 4 π (ρ + ρ^{ion}),

(94a)

and Ampère's law

\frac{\partial^{2} A}{\partial t^{2}} - c^{2} \nabla^{2} A + c^{2} \nabla (\nabla \cdot A) + c \nabla \frac{\partial φ}{\partial t} = 4 π c J,

(94b)

as the equations of motion. The definitions of the electric and magnetic fields in terms of potentials give the other two Maxwell equations. Subtracting the time derivative of from the divergence of gives the continuity equation consistent with our definitions of

ρ

and

J

⁠.

We now explicitly decompose the general variations

δ φ

and

δ A

into a part that can be expressed as a gauge transformation and the remainder. The simplest representation is

δ φ = - \frac{1}{c} \frac{\partial δ λ}{\partial t}

(95a)

and

δ A = \nabla δ λ + δ C,

(95b)

where

δ λ

and

δ C

are arbitrary. We can, thus, rearrange as

\begin{matrix} δ I = - \frac{1}{4 π c^{2}} \int d t \int d^{3} r {4 π c δ λ (\frac{\partial ρ}{\partial t} + \nabla \cdot J) \\ + δ C \cdot [\frac{\partial^{2} A}{\partial t^{2}} - c^{2} \nabla^{2} A + c^{2} \nabla (\nabla \cdot A) \\ + c \nabla \frac{\partial φ}{\partial t} - 4 π c J]} . \end{matrix}

(96)

The equations of motion for the above variation are

\frac{\partial ρ}{\partial t} + \nabla \cdot J = 0

(97a)

and

\frac{\partial^{2} A}{\partial t^{2}} - c^{2} \nabla^{2} A + c^{2} \nabla (\nabla \cdot A) + c \nabla \frac{\partial φ}{\partial t} = 4 π c J,

(97b)

and, similar to before, we can take the divergence of Ampère's Law and substitute in the other equation of motion to find

\frac{\partial}{\partial t} (\nabla^{2} φ + \frac{1}{c} \nabla \cdot \frac{\partial A}{\partial t}) = - 4 π \frac{\partial ρ}{\partial t},

(98)

which is the time derivative of Gauss's Law. Since we only have the time derivative of Gauss's Law in this formulation, the stationary ions only enter the dynamics via initial conditions. Taking

δ C = 0

results in the variations corresponding to a gauge transformation and so we see that the continuity equation holding is a prerequisite for the action to be invariant under gauge transformations. By this splitting of variations, we see that either the continuity equation must be satisfied or the variations must be made in a way that keeps

δ λ = 0

⁠. When the method of discretization results in the discrete continuity equation not holding exactly (and, thus, breaks gauge invariance), we can no longer vary

φ

without introducing contradictions into the equations of motion; since the equations of motion can be combined to yield the continuity equation this is an obvious contradiction when the continuity equation is no longer true. It has previously been shown⁶ that, in one-dimension, a contradiction existed between Ampère's Law and Gauss's Law unless a discrete version of the continuity equation held. Here, we have now shown that this contradiction arises from the gauge transformation “part” of the variation and setting

δ λ = 0

yields the correct dynamics.

Gauss's Law and charge conservation are effectively tied to one another in the continuous system. As was seen in the Lagrangian (74), the multiplier that enters the system to produce Gauss's Law in the Hamiltonian picture transforms into the structure needed to allow gauge transformations in the Lagrangian picture. Thus, when the discretization violates local charge conservation, the dynamical equations no longer imply Gauss's Law.

B. Enforcing Gauss's Law via Lagrange multipliers

Knowing that the method of discretization can change the variational principle such that Gauss's Law is no longer an equation of motion, we seek to keep it in the system by adding it to the Lagrangian as a constraint

\begin{matrix} L = - m c^{2} \sum_{α} w_{α} \sqrt{1 - {| \frac{{\dot{ξ}}^{α}}{c} |}^{2}} + q \sum_{α} w_{α} \int d^{3} r S (r - ξ^{α}) (\frac{{\dot{ξ}}^{α}}{c} \cdot A - φ) - \int d^{3} r ρ^{ion} φ \\ + \frac{1}{8 π} \int d^{3} r [\frac{1}{c^{2}} {| \frac{\partial A}{\partial t} |}^{2} + \frac{2}{c} \nabla φ \cdot \frac{\partial A}{\partial t} + {| \nabla φ |}^{2} - {| \nabla \times A |}^{2}] - \frac{1}{4 π} \int d^{3} r μ [\nabla^{2} φ + \frac{1}{c} \nabla \cdot \frac{\partial A}{\partial t} + 4 π (ρ + ρ^{ion})] . \end{matrix}

(99)

We remove gauge transformations from the variations and, thus, vary with respect only to

C

and

μ

(and again ignore the variations in

ξ^{α}

as they are not important to the discussion) giving

\begin{matrix} δ I = \frac{1}{4 π c^{2}} \int d t \int d^{3} r {c^{2} δ μ [\nabla^{2} φ + \frac{1}{c} \nabla \cdot \frac{\partial A}{\partial t} \\ + 4 π (ρ + ρ^{ion})] - δ C \cdot [\frac{\partial^{2} A}{\partial t^{2}} - c^{2} \nabla^{2} A + c^{2} \nabla (\nabla \cdot A) \\ + c \nabla (\frac{\partial φ}{\partial t} + \frac{\partial μ}{\partial t}) - 4 π c J]} . \end{matrix}

(100)

This obviously re-introduces Gauss's Law as an equation of motion with the consequence of the addition of

\partial μ / \partial t

in Ampère's Law. We combine the two equations of motion as before to find

\nabla^{2} \dot{μ} = 4 π (\frac{\partial ρ}{\partial t} + \nabla \cdot J) .

(101)

Of course for a gauge invariant system that preserves the continuity equation, the right hand side of the above is zero and the only reasonable solution is

μ = 0

⁠. This is to be expected since a gauge invariant system already has Gauss's Law as an equation of motion so the Lagrange multiplier enforcing it is unnecessary. For a non-gauge invariant system, however, we see that the error in the continuity equation enters as a source for the Lagrange multiplier's equation of motion and modifies Ampère's Law.

VI. DISCRETIZING ON A GRID

We consider a discretization of by representing the potentials on a grid. When projected onto the grid, the charge and current densities, , and will no longer satisfy the discrete form of the continuity equation⁶ and the discrete Lagrangian will not be gauge invariant. As we have seen above, we must restrict variations of the potentials to exclude gauge transformations. The simplest way to accomplish this is to choose a gauge that does not impose constraints amongst the potentials; here we will use the Weyl gauge,

φ = 0

and enforce Gauss's Law using a Lagrange multiplier; see Sec. V B. We take a uniform grid with

N_{x}

⁠,

N_{y}

⁠, and

N_{z}

grid points in the x, y, and z directions, respectively, grid spacing

Δ x

⁠,

Δ y

⁠, and

Δ z

⁠, and extent

L_{x}

⁠,

L_{y}

⁠, and

L_{z}

⁠. We enumerate the grid points as

x_{i}

⁠,

0 \leq i < N_{x}

⁠,

y_{j}

⁠,

0 \leq j < N_{y}

⁠, and

z_{k}

⁠,

0 \leq k < N_{z}

⁠. For any function of

F (r)

⁠, we denote the approximate value at a grid point as

F_{ijk} \approx F (x_{i}, y_{j}, z_{k})

⁠. As much as possible we will suppress indices and sums to reduce clutter in expressions. The relevant finite difference operators needed are

\nabla_{N}

and

\nabla_{N}^{2}

⁠, the discrete approximations of

\nabla

and

\nabla^{2}

⁠, respectively. We define

{(\nabla_{N} F)}_{ijk} = {(D^{x} F)}_{ijk} \hat{x} + {(D^{y} F)}_{ijk} \hat{y} + {(D^{z} F)}_{ijk} \hat{z},

(102)

where

{(D^{x} F)}_{ijk} = \sum_{i^{'}} D_{i i^{'}}^{x} F_{i^{'} j k},

(103a)

{(D^{y} F)}_{ijk} = \sum_{i^{'}} D_{j j^{'}}^{y} F_{i j^{'} k},

(103b)

{(D^{z} F)}_{ijk} = \sum_{k^{'}} D_{k k^{'}}^{z} F_{i j k^{'}}

(103c)

are central difference approximations to the first derivative in x, y, and z, respectively, and, consequently, are antisymmetric:

D_{i j}^{a} = - D_{j i}^{a}

⁠. Similarly, we define

{(\nabla_{N}^{2} F)}_{ijk} = {(K^{x} F)}_{ijk} + {(K^{y} F)}_{ijk} + {(K^{z} F)}_{ijk},

(104)

where

{(K^{x} F)}_{ijk} = \sum_{i^{'}} K_{i i^{'}}^{x} F_{i^{'} j k},

(105a)

{(K^{y} F)}_{ijk} = \sum_{i^{'}} K_{j j^{'}}^{y} F_{i j^{'} k},

(105b)

{(K^{z} F)}_{ijk} = \sum_{k^{'}} K_{k k^{'}}^{z} F_{i j k^{'}}

(105c)

are central difference approximations to the second derivative in x, y, and z, respectively, and, consequently, are symmetric:

K_{i j}^{a} = K_{j i}^{a}

for

a \in x, y, z

⁠. For example, a second order approximation yields

D_{i i^{'}}^{x} = \frac{δ_{i + 1, i^{'}} - δ_{i - 1, i^{'}}}{2 Δ x}

(106)

and

K_{i i^{'}}^{x} = \frac{δ_{i + 1, i^{'}} - 2 δ_{i, i^{'}} + δ_{i - 1, i^{'}}}{Δ x^{2}} .

(107)

The divergence and curl have a natural approximation in terms of ,

{(\nabla_{N} \cdot F)}_{ijk} = {(D^{a} F^{a})}_{ijk},

(108a)

{(\nabla_{N} \times F)}_{ijk} = ε_{abc} e_{a} {(D^{b} F^{c})}_{ijk},

(108b)

where

e_{a}

are the Cartesian unit vectors. A consequence of these definitions is that important vector calculus identities, which are geometric in nature, are respected by these finite difference operators, namely,

{[\nabla_{N} \times (\nabla_{N} F)]}_{ijk} = 0,

(109a)

{[\nabla_{N} \cdot (\nabla_{N} \times F)]}_{ijk} = 0,

(109b)

and

\begin{matrix} {[\nabla_{N} \times (\nabla_{N} \times F)]}_{ijk} = {[\nabla_{N} (\nabla_{N} \cdot F)]}_{ijk} - {[(\nabla_{N} \cdot \nabla_{N}) F]}_{ijk} . \end{matrix}

(109c)

Whether

[(\nabla_{N} \cdot \nabla_{N}) F]

reduces to

(\nabla_{N}^{2} F)

depends on the details of

D_{i j}^{a}

and

K_{i j}^{a}

but does not hold for commonly used compact approximations for the first and second derivative regardless of order. The symmetry properties of D and K lead to the identities

\sum_{ijk} {(F \nabla_{N} G)}_{ijk} = - \sum_{ijk} {(G \nabla_{N} F)}_{ijk}

(110a)

and

\sum_{ijk} {(F \nabla_{N}^{2} G)}_{ijk} = \sum_{ijk} {(G \nabla_{N}^{2} F)}_{ijk},

(110b)

where the sum extends over all grid points and we take either periodic boundary conditions or assume that F and G vanish at the edge of the domain. An immediate consequence of is

\begin{matrix} \sum_{ijk} {(F \nabla_{N} F)}_{ijk} = 0 . \end{matrix}

(110c)

In general, there is no finite difference analog of the product rule.

A. Lagrangian formulation

We now discretize in the Weyl gauge and add a constraint to enforce Gauss's Law. Furthermore, we replace the integration by a sum and expand

{| \nabla \times A |}^{2}

before approximating. This gives

\begin{matrix} L = - m c^{2} \sum_{α} w_{α} \sqrt{1 - {| \frac{{\dot{ξ}}^{α}}{c} |}^{2}} + \frac{q}{c} \sum_{α} w_{α} \sum_{ijk} ϱ_{ijk}^{α} {\dot{ξ}}^{α} \cdot A_{ijk} + \frac{Δ}{8 π} \sum_{ijk} {\frac{1}{c^{2}} {| {\dot{A}}_{ijk} |}^{2} \\ + {(A \cdot \nabla_{N}^{2} A)}_{ijk} + {[{(\nabla_{N} \cdot A)}_{ijk}]}^{2}} - Δ \sum_{ijk} μ_{ijk} [\frac{1}{4 π c} {(\nabla_{N} \cdot \dot{A})}_{ijk} + ρ_{ijk} + ρ_{ijk}^{ion}], \end{matrix}

(111)

where the sums extend over all grid points,

Δ = Δ x Δ y Δ z

⁠,

ρ_{ijk} = \frac{q}{Δ} \sum_{α} w_{α} ϱ_{ijk}^{α}

(112)

and

ϱ_{ijk}^{α}

is the interpolated particle shape evaluated at

ξ^{α}

and is assumed to be product of one-dimensional shapes (see Ref. 5). The equations of motion are

\begin{matrix} \frac{d}{d t} (γ_{α} m {\dot{ξ}}^{α}) = - q \sum_{ijk} {μ_{ijk} \frac{\partial ϱ_{ijk}^{α}}{\partial ξ^{α}} + \frac{1}{c} ϱ_{ijk}^{α} {\dot{A}}_{ijk} \\ + \frac{{\dot{ξ}}^{α}}{c} \times [\nabla_{ξ} \times (ϱ_{ijk}^{α} A_{ijk})]}, \end{matrix}

(113a)

and

\begin{matrix} {\ddot{A}}_{ijk} + c {(\nabla_{N} \dot{μ})}_{ijk} - c^{2} {(\nabla_{N}^{2} A)}_{ijk} + c^{2} {[\nabla_{N} (\nabla_{N} \cdot A)]}_{ijk} = 4 π c J_{ijk}, \end{matrix}

(113b)

where

J_{ijk} = \frac{q}{Δ} \sum_{α} w_{α} {\dot{ξ}}^{α} ϱ_{ijk}^{α} .

(114)

The Lagrange multiplier satisfies

\begin{matrix} {[(\nabla_{N} \cdot \nabla_{N}) \dot{μ}]}_{ijk} = c {\nabla_{N} \cdot [\nabla_{N}^{2} A - \nabla_{N} (\nabla_{N} \cdot A)]}_{ijk} \\ + \frac{4 π q}{Δ} \sum_{α} w_{α} {\dot{ξ}}^{α} \cdot [\frac{\partial ϱ_{ijk}^{α}}{\partial ξ^{α}} + {(\nabla_{N} ϱ^{α})}_{ijk}] . \end{matrix}

(115)

The first term on the right of is a measure of the discrepancy between

\nabla_{N}^{2}

and

\nabla_{N} . \nabla_{N}

and vanishes as the grid spacing is reduced in accordance with the truncation errors of the finite-difference approximations. The second term is measure of the error in the discrete version of the continuity equation and does not generally vanish as the grid is refined. This is a consequence of the usual practice of tying the particle shape to the grid spacing. This term will tend to zero as the number of grid points within the support of the particle shape function grows.

B. Canonical formulation

When discretizing the canonical Hamiltonian formulation on a grid, we expect the Poisson bracket to become

\begin{matrix} {F, G} = \sum_{α} {[F, G]}_{(ξ^{α}, p^{α})} + \sum_{ijk} (\frac{\partial F}{\partial A_{ijk}} \cdot \frac{\partial G}{\partial Π_{ijk}} - \frac{\partial G}{\partial A_{ijk}} \cdot \frac{\partial F}{\partial Π_{ijk}}), \end{matrix}

(116)

as we have a finite degree-of-freedom system. If we examine the unconstrained version of , we see that the momentum conjugate to

A

becomes

Π_{ijk} = \frac{Δ}{4 π c^{2}} {\dot{A}}_{ijk} = - \frac{Δ}{4 π c} E_{ijk},

(117)

which differs from by a factor of the discrete volume element. Working in the Weyl gauge eliminates the constraint associated with

Π_{φ}

⁠, but the Gauss's Law constraint still needs to be considered. The primary Hamiltonian is now

\begin{matrix} H_{p} = m c^{2} \sum_{α} w_{α} γ_{α} + \frac{Δ}{8 π} \sum_{ijk} {\frac{16 π^{2} c^{2}}{Δ^{2}} {| Π_{ijk} |}^{2} \\ - {(A \cdot \nabla_{N}^{2} A)}_{ijk} - {[{(\nabla_{N} \cdot A)}_{ijk}]}^{2}} + \sum_{ijk} μ_{ijk} [c {(\nabla_{N} \cdot Π)}_{ijk} + Δ (ρ_{ijk} + ρ_{ijk}^{ion})], \end{matrix}

(118)

where

γ_{α}

is given by and

π^{α} = p^{α} - \frac{w_{α} q}{c} \sum_{ijk} ϱ_{ijk}^{α} A_{ijk} .

(119)

The equations of motion for the fields are

{\dot{A}}_{ijk} = \frac{4 π c^{2}}{Δ} Π_{ijk} - c {(\nabla_{N} μ)}_{ijk},

(120a)

{\dot{Π}}_{ijk} = \frac{Δ}{4 π} {{(\nabla_{N}^{2} A)}_{ijk} - {[\nabla_{N} (\nabla_{N} \cdot A)]}_{ijk}} + \frac{Δ}{c} J_{ijk},

(120b)

which can be combined to give

\begin{matrix} {\ddot{A}}_{ijk} = - c {(\nabla_{N} \dot{μ})}_{ijk} + c^{2} {(\nabla_{N}^{2} A)}_{ijk} - c^{2} {[\nabla_{N} (\nabla_{N} \cdot A)]}_{ijk} + 4 π c J_{ijk}, \end{matrix}

(121)

which is identical to . The constraint is evaluated using to give

\begin{matrix} {[(\nabla_{N} \cdot \nabla_{N}) \dot{μ}]}_{ijk} = c {\nabla_{N} \cdot [\nabla_{N}^{2} A - \nabla_{N} (\nabla_{N} \cdot A)]}_{ijk} \\ + 4 π [{(\nabla_{N} \cdot J)}_{ijk} + \frac{\partial ρ_{ijk}}{\partial t}], \end{matrix}

(122)

in agreement with . Again, the equation of motion for the macro-particle quantities are more involved. We find

{\dot{ξ}}^{α} = \frac{π^{α}}{γ_{α} m w_{α}}

(123)

and

{\dot{p}}^{α} = q w_{α} \sum_{ijk} (- μ_{ijk} + \frac{{\dot{ξ}}_{b}^{α}}{c} A_{b, ijk}) \frac{\partial ϱ_{ijk}^{α}}{\partial ξ^{α}},

(124)

while from , we have

{\dot{p}}^{α} = {\dot{π}}^{α} + \frac{q w_{α}}{c} \sum_{ijk} (ϱ_{ijk}^{α} {\dot{A}}_{ijk} + \frac{\partial ϱ_{ijk}^{α}}{\partial ξ_{a}^{α}} {\dot{ξ}}_{a}^{α} A_{ijk}) .

(125)

Combining these two expressions, we arrive at (113a). While the canonical and Lagrangian formulations give identical dynamics, the canonical formulation allow for the use of symplectic integrators¹³ on the combined macro-particle and field phase space.^6,14

VII. NUMERICAL EXAMPLES

To illustrate our formulation, we consider a combination of an electron electromagnetic mode and a Langmuir wave in one spatial dimension with stationary ions in the nonrelativistic approximation. We solve taking

μ \equiv 0

to understand the errors in Gauss's Law due to the lack of an exact continuity equation. We take the domain size

k_{p} L = 2 π

⁠, where

k_{p} c = ω_{p}

and

ω_{p} = \sqrt{4 π q^{2} n_{0} / m}

is the plasma frequency and

n_{0}

is the equilibrium density. We examine three choices for

ϱ_{i}

⁠, which correspond to standard particle shapes¹⁷ used with the PIC algorithm, that are piecewise quadratic, cubic, and quartic, extending over three, four, and five cells, respectively. The non-zero values of

ϱ_{k}

are given by

\begin{matrix} ϱ_{k - 1} = \frac{1}{2} {(δ - \frac{1}{2})}^{2}, \\ ϱ_{k} = \frac{3}{4} - δ^{2}, \\ ϱ_{k + 1} = \frac{1}{2} {(δ + \frac{1}{2})}^{2}, \end{matrix}

(126)

quadratic;

\begin{matrix} \begin{matrix} ϱ_{k - 1} = \frac{1}{6} - \frac{1}{2} δ + \frac{1}{2} δ^{2} - \frac{1}{6} δ^{3}, \\ ϱ_{k} = \frac{2}{3} - δ^{2} + \frac{1}{2} δ^{3}, \\ ϱ_{k + 1} = \frac{1}{6} + \frac{1}{2} δ + \frac{1}{2} δ^{2} - \frac{1}{2} δ^{3}, \\ ϱ_{k + 2} = \frac{1}{6} δ^{3}, \end{matrix}} δ > 0, \\ \begin{matrix} ϱ_{k - 2} = - \frac{1}{6} δ^{3}, \\ ϱ_{k - 1} = \frac{1}{6} - \frac{1}{2} δ + \frac{1}{2} δ^{2} + \frac{1}{2} δ^{3}, \\ ϱ_{k} = \frac{2}{3} - δ^{2} - \frac{1}{2} δ^{3}, \\ ϱ_{k + 1} = \frac{1}{6} + \frac{1}{2} δ + \frac{1}{2} δ^{2} + \frac{1}{6} δ^{3} \end{matrix}} δ < 0, \end{matrix}

(127)

cubic; and

\begin{matrix} ϱ_{k - 2} = \frac{1}{24} {(δ - \frac{1}{2})}^{4}, \\ ϱ_{k - 1} = \frac{19}{96} - \frac{11}{24} δ + \frac{1}{4} δ^{2} + \frac{1}{6} δ^{3} - \frac{1}{6} δ^{4}, \\ ϱ_{k} = \frac{115}{192} - \frac{5}{8} δ^{2} + \frac{1}{4} δ^{4}, \\ ϱ_{k + 1} = \frac{19}{96} + \frac{11}{24} δ + \frac{1}{4} δ^{2} - \frac{1}{6} δ^{3} - \frac{1}{6} δ^{4}, \\ ϱ_{k + 2} = \frac{1}{24} {(δ + \frac{1}{2})}^{4}, \end{matrix}

(128)

quartic, where

δ = (ξ_{z}^{α} - z_{k}) / Δ z

and

z_{k}

is the grid point closest to the macro-particle location

ξ_{z}^{α}

⁠.

The macro-particles are initially uniformly spaced as

ξ_{z}^{α} = \frac{Δ z}{N_{ppc}} (j + \frac{1}{2}), 0 \leq j < N_{ppc} N_{z},

(129)

where

N_{ppc}

is the number of macro-particles loaded in each cell. Initially,

{\dot{ξ}}_{x}^{α} = 0

and

{\dot{ξ}}_{z}^{α} = 0.01 cos (2 k_{p} ξ_{z}^{α}) c .

(130)

Thus, the initial charge density is zero and hence

{\dot{A}}_{z, i} = 0

exactly satisfies Gauss's Law. The initial transverse vector potential is taken to be

A_{x, i} = 0.1 \frac{m c^{2}}{q} cos (k_{p} z_{i}),

(131)

with

{\dot{A}}_{x, i} = 0

⁠. We solve the equations of motion using the method of lines²⁷ with a second-order central differencing in space and fourth-order Runge–Kutta time stepping which has the stability constraint

c Δ t / Δ z \leq \sqrt{2}

⁠.²⁸ We take

N_{z} \in {32 \times 2^{j} | j \in ℕ, j \leq 8}

and

N_{ppc} \in {1, 2, 4, 8, 16, 32}

⁠. For

N_{z} = 32

⁠, we set

ω_{p} Δ t = 0.1

(⁠

c Δ t / Δ z \approx 0.509

⁠) while for

N_{z} \geq 64

⁠, we set

c Δ t = (16 / 5 π) Δ z \approx 1.019 Δ z

⁠.

Figure 1 shows $E_{x} / E_{0}$ (top row), $E_{z} / E_{0}$ (second row), $N_{e} / n_{o} = ρ / (q n_{0})$ (third row), and the error in Gauss's Law ${(D^{x} E_{z})}_{i} - k_{p} E_{0} {(N_{e} - n_{}^{ion})}_{i}$ (bottom row), where $E_{0} = m c ω_{p} / q$ ⁠. We compare $N_{z} = 32$ ⁠, $N_{ppc} = 1$ ⁠, and quadratic $ϱ_{k}$ (first and third column) with $N_{z} = 1024$ ⁠, $N_{ppc} = 32$ ⁠, and quartic $ϱ_{k}$ (second and fourth column). While $N_{z} = 32$ is insufficient resolution to capture the dynamics, the qualitative features are plausibly correct. The differences in the two cases is largely due to spatial resolution; taking $N_{z} = 32$ ⁠, $N_{ppc} = 32$ ⁠, quartic $ϱ_{k}$ and the same time step as the high-resolution gives essentially the same results as the $N_{z} = 32$ case in Fig. 1.

FIG. 1.

View large Download slide

Electric fields, density, and the error in Gauss's Law for $N_{z} = 32$ ⁠, $N_{ppc} = 1$ ⁠, and quadratic $ϱ_{K}$ (126), (a)–(d) (early time) and (i)–(l) (late time) and for $N_{z} = 1024$ ⁠, $N_{ppc} = 32$ and quartic $ϱ_{k}$ (128), (e)–(h) (early time) and (m)–(p) (late time): $E_{x} / E_{0}$ (a), (e), (i), and (m); $E_{z} / E_{0}$ (b), (f), (j), and (n); $N_{e} / n_{0}$ (c), (g), (k), and (o); and ${(D^{x} E_{z})}_{i} - k_{p} E_{0} {(N_{e} - n^{ion})}_{i}$ (d), (h), (l), and (p).

Figure 2 shows the RMS error (in z) in Gauss's Law for quadratic (top row) and quartic (bottom row) $ϱ_{k}$ with $N_{z} = 32$ (left column), $N_{z} = 512$ (center column), and $N_{z} = 8192$ (right column) for various values of $N_{ppc}$ ⁠. For quadratic $ϱ_{k}$ ⁠, the error in Gauss's Law exhibits significant dependence on $N_{ppc}$ as the resolution increases; for the range of $N_{ppc}$ ⁠, the error plateaus for $N_{z} \geq 512$ (see Fig. 4). Quartic $ρ_{k}$ yields errors that are essentially independent of $N_{ppc}$ for all but the highest resolution considered. The error in Gauss's Law arises from the difference between $\partial ϱ_{k} / \partial ξ_{x}^{α}$ and ${(D^{x} ϱ)}_{k}$ ⁠, thus, we expect the error to be smaller for $ϱ_{k}$ that span more cells and to scale with $Δ z^{2}$ due to second-order approximations of spatial derivatives. Both of these expectations are borne out in Fig. 2.

FIG. 2.

View large Download slide

The RMS error in Gauss's Law for quadratic $ϱ_{k}$ (126) (top row) and quartic $ϱ_{k}$ (128) (bottom row) and a range of $N_{ppc}$ ⁠: $N_{z} = 32$ (a) and (d); $N_{z} = 512$ (b) and (e); and $N_{z} = 8192$ (c) and (f).

Figure 3 shows the absolute value $P_{z} / m c$ for the same parameters as in Fig. 2. The departure of $P_{z}$ from its initial value of 0 is the consequence of the lack of translation invariance of the discretization. The error in $P_{z}$ is seen to be largely insensitive to both $N_{ppc}$ and to the nature of $ϱ_{k}$ ⁠. For quadratic $ϱ_{k}$ ⁠, the behavior for $N_{ppc} = 1$ and $N_{ppc} = 2$ is inconsistent with the other cases. We do not have an explanation for this observation.

FIG. 3.

View large Download slide

The absolute value of $P_{z} / m c$ for quadratic $ϱ_{k}$ (126) (top row) and quartic, $ϱ_{k}$ (128), (bottom row) and a range of $N_{ppc}$ ⁠: $N_{z} = 32$ (a) and (d); $N_{z} = 512$ (b) and (e); and $N_{z} = 8192$ (c) and (f).

Figures 4 and 5 show the RMS (in z and t) error in Gauss's Law and the RMS value of $P_{z} / m c$ ⁠, respectively, for quadratic (top), cubic (center), and quartic (bottom) $ϱ_{k}$ ⁠. The dashed gray line is proportional to $Δ z^{2}$ from which we see that both quantities scale consistent with second-order spatial differencing. At fixed $N_{ppc}$ ⁠, the error in Gauss's Law will become independent of $Δ z$ for small $k_{p} Δ z$ ⁠. As the macro-particle becomes broader and smoother, the quadratic scaling persists to smaller $k_{p} Δ z$ ⁠. We have verified that $Δ t$ does not influence this behavior: we repeated the quadratic $ϱ_{k}$ cases with the smallest value of $Δ t$ used in our computations and we obtain virtually identical results to those of Fig. 4(a).

FIG. 4.

View large Download slide

The RMS error in Gauss's Law for quadratic $ϱ_{k}$ (126), (a); cubic $ϱ_{k}$ (127), (b); and quartic $ϱ_{k}$ (128), (c).

FIG. 5.

View large Download slide

The RMS value of $P_{z} / m c$ for quadratic $ϱ_{k}$ (126), (a); cubic $ϱ_{k}$ (127), (b); and quartic $ϱ_{k}$ (128), (c).

Again we see that the behavior of $P_{z}$ is nearly independent of the details of $ϱ_{k}$ and $N_{ppc}$ ⁠; $P_{z}$ exhibits nearly perfect scaling with $Δ z$ ⁠. As with Fig. 3, for quadratic $ϱ_{k}$ ⁠, the behavior for $N_{ppc} = 1$ and $N_{ppc} = 2$ is inconsistent with the other cases. There is an indication of this behavior in the cubic $ϱ_{k}$ case as the $N_{ppc} = 1$ curve appears to be slightly above the other lines for large $k_{p} Δ z$ ⁠. Understanding the interplay of spatial resolution and macro-particle shape will be the subject of future work.

The errors in Gauss's Law and momentum generally do not grow in time and scale as expected with spatial resolution. The error in Gauss's Law can be reduced to near machine precision by using the Lagrange multiplier however doing so involves a significant increase in computational cost as an elliptic equation, (115), would have to be solved at each stage of the time stepper. We do not expect imposing the Gauss's Law constraint to effect the error in momentum. Exploring the advantages and trade-off of imposing Gauss's Law as a constraint will be explored in a future publication.

VIII. CONCLUSIONS

We have shown that a macro-particle reduction applied to the Low Lagrangian or to the noncanonical Vlasov–Maxwell system results in equivalent dynamical systems with continuous space. In the case of the noncanonical system, a complicated coordinate transformation yielded a canonical Hamiltonian and bracket. Using the Dirac–Bergmann procedure, we constructed a canonical Hamiltonian system on the extended phase-space from the macro-particle reduced Lagrangian, which differs from the transform noncanonical system as the latter yields the constrained Hamiltonian. We considered two approaches to reducing these systems to a finite degree-of-freedom approximations suitable for computational implementation. We illustrated a general truncated basis function expansion approach that preserves gauge invariance (and, thus, Gauss's Law) and showed that the special case of a truncated Fourier representation of the potentials also conserves total momentum.

The use of a Fourier basis seems unavoidable if exact momentum conservation is desired. The basis must be an eigenfunction of the translation operator which seems to suggest that a Fourier representation is required.^5,6 Provided the domain is large enough, the influence of periodic boundary conditions can be minimized. Of more concern is the all-to-all communications necessary a parallel implementation since all of the macro-particle couple to all field modes (in contrast to the localized coupling in a gridded mode). As a result, there currently appears to be no path to obtaining sufficient computational performance to allow use of a Fourier representation for production work. Nonetheless, the Fourier representation is valuable for assessing the quantitative errors introduced by the lack of exact momentum conservation.

We presented a grid-based discretization of both the Lagrangian and canonical Hamiltonian formulations in the Weyl gauge. The spatial grid results in the continuity equation only holding on average and, thus, Gauss's Law must be enforced with a constraint. Of course, all of these reductions have exact energy conservation with continuous time. We have presented numerical examples in one spatial dimension of a combined electromagnetic mode and Langmuir wave. The error in Gauss's Law does not grow secularly and scales as expected with spatial resolution. The equation of motion in canonical Hamiltonian formulation can be numerically integrated using standard symplectic methods which should result in excellent (but not exact) energy conservation.

It has been shown that variational methods can have low level of noise.^6,7 Here we have shown that the Weyl gauge, without imposing Gauss's Law with a Lagrange multiplier results in reasonable conservation properties. Thus, it is possible to use the variational formulation without having to solve an elliptic equation, thereby having performance approximately equivalent (in terms of arithmetic operations per time step) to traditional PIC algorithms but with much less noise so fewer particles should be necessary. As our examples show, the departure from exact conservation is rather benign: errors do not grow secularly and scales quadratically with grid resolution. It seems likely that with higher-order spatial approximations, these errors could become negligible.

ACKNOWLEDGMENTS

The authors gratefully acknowledge helpful discussions with P. J. Morrison, John Finn, and Evstati Evstatiev. This work was supported by NSF (Grant No. PHY-2108788). A.J.H. was supported by the University of Nebraska Program of Excellence in Atomic, Molecular, Optical, and Plasma Physics.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Adam J. Higuet: Conceptualization (supporting); Formal analysis (supporting); Writing – original draft (lead); Writing – review & editing (lead). B. A. Shadwick: Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Supervision (lead); Writing – original draft (supporting); Writing – review & editing (supporting).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX: GENERAL GAUGE INVARIANT BASIS EXPANSION

Inspired by Kraus et al.,⁹ which seeks to discretize the noncanonical formulation in a way that preserves the entire de Rahm complex present in vector calculus, we seek to construct a simpler restriction on basis elements while retaining automatic gauge invariance. The idea we keep from the exploration into discrete exterior calculus is that 0-forms and 1-forms (scalars and polar vectors) are geometrically different objects, but related to one another via a gradient operator. We, thus, choose to expand the scalar potential with some set of basis functions and the vector potential with a new basis made from derivatives of the scalar basis

φ (r, t) = \sum_{I} φ_{I} (t) u_{I} (r)

(A1a)

and

A (r, t) = \sum_{I} \sum_{a} A_{I}^{a} (t) u_{I}^{a} (r),

(A1b)

where

u_{I} (r) = X_{i} (x) Y_{j} (y) Z_{k} (z),

(A2a)

u_{I}^{a} (r) = \sum_{a} e_{a} \partial_{a} u_{I} (r),

(A2b)

capital roman letters are used as a shorthand for three lowercase indices

{i, j, k}

denoting individual Cartesian directions' basis indices, lowercase roman letters at the start of the alphabet are used for the Cartesian directions

{x, y, z}

⁠,

e_{a}

are Cartesian unit vectors, and the summation convention is not used. For example,

u_{I}^{x} (r) = u_{i, j, k}^{x} (x, y, z) = e_{x} \frac{d X_{i} (x)}{d x} Y_{j} (y) Z_{k} (z) .

(A3)

The gradient operator can be written as

\nabla = \sum_{a} e_{a} \partial_{a}

so we have that

\nabla u_{I} = \sum_{a} u_{I}^{a}

⁠, which is the only carryover from discrete exterior calculus that is necessary to maintain gauge invariance. Using this expansion, we find

\begin{matrix} \int d^{3} r S (r - ξ^{α}) {\dot{ξ}}^{α} \cdot A \\ = \sum_{I} \sum_{a} A_{I}^{a} \int d^{3} r S (r - ξ^{α}) {\dot{ξ}}^{α} \cdot u_{I}^{a} (r) \\ = \sum_{I} \sum_{a} A_{I}^{a} \int d^{3} r S (r - ξ^{α}) {\dot{ξ}}^{α} \cdot e_{a} \partial_{a} u_{I} (r) \\ = \sum_{I} \sum_{a} A_{I}^{a} \int d^{3} r S (r - ξ^{α}) {\dot{ξ}}_{a}^{α} \partial_{a} u_{I} (r) \\ = \sum_{I} \sum_{a} A_{I}^{a} {\dot{ξ}}_{a}^{α} \frac{\partial ϱ_{I} (ξ^{α})}{\partial ξ_{a}^{α}}, \end{matrix}

(A4)

where

ϱ_{I} = \int d^{3} r S (r - ξ^{α}) u_{I} (r)

(A5)

is the projection of the particle shape onto the basis. A similar calculation gives

\int d^{3} r S (r - ξ^{α}) φ = \sum_{I} φ_{I} ϱ_{I} (ξ^{α}) .

(A6)

The discretized Lagrangian becomes

\begin{matrix} L = - m c^{2} \sum_{α} w_{α} \sqrt{1 - {| \frac{{\dot{ξ}}^{α}}{c} |}^{2}} \\ + \frac{1}{8 π} \sum_{I, J} \sum_{a} [\frac{1}{c^{2}} {\dot{A}}_{I}^{a} {\dot{A}}_{J}^{a} + \frac{2}{c} φ_{I} {\dot{A}}_{J}^{a} + φ_{I} φ_{J}] M_{I J}^{a} \\ - \frac{1}{8 π} \sum_{I, J} \sum_{a, b} A_{I}^{a} A_{J}^{b} M_{I J}^{a b} - \sum_{I} φ_{I} \int d^{3} r u_{I} (r) ρ^{ion} (r) \\ + q \sum_{α} w_{α} \sum_{I} [\frac{1}{c} \sum_{a} A_{I}^{a} {\dot{ξ}}_{a}^{α} \frac{\partial ϱ_{I} (ξ^{α})}{\partial ξ_{a}^{α}} - φ_{I} ϱ_{I} (ξ^{α})], \end{matrix}

(A7)

where the so-called mass matrices are defined by

M_{I J}^{a} = \int d^{3} r \partial_{a} u_{I} (r) \partial_{a} u_{J} (r)

(A8)

and

\begin{matrix} M_{I J}^{a b} = \int d^{3} r u_{I}^{a} (r) \times u_{J}^{b} (r) \\ = \sum_{c} e_{c} ε_{cab} \int d^{3} r \partial_{a} u_{I} (r) \partial_{b} u_{J} (r) . \end{matrix}

(A9)

The mass matrix

M_{I J}^{a}

is generally expected to be sparse but have large bandwidth; its inverse (which will be dense) appears in the evolution equations for

φ_{I}

and

A_{I}^{a}

⁠. Except for very small numbers of basis functions, iterative solution of the corresponding linear systems will be necessary. Furthermore, the force terms in the macro-particle evolution equations couple to all

φ_{I}

and

A_{I}^{a}

⁠, thus, necessitating all-to-all communication in a parallel implementation.

Consider a general gauge transformation,

φ \to φ - \frac{1}{c} \frac{\partial λ}{\partial t}

(A10a)

and

A \to A + \nabla λ .

(A10b)

We expand the gauge function as

λ (r, t) = \sum_{I} λ_{I} (t) u_{I} (r),

(A11)

which gives us

\frac{\partial λ}{\partial t} = \sum_{I} {\dot{λ}}_{I} u_{I} (r)

(A12)

and

\nabla λ = \sum_{I} λ_{I} (t) \nabla u_{I} (r) = \sum_{I} \sum_{a} λ_{I} (t) u_{I}^{a} (r) .

(A13)

This implies the potentials transform mode-by-mode as

\sum_{I} φ_{I} (t) u_{I} (r) \to \sum_{I} [φ_{I} - \frac{1}{c} {\dot{λ}}_{I}] u_{I} (r)

(A14)

and

\sum_{I} \sum_{a} A_{I}^{a} (t) u_{I}^{a} (r) \to \sum_{I} \sum_{a} [A_{I}^{a} (t) + λ_{I} (t)] u_{I}^{a} (r) .

(A15)

It is easy show that the field terms in (A7) are exactly invariant under this transformation. The interaction terms in the Lagrangian are found to transform as

\begin{matrix} L_{Int} \to L_{Int} + \frac{q}{c} \sum_{α} w_{α} \sum_{I} [{\dot{λ}}_{I} ϱ_{I} (ξ^{α}) \\ + λ_{I} {\dot{ξ}}_{a}^{α} \frac{\partial ϱ_{I} (ξ^{α})}{\partial ξ_{a}^{α}}] + \frac{1}{c} \sum_{I} {\dot{λ}}_{I} \int d^{3} r u_{I} (r) ρ^{ion} (r) \\ = L_{Int} + \frac{d}{d t} [\frac{q}{c} \sum_{α} w_{α} \sum_{I} λ_{I} ϱ_{I} (ξ^{α}) \\ + \frac{1}{c} \sum_{I} λ_{I} \int d^{3} r u_{I} (r) ρ^{ion} (r)], \end{matrix}

(A16)

so we find that the transformation does not effect the dynamics and the structure of the basis expansion preserves gauge symmetry. While this expansion preserves gauge invariance, the basis functions, in general, are not translation invariant (i.e., they are not eigenfunctions of the translation operator). By the same analysis in the previous work,^5,6 we determine that translational invariance and momentum conservation are only maintained by a discrete Fourier basis, which can be viewed as a special case of this general formulation. By virtue of the analysis of Sec. III, this same expansion can be used directly in the noncanonical system and will yield a good Poisson bracket.

REFERENCES

1.

H. R.

Lewis

, “

Energy-conserving numerical approximations for Vlasov plasmas

,”

J. Comput. Phys.

6

,

136

–

141

(

1970

).

https://doi.org/10.1016/0021-9991(70)90012-4

Google Scholar

Crossref

2.

J. W.

Eastwood

, “

The virtual particle electromagnetic particle-mesh method

,”

Comput. Phys. Commun.

64

,

252

–

266

(

1991

).

https://doi.org/10.1016/0010-4655(91)90036-K

Google Scholar

Crossref

3.

F. E.

Low

, “

A Lagrangian formulation of the Boltzmann–Vlasov equation for plasmas

,”

Proc. R. Soc. London, Ser. A

248

,

282

–

287

(

1958

).

https://doi.org/10.1098/rspa.1958.0244

Google Scholar

Crossref

4.

J.

Squire

,

H.

Qin

, and

W. M.

Tang

, “

Geometric integration of the Vlasov–Maxwell system with a variational particle-in-cell scheme

,”

Phys. Plasmas

19

,

084501

(

2012

).

https://doi.org/10.1063/1.4742985

Google Scholar

Crossref

5.

E. G.

Evstatiev

and

B. A.

Shadwick

, “

Variational formulation of particle algorithms for kinetic plasma simulations

,”

J. Comput. Phys.

245

,

376

–

398

(

2013

).

https://doi.org/10.1016/j.jcp.2013.03.006

Google Scholar

Crossref

6.

B. A.

Shadwick

,

A. B.

Stamm

, and

E. G.

Evstatiev

, “

Variational formulation of macro-particle plasma simulation algorithms

,”

Phys. Plasmas

21

,

055708

(

2014

).

https://doi.org/10.1063/1.4874338

Google Scholar

Crossref

7.

A. B.

Stamm

,

B. A.

Shadwick

, and

E. G.

Evstatiev

, “

Variational formulation of macroparticle models for electromagnetic plasma simulations

,”

IEEE Trans. Plasma Sci.

42

,

1747

–

1758

(

2014

).

https://doi.org/10.1109/TPS.2014.2320461

Google Scholar

Crossref

8.

Y.

He

,

Y.

Sun

,

H.

Qin

, and

J.

Liu

, “

Hamiltonian particle-in-cell methods for Vlasov–Maxwell equations

,”

Phys. Plasmas

23

,

092108

(

2016

).

https://doi.org/10.1063/1.4962573

Google Scholar

Crossref

9.

M.

Kraus

,

K.

Kormann

,

P. J.

Morrison

, and

E.

Sonnendrücker

, “

GEMPIC: geometric electromagnetic particle-in-cell methods

,”

J. Plasma Phys.

83

,

905830401

(

2017

).

https://doi.org/10.1017/S002237781700040X

Google Scholar

Crossref

10.

A. B.

Langdon

, “

Effects of the spatial grid in simulation plasmas

,”

J. Comput. Phys.

6

,

247

–

267

(

1970

).

https://doi.org/10.1016/0021-9991(70)90024-0

Google Scholar

Crossref

11.

H.

Okuda

, “

Nonphysical noises and instabilities in plasma simulation due to a spatial grid

,”

J. Comput. Phys.

10

,

475

–

486

(

1972

).

https://doi.org/10.1016/0021-9991(72)90048-4

Google Scholar

Crossref

12.

E.

Cormier-Michel

,

B. A.

Shadwick

,

C. G. R.

Geddes

,

E.

Esarey

,

C. B.

Schroeder

, and

W. P.

Leemans

, “

Unphysical kinetic effects in particle-in-cell modeling of laser wakefield accelerators

,”

Phys. Rev. E

78

,

016404

(

2008

).

https://doi.org/10.1103/PhysRevE.78.016404

Google Scholar

Crossref

13.

P. J.

Channell

and

C.

Scovel

, “

Symplectic integration of Hamiltonian systems

,”

Nonlinearity

3

,

231

(

1990

).

https://doi.org/10.1088/0951-7715/3/2/001

Google Scholar

Crossref

14.

A. B.

Stamm

and

B. A.

Shadwick

, “

Comparison of conservation behavior in several variationally derived implementations for simulating electromagnetic plasmas

,” in

2018 IEEE Advanced Accelerator Concepts Workshop (AAC)

, edited by

E. I.

Simakov

,

N.

Yampolsky

, and

K. P.

Wootton

(

IEEE

,

2018

), pp.

1

–

5

.

Google Scholar

15.

Y. L.

Klimontovich

,

The Statistical Theory of Non-Equilibrium Processes in a Plasma

(

MIT Press

,

Cambridge, MA

,

1967

).

Google Scholar

16.

C. K.

Birdsall

and

A. B.

Langdon

,

Plasma Physics via Computer Simulations

, Series in Plasma Physics No. IP586 (

Institute of Physics Publishing

,

Bristol

,

1991

).

Google Scholar

17.

R. W.

Hockney

and

J. W.

Eastwood

,

Computer Simulation Using Particles

(

Taylor & Francis Group

,

New York

,

1988

).

Google Scholar

18.

P. J.

Morrison

, “

The Maxwell–Vlasov equations as a continuous Hamiltonian system

,”

Phys. Lett. A

80

,

383

–

386

(

1980

).

https://doi.org/10.1016/0375-9601(80)90776-8

Google Scholar

Crossref

19.

A.

Weinstein

and

P. J.

Morrison

, “

Comments on: The Maxwell–Vlasov equations as a continuous hamiltonian system

,”

Phys. Lett. A

86

,

235

–

236

(

1981

).

https://doi.org/10.1016/0375-9601(81)90496-5

Google Scholar

Crossref

20.

P. J.

Morrison

, “

Poisson brackets for fluids and plasmas

,”

AIP Conf. Proc.

88

,

13–46

(

1982

).

https://doi.org/10.1063/1.33633

Google Scholar

Crossref

21.

J. E.

Marsden

and

A.

Weinstein

, “

The Hamiltonian structure of the Maxwell–Vlasov equations

,”

Phys. D

4

,

394

(

1982

).

https://doi.org/10.1016/0167-2789(82)90043-4

Google Scholar

Crossref

22.

K.

Sundermeyer

, Constrained Dynamics: With Applications to Yang–Mills Theory, General Relativity, Classical Spin, Dual String Model, Lecture Notes in Physics, Vol. 169 (

Springer-Verlag

,

Berlin

,

1982

).

23.

P. A. M.

Dirac

, “

Generalized Hamiltonian dynamics

,”

Can. J. Math.

2

,

129

–

148

(

1950

).

https://doi.org/10.4153/CJM-1950-012-1

Google Scholar

Crossref

24.

J. L.

Anderson

and

P. G.

Bergmann

, “

Constraints in covariant field theories

,”

Phys. Rev.

83

,

1018

–

1025

(

1951

).

https://doi.org/10.1103/PhysRev.83.1018

Google Scholar

Crossref

25.

P. G.

Bergmann

and

I.

Goldberg

, “

Dirac bracket transformations in phase space

,”

Phys. Rev.

98

,

531

–

538

(

1955

).

https://doi.org/10.1103/PhysRev.98.531

Google Scholar

Crossref

26.

A. B.

Stamm

, “

Variational formulation of Macro-Particle algorithms for studying electromagnetic plasmas

,” Ph.D. thesis (

University of Nebraska

,

Lincoln

,

2015

).

Google Scholar

27.

W. E.

Schiesser

,

The Numerical Method of Lines: Integration of Partial Differential Equations

(

Academic Press

,

San Diego

,

1991

).

Google Scholar

28.

J.

Reyes

and

B.

Shadwick

, “

Higher-order explicit methods for laser-plasma interactions

,” in

2013 19th IEEE Pulsed Power Conference

(

IEEE

,

2013

), pp.

1

–

5

.

Google Scholar

© 2025 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International (CC BY-NC-ND) license (https://creativecommons.org/licenses/by-nc-nd/4.0/).

This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License .

Gauge invariance and conservation laws in the variational formulation of macro-particle plasma models

I. INTRODUCTION

II. LAGRANGIAN MACRO-PARTICLE REDUCTION

III. REDUCTION OF THE NONCANONICAL SYSTEM

A. Macro-particle reduction

B. Transformation to canonical coordinates

C. Connection to the Lagrangian formulation

D. Constructing the primary Hamiltonian

IV. GAUGE INVARIANCE DISCRETIZATION: TRUNCATED FOURIER BASIS

V. CHARGE CONSERVATION AND GAUGE SELECTION IN NON-GAUGE INVARIANT SYSTEMS

A. Varying a non-gauge invariant Lagrangian

B. Enforcing Gauss's Law via Lagrange multipliers

VI. DISCRETIZING ON A GRID

A. Lagrangian formulation

B. Canonical formulation

VII. NUMERICAL EXAMPLES

VIII. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX: GENERAL GAUGE INVARIANT BASIS EXPANSION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Gauge invariance and conservation laws in the variational formulation of macro-particle plasma models

I. INTRODUCTION

II. LAGRANGIAN MACRO-PARTICLE REDUCTION

III. REDUCTION OF THE NONCANONICAL SYSTEM

A. Macro-particle reduction

B. Transformation to canonical coordinates

C. Connection to the Lagrangian formulation

D. Constructing the primary Hamiltonian

IV. GAUGE INVARIANCE DISCRETIZATION: TRUNCATED FOURIER BASIS

V. CHARGE CONSERVATION AND GAUGE SELECTION IN NON-GAUGE INVARIANT SYSTEMS

A. Varying a non-gauge invariant Lagrangian

B. Enforcing Gauss's Law via Lagrange multipliers

VI. DISCRETIZING ON A GRID

A. Lagrangian formulation

B. Canonical formulation

VII. NUMERICAL EXAMPLES

VIII. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX: GENERAL GAUGE INVARIANT BASIS EXPANSION

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only