A self-consistent Hamiltonian model of the ponderomotive force and its structure preserving discretization

Barham, William; Güçlü, Yaman; Morrison, Philip J.; Sonnendrücker, Eric

doi:10.1063/5.0178935

In the presence of an inhomogeneous oscillatory electric field, charged particles experience a net force, averaged over the oscillatory timescale, known as the ponderomotive force. We derive a one-dimensional Hamiltonian model which self-consistently couples the electromagnetic field to a plasma which experiences the ponderomotive force. We derive a family of structure preserving discretizations of the model of varying order in space and time using conforming and broken finite element exterior calculus spectral element methods. In all variants of our discretization framework, the method is found to conserve the Casimir invariants of the continuous model to machine precision and the energy to the order of the splitting method used.

I. INTRODUCTION

A well-known laser–plasma interaction model gives rise to a nonlinear polarization of the electromagnetic medium. In the presence of an inhomogeneous oscillatory electric field, charged particles experience a net force, averaged over the oscillatory timescale, known as the ponderomotive force. Starting from the two-fluid equations, we derive a Hamiltonian model in which the ponderomotive force provides a current source for Maxwell's equations. The asymptotic assumptions under which this model is derived correspond to what is frequently called the “linear regime” in laser plasma interaction models.¹ Our model is possibly the simplest self-consistent model of the ponderomotive force and closely resembles other fluid models of the ponderomotive force in the linear regime,² although the model therein is stated in terms of the vector potential whereas our model is stated in terms of the electromagnetic fields. Our model is perhaps the first to elucidate the Hamiltonian structure of the 1D linear regime and provides an alternative approach to its derivation based on applying the asymptotics directly to the energy functional and Poisson bracket rather than the equations of motion. This self-consistent model of the ponderomotive force leads to a nonlinear polarization of the electromagnetic medium; we adapted ideas from a modeling framework for kinetic theories in nonlinear media to build the model studied herein.³ We call our model the 1D ponderomotive Maxwell system.

The simple Poisson bracket one finds for the 1D ponderomotive Maxwell system makes this model straightforward to discretize in a structure preserving manner: its Poisson bracket is a direct sum of the Poisson brackets for Maxwell and acoustic wave equations, as well as a coupling bracket. Each of these brackets is field-free, i.e., the Poisson bivector is independent of the dynamical fields. Hence, antisymmetry is a sufficient condition for the Jacobi identity to be satisfied,⁴ and a structure preserving discretization is easily accomplished using finite element exterior calculus (FEEC).⁵ In particular, as the model studied in this paper is a nonlinear wave equation, we study the model using a variety of FEEC which is built from broken elements (i.e., the method uses discontinuous shape functions which are local to each element).^6–10 These are referred to variously as CONGA (COnforming/Non-conforming GAlerkin) or broken-FEEC methods. We consider both standard and broken-FEEC methods in our numerical studies. In general, self-consistent models of the interaction of a plasma with the electromagnetic field come from field-dependent Poisson brackets making structure preserving discretizations rare. For example, Hamiltonian structure preserving discretizations of the Vlasov–Maxwell system and 2D vorticity dynamics rely on particle based discretizations: Hamiltonian structure preserving discretizations of the Vlasov–Maxwell system have been derived based on a particle based representation of the distribution function,^11,12 and the reduction of 2D vorticity to finite dimensional Hamiltonian mechanics of point vorticies is classical.¹³ The nonlinearities in the ponderomotive Maxwell model are all contained in the Hamiltonian leaving the Poisson bracket extremely simple. The amenability of this model to structure preserving discretization points to a broad class of nonlinear-wave models likewise admitting a similar numerical treatment.

This paper is organized as follows: In Sec. II, the 1D ponderomotive Maxwell model is derived from the two-fluid Maxwell–Euler system via a traditional order-by-order asymptotic approach. In Sec. III, the Hamiltonian structure of the 1D ponderomotive Maxwell system is derived by applying the same asymptotic considerations as before to the Hamiltonian structure of the two-fluid Maxwell–Euler system. In Sec. IV, a Hamiltonian structure preserving discretization of the 1D ponderomotive Maxwell system is derived. Finally, in Sec. V, the numerical results from simulating the 1D ponderomotive Maxwell system using the algorithm derived in Sec. IV are presented. Appendix A elucidates the Hamiltonian structure of Maxwell's equations in nonlinear media. Appendix B provides the details of the 1D FEEC spectral element method used in Sec. IV. Appendix C provides a technical result used in the discretization procedure.

II. DIRECT DERIVATION FROM THE TWO-FLUID MAXWELL–EULER MODEL

In the presence of a inhomogeneous oscillatory electric field, charged particles experience a net force, averaged over the oscillatory timescale known as the ponderomotive force given by

F_{p} = - \frac{q^{2}}{4 m ω_{0}^{2}} \nabla | E |^{2},

(1)

where m and q are the particles' mass and charge, respectively, and ω₀ and E are the frequency and amplitude of the oscillatory field, respectively. In an appropriate asymptotic regime, we shall find that the envelope of the oscillatory electromagnetic fields evolve according to the macroscopic Maxwell's equations

\partial_{t} D = c \nabla \times B and \partial_{t} B = - c \nabla \times E,

(2)

where

D = E + P

with some to be determined self-consistent polarization. In this section, we elucidate the asymptotic derivation of such a model which we call the ponderomotive Maxwell equation beginning from a two-fluid model coupled to Maxwell's equations and then provide a self-consistent Hamiltonian formulation of the same model.

The parent model from which we derive our asymptotic model, which we call the ponderomotive Maxwell system, is the two-fluid model. In Gaussian units, this model may be written as

\begin{matrix} \partial_{t} E - c \nabla \times B = - 4 π (q_{e} n_{e} v_{e} + q_{i} n_{i} v_{i}), \\ \partial_{t} B + c \nabla \times E = 0, \\ \partial_{t} v_{α} + \frac{1}{2} \nabla | v_{α} |^{2} - v_{α} \times \nabla \times v_{α} = \frac{q_{α}}{m_{α}} (E + \frac{1}{c} v_{α} \times B) + \frac{\nabla p_{α}}{n_{α}}, \\ \partial_{t} n_{α} + \nabla \cdot (n_{α} v_{α}) = 0, \end{matrix}

(3)

where

α = e, i

⁠, and the thermodynamic pressure,

p_{α} = n_{α}^{2} {(U_{α})}_{n_{α}}

⁠, is prescribed by an internal energy,

U_{α} (n_{α})

⁠, which depends only on density. As a matter of notation convention,

{(U_{α})}_{n_{α}} = \partial U_{α} / \partial n_{α}

⁠. Using the chain rule, one may show that

\frac{\nabla p_{α}}{n_{α}} = \nabla ({(n_{α} U_{α})}_{n_{α}}) .

(4)

We assume

q_{i} = - Z q_{e}

where q_e is the charge of an electron. Such an isothermal two-fluid model describes a plasma which has reached local thermodynamic equilibrium and which possesses a single temperature. We shall assume that all processes occur on a timescale faster than that at which the temperature changes so that the assumptions of this model are not violated. A more complete model would include an entropy advection law and an entropy dependent internal energy. While we retain a thermodynamic pressure for the purpose of deriving the 1D ponderomotive Maxwell system in this section And Sec. III, we will neglect it for simplicity in our discrete model in Secs. IV and V.

A. The two-fluid system: Longitudinal/transverse split

The problem studied in this paper possesses a distinguished direction, namely, the direction of laser propagation, which we call

\hat{z}

⁠. The plane transverse to the direction of laser propagation is denoted with the symbol “

⊥

⁠.” The system written in this longitudinal/transverse split takes the following form:

\begin{matrix} \partial_{t} E_{⊥} - c \hat{z} \times (\partial_{z} B_{⊥} - \nabla_{⊥} B_{z}) = - 4 π \sum_{α} q_{α} n_{α} v_{α, ⊥}, \\ \partial_{t} B_{⊥} + c \hat{z} \times (\partial_{z} E_{⊥} - \nabla_{⊥} E_{z}) = 0, \\ \partial_{t} E_{z} - c {curl}_{⊥} (B_{⊥}) = - 4 π \sum_{α} q_{α} n_{α} v_{α, z}, \\ \partial_{t} B_{z} + c {curl}_{⊥} (E_{⊥}) = 0, \end{matrix}

(5)

for Maxwell's equations,

\partial_{t} n_{α} + \nabla_{⊥} \cdot (n_{α} v_{α, ⊥}) + \partial_{z} (n_{α} v_{α, z}) = 0

(6)

for the continuity equations, and

\begin{matrix} \partial_{t} v_{α, ⊥} & + \frac{1}{2} \nabla_{⊥} (| v_{α, ⊥} |^{2} + v_{α, z}^{2} + {(n_{α} U_{α} (n_{α}))}_{n_{α}}) \\ = v_{α, z} (\nabla_{⊥} v_{α, z} - \partial_{z} v_{α, ⊥} + \hat{z} \times \frac{q_{α}}{m_{α} c} B_{⊥}) \\ + \hat{z} \times v_{α, ⊥} ({curl}_{⊥} (v_{α, ⊥}) - \frac{q_{α}}{m_{α} c} B_{z}) + \frac{q_{α}}{m_{α}} E_{⊥} \\ \partial_{t} v_{α, z} & + \frac{1}{2} \partial_{z} (| v_{α, ⊥} |^{2} + v_{α, z}^{2} + {(n_{α} U_{α} (n_{α}))}_{n_{α}}) - \frac{q_{α}}{m_{α}} E_{z} \\ = - v_{α, ⊥} \cdot (\nabla_{⊥} v_{α, z} - \partial_{z} v_{α, ⊥} + \hat{z} \times \frac{q_{α}}{m_{α} c} B_{⊥}) \end{matrix}

(7)

for the momentum equations, where we define

\nabla_{⊥} f = (\begin{matrix} \partial_{x} f \\ \partial_{y} f \\ 0 \end{matrix}) and {curl}_{⊥} (A) = \partial_{x} A_{y} - \partial_{y} A_{x} .

(8)

B. Laser–plasma interaction—1D linear regime

We consider here what is frequently called the linear regime¹ in laser-plasma literature. In this regime, we assume the transverse electromagnetic fields scale as

E_{⊥}, B_{⊥} \sim ϵ ≪ 1

⁠, the longitudinal electromagnetic fields scale as

E_{z}, B_{z} \sim ϵ^{2}

⁠, gradients are small in the transverse direction

\nabla_{⊥} \sim ϵ^{2}

while

\partial_{z} \sim O (1)

⁠, and we expand the fluid variables in powers of ϵ: e.g.,

n_{α} = \sum_{k} n_{α}^{(k)}

where

n_{α}^{(k)} \sim ϵ^{k}

⁠. We further assume that all transverse fields are such that

A_{⊥} (x, t) = A_{⊥ 0} (x, t) cos (ω_{0} t + k_{0} z),

(9)

where

ω_{0} \sim ϵ^{- 1}

and

\partial_{t} A_{⊥ 0} \sim ϵ

⁠.

We first consider the continuity and momentum equations order by order. At leading order,

O (ϵ^{0})

⁠, we find that the fluid variables become decoupled from the electromagnetic fields as the electromagnetic fields are asymptotically smaller and we find that

n_{α}^{(0)} = n_{0 α}, v_{α, ⊥}^{(0)} = 0, and v_{α, z}^{(0)} = 0.

(10)

At

O (ϵ^{1})

⁠, we have

\partial_{t} v_{α, ⊥}^{(1)} = \frac{q_{α}}{m_{α}} E_{⊥}, v_{α, z}^{(1)} = 0, and n_{α}^{(1)} = 0.

(11)

Hence, due to the oscillatory characteristic of the transverse electromagnetic fields, we find

v_{α, ⊥}^{(1)} = \frac{q_{α}}{m_{α} ω_{0}} E_{0} sin (ω_{0} t + k_{0} z) .

(12)

The transverse momentum equation along with Faraday's equation at

O (ϵ^{1})

also tells us

\hat{z} \times \partial_{z} v_{α, ⊥}^{(1)} - \frac{q_{α}}{m_{α} c} B_{⊥} = O (ϵ^{2}) .

(13)

Finally, at

O (ϵ^{2})

⁠, we find that

\begin{matrix} \partial_{t} v_{α, ⊥}^{(2)} = 0, \\ \partial_{t} v_{α, z}^{(2)} + \frac{1}{2} \partial_{z} (| v_{α, ⊥}^{(1)} |^{2} + (n_{α} {(U_{α} (n_{α}))}_{n_{α}})) = \frac{q_{α}}{m_{α}} E_{z}, \\ \partial_{t} n_{α}^{(2)} + \partial_{z} (n_{α}^{(0)} v_{α, z}^{(2)}) = 0. \end{matrix}

(14)

Because of the asymptotic ordering, Maxwell's equations in the transverse direction become

\begin{matrix} \partial_{t} E_{0 ⊥} cos (ω_{0} t + k_{0} z) - c \hat{z} \times \partial_{z} B_{0 ⊥} cos (ω_{0} t + k_{0} z) \\ = - 4 π \sum_{α} q_{α} n_{α}^{(0)} v_{α, ⊥}^{(1)} \\ \partial_{t} B_{0 ⊥} + c \hat{z} \times \partial_{z} E_{0 ⊥} = 0. \end{matrix}

(15)

We can see that we have simply obtained evolution equations for the envelopes of the transverse fields. Additional simplification is needed in order to address the current source in the transverse Ampère's law; however, this is deferred to later when the model is described as a Hamiltonian system. What we will find is that the current source in the transverse equation yields a nonlinear polarization of the medium. To leading order, the longitudinal Maxwell equations become

\begin{matrix} \partial_{t} E_{z} = - 4 π \sum_{α} q_{α} n_{α}^{(0)} v_{α, z}^{(2)}, \\ \partial_{t} B_{z} = 0. \end{matrix}

(16)

Hence, the longitudinal field is electrostatic. The longitudinal and transverse fields have decoupled (except through the current sources in Ampère's law).

While this direct asymptotic derivation reveals many of the characteristics of the desired model, the full details of the reduced model are more easily obtained through consideration of the Hamiltonian formulation of the model. This approach ensures energy conservation in the reduced model and better motivates the intricate manner in which the plasma and the field couple in this asymptotic regime.

III. HAMILTONIAN FORMULATION OF THE 1D PONDEROMOTIVE MAXWELL SYSTEM

As we saw in Sec. II, the model truncates at $O (ϵ^{2})$ yielding a complete set of equations. We now derive the Hamiltonian structure of this model beginning from the two-fluid Maxwell–Euler equations, which is known to be Hamiltonian.¹⁴

A. Deriving the Hamiltonian

In order to add special relativistic effects to the two-fluid Maxwell–Euler model, only the kinetic energy needs to be modified

\begin{matrix} H [p_{α, ⊥}, E_{⊥}, B_{⊥}, E_{z}, B_{z}, p_{α, z}, n_{α}] \\ = \sum_{α} \int n_{α} (m_{α} (γ_{α} - 1) c^{2} + U_{α} (n_{α})) d^{3} x + \frac{1}{8 π} \int (| E_{⊥} |^{2} + E_{z}^{2} + | B_{⊥} |^{2} + B_{z}^{2}) d^{3} x, \end{matrix}

(17)

where

p_{α} = γ_{α} m_{α} v_{α}

is the relativistic momentum,

(p_{α, ⊥}, p_{α, z})

are its transverse and longitudinal components, and the Lorentz factor is defined as

γ_{α} = {(1 + \frac{| p_{α, ⊥} |^{2} + p_{α, z}^{2}}{{(m_{α} c)}^{2}})}^{1 / 2} .

(18)

For the purpose of this paper, we need only weakly relativistic effects. Assuming

| p_{α} | ≪ m_{α} c

⁠, we obtain the lowest order relativistic correction to the classical Hamiltonian

\begin{matrix} H [p_{α, ⊥}, E_{⊥}, B_{⊥}, E_{z}, B_{z}, p_{α, z}, n_{α}] \\ = \sum_{α} \int n_{α} [U_{α} (n_{α}) + \frac{1}{2 m_{α}} (| p_{α, ⊥} |^{2} + p_{α, z}^{2} \\ - \frac{1}{4 m_{α}^{2} c^{2}} (| p_{α, ⊥} |^{4} + 2 | p_{α, ⊥} |^{2} p_{α, z}^{2} + p_{z}^{4}))] d^{3} x \\ + \frac{1}{8 π} \int (| E_{⊥} |^{2} + E_{z}^{2} + | B_{⊥} |^{2} + B_{z}^{2}) d^{3} x . \end{matrix}

(19)

We will find that the inclusion of this leading order relativistic correction of the Hamiltonian results in the electromagnetic medium having an intensity dependent index of refraction.

In order to obtain the Hamiltonian for the reduced model, we make the following substitutions:

\begin{matrix} E_{⊥} = E_{0 ⊥} cos (θ), B_{⊥} = B_{0 ⊥} cos (θ), B_{z} = O (ϵ^{3}), \\ p_{α, ⊥} = \frac{q_{α}}{ω_{0}} E_{0} sin (θ) + O (ϵ^{3}), p_{z} = p_{z}^{(2)} + O (ϵ^{3}), \\ n_{α} = n_{α}^{(0)} + n_{α}^{(2)} + O (ϵ^{3}), \end{matrix}

(20)

where

θ = ω_{0} t + k_{0} z

⁠. At leading order, the relativistic transverse momentum equation is

\partial_{t} p_{α, ⊥} = q_{α} E_{⊥},

(21)

which implies the expression for

p_{α, ⊥}

in Eq. . Furthermore, by assuming gradients in the transverse plane are small, we effectively assume the fields are functions of the longitudinal coordinate and time alone.

Retaining only up to

O (ϵ^{4})

and averaging over θ, the kinetic energy becomes

\begin{matrix} H_{K E} = \sum_{α} \int [\frac{n_{α}^{(0)}}{2 m_{α}} p_{α, z}^{2} + \frac{q_{α}^{2} n_{α}^{(0)}}{4 m_{α} ω_{0}^{2}} \\ (1 + \frac{n_{α}^{(2)}}{n_{α}^{(0)}} - \frac{3 q_{α}^{2}}{16 c^{2} m_{α}^{2} ω_{0}^{2}} | E_{0 ⊥} |^{2}) | E_{0 ⊥} |^{2}] d z . \end{matrix}

(22)

We can see that the longitudinal momentum reduces to the purely classical expression (so that we may let

p_{α, z} = m_{α} v_{α, z}

⁠) and the transverse momentum contributes a polarization to the transverse electromagnetic field. The longitudinal electromagnetic energy becomes

H_{E M, z} = \frac{1}{8 π} \int E_{z}^{2} d z .

(23)

The longitudinal magnetic energy is a lower order contribution. Finally, the transverse electromagnetic energy becomes

H_{E M, ⊥} = \frac{1}{16 π} \int (| E_{0 ⊥} |^{2} + | B_{0 ⊥} |^{2}) d z .

(24)

In averaging over θ in the transverse electromagnetic energy, we have picked up an additional factor of 1/2. This seemingly spurious factor of 1/2 in the energy cancels out with a compensating factor which appears in the Poisson bracket. This is most easily seen by considering the phase space Lagrangian formulation of the Maxwell subsystem. Explicitly, if

E_{⊥} = E_{0 ⊥} sin (θ)

and

A_{⊥} = A_{0 ⊥} sin (θ)

⁠, where

A_{⊥}

is the vector potential for the transverse field and

\dot{θ} = ω_{0} = const

⁠, then we find that

\frac{1}{2 π} \int_{0}^{2 π} \int_{Ω} E_{⊥} \cdot {\dot{A}}_{⊥} d z d θ = \frac{1}{2 π} \int_{0}^{2 π} \int_{Ω} (E_{0 ⊥} \cdot {\dot{A}}_{0 ⊥} {cos}^{2} (θ) - ω_{0} E_{0 ⊥} \cdot A_{0 ⊥} sin (θ) cos (θ)) d z d θ = \frac{1}{2} \int_{Ω} E_{0 ⊥} \cdot {\dot{A}}_{0 ⊥} d z .

(25)

Hence, the canonical symplectic form picks up a factor of 1/2. The canonical Poisson bracket, being the inverse of the symplectic form, is then multiplied by two.

It is awkward to perform this averaging directly in the Poisson bracket. It would be interesting to study a more complex model with a nonlinear WKB ansatz, which retains the dynamics of the phase variable, but this is beyond the scope of this paper. Previous work has been done in this vein starting from the phase-space Lagrangian formulation in order to study acoustic waves.¹⁵

Hence, we find that the Hamiltonian for the 1D ponderomotive Maxwell system is

\begin{matrix} H [v_{α, z}, n_{α}^{(2)}, E_{z}, E_{0 ⊥}, B_{0 ⊥}] = \frac{1}{2} \sum_{α} \int (m_{α} n_{α}^{(0)} v_{α, z}^{2} + n_{α}^{(0)} {\bar{U}}_{α} (n_{α}^{(2)})) d z + \frac{1}{8 π} \int E_{z}^{2} d z \\ + \frac{1}{16 π} \int [(1 + \frac{4 π q_{α}^{2} n_{α}^{(0)}}{m_{α} ω_{0}^{2}} (1 + \frac{n_{α}^{(2)}}{n_{α}^{(0)}} - \frac{3 q_{α}^{2}}{16 c^{2} m_{α}^{2} ω_{0}^{2}} | E_{0 ⊥} |^{2})) | E_{0 ⊥} |^{2} + | B_{0 ⊥} |^{2}] d z, \end{matrix}

(26)

where

{\bar{U}}_{α}

is the modified internal energy we obtain upon truncating

U_{α}

to the appropriate asymptotic order.

B. Derivatives of the Hamiltonian

In order to conform to an existing Hamiltonian modeling paradigm for nonlinear polarization of the electromagnetic field,³ we seek a functional K such that

\begin{matrix} H = K - \int E_{0 ⊥} \cdot \frac{δ K}{δ E_{0 ⊥}} d z + \frac{1}{16 π} \int (| E_{0 ⊥} |^{2} + | B_{0 ⊥} |^{2}) d z, \end{matrix}

(27)

where

δ K / δ E_{⊥}

is the functional derivative, a concept of fundamental importance in Hamiltonian field theories.⁴ See Eq. for a definition of the functional derivative. Small modifications must be made to the Hamiltonian modeling framework for nonlinear polarization³ due to the additional factor of 1/2 appearing in the electromagnetic energy from the averaging procedure.

The functional K satisfying Eq. is

K = \frac{1}{2} \sum_{α} \int [m_{α} n_{α}^{(0)} {(v_{α, z}^{(2)})}^{2} + E_{z}^{2} + n_{α}^{(0)} {\bar{U}}_{α} (n_{α}^{(2)}) - \frac{q_{α}^{2} n_{α}^{(0)}}{4 m_{α} ω_{0}^{2}} (1 + \frac{n_{α}^{(2)}}{n_{α}^{(0)}} - \frac{q_{α}^{2}}{16 c^{2} m_{α}^{2} ω_{0}^{2}} | E_{0 ⊥} |^{2}) | E_{0 ⊥} |^{2}] d z .

(28)

Hence, we define the displacement field via

D_{0 ⊥} = E_{0 ⊥} - 8 π \frac{δ K}{δ E_{0 ⊥}} = [1 + \sum_{α} \frac{4 π q_{α}^{2} n_{α}^{(0)}}{m_{α} ω_{0}^{2}} (1 + \frac{n_{α}^{(2)}}{n_{α}^{(0)}} - \frac{q_{α}^{2}}{8 c^{2} m_{α}^{2} ω_{0}^{2}} | E_{0 ⊥} |^{2})] E_{0 ⊥} .

(29)

We have defined

D_{0 ⊥} = E_{0 ⊥} - 8 π δ K / δ E_{0 ⊥}

rather than the expected expression

D_{0 ⊥} = E_{0 ⊥} - 4 π δ K / δ E_{0 ⊥}

to compensate for the additional factor of 1/2 in the electromagnetic Hamiltonian. The additional factor of 1/2 in the energy is canceled by a corresponding multiplicative factor of 1/2 in the symplectic form. Since

D_{0 ⊥}

is defined via a Legendre transformation, see Appendix A, it likewise must be scaled accordingly.

It will be convenient to perform a change of variables so that

D_{0 ⊥}

is the evolving quantity and not

E_{0 ⊥}

⁠. Consider the coordinate change

\begin{matrix} Φ [E_{0 ⊥}, B_{0 ⊥}, E_{z}, v_{α, z}^{(2)}, n_{α}^{(2)}] = (D_{0 ⊥}, B_{0 ⊥}, E_{z}, v_{α, z}^{(2)}, n_{α}^{(2)}) . \end{matrix}

(30)

Then

\begin{matrix} D Φ [E_{0 ⊥}, B_{0 ⊥}, E_{z}, v_{α, z}^{(2)}, n_{α}^{(2)}] = (\begin{matrix} E & 0 & 0 & 0 & \frac{4 π q_{α}^{2}}{m_{α} ω_{0}^{2}} (E_{0 ⊥} \cdot) \\ I \\ I \\ I \\ I \end{matrix}), \end{matrix}

(31)

where

E = I - 8 π \frac{δ^{2} K}{δ E_{0 ⊥} δ E_{0 ⊥}}

and

(E_{0 ⊥} \cdot)

means the operator which applies the dot product of

E_{0 ⊥}

to a given vector. The derivative of the inverse transformation is then

\begin{matrix} D Φ^{- 1} [E_{0 ⊥}, B_{0 ⊥}, E_{z}, v_{α, z}^{(2)}, n_{α}^{(2)}] = (\begin{matrix} E^{- 1} & 0 & 0 & 0 & - \frac{4 π q_{α}^{2}}{m_{α} ω_{0}^{2}} E^{- 1} (E_{0 ⊥} \cdot) \\ I \\ I \\ I \\ I \end{matrix}) . \end{matrix}

(32)

Therefore, if we define

\begin{matrix} \bar{H} [D_{0 ⊥}, B_{0 ⊥}, E_{z}, v_{α, z}^{(2)}, n_{α}^{(2)}] = H [E_{0 ⊥}, B_{0 ⊥}, E_{z}, v_{α, z}^{(2)}, n_{α}^{(2)}], \end{matrix}

(33)

then (due to the self-adjointness of

E

⁠)

\begin{matrix} \frac{δ \bar{H}}{δ n_{α}^{(2)}} = \frac{δ H}{δ n_{α}^{(2)}} + {(\frac{δ E_{0 ⊥}}{δ n_{α}^{(2)}})}^{*} \frac{δ H}{δ E_{0 ⊥}} \\ = \frac{q_{α}^{2}}{4 m_{α} ω_{0}^{2}} | E_{0 ⊥} |^{2} + {(n_{α}^{(0)} {\bar{U}}_{α} (n_{α}^{(2)}))}_{n_{α}^{(2)}} + {(\frac{δ E_{0 ⊥}}{δ n_{α}^{(2)}})}^{*} E \frac{E_{0 ⊥}}{8 π} \\ = - \frac{q_{α}^{2}}{4 m_{α} ω_{0}^{2}} | E_{0 ⊥} |^{2} + {(n_{α}^{(0)} {\bar{U}}_{α} (n_{α}^{(2)}))}_{n_{α}^{(2)}}, \end{matrix}

(34)

where the asterisks denote the L² adjoint of an operator. Similarly, one finds that

\frac{δ \bar{H}}{δ D_{0 ⊥}} = {(\frac{δ E_{0 ⊥}}{δ D_{0 ⊥}})}^{*} \frac{δ H}{δ E_{0 ⊥}} = \frac{E_{0 ⊥}}{8 π} .

(35)

The functional derivatives of the remaining variables do not change due to the coordinate transformation. Hence, we find that

\begin{matrix} \frac{δ \bar{H}}{δ B_{0 ⊥}} = \frac{B_{0 ⊥}}{8 π}, \frac{δ \bar{H}}{δ E_{z}} = \frac{E_{z}}{4 π}, and \frac{δ \bar{H}}{δ v_{α, z}^{(2)}} = m_{α} n_{α}^{(0)} v_{α, z}^{(2)} . \end{matrix}

(36)

These along with the Poisson bracket generate the flow for the 1D ponderomotive Maxwell system.

C. Deriving the Poisson bracket

The derivation of the Poisson bracket for the reduced model is heuristic and lacks a complete justification. This is because it is awkward to handle the inclusion of the WKB ansatz for the rapidly oscillating fields directly in the Poisson bracket. Instead, one should consider the WKB ansatz beginning from the phase-space Lagrangian for the two-fluid Maxwell–Euler and then derive the Poisson bracket after performing the asymptotic reductions in the phase-space Lagrangian.¹⁵

The two-fluid Maxwell–Euler system was shown to be Hamiltonian.¹⁴ The bracket given therein uses momentum coordinates,

M_{α} = ρ_{α} v_{α}

⁠. One may change variables to formulate the model in terms of velocity⁴ to obtain the Poisson bracket

\begin{matrix} {F, G} = \sum_{α} [\frac{1}{m_{α}} \int (\frac{\nabla \times v_{α}}{n_{α}} \cdot \frac{δ F}{δ v_{α}} \times \frac{δ G}{δ v_{α}}) d^{3} x \\ - \frac{1}{m_{α}} \int (\frac{δ F}{δ v_{α}} \cdot \nabla \frac{δ G}{δ n_{α}} - \frac{δ G}{δ v_{α}} \cdot \nabla \frac{δ F}{δ n_{α}}) d^{3} x \\ + \frac{q_{α}}{c m_{α}^{2}} \int (\frac{B}{n_{α}} \cdot \frac{δ F}{δ v_{α}} \times \frac{δ G}{δ v_{α}}) d^{3} x \\ - \frac{4 π q_{α}}{m_{α}} \int (\frac{δ F}{δ D} \cdot \frac{δ G}{δ v_{α}} - \frac{δ G}{δ D} \cdot \frac{δ F}{δ v_{α}}) d^{3} x] \\ + 4 π c \int (\frac{δ F}{δ D} \cdot \nabla \times \frac{δ G}{δ B} - \frac{δ G}{δ D} \cdot \nabla \times \frac{δ F}{δ B}) d^{3} x, \end{matrix}

(37)

where we have made a further alteration to the model so that

D = E - 4 π δ K / δ E

⁠. This notational change is made to better accommodate a self-consistent treatment of nonlinear polarization,³ which arises in the reduced model. Because in this parent two-fluid model,

K = K [n_{α}, v_{α}]

⁠, it follows that

D = E

⁠. In passing to the reduced model, this will no longer be the case.

The Poisson bracket for the full longitudinal/transverse split Maxwell–Euler system may be written as

\begin{matrix} {F, G} = \sum_{α} [- \frac{1}{m_{α}} \int (\frac{δ F}{δ v_{α, z}} \partial_{z} \frac{δ G}{δ n_{α}} - \frac{δ G}{δ v_{α, z}} \partial_{z} \frac{δ F}{δ n_{α}}) d^{3} x + \int (\frac{δ F}{δ v_{α, ⊥}} \cdot \nabla_{⊥} \frac{δ G}{δ n_{α}} - \frac{δ G}{δ v_{α, ⊥}} \cdot \nabla_{⊥} \frac{δ F}{δ n_{α}}) d^{3} x] \\ + \sum_{α} \frac{1}{m_{α}} [\int \frac{\nabla_{⊥} v_{α, z} - \partial_{z} v_{α, ⊥}}{n_{α}} \cdot (\frac{δ F}{δ v_{α, z}} \frac{δ G}{δ v_{⊥, α}} - \frac{δ G}{δ v_{α, z}} \frac{δ F}{δ v_{⊥, α}}) d^{3} x + \int \frac{{curl}_{⊥} (v_{α, ⊥})}{n_{α}} (\frac{δ F}{δ v_{α, ⊥}} \cdot \hat{z} \times \frac{δ G}{δ v_{α, ⊥}}) d^{3} x] \\ + \sum_{α} \frac{q_{α}}{m_{α}^{2} c} [\int \frac{\hat{z} \times B_{⊥}}{n_{α}} \cdot (\frac{δ F}{δ v_{α, z}} \frac{δ G}{δ v_{⊥, α}} - \frac{δ G}{δ v_{α, z}} \frac{δ F}{δ v_{⊥, α}}) d^{3} x - \int \frac{B_{z}}{n_{α}} (\frac{δ F}{δ v_{α, ⊥}} \cdot \hat{z} \times \frac{δ G}{δ v_{α, ⊥}}) d^{3} x] \\ - \sum_{α} \frac{4 π q_{α}}{m_{α}} [\int (\frac{δ F}{δ D_{z}} \frac{δ G}{δ v_{α, z}} - \frac{δ G}{δ D_{z}} \frac{δ F}{δ v_{α, z}}) d^{3} x + \int (\frac{δ F}{δ D_{⊥}} \cdot \frac{δ G}{δ v_{α, ⊥}} - \frac{δ G}{δ D_{⊥}} \cdot \frac{δ F}{δ v_{α, ⊥}}) d^{3} x] \\ + 4 π c \int [\frac{δ F}{δ D_{⊥}} \times \nabla_{⊥} \frac{δ G}{δ B_{z}} - \frac{δ G}{δ D_{⊥}} \times \nabla_{⊥} \frac{δ F}{δ B_{z}} + \frac{δ F}{δ B_{⊥}} \times \nabla_{⊥} \frac{δ G}{δ D_{z}} - \frac{δ G}{δ B_{⊥}} \times \nabla_{⊥} \frac{δ F}{δ D_{z}} \\ - \hat{z} \cdot (\frac{δ F}{δ D_{⊥}} \times \partial_{z} \frac{δ G}{δ B_{⊥}} - \frac{δ G}{δ D_{⊥}} \times \partial_{z} \frac{δ F}{δ B_{⊥}})] d^{3} x . \end{matrix}

(38)

Using the asymptotic scaling assumptions and the results from Sec. II, we find that, to leading order, all transverse gradients vanish, the bracket generating the Lorentz force (the fifth line of the bracket) cancels with the vorticity bracket (the third and fourth lines of the bracket), the longitudinal magnetic field vanishes, and the current in the transverse plane vanishes (having been instead modeled as a nonlinear dielectric function rather than an independent field). Hence, we obtain the following Poisson bracket for the evolution of the fields

(D_{0 ⊥}, B_{0 ⊥}, E_{z}, v_{α, z}^{(2)}, n_{α}^{(2)})

⁠:

\begin{matrix} {F, G} = - \sum_{α} \frac{1}{m_{α}} \int (\frac{δ F}{δ v_{α, z}^{(2)}} \partial_{z} \frac{δ G}{δ n_{α}^{(2)}} - \frac{δ G}{δ v_{α, z}^{(2)}} \partial_{z} \frac{δ F}{δ n_{α}^{(2)}}) d z - \sum_{α} \frac{4 π q_{α}}{m_{α}} \int (\frac{δ F}{δ E_{z}} \frac{δ G}{δ v_{α, z}^{(2)}} - \frac{δ G}{δ E_{z}} \frac{δ F}{δ v_{α, z}^{(2)}}) d z \\ + 8 π c \int \hat{z} \cdot (\frac{δ F}{δ D_{0 ⊥}} \times \partial_{z} \frac{δ G}{δ B_{0 ⊥}} - \frac{δ G}{δ D_{0 ⊥}} \times \partial_{z} \frac{δ F}{δ B_{0 ⊥}}) d z . \end{matrix}

(39)

As elucidated in Eq. , the averaging procedure causes the electromagnetic symplectic form to be halved. Hence, as the (canonical) Poisson bracket is the inverse of the symplectic form, the electromagnetic Poisson bracket is multiplied by two.

It is generally a delicate procedure to freely eliminate terms from the Poisson bracket based on asymptotic arguments as the resulting Poisson bracket may fail to satisfy the Jacobi identity. However, in this case, there is no cause for concern since the resulting bracket is field-free and anti-symmetric thus automatically satisfying the Jacobi identity.⁴

D. 1D ponderomotive Maxwell equations of motion

We may now obtain the equations of motion for the 1D ponderomotive Maxwell system via the usual means, i.e., for any observable F,

\dot{F} = {F, H}

⁠. We find that

\begin{matrix} \dot{F} = {F, H} = - \sum_{α} \frac{1}{m_{α}} \int (\frac{δ F}{δ v_{α, z}^{(2)}} \partial_{z} (- \frac{q_{α}^{2}}{4 m_{α} ω_{0}^{2}} | E_{0 ⊥} |^{2} + {(n_{α}^{(0)} {\bar{U}}_{α} (n_{α}^{(2)}))}_{n_{α}^{(2)}}) - (m_{α} n_{α}^{(0)} v_{α, z}^{(2)}) \partial_{z} \frac{δ F}{δ n_{α}^{(2)}}) d z \\ - \sum_{α} \frac{4 π q_{α}}{m_{α}} \int (\frac{δ F}{δ E_{z}} (m_{α} n_{α}^{(0)} v_{α, z}^{(2)}) - \frac{E_{z}}{4 π} \frac{δ F}{δ v_{α, z}^{(2)}}) d z + c \int \hat{z} \cdot (\frac{δ F}{δ D_{0 ⊥}} \times \partial_{z} B_{0 ⊥} - E_{0 ⊥} \times \partial_{z} \frac{δ F}{δ B_{0 ⊥}}) d z . \end{matrix}

Hence, we find that

\begin{matrix} \partial_{t} D_{0 ⊥} - c \hat{z} \times \partial_{z} B_{0 ⊥} = 0, \\ \partial_{t} B_{0 ⊥} + c \hat{z} \times \partial_{z} E_{0 ⊥} = 0, \end{matrix} \partial_{t} E_{z} + 4 π q_{α} n_{α}^{(0)} v_{α, z}^{(2)} = 0

(40)

for Maxwell's equations, and

\begin{matrix} m_{α} \partial_{t} v_{α, z}^{(2)} - q_{α} E_{z} = \partial_{z} (- \frac{q_{α}^{2} | E_{0 ⊥} |^{2}}{4 m_{α} ω_{0}^{2}} + {(n_{α}^{(0)} {\bar{U}}_{α} (n_{α}^{(2)}))}_{n_{α}^{(2)}}) \\ \partial_{t} n_{α}^{(2)} + n_{α}^{(0)} \partial_{z} v_{α, z}^{(2)} = 0 \end{matrix}

(41)

for the fluid equations.

If we align the x-axis with the electric field and the y-axis with the magnetic field, we obtain

\begin{matrix} \partial_{t} D_{x} + c \partial_{z} B_{y} = 0, \\ \partial_{t} B_{y} + c \partial_{z} E_{x} = 0, \end{matrix} \partial_{t} E_{z} + 4 π q_{α} n_{α}^{(0)} v_{α, z}^{(2)} = 0

(42)

for Maxwell's equations and

\begin{matrix} m_{α} \partial_{t} v_{α, z}^{(2)} - q_{α} E_{z} = \partial_{z} (- \frac{q_{α}^{2}}{4 m_{α} ω_{0}^{2}} E_{x}^{2} + {(n_{α}^{(0)} {\bar{U}}_{α} (n_{α}^{(2)}))}_{n_{α}^{(2)}}) \\ \partial_{t} n_{α}^{(2)} + n_{α}^{(0)} \partial_{z} v_{α, z}^{(2)} = 0 \end{matrix}

(43)

for the fluid equations. These evolution equations are completed by the constitutive relation

D_{x} = [1 + \sum_{α} \frac{4 π q_{α}^{2} n_{α}^{(0)}}{m_{α} ω_{0}^{2}} (1 + \frac{n_{α}^{(2)}}{n_{α}^{(0)}} - \frac{q_{α}^{2}}{8 c^{2} m_{α}^{2} ω_{0}^{2}} E_{x}^{2})] E_{x} .

(44)

The system of first-order equations given in Eqs. and may be rewritten as the following system of second-order nonlinear wave equations:

\begin{matrix} \partial_{t}^{2} D_{x} = c^{2} \partial_{z}^{2} E_{x} \end{matrix}

and

\begin{matrix} \partial_{t}^{2} n_{α}^{(2)} + \sum_{β} \frac{4 π n_{α}^{(0)} q_{α} q_{β}}{m_{α}} (n_{β}^{(0)} + n_{β}^{(2)}) \\ = \frac{1}{m_{α}} \partial_{z}^{2} (\frac{q_{α}^{2}}{4 m_{α} ω_{0}^{2}} E_{x}^{2} - {(n_{α}^{(0)} {\bar{U}}_{α} (n_{α}^{(2)}))}_{n_{α}^{(2)}}), \end{matrix}

(45)

where we have used the fact that

\partial_{z} E_{z} = \sum_{β} 4 π q_{β} (n_{β}^{(0)} + n_{β}^{(2)}),

(46)

which is conserved as a Casimir invariant of the bracket.

Note, if we model only the electron dynamics, i.e., let

n_{i}^{(2)} = 0

with a neutralizing background,

n_{e}^{(0)} = Z n_{i}^{(0)}

⁠, we find

\begin{matrix} \partial_{t}^{2} D_{x} = c^{2} \partial_{z}^{2} E_{x} \end{matrix}

and

\begin{matrix} \partial_{t}^{2} n_{e}^{(2)} + \frac{4 π n_{e}^{(0)} q_{e}^{2}}{m_{e}} n_{e}^{(2)} = \frac{1}{m_{e}} \partial_{z}^{2} (\frac{q_{e}^{2}}{4 m_{e} ω_{0}^{2}} E_{x}^{2} - {(n_{α}^{(0)} {\bar{U}}_{α} (n_{e}^{(2)}))}_{n_{e}^{(2)}}), \end{matrix}

(47)

where

D_{x} = [1 + \frac{4 π q_{e}^{2} n_{e}^{(0)}}{m_{e} ω_{0}^{2}} (1 + \frac{n_{e}^{(2)}}{n_{e}^{(0)}} - \frac{q_{e}^{2}}{8 c^{2} m_{e}^{2} ω_{0}^{2}} E_{x}^{2})] E_{x} .

(48)

This closely resembles previously studied models² except that the transverse electromagnetic field is described by a wave equation for the electric field rather than one for the vector potential.

E. Dimensionless 1D ponderomotive Maxwell system

Suppose we consider the electron dynamics in which

n_{i}^{(2)} = 0

and

n_{e}^{(0)} = Z n_{i}^{(0)}

neutralizes the background ion density. For notational simplicity, we denote the dynamical fields by

(E_{x}, B_{y}, E_{z}, v_{z}, n)

since there is no longer possibility for confusion of the plasma species. The background electron density is assumed constant and denoted n₀. Suppose we nondimensionalize the fields such that

\begin{matrix} v_{z} = c {\tilde{v}}_{z}, n = n_{0} \tilde{n}, B_{y} = B_{0} {\tilde{B}}_{y}, \\ E_{x} = B_{0} {\tilde{E}}_{x}, and E_{z} = B_{0} {\tilde{E}}_{z} . \end{matrix}

(49)

Recall, the (electron) plasma and cyclotron frequencies are given by

ω_{p} = \sqrt{\frac{4 π n_{0} q_{e}^{2}}{m_{e}}}, and ω_{c} = \frac{q_{e} B_{0}}{m_{e} c} .

(50)

We nondimensionalize space and time via

t = ω_{p}^{- 1} \tilde{t}

and

z = c / ω_{p} \tilde{z}

⁠.

Dropping the tildes for notational convenience, the dimensionless equations of motion may be written as

\begin{matrix} \partial_{t} D_{x} + \partial_{z} B_{y} = 0 \\ \partial_{t} B_{y} + \partial_{z} E_{x} = 0 \end{matrix} \partial_{t} E_{z} + \frac{ω_{p}}{ω_{c}} v_{z} = 0

(51)

for Maxwell's equations, and

\begin{matrix} \partial_{t} v_{z} - \frac{ω_{c}}{ω_{p}} E_{z} + \partial_{z} (\frac{ω_{c}^{2}}{4 ω_{0}^{2}} E_{x}^{2} - {(\bar{U} (n))}_{n}) = 0 \\ \partial_{t} n + \partial_{z} v_{z} = 0 \end{matrix}

(52)

for the fluid equations where

D_{x} = [1 + \frac{ω_{p}^{2}}{ω_{0}^{2}} (1 + n - \frac{ω_{c}^{2}}{8 ω_{0}^{2}} E_{x}^{2})] E_{x} .

(53)

We have implicitly rescaled the internal energy so as to make it dimensionless. There are three dimensionless parameters of interest in this system:

ω_{p} / ω_{c}, ω_{p} / ω_{0}

⁠, and

ω_{c} / ω_{0}

⁠. As a nonlinear wave equation, this system may be written as

\begin{matrix} \partial_{t}^{2} D_{x} = \partial_{z}^{2} E_{x} \\ \partial_{t}^{2} n + n = \partial_{z}^{2} (\frac{ω_{c}^{2}}{4 ω_{0}^{2}} E_{x}^{2} - {(\bar{U} (n))}_{n}) . \end{matrix}

(54)

In this form, it is clear that this model is simply the Maxwell wave equation (for a plane wave) nonlinearly coupled to an acoustic wave equation.

The dimensionless formulation of the model in Eqs. (53) and (54) as a second-order nonlinear wave equation shows that there are only two independent parameters: $(ω_{p} / ω_{0}, ω_{c} / ω_{0})$ ⁠. The magnitude of the cyclotron frequency controls the strength of the ponderomotive force (which is intuitive since the ponderomotive force arises from the fast periodic cyclotron motion in the transverse plane), while the overall strength of the nonlinear index of refraction is controlled by the plasma frequency. There is the additional effect of the relativistic correction whose strength is proportional to the product of the squares of the two frequencies, but this term is smaller (assuming the dimensionless parameters are small). Hence, for the most part, these two parameters independently control the strength of the two nonlinear features of the model.

This system may be written as a Hamiltonian system by letting

\begin{matrix} H = \frac{1}{2} \int [\frac{ω_{p}^{2}}{ω_{c}^{2}} (v_{z}^{2} + \bar{U} (n)) + E_{z}^{2} + \frac{1}{2} (B_{y}^{2} \\ + (1 + \frac{ω_{p}^{2}}{ω_{0}^{2}} (1 + n - \frac{3 ω_{c}^{2}}{16 ω_{0}^{2}} E_{x}^{2})) E_{x}^{2})] d z . \end{matrix}

(55)

One can show that

\begin{matrix} \frac{δ H}{δ v_{z}} = \frac{ω_{p}^{2}}{ω_{c}^{2}} v_{z}, \frac{δ H}{δ n} = - \frac{ω_{p}^{2}}{4 ω_{0}^{2}} E_{x}^{2} + \frac{ω_{p}^{2}}{ω_{c}^{2}} {(\bar{U} (n))}_{n}, \\ \frac{δ H}{δ E_{z}} = E_{z}, \frac{δ H}{δ D_{x}} = \frac{E_{x}}{2}, and \frac{δ H}{δ B_{y}} = \frac{B_{y}}{2} . \end{matrix}

(56)

The Poisson bracket is given by

\begin{matrix} {F, G} = - \frac{ω_{c}^{2}}{ω_{p}^{2}} \int (\frac{δ F}{δ v_{z}} \partial_{z} \frac{δ G}{δ n} - \frac{δ G}{δ v_{z}} \partial_{z} \frac{δ F}{δ n}) d z \\ - \frac{ω_{c}}{ω_{p}} \int (\frac{δ F}{δ E_{z}} \frac{δ G}{δ v_{z}} - \frac{δ G}{δ E_{z}} \frac{δ F}{δ v_{z}}) d z \\ - 2 \int (\frac{δ F}{δ D_{x}} \partial_{z} \frac{δ G}{δ B_{y}} - \frac{δ G}{δ D_{x}} \partial_{z} \frac{δ F}{δ B_{y}}) d z . \end{matrix}

(57)

One can show that this bracket and Hamiltonian generate the desired equations of motion. This Poisson bracket possesses several Casimir invariants (assuming periodic or homogeneous boundary conditions). Three are the consequence of the fact that the partial derivative of a constant is zero

C_{D_{x}} = \int D_{x} d z, C_{B_{y}} = \int B_{y} d z, and C_{n} = \int n d z .

(58)

The final of these represents conservation of the total number of particles whereas the other two conservation laws arise from the symmetric configuration of plane electromagnetic waves (analogous conservation laws do not exist in the full 3D Maxwell's equations). Finally, we also see that

C_{E_{z}} = \int C (\partial_{z} E_{z} + \frac{ω_{p}}{ω_{c}} n) d z \forall C \in C^{\infty}

(59)

is a Casimir invariant. This invariant corresponds to conservation of charge.

IV. A STRUCTURE PRESERVING DISCRETIZATION OF THE 1D PONDEROMOTIVE MAXWELL SYSTEM

We use a finite element exterior calculus (FEEC) spectral element method to discretize the 1D ponderomotive Maxwell system. Moreover, we consider both conforming and broken FEEC methods (broken FEEC methods are discontinuous Galerkin methods which nonetheless preserve the conforming de Rham complex structure^6–10). Here, we will only provide the briefest possible exposition of notation used in this section. See Appendix B for a more detailed explanation of the spectral element FEEC method and the definitions of the notation used herein.

Briefly, in 1D, the FEEC method models fields as expansions in one of two different bases

ϕ_{h} = \sum_{j = 1}^{N_{0}} ϕ_{j} Λ_{0, j} (z) and ρ_{h} = \sum_{j = 1}^{N_{1}} ρ_{j} Λ_{1, j} (z),

(60)

where

{Λ_{ℓ, j}}_{j = 1}^{N_{ℓ}}

are the finite element basis functions and

ϕ

and

ρ

are the Galerkin coefficients (the degrees of freedom). The reason for using two different bases of finite element expansions is that this allows us to capture physically meaningful qualities of our fields by specializing the bases used to represent them. The basis functions,

{Λ_{ℓ, j}}_{j = 1}^{N_{ℓ}}

⁠, provide bases for the finite element spaces

V_{h}^{0}

and

V_{h}^{1}

which in turn approximate the infinite dimensional function spaces

V^{0} : = H^{1} Λ^{0} (Ω)

and

V^{1} : = L^{2} Λ^{1} (Ω)

⁠. The former space models fields which transform as scalar 0-forms while the latter models fields which transform as scalar 1-forms (see Appendix B).

The interpolation operators $I^{ℓ} : C_{ℓ} \to V_{h}^{ℓ}$ are simply those operators which take a given vector of Galerkin coefficients and returns the interpolated field: e.g., $I^{0} ϕ = ϕ_{h}$ ⁠. We also define degrees of freedom operators $σ^{ℓ} : V^{ℓ} \to C_{ℓ}$ which simply take a field and return the appropriate degrees of freedom associated with that field, e.g., $σ^{0} (ϕ_{h}) = ϕ$ ⁠.

Finally, without going too much into the technical details, there are also dual spaces, ${(V_{h}^{ℓ})}^{*}$ ⁠, dual interpolation, $I_{ℓ}^{*}$ ⁠, and dual degrees of freedom operators, $σ_{ℓ}^{*}$ ⁠, which are introduced to discretize the adjoint of the derivative. We represent some variables as coming from these dual spaces and thus have coefficients in spaces which we denote $C_{ℓ}^{*}$ ⁠. Roughly speaking, $σ_{ℓ}^{*} : {(V^{ℓ})}^{*} \to C_{ℓ}^{*}$ are L²-projections. Any equation which contains an adjoint derivative must be discretized by applying such an L²-projection in order to integrate by parts and discretize the adjoint derivative weakly.

The finite element mass matrices,

M_{ℓ} : C_{ℓ} \to C_{ℓ}^{*}

⁠, are defined

{(M_{ℓ})}_{i j} = {(Λ_{ℓ, i}, Λ_{ℓ, j})}_{L^{2} (Ω)}

⁠. By Proposition 3 in Appendix B, a variable on the dual complex is related to its primal complex counterpart by the mass matrix

\tilde{A} \in C_{ℓ}^{*} \Leftrightarrow A \in C^{ℓ} and \tilde{A} = M_{ℓ} A .

(61)

A. Direct projection of the equations of motion

In this section, we discretize the dimensionless 1D ponderomotive Maxwell equations derived in Sec. III. For simplicity, we have taken the internal energy to be $\bar{U} (n) = 0$ so that there is no thermodynamic pressure. Before considering the Hamiltonian structure preserving discretization in detail, we first look at a discretization achieved via a direct application of the degrees of freedom operators to the equations of motion.

The first step to discretize the model is to decide which variables belong in which vector space. This is achieved by inspecting the equations and the Poisson bracket. We can see from the Poisson bracket that: D_x and B_y are dual to each other, n and v_z are dual to each other, and E_z and v_z are dual to each other. One option (there are others) is to let

(E_{z}, B_{y}, n)

be modeled on the primal sequence while (D_x, v_z) are on the dual sequence. We chose to model v_z on the dual sequence so that the derivative on the electric field squared might be weakened via integration by parts. Hence,

E_{z} \in V^{0}, B_{y}, n \in V^{1}, and {\tilde{D}}_{x}, {\tilde{v}}_{z} \in V_{0}^{*},

(62)

where the tildes indicate that a variable is an element of the dual space. We shall occasionally need to make use of representations of a same field in both the primal and dual space. The presence or absence of a tilde indicates which space is meant: e.g.,

B_{y} \in V^{1}

while

{\tilde{B}}_{y} \in {(V^{1})}^{*}

⁠. The two representations are related through the L² duality map. A subscript “h” denotes a field's FEM representation: e.g.,

A \in V^{ℓ} \Rightarrow A_{h} = Π^{ℓ} A

⁠.

We may directly discretize the equations by applying the appropriate degrees of freedom operators to each equation

\begin{matrix} \partial_{t} σ_{0}^{*} ({\tilde{D}}_{x}) + d_{0}^{*} σ_{1}^{*} ({\tilde{B}}_{y}) & = 0 \\ \partial_{t} σ^{1} (B_{y}) + d_{0} σ^{0} (E_{x}) & = 0 \\ \partial_{t} σ^{0} (E_{z}) + \frac{ω_{p}}{ω_{c}} σ^{0} (v_{z}) & = 0 \end{matrix} \Rightarrow \begin{matrix} \partial_{t} {\tilde{D}}_{x} + d_{0}^{*} {\tilde{B}}_{y} & = 0, \\ \partial_{t} B_{y} + d_{0} E_{x} & = 0, \\ \partial_{t} E_{z} + \frac{ω_{p}}{ω_{c}} v_{z} & = 0, \end{matrix}

(63)

for Maxwell's equations, and the fluid equations become

\begin{matrix} \begin{matrix} \partial_{t} σ_{0}^{*} ({\tilde{v}}_{z}) - \frac{ω_{c}}{ω_{p}} σ_{0}^{*} ({\tilde{E}}_{z}) + \frac{ω_{c}^{2}}{4 ω_{0}^{2}} d_{0}^{*} σ_{1}^{*} (\tilde{E_{x}^{2}}) = 0 \\ \partial_{t} σ^{1} (n) + d_{0} σ^{1} (v_{z}) = 0 \end{matrix} \\ \Rightarrow \begin{matrix} \partial_{t} {\tilde{v}}_{z} - \frac{ω_{c}}{ω_{p}} {\tilde{E}}_{z} + \frac{ω_{c}^{2}}{4 ω_{0}^{2}} d_{0}^{*} \tilde{I} = 0 \\ \partial_{t} n + d_{0} v_{z} = 0, \end{matrix} \end{matrix}

(64)

where

\tilde{I} = σ_{1}^{*} (E_{x, h}^{2}), E_{x, h} = Π_{0} E_{x}

⁠, and the adjoint discrete exterior derivative is defined

d_{0}^{*} = M_{0}^{- 1} d_{0}^{T} M_{1}

⁠. The constitutive relation is given by

\begin{matrix} {({\tilde{D}}_{x})}_{i} & = σ_{0, i}^{*} ({\tilde{D}}_{x, h}) \\ = σ_{0, i}^{*} ([1 + \frac{ω_{p}^{2}}{ω_{0}^{2}} (1 + n_{h} - \frac{ω_{c}^{2}}{8 ω_{0}^{2}} E_{x, h}^{2})] E_{x, h}) \\ = {(Λ_{0, i}, [1 + \frac{ω_{p}^{2}}{ω_{0}^{2}} (1 + n_{h} - \frac{ω_{c}^{2}}{8 ω_{0}^{2}} E_{x, h}^{2})] E_{x, h})}_{L^{2} (Ω)} . \end{matrix}

(65)

We shall find that precisely the same semi-discrete equations are obtained by directly discretizing the Hamiltonian structure. Hence, this system of ODEs is Hamiltonian.

B. A Hamiltonian structure preserving discretization of the 1D ponderomotive Maxwell system

In order to derive the spatially discrete Hamiltonian structure, we simply make the ansatz

E_{z, h} \in V_{h}^{0}, B_{y, h}, n_{h} \in V_{h}^{1}, and {\tilde{D}}_{x, h}, {\tilde{v}}_{z, h} \in V_{0, h}^{*} .

(66)

We then restrict the Hamiltonian and Poisson bracket to act on functionals of the form

F [{\tilde{D}}_{x, h}, B_{y, h}, {\tilde{v}}_{z, h}, E_{z, h}, n_{h}] = F [{\tilde{D}}_{x}, B_{y}, {\tilde{v}}_{z}, E_{z}, n] .

(67)

To compute the discretized bracket, we use the functional chain rule, the fact that

{({σ^{ℓ} |}_{V_{h}^{ℓ}})}^{†} = I_{ℓ}^{*}

and

{({σ_{ℓ}^{*} |}_{{(V_{h}^{ℓ})}^{*}})}^{†} = I^{ℓ}

⁠, where

†

denotes the adjoint operator (see Appendix C), and

\int Λ_{ℓ, i} {\tilde{Λ}}_{j}^{ℓ} d z = δ_{i j}

⁠. We then find that

\begin{matrix} {F, G} & = [F, G] \\ = - \frac{ω_{c}^{2}}{ω_{p}^{2}} \int ({(σ^{1})}^{†} (\frac{\partial F}{\partial n}) \partial_{z} {(σ_{0}^{*})}^{†} (\frac{\partial G}{\partial {\tilde{v}}_{z}}) - {(σ^{1})}^{†} (\frac{\partial G}{\partial n}) \partial_{z} {(σ_{0}^{*})}^{†} (\frac{\partial F}{\partial {\tilde{v}}_{z}})) d z \\ - \frac{ω_{c}}{ω_{p}} \int ({(σ^{0})}^{†} (\frac{\partial F}{\partial E_{z}}) {(σ_{0}^{*})}^{†} (\frac{\partial G}{\partial {\tilde{v}}_{z}}) - {(σ^{0})}^{†} (\frac{\partial G}{\partial E_{z}}) {(σ_{0}^{*})}^{†} (\frac{\partial F}{\partial {\tilde{v}}_{z}})) d z \\ - 2 \int ({(σ^{1})}^{†} (\frac{\partial F}{\partial B_{y}}) \partial_{z} {(σ_{0}^{*})}^{†} (\frac{\partial G}{\partial {\tilde{D}}_{x}}) - {(σ^{1})}^{†} (\frac{\partial G}{\partial B_{y}}) \partial_{z} {(σ_{0}^{*})}^{†} (\frac{\partial F}{\partial {\tilde{D}}_{x}})) d z, \end{matrix}

(68)

which in turn simplifies to become

\begin{matrix} [F, G] = - \frac{ω_{c}^{2}}{ω_{p}^{2}} [{(\frac{\partial F}{\partial n})}_{i} {(d_{0})}_{i j} {(\frac{\partial G}{\partial {\tilde{v}}_{z}})}_{j} - {(\frac{\partial G}{\partial n})}_{i} {(d_{0})}_{i j} {(\frac{\partial F}{\partial {\tilde{v}}_{z}})}_{j}] - \frac{ω_{c}}{ω_{p}} [{(\frac{\partial F}{\partial E_{z}})}_{i} {(\frac{\partial G}{\partial {\tilde{v}}_{z}})}_{i} - {(\frac{\partial G}{\partial E_{z}})}_{i} {(\frac{\partial F}{\partial {\tilde{v}}_{z}})}_{i}] \\ - 2 [{(\frac{\partial F}{\partial B_{y}})}_{i} {(d_{0})}_{i j} {(\frac{\partial G}{\partial {\tilde{D}}_{x}})}_{j} - {(\frac{\partial G}{\partial B_{y}})}_{i} {(d_{0})}_{i j} {(\frac{\partial F}{\partial {\tilde{D}}_{x}})}_{j}] . \end{matrix}

(69)

Hence, simply by restricting the functional derivatives to act on fields which are Galerkin expansions in the appropriate finite element spaces, we obtain a Poisson bracket for the finite dimensional dynamical system with the coefficients as the dynamical quantities. Because the Poisson bivector is field independent and antisymmetric, this bracket satisfies the Jacobi identity.⁴

The discrete Hamiltonian is likewise obtained by restricting the Hamiltonian functional to act only on fields in the finite element space. Let

E_{x, h} \in V^{0}

and the other fields are as before

\begin{matrix} H [E_{x}, B_{y}, {\tilde{v}}_{z}, E_{z}, n] = H [E_{x, h}, B_{y, h}, {\tilde{v}}_{z, h}, E_{z, h}, n_{h}] \\ = \frac{1}{2} \int [\frac{ω_{p}^{2}}{ω_{c}^{2}} {\tilde{v}}_{z, h}^{2} + E_{z, h}^{2} + \frac{1}{2} ((1 + \frac{ω_{p}^{2}}{ω_{0}^{2}} (1 + n_{h} - \frac{3 ω_{c}^{2}}{16 ω_{0}^{2}} E_{x, h}^{2})) E_{x, h}^{2} + B_{y, h}^{2})] d z \\ = \frac{1}{2} [\frac{ω_{p}^{2}}{ω_{c}^{2}} {\tilde{v}}_{z}^{T} M_{0}^{- 1} {\tilde{v}}_{z} + E_{z}^{T} M_{0} E_{z} + \frac{1}{2} ({\tilde{E}}_{*}^{T} E_{x} + B_{y}^{T} M_{1} B_{y})], \end{matrix}

(70)

where we have the Galerkin projection

{({\tilde{E}}_{*})}_{i} = {(Λ_{0, i}, [1 + \frac{ω_{p}^{2}}{ω_{0}^{2}} (1 + n_{h} - \frac{3 ω_{c}^{2}}{16 ω_{0}^{2}} E_{x, h}^{2})] E_{x, h})}_{L^{2} (Ω)} .

(71)

As in the continuous setting, we must define a constitutive relation:

{\tilde{D}}_{x} = {\tilde{D}}_{x} (E_{x}, n)

⁠. This is accomplished through a direct Galerkin projection of the continuous constitutive relation

{({\tilde{D}}_{x})}_{i} = {(Λ_{0, i}, [1 + \frac{ω_{p}^{2}}{ω_{0}^{2}} (1 + n_{h} - \frac{ω_{c}^{2}}{8 ω_{0}^{2}} E_{x, h}^{2})] E_{x, h})}_{L^{2} (Ω)} .

(72)

Notice the difference between the coefficients on the cubic nonlinearity in the definition of

{\tilde{E}}_{*}

and

{\tilde{D}}_{x}

⁠. Recall, we reconstruct the finite element fields from the coefficients via a Galerkin expansion

E_{x, h} = \sum_{i} {(E_{x})}_{i} Λ_{0, i} (z) and n_{h} = \sum_{i} n_{i} Λ_{1, i} (z) .

(73)

As in the continuous setting, the Hamiltonian and Poisson bracket are defined in terms of different variables. Hence, we perform a change of variables in order to define the Hamiltonian in terms of the correct fields

\bar{H} [{\tilde{D}}_{x}, B_{y}, {\tilde{v}}_{z}, E_{z}, n] = H [E_{x}, B_{y}, {\tilde{v}}_{z}, E_{z}, n],

(74)

where we use the implicit function theorem to find the coordinates

(E_{x}, B_{y}, {\tilde{v}}_{z}, E_{z}, n)

in terms of

({\tilde{D}}_{x}, B_{y}, {\tilde{v}}_{z}, E_{z}, n)

⁠.

C. Derivatives of the discrete Hamiltonian

Just as in the continuous setting, we take derivatives of the Hamiltonian using the chain rule. This directly mirrors the process in the continuous setting. Define a functional

K [E_{x}, {\tilde{v}}_{z}, n]

such that

H = K - {(\frac{\partial K}{\partial E_{x}})}^{T} E_{x} + \frac{1}{4} (E_{x}^{T} M_{0} E_{x} + B_{y}^{T} M_{1} B_{y}) .

(75)

Evidently, this is accomplished by letting

K = \frac{1}{2} [\frac{ω_{p}^{2}}{ω_{c}^{2}} {\tilde{v}}_{z}^{T} M_{0}^{- 1} {\tilde{v}}_{z} + E_{z}^{T} M_{0} E_{z} + N [E_{x}, n]],

(76)

where

N [E_{x}, n] = - \frac{ω_{p}^{2}}{ω_{0}^{2}} \int (1 + n_{h} - \frac{ω_{c}^{2}}{16 ω_{0}^{2}} E_{x, h}^{2}) E_{x, h}^{2} d z .

(77)

It then follows that:

{\tilde{D}}_{x} = M_{0} E_{x} - \frac{\partial K}{\partial E_{x}}

(78)

corresponding with the expression in Eq. . Hence, the polarization is characterized by this

K

functional. This will prove a convenient tool in simplifying the chain rule calculation.

We begin by taking partial derivatives of

H

⁠,

\frac{\partial H}{\partial {\tilde{v}}_{z}} = \frac{ω_{p}^{2}}{ω_{c}^{2}} M_{0}^{- 1} {\tilde{v}}_{z}, \frac{\partial H}{\partial E_{z}} = M_{0} E_{z}, and \frac{\partial H}{\partial B_{y}} = \frac{M_{1} B_{y}}{2} .

(79)

If we define the following matrix:

E [E_{x}, n] = M_{0} - \frac{\partial^{2} K}{\partial E_{x} \partial E_{x}},

(80)

then we may write

\frac{\partial H}{\partial E_{x}} = E [E_{x}, n] E_{x} .

(81)

Finally, the derivative with respect to the density is

{(\frac{\partial H}{\partial n})}_{i} = \frac{ω_{p}^{2}}{4 ω_{0}^{2}} {(Λ_{1, i}, E_{x, h}^{2})}_{L^{2} (Ω)} .

(82)

We find that if we define the coordinate change

Φ [E_{x}, B_{y}, E_{z}, {\tilde{v}}_{z}, n] = ({\tilde{D}}_{x}, B_{y}, E_{z}, {\tilde{v}}_{z}, n),

(83)

then

\begin{matrix} D Φ [E_{x}, B_{y}, E_{z}, {\tilde{v}}_{z}, n] (δ E_{x}, δ B_{y}, δ E_{z}, δ {\tilde{v}}_{z}, δ n) \\ = (E [E_{x}, n] δ E_{x}, δ B_{y}, δ E_{z}, δ {\tilde{v}}_{z}, O [E_{x}] δ n), \end{matrix}

(84)

where

{(O [E_{x}] δ n)}_{i} = \frac{ω_{p}^{2}}{ω_{0}^{2}} \sum_{j} \int Λ_{0, i} (z) Λ_{j}^{1} (z) E_{x, h} (z) δ n_{j} d z .

(85)

Hence, we find that

D Φ^{- 1} [E_{x}, B_{y}, E_{z}, {\tilde{v}}_{z}, n] = (\begin{matrix} E^{- 1} & 0 & 0 & 0 & - E^{- 1} O [E_{x}] \\ I \\ I \\ I \\ I \end{matrix}),

(86)

where

I

is the identity matrix.

We define

\bar{H} [{\tilde{D}}_{x}, B_{y}, E_{z}, {\tilde{v}}_{z}, n] = H [E_{x}, B_{y}, E_{z}, {\tilde{v}}_{z}, n]

⁠. Due to the self-adjointness of

E

⁠, we find

\begin{matrix} \frac{\partial \bar{H}}{\partial n} & = \frac{\partial H}{\partial n} + {(\frac{\partial E_{x}}{\partial n})}^{†} \frac{\partial H}{\partial E_{x}} \\ = \frac{ω_{p}^{2}}{4 ω_{0}^{2}} {(Λ^{1}, E_{x, h}^{2})}_{L^{2} (Ω)} - \frac{1}{2} O^{†} [e_{x}] E_{x} \\ = - \frac{ω_{p}^{2}}{4 ω_{0}^{2}} {(Λ^{1}, E_{x, h}^{2})}_{L^{2} (Ω)}, \end{matrix}

(87)

since

\begin{matrix} {(O^{†} [E_{x}] E_{x})}_{i} = \frac{ω_{p}^{2}}{ω_{0}^{2}} \sum_{j} \int Λ_{j}^{0} (z) Λ_{1, i} (z) E_{x, h} (z) {(E_{x})}_{j} d z \\ = \frac{ω_{p}^{2}}{ω_{0}^{2}} {(Λ_{1, i}, E_{x, h}^{2})}_{L^{2} (Ω)} . \end{matrix}

(88)

Similarly, one finds that

\frac{δ \bar{H}}{δ {\tilde{D}}_{x}} = {(\frac{\partial E_{x}}{\partial {\tilde{D}}_{x}})}^{†} \frac{\partial H}{\partial E_{x}} = \frac{E_{x}}{2} .

(89)

The derivatives with respect to the remaining variables are unchanged by the coordinate transformation

\begin{matrix} \frac{\partial \bar{H}}{\partial {\tilde{v}}_{z}} = \frac{\partial H}{\partial {\tilde{v}}_{z}} = \frac{ω_{p}^{2}}{ω_{c}^{2}} M_{0}^{- 1} {\tilde{v}}_{z}, \frac{\partial \bar{H}}{\partial E_{z}} = \frac{\partial H}{\partial E_{z}} = M_{0} E_{z}, \end{matrix}

and

\begin{matrix} \frac{\partial \bar{H}}{\partial B_{y}} = \frac{\partial H}{\partial B_{y}} = \frac{M_{1} B_{y}}{2} . \end{matrix}

(90)

These along with the spatially discrete Poisson bracket generate the flow for the discrete 1D ponderomotive Maxwell system.

D. Spatially discrete equations of motion

The spatially discrete 1D ponderomotive Maxwell system may be obtained via the usual means for a Hamiltonian system

\dot{F} = [F, H]

⁠. Doing so, we find the discrete equations of motion

\begin{matrix} {\dot{B}}_{y} & = - d_{0} E_{x}, \\ {\dot{\tilde{D}}}_{x} & = d_{0}^{T} M_{1} B_{y}, \\ {\dot{E}}_{z} & = - \frac{ω_{p}}{ω_{c}} M_{0}^{- 1} {\tilde{v}}_{z}, \end{matrix} and \begin{matrix} {\dot{\tilde{v}}}_{z} & = - \frac{ω_{c}^{2}}{4 ω_{0}^{2}} d_{0}^{T} \tilde{I} + \frac{ω_{c}}{ω_{p}} M_{0} E_{z}, \\ \dot{n} & = - d_{0} M_{0}^{- 1} {\tilde{v}}_{z}, \end{matrix}

(91)

where

{\tilde{D}}_{x}

is given by Eq. and

\tilde{I} = {(Λ^{1}, E_{x, h}^{2})}_{L^{2} (Ω)} .

(92)

Notice, this is identical to the system obtained via the direct application of the appropriate degrees of freedom operators. It is possible to write a discrete analog of the nonlinear wave equations

\frac{d^{2} {\tilde{D}}_{x}}{d t^{2}} = d_{0}^{T} M_{1} d_{0} E_{x} and \frac{d^{2} n}{d t^{2}} + n = \frac{ω_{c}^{2}}{4 ω_{0}^{2}} d_{0} M_{0}^{- 1} d_{0}^{T} \tilde{I} .

(93)

What is gained from this discretization approach based on the Hamiltonian and Poisson bracket is a clear motivation for the choice of finite dimensional spaces used for each variable and elucidation of the discrete Hamiltonian structure.

E. Discrete conservation laws

From the discrete Hamiltonian structure, we immediately know the conserved energy of the system as well as several discrete conservation laws associated with Casimir invariants of the bracket. Conservation of the Hamiltonian comes from the antisymmetry of the Poisson bracket

\dot{H} = [H, H] = 0.

(94)

From Subsection 2 of Appendix B, we know that the discrete Casimir invariants are given by

\begin{matrix} C_{D_{x}} = \int {\tilde{D}}_{x, h} d z = 1^{T} M_{0} D_{x} = 1^{T} {\tilde{D}}_{x}, \\ C_{B_{y}} = \int B_{y, h} d z = 1^{T} B_{y}, and C_{n} = \int n_{h} d z = 1^{T} n, \end{matrix}

(95)

where

1 = {(1, 1, \dots, 1)}^{T}

⁠. Hence,

\frac{\partial C_{D}}{\partial {\tilde{D}}_{x}} = 1, \frac{\partial C_{B}}{\partial B_{y}} = 1, and \frac{\partial C_{n}}{\partial n} = 1 .

(96)

From this, and the fact that

1 \in null (d_{0})

and

1 \in null (d_{0}^{T})

⁠, we conclude that these are Casimir invariants of the discrete Poisson bracket.

A discrete charge conservation law arises as another Casimir invariant of the discrete bracket

\frac{d}{d t} (d_{0} E_{z} - \frac{ω_{p}}{ω_{c}} n) = 0.

(97)

Hence, if we set

d_{0} E_{z} = \frac{ω_{p}}{ω_{c}} n

as an initial condition, which is the discrete dimensionless Gauss's law, it remains conserved for all of time. For the purpose of numerically studying the charge conservation law, we define the scalar quantity

\begin{matrix} C_{E_{z}} = {| | \partial_{z} E_{z, h} - \frac{ω_{p}}{ω_{c}} n_{h} | |}_{L^{2}} \\ = \sqrt{{(d_{0} E_{z} - \frac{ω_{p}}{ω_{c}} n)}^{T} M_{1} (d_{0} E_{z} - \frac{ω_{p}}{ω_{c}} n)} . \end{matrix}

(98)

F. Temporal discretization via Hamiltonian splitting

The Hamiltonian splits into two pieces giving us two exactly integrable subsystems. Define

\begin{matrix} H_{E} = \frac{1}{2} [E_{z}^{T} M_{0} E_{z} + \frac{1}{2} {\tilde{E}}_{*}^{T} E_{x}] \end{matrix}

and

\begin{matrix} H_{B v} = \frac{1}{2} [\frac{ω_{p}^{2}}{ω_{c}^{2}} {\tilde{v}}_{z}^{T} M_{0}^{- 1} {\tilde{v}}_{z} + \frac{1}{2} B_{y}^{T} M_{1} B_{y}] . \end{matrix}

(99)

It is straightforward to show that these give rise to two subsystems

\begin{matrix} {\dot{B}}_{y} & = - d_{0} E_{x}, \\ {\dot{\tilde{D}}}_{x} & = 0 \\ {\dot{E}}_{z} & = 0, \end{matrix} and \begin{matrix} {\dot{\tilde{v}}}_{z} & = - \frac{ω_{c}^{2}}{4 ω_{0}^{2}} d_{0}^{T} \tilde{I} + \frac{ω_{c}}{ω_{p}} M_{0} E_{z}, \\ \dot{n} & = 0 \end{matrix}

(100)

for

H_{E}

and

\begin{matrix} {\dot{B}}_{y} & = 0, \\ {\dot{\tilde{D}}}_{x} & = d_{0}^{T} M_{1} B_{y} \\ {\dot{E}}_{z} & = - \frac{ω_{p}}{ω_{c}} M_{0}^{- 1} {\tilde{v}}_{z}, \end{matrix} and \begin{matrix} {\dot{\tilde{v}}}_{z} & = 0, \\ \dot{n} & = - d_{0} M_{0}^{- 1} {\tilde{v}}_{z} \end{matrix}

(101)

for

H_{B v}

⁠. It is clear that these two systems are exactly integrable with update rules given by

\begin{matrix} {\tilde{v}}_{z} (t + Δ t) = {\tilde{v}}_{z} (t) + Δ t (- \frac{ω_{c}^{2}}{4 ω_{0}^{2}} d_{0}^{T} \tilde{I} (t) + \frac{ω_{c}}{ω_{p}} M_{0} E_{z} (t)), \end{matrix}

and

\begin{matrix} B_{y} (t + Δ t) = B_{y} (t) - Δ t d_{0} E_{x} (t) \end{matrix}

(102)

for the

H_{E}

subsystem, and

\begin{matrix} n (t + Δ t) = n (t) - Δ t d_{0} M_{0}^{- 1} {\tilde{v}}_{z} (t), \\ E_{z} (t + Δ t) = E_{z} (t) - \frac{ω_{p}}{ω_{c}} Δ t M_{0}^{- 1} {\tilde{v}}_{z} (t), \end{matrix}

and

\begin{matrix} {\tilde{D}}_{x} (t + Δ t) = {\tilde{D}}_{x} (t) + Δ t d_{0}^{T} M_{1} B_{y} (t) \end{matrix}

(103)

for the

H_{B v}

subsystem. Note, each time we solve the

H_{B v}

subsystem, we must use the constitutive relation in Eq. to solve for

E_{x}

and perform the Galerkin projection in Eq. to get

\tilde{I}

before we can advance the

H_{E}

subsystem again.

In order to approximate the flow for the full system we use operator splitting methods.¹⁶ Strang splitting¹⁷ is a popular symmetric second order method

\begin{matrix} e^{Δ t {\cdot, H}} & = e^{Δ t {\cdot, H_{E}} + Δ t {\cdot, H_{B}}} \\ = e^{Δ t / 2 {\cdot H_{E}}} e^{Δ t {\cdot, H_{B}}} e^{Δ t / 2 {\cdot H_{E}}} + O (Δ t^{3}) . \end{matrix}

(104)

This is a favorable approach at second order because only one nonlinear solve is needed per time step (since the

H_{B v}

flow is only computed once per time step). There are also higher order splitting methods which achieve a higher temporal order by having a larger number of substeps for each time step.¹⁸ The derivation of symplectic integrators based on splitting methods for exponential maps is a well-established theory.¹⁹ We only consider second and sixth order symmetric splitting methods.

The symplecticity of the algorithm is only guaranteed up to the precision of the nonlinear solver. Hence, while in the idealized world of infinite precision arithmetic this algorithm is a Poisson integrator, in reality, the flow is only approximately so to the precision of the nonlinear solver. As the nonlinear system is a perturbation of a linear symmetric positive definite system by a nonlinear operator with a small coefficient, it is sufficient to use Picard iteration with Anderson acceleration²⁰ to solve for $E_{x}$ at each time step.

G. Full algorithm to advance 1D ponderomotive Maxwell system

In the following, we shall assume that we may treat the subroutines associated with the FEEC method as a black box. Only five principle routines are needed to accomplish each time step: (1) a routine to perform the nonlinear solve for $E_{x}$ (2) a routine to compute $\tilde{I}$ via an L² projection, (3) a routine to advance the $H_{E}$ subsystem, (4) a routine to advance the $H_{B v}$ subsystem, and (5) a routine which composes together the partial flows to get the flow of the entire system to some specified temporal order.

We assume the availability of some Gaussian quadrature scheme to perform Gauss-Lobatto quadrature at a given number of quadrature points:

\tilde{g} \leftarrow gaussQuadrature (g, N_{Q})

⁠. Moreover, for algorithmic convenience, we consider an expanded state vector

stateVec : = ({\tilde{D}}_{x}, E_{x}, \tilde{I}, E_{z}, B_{y}, {\tilde{v}}_{z}, n) .

(105)

While

E_{x}

and

\tilde{I}

are technically redundant information which might be obtained from the other fields, it is more convenient to store this information rather than recompute as needed. We also have a helper routine given in Algorithm 1 which simply projects a field onto a given FEM basis.

Algorithm 1.

$assembl e_{ℓ} (f; N_{Q})$ ⁠.

\tilde{f} \leftarrow gaussQuadrature ({(Λ_{ℓ, i}, f)}_{L^{2}}, N_{Q})

▹ Quadrature with N_Q points

return

\tilde{f}

We compute $E_{x}$ via fixed point iteration. The initial guess for $E_{x}$ is its value at the previous time step. Assuming a sufficiently small time step, the fixed point iteration converges relatively quickly. In order to further speed up the fixed point iteration, we employed Anderson acceleration, although this is not reflected in the basic algorithm outlined in Algorithm 2. The algorithm to compute $\tilde{I}$ is fairly trivial and given in Algorithm 3.

Algorithm 2.

picardSolveE(⁠ $stateVec; ATOL, RTOL, N_{Q}$ ⁠).

Input:

{\tilde{D}}_{x}, E_{x}, n

Output:

E_{x}

field

(x; E_{x}, n) : = \frac{ω_{p}^{2}}{ω_{0}^{2}} [1 + \sum_{i} n_{i} Λ_{1, i} (x)

\begin{array}{l} - \frac{ω_{c}^{2}}{8 ω_{0}^{2}} {| \sum_{i} {(E_{x})}_{i} Λ_{0, i} (x) |}^{2}] \sum_{i} {(E_{x})}_{i} Λ_{0, i} (x) \\ D_{x} \leftarrow M_{1}^{- 1} {\tilde{D}}_{x} \\ res (E_{x}, P_{x}) : = E_{x} - P_{x} - D_{x} \end{array}

do

\begin{matrix} {\tilde{P}}_{x} \leftarrow assembl e_{1} (field (x; E_{x}, n), N_{Q}) \\ P_{x} \leftarrow M_{1}^{- 1} {\tilde{P}}_{x} \\ R \leftarrow res (E_{x}, P_{x}) \\ E_{x} \leftarrow E_{x} - R \end{matrix}

while (⁠

| R_{i} | < RTOL | D_{i} | + ATOL, \forall i

⁠)

Algorithm 3.

computeIntensity(⁠ $stateVec; N_{Q}$ ⁠).

Input:

E_{x}

Output:

\tilde{I}

\begin{matrix} field (x; E_{x}) : = {| \sum_{i} {(E_{x})}_{i} Λ_{0, i} (x) |}^{2} \\ \tilde{I} \leftarrow assembl e_{1} (field (x; E_{x}), N_{Q}) \end{matrix}

The remaining three routines are straightforward. We simply use the update rules given in Eqs. (102) and (103) to define the two partial flows given by $H_{E}$ and $H_{B v}$ ⁠, and then compose them together using Strang splitting in Algorithm 6. In several tests, we also use higher order splitting methods.¹⁸ We do not write down an algorithm for these methods as the idea is analogous to Strang splitting, but with more sub-steps. A single call of $adv2ndOrd$ advances the system a single time step. In order to ensure stability of the scheme, we require $Δ t / Δ x < 1$ where $Δ x$ is the minimum Gauss-Lobatto grid-spacing among all of the mapped elements. In practice, we let $Δ t / Δ x = 1 / 2$ ⁠.

Algorithm 4.

updateBV(⁠ $Δ t$ ⁠; stateVec).

Input: stateVec

Output:

{\tilde{v}}_{z}, B_{y}

\begin{matrix} {\tilde{v}}_{z} \leftarrow {\tilde{v}}_{z} + Δ t (- \frac{ω_{c}^{2}}{4 ω_{0}^{2}} d_{0}^{T} \tilde{I} + \frac{ω_{c}}{ω_{p}} M_{0} E_{z}) \\ B_{y} \leftarrow B_{y} - Δ t d_{0} E_{x} \end{matrix}

Algorithm 5.

updateEN(⁠ $Δ t; stateVec, ATOL, RTOL, N_{Q}$ ⁠).

Input: stateVec

Output:

n, E_{z}, {\tilde{D}}_{x}

\begin{array}{l} n \leftarrow n - Δ t d_{0} M_{0}^{- 1} {\tilde{v}}_{z} \\ E_{z} \leftarrow E_{z} - \frac{ω_{p}}{ω_{c}} Δ t M_{0}^{- 1} {\tilde{v}}_{z} \\ {\tilde{D}}_{x} \leftarrow {\tilde{D}}_{x} + Δ t d_{0}^{T} M_{1} B_{y} \end{array}

picardSolveE(⁠

stateVec; ATOL, RTOL, N_{Q}

⁠)

computeIntensity(⁠

stateVec; N_{Q}

⁠)

Algorithm 6.

adv2ndOrd(⁠ $Δ t; stateVec; ATOL, RTOL, N_{Q}$ ⁠).

Input: stateVec

Output: stateVec

updateBV(⁠

Δ t / 2

⁠; stateVec)

updateEN(⁠

Δ t; stateVec, ATOL, RTOL, N_{Q}

⁠)

updateBV(⁠

Δ t / 2

⁠; stateVec)

V. NUMERICAL RESULTS FOR 1D PONDEROMOTIVE MAXWELL SYSTEM

In this section, we examine the behavior of the algorithm derived in Sec. IV. All quantities in these tests are dimensionless, see Sec. III E. We take

t \in (0, 15)

⁠, and let the spatial domain,

Ω = (0, 40)

⁠, be periodic. The domain is made long enough that any deflected waves propagating in the reverse direction do not pollute the forward propagating wave. All fields except the transverse electromagnetic fields are initially zero. The transverse electromagnetic fields are Gaussian

E_{0} (z) = B_{0} (z) = exp (- {(z - 4)}^{2}) .

(106)

In a vacuum, this wave would propagate unchanged with unit speed. Hence, we will plot the solutions only on the domain

z \in (0, 25)

since the forward propagating wave remains in this portion of the domain.

Let

ξ_{i, k}

denote the mapped Gauss–Lobatto points in each element and define

Δ x = min_{i, k} | ξ_{i, k} - ξ_{i + 1, k} | .

(107)

The CFL condition requires that

u Δ t / Δ x < C_{max}

where u is the maximum wave speed in our system and

C_{max}

is the maximum Courant parameter allowed by the scheme. It is difficult to determine u exactly, but we make the following heuristic argument. In linear media, the speed of light is

c / \sqrt{ϵ}

where ϵ is the dielectric constant. Hence, assuming our dielectric function remains greater than one, i.e.,

1 + n - \frac{ω_{c}^{2}}{8 ω_{0}^{2}} E_{x}^{2} > 0,

(108)

light waves travel slower than the speed of light in vacuum which, in our dimensionless units, has value unity. One should not blindly assume this inequality remains satisfied, but must consider whether the initial conditions and parameters of the system justify this assumption and whether it continues to remain satisfied through the whole simulation. We found that the inequality was satisfied in all of our simulations as we choose relatively small physical parameters and the density perturbation and electric field remained small. Likewise, determining

C_{max}

exactly would require rigorous linear stability analysis which is beyond the scope of this paper. We empirically observe that the solution is stable when

Δ t = Δ x / 2

for simulations in which the nonlinearities of the system are relatively weak, and when

Δ t = Δ x / 4

for simulations in which the nonlinearity plays a more dominant role (the parameters associated with these two regimes will be made definite subsequently). The temporal resolution is, therefore, set by the spatial resolution of the method and also by the state of the system and the parametric regime in a manner which requires further investigation.

A. Numerical results

In order to understand the properties of the numerical solver, we vary several physical and numerical parameters. Regarding the physical parameters, we consider a weakly nonlinear regime in which $(ω_{p} / ω_{0}, ω_{c} / ω_{0}) = (1 / 5, - 1 / 5)$ ⁠, and a strongly nonlinear regime in which $(ω_{p} / ω_{0}, ω_{c} / ω_{0}) = (2 / 5, - 9 / 10)$ ⁠. The two dimensionless parameters independently control the strength of the two nonlinear features of the model: $ω_{c} / ω_{0}$ controls the strength of the ponderomotive force, and $ω_{p} / ω_{0}$ controls the strength of the nonlinear refractive index. Recall, it is easiest to see the effects of these parameters by looking at the nondimensional system formulated as a coupled system of second-order wave equations in Eq. (54) and the nondimensional constitutive relation (53).

In the weakly nonlinear case, we let the plasma and cyclotron frequencies be the same while in the latter case, the cyclotron frequency is double the plasma frequency. We do not perform an exhaustive parameter sweep in this work to understand the physical effects of the relative strength of these two parameters, but one should expect different regimes in which the physics is dominated by a strong ponderomotive force or a strong nonlinear refractive index. However, as both of these phenomena lead to wave steepening, their effects may be qualitatively similar in this one-dimensional model.

Regarding numerical parameters, we consider both high and low spatial and temporal order solutions for both conforming and broken FEEC discretizations. Recall, from the construction of the spectral element method in Appendix B, the polynomials interpolating the 0-forms are one degree higher than those interpolating the 1-forms. For the low order tests, we use a second order splitting method in time, N = 3 (which yields cubic polynomials for the 0-forms and quadratic for the 1-forms), and K = 200 cells. In the conforming FEEC case, the 0-forms have 600 degrees of freedom while in broken FEEC they have 800. In both cases, the 1-forms have 600 degrees of freedom. For the high order tests, we use a sixth order splitting method in time, N = 6 (which yields sixth-order polynomials for the 0-forms and fifth-order polynomials for the 1-forms), and K = 100 cells. In the conforming FEEC case, the 0-forms have 600 degrees of freedom while in broken FEEC they have 700. In both cases, the 1-forms have 600 degrees of freedom. The time step is implicitly determined by the spatial resolution and was empirically found to require greater resolution in the regimes where the nonlinearities of the system are more pronounced. The low spatial order test cases used $Δ t \approx 0.028$ for the weakly nonlinear case, and $Δ t \approx 0.014$ for the strongly nonlinear case. The high spatial order test cases used $Δ t \approx 0.017$ in the weakly nonlinear case, and $Δ t \approx 0.0085$ in the strongly nonlinear case. In all cases, the nonlinear solver tolerance has been set to $10^{- 13}$ ⁠.

1. Weakly nonlinear regime

We first consider a weakly nonlinear regime in which $(ω_{p} / ω_{0}, ω_{c} / ω_{0}) = (1 / 5, - 1 / 5)$ ⁠. This choice of parameters is within the valid bounds for which the asymptotic assumptions in the derivation of the 1D ponderomotive Maxwell model hold.

Figure 1 summarizes the results of the numerical experiments in this weakly nonlinear regime. The first column plots the error in the energy and Casimir invariants as a function of time. One can see that in each case, the Casimir invariants are conserved up to machine precision while the error in the energy does not grow and has magnitude determined by the order of the splitting method. Symplectic integrators do not exactly conserve the Hamiltonian of the continuous system, but rather that of another system which is a perturbation of the original system. The relationship between the true and perturbed systems may be determined with a technique known as backward error analysis.¹⁹ As a result of this property of symplectic integrators, the computed value of the Hamiltonian oscillates with finite amplitude about the initial value of the energy (the amplitude of this oscillation decreases as the order of the method increases). With second order splitting, the error in the Hamiltonian is found to be $\sim 10^{- 7}$ ⁠. With sixth order splitting, the Hamiltonian is conserved very near machine precision likely because the amplitude of the oscillation is smaller than machine precision.

FIG. 1.

View large Download slide

Weakly nonlinear regime: $(ω_{p} / ω_{0}, ω_{c} / ω_{0}) = (1 / 5, - 1 / 5)$ ⁠. The left and right columns plot the conservation law errors and the solution, respectively. Section IV E defines the discrete Casimir invariants. (a) and (b) are the conforming low and high order methods, respectively; (c) and (d) are the low and high order broken methods, respectively.

The second column shows the solution at three consecutive times. One can see that each solution qualitatively agrees in the behavior of the system, but that the low-order broken-FEEC scheme exhibits spurious oscillations. In general, it was found that the broken-FEEC scheme required higher resolution to suppress these oscillations. The solution behaves as one might anticipate. At t = 7.5 in the second column of Fig. 1, one can see that the transverse field has induced a ponderomotive force pointing in the longitudinal direction and that this force induced a charge separation in the longitudinal direction which in turn excited a longitudinal electrostatic field. Notice, in the first column of the same figures, this longitudinal field satisfies an exact charge conservation law due to a Casimir invariant of the discrete system. The longitudinal field is relatively low amplitude because of the weakness of the nonlinearity in this system. Moreover, the transverse fields remain largely undeformed because the nonlinear polarization is so weak.

2. Strongly nonlinear regime

We now consider a strongly nonlinear regime in which $(ω_{p} / ω_{0}, ω_{c} / ω_{0}) = (2 / 5, - 9 / 10)$ ⁠. These parameters are somewhat larger than the asymptotic derivation allows, however we include these parameters in this numerical study because they more strongly activate the nonlinearities in the model. This serves as a stiffer test for the numerical method and allows us to examine its conservative properties in a more interesting context.

Figure 2 summarizes the results of the numerical experiments in this strongly nonlinear regime. Once again, the error in the conservation laws is in the first column while the computed solutions are shown in the second column. One can see that in each case, all Casimir invariants are conserved up to machine precision. The charge conservation Casimir invariant $C_{E_{z}}$ very slowly grows in magnitude, possibly due to the round-off errors incurred by the nonlinear solver in each step (which has been set to $10^{- 13}$ ⁠). Despite this slow drift, the charge conservation error remains near machine precision through the whole simulation. The Hamiltonian remains well conserved (the magnitude of the error being determined by the order of the splitting method) until a shock begins forming, discussed below, around t = 9 causing spurious oscillations. After this time, the error in the Hamiltonian begins growing presumably because of increased roundoff error as the regularity of the solution breaks down.

FIG. 2.

View large Download slide

Strongly nonlinear regime: $(ω_{p} / ω_{0}, ω_{c} / ω_{0}) = (2 / 5, - 9 / 10)$ ⁠. The left and right columns plot the conservation law errors and the solution respectively. Section IV E defines the discrete Casimir invariants. (a) and (b) The conforming low and high order methods, respectively; (c) and (d) are the low and high order broken methods, respectively.

In the strongly nonlinear regime, one can see that the transverse wave steepens over time, see the second column of Fig. 2. This eventually leads to the formation of a shock which our scheme does not adequately capture. This is because high order polynomial interpolation exhibits Runge's phenomenon when gradients become too large. Hence, in this strongly nonlinear regime, a lower order method with a larger number of cells yields better results, see Sec. V B. While the solution does exhibit spurious oscillations, it continues to conserve the invariants of motion to high precision. This issue serves as a reminder that conservative integrators, while beneficial in their ability to preserve the invariants of motion, are not immune to the challenges and failures of standard numerical methods. In addition to increasing grid resolution in regions where the gradient is large, for example by using an adaptive mesh, there are various remedies one might try in order to better resolve the shocks, i.e., shock capturing methods, but how such techniques might be integrated into a structure preserving framework is beyond the scope of this paper.

3. Sensitivity to tolerance of the nonlinear solver

Each time the field $D_{x}$ is updated (in the $H_{B v}$ sub-step when one calls $updateEN$ ⁠), a new value of $E_{x}$ must be computed using an iterative nonlinear solver. With infinite precision arithmetic, the time advance of the entire system would be a Poisson integrator exactly conserving the Casimir invariants by keeping the solution on the proper symplectic leaf and exactly conserving the symplectic two-form on that leaf.¹⁹ However, the need to perform a nonlinear solve after each $H_{B v}$ sub-step incurs an error which takes the solution off of the symplectic leaf. Therefore, the algorithm is not truly a Poisson integrator, but only approximately so. In this section, we examine the sensitivity of the conservation of the Casimir invariants to the specified tolerance of the nonlinear solver.

In order to see the effects of the tolerance of the nonlinear solver in isolation, we run high-resolution (N = 6, K = 80 cells, and $6 th$ -order splitting), conforming and broken FEEC tests in the strongly nonlinear regime with the same parameters as before, but only up to t = 7.5 to avoid pollution of the results from the known issue the solvers face when a shock forms. In both the conforming and broken FEEC cases, $d t \approx 0.021$ ⁠. In the conforming case, both the 0-forms and 1-forms have about 480 DOFs while in the broken FEEC case, the 0-forms have 560 DOFs (the 1-forms are the same in conforming and broken FEEC). The results for the weakly nonlinear regime aremittedd because conservation in the strongly nonlinear regime is a much stiffer test. See Fig. 3 for the results. One can see that the Casimirs are conserved up to machine precision regardless of the precision of the nonlinear solver. However, conservation of the Hamiltonian is dependent on the tolerance of the nonlinear solver.

FIG. 3.

View large Download slide

Numerical test of the sensitivity of discrete conservation laws to nonlinear solver tolerance: conforming-FEEC on the left and broken-FEEC on the right. See Sec. IV E for definition of the discrete Casimir invariants.

B. Hamiltonian finite difference/lowest order conforming-FEEC scheme

As a point of comparison with the FEEC spectral element methods previously discussed, it is worth briefly considering how a similar Hamiltonian finite difference scheme might be derived. We compare the finite element approximation with a Hamiltonian structure preserving finite difference approximation because the obvious finite difference scheme obtained from replacing all derivatives in Eqs. (51)–(53) with first order finite difference stencils results in an unstable approximation. We find that the Hamiltonian structure preserving finite difference approximation is a special case of finite element method in which one uses uniform grids and lowest order interpolating polynomials. Despite this, it is still instructive to show how the method might be directly derived as a finite difference approximation.

Suppose the domain is subdivided into an equispaced grid

{z_{i}}_{i = 1}^{N}

⁠. The Hamiltonian is discretized via trapezoidal rule which, on a periodic domain, is simply

H_{F D} = \frac{Δ z}{2} \sum_{i = 1}^{N} [\frac{ω_{p}^{2}}{ω_{c}^{2}} v_{z, i}^{2} + E_{z, i}^{2} + \frac{1}{2} ((1 + \frac{ω_{p}^{2}}{ω_{0}^{2}} (1 + n_{i} - \frac{3 ω_{c}^{2}}{16 ω_{0}^{2}} E_{x, i}^{2})) E_{x, i}^{2} + B_{y, i}^{2})],

(109)

where, if A(z) is a continuous field,

A_{i} = A (z_{i})

⁠. We define the relationship between

D_{x}

and

E_{x}

by identifying values at the collocation points

D_{x, i} = [1 + \frac{ω_{p}^{2}}{ω_{0}^{2}} (1 + n_{i} - \frac{ω_{c}^{2}}{8 ω_{0}^{2}} {(E_{x, i})}^{2})] E_{x, i} .

(110)

Finally, taking derivatives of the Hamiltonian proceeds much as in the case of the finite element method. Letting

\bar{H} (D_{x}, B_{y}, v_{z}, n, E_{z}) = H (E_{x}, B_{y}, v_{z}, n, E_{z})

⁠, we obtain

\begin{matrix} \frac{\partial {\bar{H}}_{F D}}{\partial D_{x, i}} = \frac{Δ z E_{x, i}}{2}, \frac{\partial {\bar{H}}_{F D}}{\partial B_{y, i}} = \frac{Δ z B_{y, i}}{2}, \\ \frac{\partial {\bar{H}}_{F D}}{\partial E_{z, i}} = Δ z E_{z, i}, \frac{\partial {\bar{H}}_{F D}}{\partial v_{z, i}} = Δ z \frac{ω_{p}^{2}}{ω_{c}^{2}} v_{z, i}, \end{matrix}

and

\begin{matrix} \frac{\partial {\bar{H}}_{F D}}{\partial n_{i}} = - Δ z \frac{ω_{p}^{2}}{4 ω_{0}^{2}} E_{x, i}^{2} . \end{matrix}

(111)

Define

d_{F D}

to be the circulant matrix with stencil

{(d_{F D})}_{i j} = {\begin{matrix} 1, & j = i + 1, \\ - 1, & j = i, \\ 0, & else . \end{matrix}

(112)

Notice, this derivative matrix is precisely that of the conforming spectral element FEEC method.

Before discretizing the Poisson bracket, we must briefly explain how functional derivatives discretize in this finite difference setting. Letting

{(δ F / δ A) |}_{z = z_{i}}

be the functional derivative evaluated at z = z_i, we have that

\int_{Ω} \frac{δ F}{δ A} δ A d z \approx \sum_{i} {(\frac{δ F}{δ A})}_{z = z_{i}} δ A_{i} Δ z = \sum_{i} \frac{\partial F}{\partial A_{i}} δ A_{i},

(113)

where the integral is approximated with trapezoidal rule. Hence, we find

{(\frac{δ F}{δ A})}_{z = z_{i}} = \frac{1}{Δ z} \frac{\partial F}{\partial A_{i}} .

(114)

The discrete Poisson bracket is then obtained by again replacing the integral with its trapezoidal rule approximation and the derivatives by first order finite differences. One finds

\begin{matrix} {[F, G]}_{F D} = - \frac{ω_{c}^{2}}{ω_{p}^{2}} \sum_{i, j = 1}^{N} [\frac{\partial F}{\partial n_{i}} {(d_{F D})}_{i j} \frac{\partial G}{\partial v_{z, j}} - \frac{\partial G}{\partial n_{i}} {(d_{F D})}_{i j} \frac{\partial F}{\partial v_{z, j}}] \frac{1}{Δ z^{2}} - \frac{ω_{c}}{ω_{p}} \sum_{i = 1}^{N} [\frac{\partial F}{\partial E_{z, i}} \frac{\partial G}{\partial v_{z, i}} - \frac{\partial G}{\partial E_{z, i}} \frac{\partial F}{\partial v_{z, i}}] \frac{1}{Δ z} \\ - 2 \sum_{i, j = 1}^{N} [\frac{\partial F}{\partial B_{y, i}} {(d_{F D})}_{i j} \frac{\partial G}{\partial D_{x, j}} - \frac{\partial G}{\partial B_{y, i}} {(d_{F D})}_{i j} \frac{\partial F}{\partial D_{x, j}}] \frac{1}{Δ z^{2}} . \end{matrix}

(115)

The scaling of the discrete Poisson bracket by

Δ z^{- 2}

comes from appropriately treatment of the discrete functional derivatives. However, this scaling may also be understood via an appropriate interpretation of the variables as living on the primal/dual de Rham complex. Scaling a vector by

Δ z^{2}

may be interpreted as applying the diagonal matrix

Δ z^{2} I

⁠. The finite difference method in fact is identical to the spectral element FEEC method on a uniform grid using lowest order interpolating polynomials. In that method, the mass matrix for the 1-forms is precisely

Δ z^{2} I

⁠. As the bracket obtained with finite differences is similar to the FEEC spectral element Poisson bracket, it possesses an analogous set of Casimir invariants

\begin{matrix} C_{D_{x}} = 1^{T} D_{x} Δ z, C_{B_{y}} = 1^{T} B_{y} Δ z, C_{n} = 1^{T} n Δ z, \end{matrix}

and

\begin{matrix} C_{E_{z}} = {| | Δ z^{- 1} d_{F D} E_{z} - \frac{ω_{p}}{ω_{c}} n | |}_{2} \end{matrix},

(116)

which are the trapezoidal rule approximations of the integrals of the continuous fields (note,

| | \cdot | |_{2}

is the vector 2-norm).

Using the fact that

\dot{F} = {[F, H_{F D}]}_{F D}

⁠, we obtain the equations of motion

\begin{matrix} {\dot{D}}_{x} & = Δ z^{- 1} d_{F D}^{T} B_{y} \\ {\dot{B}}_{y} & = - Δ z^{- 1} d_{F D} E_{x} \\ {\dot{E}}_{z} & = - \frac{ω_{p}}{ω_{c}} v_{z} \end{matrix} \begin{matrix} {\dot{v}}_{z} & = \frac{ω_{c}}{ω_{p}} E_{z} - Δ z^{- 1} \frac{ω_{c}^{2}}{4 ω_{0}^{2}} d_{F D}^{T} (E_{x}^{2}) \\ \dot{n} & = - Δ z^{- 1} d_{F D} v_{z}, \end{matrix}

(117)

which are completed with the constitutive relation for

D_{x} (E_{x}, n)

in Eq. . Notice, this is precisely what is obtained from simply applying a finite difference approximation to the spatial derivatives in Eqs. .

Hence, we found that the Hamiltonian structure preserving spatial discretization using finite differences or lowest order finite elements is simply the finite difference method one might naively guess without all of this extra theoretical machinery. However, in addition to a deeper understanding of the dynamics of the discrete system and its relationship to the continuous system, what is gained from this Hamiltonian perspective is an appropriate operator splitting method to approximate temporal flow. A simple finite difference approximation to the time derivatives yields an unstable method. However, the Hamiltonian is separable and, when appropriately split, we obtain two exactly integrable subsystems

\begin{matrix} {\begin{matrix} {\dot{D}}_{x} = Δ z^{- 1} d_{F D}^{T} B_{y}, \\ {\dot{E}}_{z} = - \frac{ω_{p}}{ω_{c}} v_{z}, \\ \dot{n} = - Δ z^{- 1} d_{F D} v_{z} \end{matrix} \\ and {\begin{matrix} {\dot{B}}_{y} = - Δ z^{- 1} d_{F D} E_{x}, \\ {\dot{v}}_{z} = \frac{ω_{c}}{ω_{p}} E_{z} - Δ z^{- 1} \frac{ω_{c}^{2}}{4 ω_{0}^{2}} d_{F D}^{T} (E_{x}^{2}), \end{matrix} \end{matrix}

(118)

which possess exact solutions

\begin{array}{l} {\begin{array}{l} D_{x} (t + Δ t) = D_{x} (t) + \frac{Δ t}{Δ z} d_{F D}^{T} B_{y} (t), \\ E_{z} (t + Δ t) = E_{z} (t) - Δ t \frac{ω_{p}}{ω_{c}} v_{z} (t), \\ n (t + Δ t) = n (t) - \frac{Δ t}{Δ z} d_{F D} v_{z} (t) \end{array} \end{array}

and

\begin{array}{l} {\begin{array}{l} B_{y} (t + Δ t) = B_{y} (t) - \frac{Δ t}{Δ z} d_{F D} E_{x} (t), \\ v_{z} (t + Δ t) = v_{z} (t) + Δ t \frac{ω_{c}}{ω_{p}} E_{z} (t), \\ - \frac{Δ t}{Δ z} \frac{ω_{c}^{2}}{4 ω_{0}^{2}} d_{F D}^{T} (E_{x}^{2} (t)) . \end{array} \end{array}

(119)

These may then be composed together as in Sec. IV F to obtain approximations of the full flow.

We run a simulation of the 1D ponderomotive Maxwell system in the strongly nonlinear regime with the same parameters as in Sec. V A 2. We use second order splitting and N = 600 as this is the same number of degrees of freedom as was used in the previous spectral element FEEC simulations. The time step is $d t \approx 0.033$ ⁠, see Fig. 4 for the results of this simulation. One can see that a qualitatively similar solution is obtained, but it does not possess the spurious oscillations of the previous simulations. This is because high order methods suffer from Runge's phenomenon when gradients are too large, and, in such situations, a larger number of cells is preferable to high order elements.

FIG. 4.

View large Download slide

Numerical tests for strongly nonlinear regime using the finite difference/lowest order FEM solver: $(ω_{p} / ω_{0}, ω_{c} / ω_{0}) = (1 / 2, - 1)$ ⁠. See Eq. (116) for the definition of the discrete Casimir invariants. The figures on the left and right are the errors in the conservation laws and a visualization of the solution, respectively.

VI. CONCLUSION

The objectives of this paper are twofold: to derive a self-consistent Hamiltonian model of the ponderomotive force and to apply the general discretization procedure to study this model. The simple structure of the 1D ponderomotive Maxwell bracket is due to the fact that the model linearizes the fluid equations about a quiescent, uniform density background. The bracket, thus, becomes extraordinarily simple with all of the nonlinearity in the model arising from the complexity of the Hamiltonian. The asymptotic nature of the 1D ponderomotive Maxwell system causes it to miss certain physically significant features of the true dynamics it approximates. Moreover, even in parameter regimes where it is valid, the validity holds only on a limited time horizon. A high resolution simulation of the physics of the full, non-asymptotic system was previously studied.²¹ One can see that the initial phase of the evolution of the 1D ponderomotive Maxwell system captures some of the relevant features. The oscillatory transverse electromagnetic field induces charge separation in the longitudinal direction exciting a longitudinal electric field. Moreover, the transverse Gaussian pulse becomes deformed as it propagates initially experiencing pulse steepening. However, the subsequent pulse elongation as well as the frequency modulation cannot be captured by our model as the asymptotic model omits phase dynamics.

We also neglected the thermodynamic pressure in the simulations. This is because this term would be $O (c_{S}^{2} / c^{2})$ ⁠, where c_S is the sound speed. In many cases, this term is asymptotically small, and its effects would be negligible over the short time horizon where this asymptotic regime remains valid. Including it would introduce sound waves to the model. These acoustic waves were eliminated in the simulations so that we could more clearly see the behavior of the nonlinear wave induced by the ponderomotive force. To include sound waves in the numerical method, one simply includes the appropriate pressure term on the righthand side of the continuity equation time step in Eq. (102) [or Eq. (119) for the finite difference implementation].

The numerical method is not without its own limitations. We saw that the methods using high order interpolation in space underperformed in the strongly nonlinear regime with respect to discretizations using an equivalent number of degrees of freedom with lower order polynomial interpolation (which may be interpreted as a low order finite difference method) because of Runge's phenomenon. Likewise, the broken-FEEC approximations typically required greater resolution than the conforming-FEEC approximations in order to suppress spurious oscillations. Despite these current limitations, we have reason to believe that the numerical techniques introduced in this paper might provide a foundation for future powerful methods for simulating electrodynamics in nonlinearly polarized media.

In order to avoid spurious oscillations in regions of large gradient while maintaining the efficiency of high order methods in smoother regions, several approaches might be taken. First, one might investigate adaptive mesh refinement in a structure preserving context. The problem of refinement in a structure preserving algorithm becomes that of finding a suitable projection operator from one mesh to another which does not significantly alter the Casimir invariants or the energy: i.e., each time the mesh is refined, the data are mapped to a new Hamiltonian system of a different dimension. The performance of the broken-FEEC method might be enhanced by investigating conforming projections other than the simplest averaging procedure used in this work. The FEEC approach is preferable to the finite differences because (1) it more easily accommodates general boundary conditions, (2) it more easily generalizes to higher dimensions and more complex geometries, and (3) it scales to higher order discretizations. While this flexibility comes at the cost of algorithmic complexity, it has greater potential to yield a high performance nonlinear electrodynamics solver. Furthermore, the broken-FEEC approach localizes the basis functions to each element thus allowing for a scalable, parallel implementation. Particularly notable is that this localizes the nonlinear solve, the largest bottleneck in each time step. As a final note, the approach taken in this paper to discretize Maxwell's equation in a nonlinearly polarized medium (the ponderomotive Maxwell system is such a model) might be adapted to apply to a broad class of models in nonlinear optics. Subsequent papers will investigate this further application of the results of this paper in detail.

ACKNOWLEDGMENTS

We gratefully acknowledge the support of U.S. Department of Energy under Contract No. DE-FG05-80ET-53088, the NSF Graduate Research Fellowship via No. DGE-1610403, and the Humboldt foundation. We would also like to thank the three anonymous referees for their attentive reading of this manuscript and for their insightful feedback.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

William Barham: Conceptualization (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (equal); Writing – original draft (lead); Writing – review & editing (equal). Yaman Güçlü: Investigation (supporting); Software (equal); Supervision (equal); Writing – review & editing (equal). Philip J. Morrison: Conceptualization (supporting); Investigation (supporting); Methodology (supporting); Supervision (equal); Writing – review & editing (equal). Eric Sonnendrücker: Investigation (supporting); Supervision (equal); Writing – review & editing (equal).

DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

APPENDIX A: THE EMERGENCE OF MACROSCOPIC FIELDS FROM THE LEGENDRE TRANSFORM IN THE ELECTROMAGNETIC LAGRANGIAN

In this appendix, we briefly consider the reason for the particular form of the nonlinear constitutive relations in the Hamiltonian formulation of electromagnetism. In Gaussian units, D and A are canonically conjugate up to the constant scaling factor

4 π c

⁠. We may define the electromagnetic “kinetic energy” and “potential energy” as

T = \frac{c^{2}}{8 π} \int | \dot{A} |^{2} d^{3} x and V = \frac{1}{8 π} \int | \nabla \times A |^{2} d^{3} x,

(A1)

respectively. The Lagrangian is then given by

L = T - V .

(A2)

We then define the canonical momentum via a Legendre transform:

D / 4 π c = δ L / δ \dot{A}

⁠, and the Hamiltonian is given by

\begin{matrix} H = \frac{1}{4 π c} \int D \cdot \dot{A} d^{3} x - L \\ = \frac{1}{4 π c} \int D \cdot \dot{A} d^{3} x - \frac{1}{8 π} \int (c^{2} | \dot{A} |^{2} - | \nabla \times A |^{2}) d^{3} x \end{matrix} .

(A3)

Letting

E = c^{- 1} \dot{A}

and

B = \nabla \times A

⁠, we find

H = \int D \cdot E d^{3} x - \frac{1}{8 π} \int (| E |^{2} - | B |^{2}) d^{3} x,

(A4)

which is the usual electromagnetic energy in a vacuum since

D = E

for this choice of T. Hence, we recover the vacuum Maxwell equations.

We now generalize to nonlinear media in which the refractive index is a perturbation from the identity. If we let

T = \frac{c^{2}}{8 π} \int | \dot{A} |^{2} d^{3} x - K,

(A5)

where K is an arbitrary functional of the fields, we find

\begin{matrix} H = \frac{1}{4 π c} \int D \cdot \dot{A} d^{3} x + K - \frac{1}{8 π} \int (c^{2} | \dot{A} |^{2} - | \nabla \times A |^{2}) d^{3} x \end{matrix},

(A6)

where, since

E = c^{- 1} \dot{A}

⁠,

\frac{D}{4 π c} = \frac{δ L}{δ \dot{A}} = \frac{c^{2} \dot{A}}{4 π} - \frac{δ K}{δ \dot{A}} = \frac{E}{4 π c} - \frac{1}{c} \frac{δ K}{δ E}

(A7)

and

\begin{matrix} H = \frac{1}{4 π} \int D \cdot E d^{3} x + K - \frac{1}{8 π} \int (| E |^{2} - | B |^{2}) d^{3} x \\ = K - 4 π \int \frac{δ K}{δ E} \cdot E d^{3} x + \frac{1}{8 π} \int (| E |^{2} + | B |^{2}) d^{3} x . \end{matrix}

(A8)

This explains the particular form of the electromagnetic energy in a perturbed nonlinear medium.

We can see that the (scaled) canonical bracket is simply

{F, G} = 4 π c \int (\frac{δ F}{δ D} \cdot \frac{δ G}{δ A} - \frac{δ G}{δ D} \cdot \frac{δ F}{δ A}) d^{3} x .

(A9)

Hence, in the variables

(D, B)

⁠, the bracket is

{F, G} = 4 π c \int (\frac{δ F}{δ D} \cdot \nabla \times \frac{δ G}{δ B} - \frac{δ G}{δ D} \cdot \nabla \times \frac{δ F}{δ B}) d^{3} x .

(A10)

Assuming the Lagrangian is hyperregular, because the Legendre transform which maps

\dot{A}

to D is involutive, we have that

δ H / δ D = {(4 π c)}^{- 1} \dot{A} = {(4 π)}^{- 1} E

⁠. Finally, we simply define

H = 4 π δ H / δ B = B + 4 π δ K / δ B

⁠. We, therefore, recover the macroscopic Maxwell equations.

Hence, an asymptotic procedure which introduces changes to the electromagnetic energy must likewise induce the appropriate compensating polarization and magnetization prescribed by the formalism discussed herein. Moreover, even in linear media, $(D, A)$ are the naturally canonically conjugate variables with D and E being related to each other in the same manner in which the momentum and velocity are related in classical mechanics. We further note that this framework facilitates building theories of nonlinear electrodynamics which are not perturbations of the linear theory.

APPENDIX B: 1D SPECTRAL ELEMENT FEEC METHOD

Here, we present a brief overview of the conforming/broken FEEC method^9,10 and urge the interested reader to consult the numerical analysis literature for the full details. The FEEC method we employ is a spectral element method which uses an interpolation/histopolation approach.²²

1. Construction of the FEM basis

In a FEEC method, different kinds of quantities are approximated with a different set of basis functions in order to better represent the physics at a discrete level. In one spatial dimension, there are only scalar fields; however, we may distinguish between two different varieties of scalars transformation rules. Suppose

F : \hat{Ω} \to Ω

is a diffeomorphism. A 0-form transforms as

ϕ (z) = ϕ (F (\hat{z})) = \hat{ϕ} (\hat{z}),

(B1)

where

\hat{ϕ} = ϕ ○ F

⁠. Alternatively, a 1-form transforms as

ρ (z) d z = \hat{ρ} (\hat{z}) d \hat{z} \Rightarrow \hat{ρ} (\hat{z}) = F' (\hat{z}) ρ (F (\hat{z})) .

(B2)

The vector space of 0-forms is defined as

V^{0} : = H^{1} Λ^{0} (Ω) : = {ϕ \in L^{2} Λ^{0} (Ω) : \partial_{z} ϕ d z \in L^{2} Λ^{1} (Ω)} .

(B3)

This is the space of square integrable functions with square integrable derivatives which transform like differential 0-forms. The vector space

V^{1} : = L^{2} Λ^{1} (Ω)

is the space of square integrable functions which transform as 1-forms. These spaces are the usual spaces H¹ and L² studied in functional analysis augmented with the transformation rules of 0-forms and 1-forms.

As alluded to before, the finite element spaces which approximate V⁰ and V¹ use different basis functions for reasons we will discuss subsequently. We denote the finite element spaces by

V_{h}^{0}

and

V_{h}^{1}

the basis functions (also called shape functions) are denoted

{Λ_{0, j}}_{j = 1}^{N_{0}}

and

{Λ_{1, j}}_{j = 1}^{N_{1}}

⁠. Hence, if

ϕ_{h} \in V_{h}^{0}

and

ρ_{h} \in V_{h}^{1}

⁠, we have

ϕ_{h} (z) = \sum_{j = 1}^{N_{0}} ϕ_{j} Λ_{0, j} and ρ_{h} (z) = \sum_{j = 1}^{N_{1}} ρ_{j} Λ_{1, j},

(B4)

where

ϕ \in C_{0}

and

ρ \in C_{1}

are the Galerkin coefficients in the finite element expansions and

C_{0} \subset R^{N_{0}}

and

C_{1} \subset R^{N_{1}}

are the vector spaces of Galerkin coefficients. We shall uniformly use upright/sans-serif notation to denote the coefficient variables.

The interpolation operators are simply the operators which take a vector of Galerkin coefficients and return the interpolated field in the finite element space

\begin{matrix} ϕ_{h} = I^{0} ϕ = \sum_{j = 1}^{N_{0}} ϕ_{j} Λ_{j, 0} (z) and ρ_{h} = I^{1} ρ = \sum_{j = 1}^{N_{1}} ρ_{j} Λ_{j, 1} (z) . \end{matrix}

(B5)

We define the degrees of freedom operators

σ^{0} : V^{0} \to C_{0}

and

σ^{1} : V^{1} \to C_{1}

such that

σ_{i}^{0} (Λ_{0, j}) = δ_{i j} and σ_{i}^{1} (Λ_{1, j}) = δ_{i j} .

(B6)

Hence, applying the degrees of freedom operator to a field returns the degrees of freedom, the Galerkin coefficients, of that field in the appropriate finite element space, e.g.,

σ^{1} (ρ_{h}) = ρ

⁠. Finally, we define the projection operators

Π^{0} : V^{0} \to V_{h}^{0}

and

Π^{1} : V^{1} \to V_{h}^{1}

to be the composition of the degrees of freedom and interpolation operators

Π^{0} ϕ = \sum_{j = 1}^{N_{0}} σ_{j}^{0} (ϕ) Λ_{j, 0} (z) and Π^{1} ρ = \sum_{j = 1}^{N_{1}} σ_{j}^{1} (ρ) Λ_{j, 1} (z) .

(B7)

We may concisely summarize the spaces and operators making up our 1D FEEC method by the following commuting diagram:

.

(B8)

We define

d_{0}

below in definition 3. Commutativity of the diagram means

d_{0} σ^{0} = σ^{1} \partial_{z}, \partial_{z} I^{0} = I^{1} d_{0}, and \partial_{z, h} Π^{0} = Π^{1} \partial_{z} .

(B9)

The reason why we utilize this FEEC approach is to ensure the consistent treatment of derivatives in our discretized models

ϕ_{h} \in V_{h}^{0} \Rightarrow \partial_{z} ϕ_{h} \in V_{h}^{1}

(B10)

and

σ^{0} (ϕ_{h}) = ϕ \in C_{0} \Rightarrow σ^{1} (\partial_{z} ϕ_{h}) = d_{0} ϕ \in C_{1} .

(B11)

The framework just described gives us these desirable properties. The remainder of this section provides a brief overview of the definition of the finite element de Rham complex in one spatial dimension.

To begin, let $Ω = [0, L]$ be the physical domain. We divide this into K cells which we denote ${Ω_{k}}_{k = 1}^{K}$ such that $\underset{k}{\cup} Ω_{k} = Ω$ ⁠. We first describe the construction on the reference element $\hat{Ω} : = [- 1, 1]$ which is then mapped to each Ω_k by a diffeomorphism $F_{k} : \hat{Ω} \to Ω_{k}$ ⁠. We shall denote all operators defined on the reference element with a hat to distinguish these from their mapped counterparts.

The finite element shape functions used to interpolate the 0-forms are simply Lagrange interpolating polynomials with the Gauss-Lobatto quadrature points as nodes.

Definition 1. We denote the Gauss-Lobatto points in

[- 1, 1]

by ${ξ_{i}}_{i = 0}^{N}$ and we denote the Lagrange interpolating polynomials over this grid by ${\hat{Λ}}_{0, i}$ ⁠. We let there be N + 1 quadrature points so that the order of the interpolating polynomials be N. The 0-form degrees of freedom operators are defined to be

{\hat{σ}}_{0, i} (\hat{v}) = \hat{v} (ξ_{i}) \Rightarrow {\hat{σ}}_{0, i} ({\hat{Λ}}_{0, j}) = δ_{i j}

(B12)

and 0-form interpolation is simply

{\hat{I}}^{0} v = \sum_{i} v_{i} {\hat{Λ}}_{0, i}

⁠.

Definition 2. The degrees of freedom for the 1-forms are cell-wise integrals

{\hat{σ}}_{1, i} (v) = \int_{ξ_{i}}^{ξ_{i + 1}} v (ξ) d ξ for i = 0, 2, \dots, N - 1.

(B13)

The 1-form interpolating polynomials are

\begin{matrix} {\hat{Λ}}_{1, i} (ξ) = - \sum_{k = 0}^{i} {\hat{Λ}}_{0, k}' (ξ) = \sum_{k = i + 1}^{N} {\hat{Λ}}_{0, k}' (ξ) \\ = \frac{1}{2} [\sum_{k = i + 1}^{N} {\hat{Λ}}_{0, k}' (ξ) - \sum_{k = 0}^{i} {\hat{Λ}}_{0, k}' (ξ)], \end{matrix}

(B14)

where the equality of these three expressions follows from the fact that

\sum_{i = 0}^{N} {\hat{Λ}}_{0, i} (ξ) = 1 \Rightarrow \sum_{i = 0}^{N} {\hat{Λ}}_{0, i}' (ξ) = 0

(B15)

since the Lagrange interpolating polynomials form a partition of unity.

The 1-form basis functions were selected such that the following is true:

\begin{matrix} {\hat{σ}}_{1, i} ({\hat{Λ}}_{1, j}) & = \int_{ξ_{i}}^{ξ_{i + 1}} {\hat{Λ}}_{1, j} (ξ) d ξ \\ = \sum_{k = j + 1}^{N} [δ_{i + 1, k} - δ_{i, k}] \\ = δ_{i, j} + \sum_{k = j + 1}^{N - 1} (δ_{i, k} - δ_{i, k}) - δ_{i, N} \\ = δ_{i, j} - δ_{i, N}, \end{matrix}

(B16)

where we made use of the identity

δ_{i + c, k + c} = δ_{i, k} \forall (i, k, c) \in Z^{3}

⁠. On a grid with N elements, we have that

i = 0, \dots, N - 1

⁠, and therefore,

δ_{i, N} = 0

always. This yields the final result

{\hat{σ}}_{1, i} ({\hat{Λ}}_{1, j}) = δ_{i, j} .

(B17)

Definition 3. We define the discrete exterior derivative matrix to be

{(d_{0})}_{i j} = {\hat{σ}}_{1, i} (\partial_{ξ} {\hat{Λ}}_{0, j} (ξ)) .

(B18)

One may show that

{\hat{σ}}_{1, i} (\partial_{ξ} v) = \int_{ξ_{i}}^{ξ_{i + 1}} \partial_{ξ} v d ξ = v (ξ_{i + 1}) - v (ξ_{i}) = {(d_{0})}_{i j} {\hat{σ}}_{0, j} (v) .

(B19)

Hence, we find that the derivative matrix is simply the adjacency matrix of the grid

d_{0} : {\hat{C}}_{0} \to {\hat{C}}_{1} such that {(d_{0})}_{i j} = {\begin{array}{l} 1, & j = i + 1 \\ - 1, & j = i \\ 0, & else . \end{array}

(B20)

It is further possible to show

\partial_{ξ, h} {\hat{Π}}^{0} = {\hat{Π}}^{1} \partial_{ξ}

⁠, which establishes the commutativity of the diagram.

This defines all of the discrete operators on the reference element. We then define maps to each element in the physical domain, $F_{k} : \hat{Ω} \to Ω_{k}$ ⁠, and pullback all quantities to the reference element in order to perform all computation. The manner by which 0-forms and 1-forms transform differ.

Definition 4. We define the push-forward operators

\begin{matrix} F_{k}^{0} : \hat{u} \mapsto u : = \hat{u} ○ F_{k}^{- 1} \end{matrix}

and

\begin{matrix} F_{k}^{1} : \hat{u} \mapsto u : = \frac{1}{\partial_{\hat{ξ}} F_{k} (\hat{ξ})} \hat{u} ○ F_{k}^{- 1} . \end{matrix}

(B21)

We see that 0-forms transform like scalars while 1-forms transform like densities.

Proposition 1. The derivative commutes with the push-forward

\partial_{x} F_{k}^{0} \hat{u} = F_{k}^{1} \partial_{ξ} \hat{u} \forall \hat{u} = \hat{u} (\hat{ξ}) .

(B22)

Proof: This is simply the chain rule. □

Definition 5. The local mapped finite element spaces are simply

V_{h, k}^{ℓ} = F_{k}^{ℓ} {\hat{V}}_{h}^{ℓ}

. We define the global finite element space by

V_{h}^{ℓ} = \sum_{k} V_{h, k}^{ℓ},

(B23)

where we have extended functions outside of their local elements by zero. We then define the global degrees of freedom and global interpolation operators on the mapped domain by

σ_{ℓ, i, k} (u) : = {\hat{σ}}_{ℓ, i} ({(F_{k}^{ℓ})}^{- 1} u) and Λ_{ℓ, i, k} = F_{k}^{ℓ} {\hat{Λ}}_{ℓ, i} .

(B24)

That the commuting diagram property is satisfied on the mapped domain follows from simple algebraic manipulation. When computing the degrees of freedom or inner products with the shape functions, one first pulls back to the reference element.

Finally, in this 1D context, the distinction between conforming- and broken-FEEC methods is especially simple. Neighboring mapped elements will have their endpoints overlapping. That is, if we have $F_{1} : [- 1, 1] \to [0, 1 / 2]$ and $F_{2} : [- 1, 1] \to [1 / 2, 1]$ map to two distinct elements, $F_{1} (1) = F_{2} (- 1) = 1 / 2$ ⁠. In a conforming FEEC method, we collapse the two overlapping endpoints into one single degree of freedom. In broken-FEEC, we keep these degrees of freedom distinct to allow for discontinuity across element boundaries. We then define a conforming projection operator which maps from the discontinuous space into the continuous space in order to perform derivatives. A simple choice is to average the degrees of freedom on the element boundaries. From the point of view of one using a broken-FEEC method, the primary difference then is that the discrete derivative operator must be composed with this (sparse) projection operator.

The presence of adjoint differential operators in our Hamiltonian formulation of the 1D ponderomotive Maxwell model requires that we consider a dual de Rham complex (the complex for the adjoint differential operators) in our discretization scheme. Once the primal de Rham complex has been defined (as above), it is straightforward to define the dual complex in a subsidiary manner entirely in terms of quantities previously defined.

Definition 6. The dual degrees of freedom

σ_{ℓ}^{*} : {(V^{ℓ})}^{*} \to C_{ℓ}^{*}

are defined such that

σ_{ℓ, i}^{*} (u) = ⟨ u, Λ_{ℓ, i} ⟩ for u \in {(V^{ℓ})}^{*},

(B25)

where

⟨ \cdot, \cdot ⟩ : {(V^{ℓ})}^{*} \times V^{ℓ} \to R

denotes the evaluation pairing of a vector space with its dual. We define the dual basis,

{Λ_{ℓ, i}^{*}}_{i = 0}^{N_{ℓ}}

, such that

(Λ_{ℓ, i}, Λ_{ℓ, j}^{*}) = δ_{i, j}

. The dual interpolation operator

I_{ℓ}^{*} : C_{ℓ}^{*} \to {(V_{h}^{ℓ})}^{*}

is defined such that

I_{ℓ}^{*} c_{*} = \sum_{i = 0}^{N_{ℓ}} {(c_{*})}_{i} (Λ_{ℓ, i}^{*}, \cdot) \in {(V_{h}^{ℓ})}^{*} .

(B26)

That is,

I_{ℓ}^{*} c_{*}

is a functional which acts on elements of

V_{h}^{ℓ}

such that

I_{ℓ}^{*} c_{*} [v_{h}] = \sum_{i = 0}^{N_{ℓ}} {(c_{*})}_{i} (Λ_{ℓ, i}^{*}, v_{h}) .

(B27)

The mass matrix is defined as ${(M_{ℓ})}_{i, j} = (Λ_{ℓ, i}, Λ_{ℓ, j}) : C_{ℓ} \to C_{ℓ}^{*}$ ⁠. The dual projection operator is defined such that

Π_{ℓ}^{*} = I_{ℓ}^{*} ° σ_{ℓ}^{*} : {(V^{ℓ})}^{*} \to {(V_{h}^{ℓ})}^{*} .

(B28)

One can see that this is simply the L² projection.

Proposition 2. If

c_{*} \in C_{ℓ}^{*}

and

v_{h} \in V_{h}^{ℓ}

such that

σ^{ℓ} (v_{h}) = v \in C_{ℓ}

, then

I_{ℓ}^{*} c_{*} [v_{h}] = I_{ℓ}^{*} c_{*} [I^{ℓ} v] = c_{*}^{T} v .

(B29)

Proof: This simply follows from the fact that

(Λ_{ℓ, i}, Λ_{ℓ, j}^{*}) = δ_{i, j}

⁠. □

Proposition 3. The mass matrix provides the change of basis from

C_{ℓ}

to

C_{ℓ}^{*}

. That is, if

u \in V_{h}^{ℓ}

and we denote its corresponding element in the dual space

u_{h}^{*} = (u_{h}, \cdot) \in {(V_{h}^{ℓ})}^{*}

, then

M_{ℓ} σ_{ℓ} (u_{h}) = σ_{ℓ}^{*} (u_{h}^{*}) .

(B30)

Proof:

\begin{matrix} σ_{ℓ}^{*} (u_{h}^{*}) = ⟨ u_{h}^{*}, Λ_{ℓ, i} ⟩ = (u_{h}, Λ_{ℓ, i}) \\ = \sum_{j} (Λ_{ℓ, i}, Λ_{ℓ, j}) u_{j} = M_{ℓ} σ_{ℓ} (u_{h}) . \end{matrix}

(B31)

□

The primal and dual de Rham complexes together are frequently called the double de Rham complex. The discrete double de Rham complex may be schematically represented by

(B32)

where

R_{ℓ} : V_{h}^{ℓ} \to {(V_{h}^{ℓ})}^{*}

is the isomorphism guaranteed by the Reisz representation theorem, or, at the coefficient level,

(B33)

where each diagram commutes. We can see that the discrete integration by parts formula is given by

{(d_{0}^{*})}^{T} M_{0} = M_{1} d_{0} .

(B34)

Finally, it is worth briefly discussing the advantages of broken-FEEC. The advantage of broken-FEEC is that the basis functions are localized to each element. This gives the mass matrix a block diagonal structure and also localizes the nonlinear solves required by this method. Hence, at the cost of increasing the number of degrees of freedom in a broken-FEEC method, one localizes many operations allowing for greater parallelism. The one-dimensional problem considered in this paper does not require such fancy computational techniques; however, similar problems in nonlinear electrodynamics in higher dimensions might benefit from these considerations. By introducing conforming projection operators^9,10 which project fields from the broken spaces into appropriate spaces in the conforming finite element de Rham complex, one can construct a commuting diagram of vector spaces like Eq. (B8) (and generalizations thereof in 2D and 3D). Hence, this version of FEEC which utilizes broken elements retains the structure preserving properties celebrated in conforming FEEC discretizations.

2. Integrals of functions in FEM space

It is useful to know how to compute integrals of functions in the FEM spaces over the entire domain. The ponderomotive Maxwell system possesses several conservation laws which are integrals of fields over the entire domain. The following lemma tells us how such integrals may be computed. In the following, we combine the indices for the DOF within an element, and the element index into one single multiindex.

Lemma 1. Suppose

v_{h} \in V_{h}^{0}

and

u_{h} \in V_{h}^{1}

. Then if we denote

v_{h} = \sum_{α} Λ_{0, α} v_{α}

and

u_{h} = \sum_{α} Λ_{1, α} u_{α}

, it follows that:

\begin{matrix} \int_{0}^{L} v_{h} d z = \sum_{α, β} {(M_{0})}_{αβ} v_{β} = 1^{T} M_{0} v \end{matrix}

and

\begin{matrix} \int_{0}^{L} u_{h} d z = \sum_{α} u_{α} = 1^{T} u, \end{matrix}

(B35)

where

1_{α} = 1 \forall α

⁠.

Proof: To begin, we consider the integral of 1-forms in

V_{h}^{1}

⁠. The proof proceeds by direct calculation and application of the definitions

\begin{matrix} \int_{0}^{L} u_{h} d z & = \sum_{i, k} \int_{Ω_{k}} Λ_{1, i, k} (z) u_{i, k} d z \\ = \sum_{i, k} (\int_{\hat{Ω}} {\hat{Λ}}_{1, i} (ξ) d ξ) u_{i, k} \\ = \sum_{i, j, k} (\int_{ξ_{j}}^{ξ_{j + 1}} {\hat{Λ}}_{1, i} (ξ) d ξ) u_{i, k} \\ = \sum_{i, j, k} δ_{i j} u_{i, k} = \sum_{i, k} u_{i, k} = \sum_{α} u_{α} \end{matrix},

(B36)

where we have used the fact that

Λ_{1, i, k} (z) d z = {\hat{Λ}}_{1, i} (ξ) d ξ

where

F_{k} (ξ) = z

and

F_{k} (\hat{Ω}) = Ω_{k}

⁠.

We proceed similarly with the 0-forms

\begin{matrix} \int_{0}^{L} v_{h} d z & = \sum_{i, k} \int_{Ω_{k}} Λ_{0, i, k} (z) v_{i, k} d z \\ = \sum_{i, k} (\int_{\hat{Ω}} {\hat{Λ}}_{0, i} (ξ) \partial_{ξ} F_{k} (ξ) d ξ) v_{i, k} . \end{matrix}

(B37)

Recall, the 0-form shape functions are a partition of unity:

\sum_{l} {\hat{Λ}}_{0, l} (ξ) = 1

⁠. Hence,

\begin{matrix} \int_{\hat{Ω}} {\hat{Λ}}_{0, i} (ξ) \partial_{ξ} F_{k} (ξ) d ξ = \sum_{l} \int_{\hat{Ω}} {\hat{Λ}}_{0, i} (ξ) {\hat{Λ}}_{0, l} (ξ) \partial_{ξ} F_{k} (ξ) d ξ \\ = \sum_{l} {(M_{0}^{(k)})}_{l i}, \end{matrix}

(B38)

where

M_{0}^{(k)}

represents the 0-form mass matrix for the kth element. The global mass matrix

M_{0}

is related to each local mass matrix according to

\sum_{i} {(M_{0}^{(k)})}_{l i} v_{i, k} = {(M_{0} v)}_{l, k},

(B39)

where care must be taken in how the degrees of freedom at element boundaries are handled (this differs between conforming and broken FEEC methods). Hence, it follows that:

\int_{0}^{L} v_{h} d z = \sum_{i, k} {(M_{0} v)}_{i, k} = 1^{T} M_{0} v,

(B40)

which establishes the result. □

APPENDIX C: ON DISCRETIZING FUNCTIONAL DERIVATIVES

This appendix considers the mathematically rigorous details in discretizing functional derivatives. Much of the content in this appendix assumes a background in functional analysis and may be skipped if one is only interested in the results of the paper and not the mathematical details.

1. Discrete functional derivatives

As Hamiltonian field theories are stated in terms of functional derivatives, in order to discretize the Hamiltonian structure directly, it is necessary to consider what happens when a functional is restricted to a finite dimensional subspace (e.g., a finite element subspace). What we find is a convenient approach to discretize the Poisson bracket and Hamiltonian.

Let

(U, | | \cdot | |_{U})

be a normed vector space. For any functional

K : U \to R

⁠, we define its Fréchet derivative at

u \in U

(if it exists) to be the linear functional

D K [u]

such that

\frac{| | K [u + δ u] - K [u] - D K [u] δ u | |_{U}}{| | δ u | |_{U}} = O (| | δ u | |_{U}) .

(C1)

We shall assume in this paper that we restrict our attention to those functionals which are Fréchet differentiable. When

U = R^{n}

⁠, this reduces to the usual definition of the Jacobian matrix. Hence, this definition generalizes the notion of derivatives to arbitrary normed vector spaces (e.g., function spaces). We shall use the two common notations for the functional derivative

D K [u]

and

δ K / δ u

interchangeably.

Derivatives with respect to functions which are Galerkin expansions in our finite element framework reduce to finite dimensional derivatives with respect to the Galerkin coefficients.

Theorem 1. Let

K : V^{ℓ} \to R

be an arbitrary functional on the continuous

ℓ

-forms and let

K : = K ○ I^{ℓ} : C_{ℓ} \to R

represent the discrete analog of the functional K. Moreover, define

u = σ^{ℓ} (u)

and

v = σ^{ℓ} (v)

. Then

\begin{matrix} {⟨ D (K ○ Π^{ℓ}) [u], v ⟩}_{{(V^{ℓ})}^{*}, V^{ℓ}} \\ = {⟨ D (K ○ I^{ℓ}) [σ^{ℓ} (u)], σ^{ℓ} (v) ⟩}_{C_{ℓ}^{*}, C_{ℓ}} \\ = {⟨ D K [u], v ⟩}_{C_{ℓ}^{*}, C_{ℓ}} = {(\frac{\partial K}{\partial u})}^{T} v . \end{matrix}

(C2)

Similarly, if we have a functional of the dual de Rham sequence,

K_{*} : {(V^{ℓ})}^{*} \to R

, defining

K_{*} = K ○ I_{ℓ}^{*}, u_{*} = σ_{ℓ}^{*} (u)

, and

v_{*} = σ_{ℓ}^{*} (v)

, we have

\begin{matrix} {⟨ D (K_{*} ○ Π_{ℓ}^{*}) [u], v ⟩}_{V^{ℓ}, {(V^{ℓ})}^{*}} \\ = {⟨ D (K_{*} ○ I_{ℓ}^{*}) [σ_{ℓ}^{*} (u)], σ_{ℓ}^{*} (v) ⟩}_{C_{ℓ}, C_{ℓ}^{*}} \\ = {⟨ D K_{*} [u_{*}], v_{*} ⟩}_{C_{ℓ}, C_{ℓ}^{*}} = {(\frac{\partial K_{*}}{\partial u_{*}})}^{T} v_{*} . \end{matrix}

(C3)

Proof: The result follows from the chain rule and the linearity of the degrees of freedom operator. By the chain rule, we have

\begin{matrix} {⟨ D (K ○ Π^{ℓ}) [u], v ⟩}_{{(V^{ℓ})}^{*}, V^{ℓ}} \\ = {⟨ D (K ○ I^{ℓ}) [σ^{ℓ} (u)] ○ D σ^{ℓ} [u], v ⟩}_{{(V^{ℓ})}^{*}, V^{ℓ}} . \end{matrix}

(C4)

However, by the linearity of

σ^{ℓ}, D σ^{ℓ} [u] = σ^{ℓ}

⁠. Hence, we find

\begin{matrix} {⟨ D (K ○ Π^{ℓ}) [u], v ⟩}_{{(V^{ℓ})}^{*}, V^{ℓ}} = {⟨ D (K ○ I^{ℓ}) [σ^{ℓ} (u)], σ^{ℓ} (v) ⟩}_{C_{ℓ}^{*}, C_{ℓ}} . \end{matrix}

(C5)

The remainder of the result follows from simple notational manipulation. The result for

K_{*}

follows from an entirely analogous argument. □

Hence, one may directly translate continuous functional derivatives into discrete functional derivatives. Moreover, by appropriate interpretation of the discrete representation of the variables on the dual sequence, we may also prescribe discrete versions of functional derivatives with respect to variables in the dual space.

2. Adjoint of degrees of freedom operator

In this section, we prove a technical result which is needed to discretize the Poisson bracket. This may safely be skipped if one is only interested in the results and not the mathematical substructure of the numerical method. In the following, we use $†$ to denote operator adjoints.

Lemma 2. Suppose

u_{*} \in C_{ℓ}^{*}

and

v \in C^{ℓ}

. Then

{⟨ u_{*}, v ⟩}_{C_{ℓ}^{*}, C^{ℓ}} = {⟨ I_{ℓ}^{*} u_{*}, I^{ℓ} v ⟩}_{{(V^{ℓ})}^{*}, V^{ℓ}} .

(C6)

Proof: The result follows by simple application of the definitions:

\begin{matrix} {⟨ I_{ℓ}^{*} u_{*}, I^{ℓ} v ⟩}_{{(V^{ℓ})}^{*}, V^{ℓ}} \\ = {⟨ \sum_{i = 1}^{N_{ℓ}} {(u_{*})}_{i} (Λ_{ℓ, i}^{*}, \cdot), \sum_{j = 1}^{N_{ℓ}} v_{j} Λ_{ℓ, j} ⟩}_{{(V^{ℓ})}^{*}, V^{ℓ}} \\ = \sum_{i, j = 1}^{N_{ℓ}} {(u_{*})}_{i} v_{j} (Λ_{ℓ, i}^{*}, Λ_{ℓ, j}) \\ = u_{*}^{T} v : = {⟨ u_{*}, v ⟩}_{C_{ℓ}^{*}, C^{ℓ}} . \end{matrix}

(C7)

□

Theorem 2. With respect to the natural evaluation duality pairing,

{({σ^{ℓ} |}_{V^{ℓ}})}^{†} = I_{ℓ}^{*}

. Moreover,

\begin{matrix} \frac{| ([{(σ^{ℓ})}^{†} - I_{ℓ}^{*}] u_{*}, v) |}{| | u_{*} | |_{C_{ℓ}^{*}} | | v | |_{V^{ℓ}}} \leq | | I - Π_{ℓ} | | | | I_{ℓ}^{*} | | \\ \forall u \in C_{ℓ}^{*} and \forall v \in V^{ℓ} . \end{matrix}

(C8)

Hence,

I_{ℓ}^{*}

approximates

σ_{ℓ}^{†}

in the following sense:

\begin{matrix} | | {(σ^{ℓ})}^{†} - I_{ℓ}^{*} | | : = sup_{\begin{matrix} | | u | |_{C_{ℓ}^{*}} \leq 1 \\ u_{*} \in C_{ℓ}^{*} \end{matrix}} sup_{} \end{matrix}