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.^{1} 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,^{2} 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.^{3} 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,^{4} and a structure preserving discretization is easily accomplished using finite element exterior calculus (FEEC).^{5} 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.^{13} 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

*m*and

*q*are the particles' mass and charge, respectively, and

*ω*

_{0}and

**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**

*E**q*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.

_{e}### A. The two-fluid system: Longitudinal/transverse split

### B. Laser–plasma interaction—1D linear regime

^{1}in laser-plasma literature. In this regime, we assume the transverse electromagnetic fields scale as $ E \u22a5 , B \u22a5 \u223c \u03f5 \u226a 1$, the longitudinal electromagnetic fields scale as $ E z , B z \u223c \u03f5 2$, gradients are small in the transverse direction $ \u2207 \u22a5 \u223c \u03f5 2$ while $ \u2202 z \u223c O ( 1 )$, and we expand the fluid variables in powers of

*ϵ*: e.g., $ n \alpha = \u2211 k n \alpha ( k )$ where $ n \alpha ( k ) \u223c \u03f5 k$. We further assume that all transverse fields are such that

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

### A. Deriving the Hamiltonian

*θ*, the kinetic energy becomes

*θ*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 \u22a5 = E 0 \u22a5 \u2009 sin ( \theta )$ and $ A \u22a5 = A 0 \u22a5 \u2009 sin ( \theta )$, where $ A \u22a5$ is the vector potential for the transverse field and $ \theta \u0307 = \omega 0 = const$, then we find that

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.^{15}

### B. Derivatives of the Hamiltonian

^{3}we seek a functional

*K*such that

^{4}See Eq. (C1) for a definition of the functional derivative. Small modifications must be made to the Hamiltonian modeling framework for nonlinear polarization

^{3}due to the additional factor of 1/2 appearing in the electromagnetic energy from the averaging procedure.

*K*satisfying Eq. (27) is

*L*

^{2}adjoint of an operator. Similarly, one finds that

### 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.^{15}

^{14}The bracket given therein uses momentum coordinates, $ M \alpha = \rho \alpha v \alpha $. One may change variables to formulate the model in terms of velocity

^{4}to obtain the Poisson bracket

^{3}which arises in the reduced model. Because in this parent two-fluid model, $ K = K [ n \alpha , v \alpha ]$, it follows that $ D = E$. In passing to the reduced model, this will no longer be the case.

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.^{4}

### D. 1D ponderomotive Maxwell equations of motion

*F*, $ F \u0307 = { F , H}$. We find that

*x*-axis with the electric field and the

*y*-axis with the magnetic field, we obtain

^{2}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

*n*

_{0}. Suppose we nondimensionalize the fields such that

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: $ ( \omega p / \omega 0 , \omega c / \omega 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.

## 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.

The interpolation operators $ I \u2113 : C \u2113 \u2192 V h \u2113$ are simply those operators which take a given vector of Galerkin coefficients and returns the interpolated field: e.g., $ I 0 \varphi = \varphi h$. We also define degrees of freedom operators $ \sigma \u2113 : V \u2113 \u2192 C \u2113$ which simply take a field and return the appropriate degrees of freedom associated with that field, e.g., $ \sigma 0 ( \varphi h ) = \varphi $.

Finally, without going too much into the technical details, there are also dual spaces, $ ( V h \u2113 ) *$, dual interpolation, $ I \u2113 *$, and dual degrees of freedom operators, $ \sigma \u2113 *$, 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 \u2113 *$. Roughly speaking, $ \sigma \u2113 * : ( V \u2113 ) * \u2192 C \u2113 *$ are *L*^{2}-projections. Any equation which contains an adjoint derivative must be discretized by applying such an *L*^{2}-projection in order to integrate by parts and discretize the adjoint derivative weakly.

### 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 $ U \xaf ( 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.

*D*and

_{x}*B*are dual to each other,

_{y}*n*and

*v*are dual to each other, and

_{z}*E*and

_{z}*v*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 (

_{z}*D*,

_{x}*v*) are on the dual sequence. We chose to model

_{z}*v*on the dual sequence so that the derivative on the electric field squared might be weakened via integration by parts. Hence,

_{z}*L*

^{2}duality map. A subscript “

*h*” denotes a field's FEM representation: e.g., $ A \u2208 V \u2113 \u21d2 A h = \Pi \u2113 A$.

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

^{4}

### C. Derivatives of the discrete Hamiltonian

### D. Spatially discrete equations of motion

### E. Discrete conservation laws

### F. Temporal discretization via Hamiltonian splitting

^{16}Strang splitting

^{17}is a popular symmetric second order method

^{18}The derivation of symplectic integrators based on splitting methods for exponential maps is a well-established theory.

^{19}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^{20} 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 $ I \u0303$ via an *L*^{2} 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.

$ f \u0303 \u2190 gaussQuadrature ( ( \Lambda \u2113 , i , f ) L 2 , N Q )$ ▹ Quadrature with N _{Q} points |

return $ f \u0303$ |

$ f \u0303 \u2190 gaussQuadrature ( ( \Lambda \u2113 , i , f ) L 2 , N Q )$ ▹ Quadrature with N _{Q} points |

return $ f \u0303$ |

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 $ I \u0303$ is fairly trivial and given in Algorithm 3.

Input: $ D \u0303 x , E x , n$ |

Output: $ E x$ |

field $ ( x ; \u2009 E x , n ) : = \omega p 2 \omega 0 2 [ 1 + \u2211 i n i \Lambda 1 , i ( x )$ |

$ \u2212 \omega c 2 8 \omega 0 2 | \u2211 i ( E x ) i \Lambda 0 , i ( x ) | 2 ] \u2211 i ( E x ) i \Lambda 0 , i ( x ) D x \u2190 M 1 \u2212 1 D \u0303 x res ( E x , P x ) : = E x \u2212 P x \u2212 D x$ |

do |

$ P \u0303 x \u2190 assembl e 1 ( field ( x ; \u2009 E x , n ) , N Q ) P x \u2190 M 1 \u2212 1 P \u0303 x R \u2190 res ( E x , P x ) E x \u2190 E x \u2212 R$ |

while ( $ | R i | < RTOL | D i | + ATOL , \u2009 \u2200 i$) |

Input: $ D \u0303 x , E x , n$ |

Output: $ E x$ |

field $ ( x ; \u2009 E x , n ) : = \omega p 2 \omega 0 2 [ 1 + \u2211 i n i \Lambda 1 , i ( x )$ |

$ \u2212 \omega c 2 8 \omega 0 2 | \u2211 i ( E x ) i \Lambda 0 , i ( x ) | 2 ] \u2211 i ( E x ) i \Lambda 0 , i ( x ) D x \u2190 M 1 \u2212 1 D \u0303 x res ( E x , P x ) : = E x \u2212 P x \u2212 D x$ |

do |

$ P \u0303 x \u2190 assembl e 1 ( field ( x ; \u2009 E x , n ) , N Q ) P x \u2190 M 1 \u2212 1 P \u0303 x R \u2190 res ( E x , P x ) E x \u2190 E x \u2212 R$ |

while ( $ | R i | < RTOL | D i | + ATOL , \u2009 \u2200 i$) |

Input: $ E x$ |

Output: $ I \u0303$ |

$ field ( x ; \u2009 E x ) : = | \u2211 i ( E x ) i \Lambda 0 , i ( x ) | 2 I \u0303 \u2190 assembl e 1 ( field ( x ; \u2009 E x ) , N Q )$ |

Input: $ E x$ |

Output: $ I \u0303$ |

$ field ( x ; \u2009 E x ) : = | \u2211 i ( E x ) i \Lambda 0 , i ( x ) | 2 I \u0303 \u2190 assembl e 1 ( field ( x ; \u2009 E x ) , N Q )$ |

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.^{18} 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 $ \Delta t / \Delta x < 1$ where $ \Delta x$ is the minimum Gauss-Lobatto grid-spacing among all of the mapped elements. In practice, we let $ \Delta t / \Delta x = 1 / 2$.

Input: stateVec |

Output: $ v \u0303 z , B y$ |

$ v \u0303 z \u2190 v \u0303 z + \Delta t ( \u2212 \omega c 2 4 \omega 0 2 d 0 T I \u0303 + \omega c \omega p M 0 E z ) B y \u2190 B y \u2212 \Delta t d 0 E x$ |

Input: stateVec |

Output: $ v \u0303 z , B y$ |

$ v \u0303 z \u2190 v \u0303 z + \Delta t ( \u2212 \omega c 2 4 \omega 0 2 d 0 T I \u0303 + \omega c \omega p M 0 E z ) B y \u2190 B y \u2212 \Delta t d 0 E x$ |

Input: stateVec |

Output: $ n , E z , D \u0303 x$ |

$ n \u2190 n \u2212 \Delta t d 0 M 0 \u2212 1 v \u0303 z E z \u2190 E z \u2212 \omega p \omega c \Delta t \u2009 M 0 \u2212 1 v \u0303 z D \u0303 x \u2190 D \u0303 x + \Delta t \u2009 d 0 T M 1 B y$ |

picardSolveE( $ stateVec ; \u2009 ATOL , RTOL , N Q$) |

computeIntensity( $ stateVec ; \u2009 N Q$) |

Input: stateVec |

Output: $ n , E z , D \u0303 x$ |

$ n \u2190 n \u2212 \Delta t d 0 M 0 \u2212 1 v \u0303 z E z \u2190 E z \u2212 \omega p \omega c \Delta t \u2009 M 0 \u2212 1 v \u0303 z D \u0303 x \u2190 D \u0303 x + \Delta t \u2009 d 0 T M 1 B y$ |

picardSolveE( $ stateVec ; \u2009 ATOL , RTOL , N Q$) |

computeIntensity( $ stateVec ; \u2009 N Q$) |

Input: stateVec |

Output: stateVec |

updateBV( $ \Delta t / 2$; stateVec) |

updateEN( $ \Delta t ; \u2009 stateVec , \u2009 ATOL , \u2009 RTOL , N Q$) |

updateBV( $ \Delta t / 2$; stateVec) |

Input: stateVec |

Output: stateVec |

updateBV( $ \Delta t / 2$; stateVec) |

updateEN( $ \Delta t ; \u2009 stateVec , \u2009 ATOL , \u2009 RTOL , N Q$) |

updateBV( $ \Delta t / 2$; stateVec) |

## V. NUMERICAL RESULTS FOR 1D PONDEROMOTIVE MAXWELL SYSTEM

*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 / \u03f5$ where

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

### 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 $ ( \omega p / \omega 0 , \omega c / \omega 0 ) = ( 1 / 5 , \u2212 1 / 5 )$, and a strongly nonlinear regime in which $ ( \omega p / \omega 0 , \omega c / \omega 0 ) = ( 2 / 5 , \u2212 9 / 10 )$. The two dimensionless parameters independently control the strength of the two nonlinear features of the model: $ \omega c / \omega 0$ controls the strength of the ponderomotive force, and $ \omega p / \omega 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 $ \Delta t \u2248 0.028$ for the weakly nonlinear case, and $ \Delta t \u2248 0.014$ for the strongly nonlinear case. The high spatial order test cases used $ \Delta t \u2248 0.017$ in the weakly nonlinear case, and $ \Delta t \u2248 0.0085$ in the strongly nonlinear case. In all cases, the nonlinear solver tolerance has been set to $ 10 \u2212 13$.

#### 1. Weakly nonlinear regime

We first consider a weakly nonlinear regime in which $ ( \omega p / \omega 0 , \omega c / \omega 0 ) = ( 1 / 5 , \u2212 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.^{19} 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 $ \u223c 10 \u2212 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.

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 $ ( \omega p / \omega 0 , \omega c / \omega 0 ) = ( 2 / 5 , \u2212 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 \u2212 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.

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.^{19} 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 \u2248 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.

### 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.

*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

*z*=

*z*, we have that

_{i}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 \u2248 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.

## 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.^{21} 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

**and**

*D***are canonically conjugate up to the constant scaling factor $ 4 \pi c$. We may define the electromagnetic “kinetic energy” and “potential energy” as**

*A**T*. Hence, we recover the vacuum Maxwell equations.

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

**is involutive, we have that $ \delta H / \delta D = ( 4 \pi c ) \u2212 1 A \u0307 = ( 4 \pi ) \u2212 1 E$. Finally, we simply define $ H = 4 \pi \delta H / \delta B = B + 4 \pi \delta K / \delta B$. We, therefore, recover the macroscopic Maxwell equations.**

*D*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

**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.**

*E*### APPENDIX B: 1D SPECTRAL ELEMENT FEEC METHOD

##### 1. Construction of the FEM basis

*H*

^{1}and

*L*

^{2}studied in functional analysis augmented with the transformation rules of 0-forms and 1-forms.

*V*

^{0}and

*V*

^{1}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 $ { \Lambda 0 , j} j = 1 N 0$ and $ { \Lambda 1 , j} j = 1 N 1$. Hence, if $ \varphi h \u2208 V h 0$ and $ \rho h \u2208 V h 1$, we have

To begin, let $ \Omega = [ 0 , L ]$ be the physical domain. We divide this into *K* cells which we denote $ { \Omega k} k = 1 K$ such that $ \u222a k \Omega k = \Omega $. We first describe the construction on the reference element $ \Omega \u0302 : = [ \u2212 1 , 1 ]$ which is then mapped to each Ω_{k} by a diffeomorphism $ F k : \Omega \u0302 \u2192 \Omega 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*$ [ \u2212 1 , 1 ]$

*by $ { \xi i} i = 0 N$ and we denote the Lagrange interpolating polynomials over this grid by $ \Lambda \u0302 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*

**is simply $ I \u0302 0 v = \u2211 i v i \Lambda \u0302 0 , i$.**

*-form interpolation***Definition 2.**The

**1**

*degrees of freedom for the***are cell-wise integrals**

*-forms*

*-form interpolating polynomials are**where the equality of these three expressions follows from the fact that*

*since the Lagrange interpolating polynomials form a partition of unity.*

*N*elements, we have that $ i = 0 , \u2026 , N \u2212 1$, and therefore, $ \delta i , N = 0$ always. This yields the final result

**Definition 3.**We define the

**to be**

*discrete exterior derivative matrix*This defines all of the discrete operators on the reference element. We then define maps to each element in the physical domain, $ F k : \Omega \u0302 \u2192 \Omega 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***Proposition 1.**

*The derivative commutes with the push-forward*

*Proof:*This is simply the chain rule. □

**Definition 5.**

*The*

*local mapped finite element spaces**are simply*$ V h , k \u2113 = F k \u2113 V \u0302 h \u2113$

*. We define the*

*global finite element space**by*

*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*

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 : [ \u2212 1 , 1 ] \u2192 [ 0 , 1 / 2 ]$ and $ F 2 : [ \u2212 1 , 1 ] \u2192 [ 1 / 2 , 1 ]$ map to two distinct elements, $ F 1 ( 1 ) = F 2 ( \u2212 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

**$ \sigma \u2113 * : ( V \u2113 ) * \u2192 C \u2113 *$ are defined such that**

*dual degrees of freedom**where*$ \u27e8 \xb7 , \xb7 \u27e9 : ( V \u2113 ) * \xd7 V \u2113 \u2192 R$

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

**, $ { \Lambda \u2113 , i *} i = 0 N \u2113$**

*dual basis**, such that*$ ( \Lambda \u2113 , i , \Lambda \u2113 , j * ) = \delta i , j$

*. The*

**$ I \u2113 * : C \u2113 * \u2192 ( V h \u2113 ) *$**

*dual interpolation operator**is defined such that*

*That is,*$ I \u2113 * c *$

*is a functional which acts on elements of*$ V h \u2113$

*such that*

*The*

**mass matrix**is defined as $ ( M \u2113 ) i , j = ( \Lambda \u2113 , i , \Lambda \u2113 , j ) : C \u2113 \u2192 C \u2113 *$. The**dual projection operator**is defined such that*One can see that this is simply the L*

^{2}

*projection.*

**Proposition 2.**

*If*$ c * \u2208 C \u2113 *$

*and*$ v h \u2208 V h \u2113$

*such that*$ \sigma \u2113 ( v h ) = v \u2208 C \u2113$

*, then*

*Proof:*This simply follows from the fact that $ ( \Lambda \u2113 , i , \Lambda \u2113 , j * ) = \delta i , j$. □

**Proposition 3.**

*The mass matrix provides the change of basis from*$ C \u2113$

*to*$ C \u2113 *$

*. That is, if*$ u \u2208 V h \u2113$

*and we denote its corresponding element in the dual space*$ u h * = ( u h , \xb7 ) \u2208 ( V h \u2113 ) *$

*, then*

*Proof:*

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 \u2208 V h 0$

*and*$ u h \u2208 V h 1$

*. Then if we denote*$ v h = \u2211 \alpha \Lambda 0 , \alpha v \alpha $

*and*$ u h = \u2211 \alpha \Lambda 1 , \alpha u \alpha $

*, it follows that:*

*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

*k*th element. The global mass matrix $ M 0$ is related to each local mass matrix according to

### 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.

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 \u2113 \u2192 R$

*be an arbitrary functional on the continuous*$\u2113$

*-forms and let*$ K : = K \u2009 \u25cb \u2009 I \u2113 : C \u2113 \u2192 R$

*represent the discrete analog of the functional K. Moreover, define*$ u = \sigma \u2113 ( u )$

*and*$ v = \sigma \u2113 ( v )$

*. Then*

*Similarly, if we have a functional of the dual de Rham sequence,*$ K * : ( V \u2113 ) * \u2192 R$

*, defining*$ K * = K \u2009 \u25cb \u2009 I \u2113 * , \u2009 u * = \sigma \u2113 * ( u )$

*, and*$ v * = \sigma \u2113 * ( v )$

*, we have*

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

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 $\u2020$ to denote operator adjoints.

**Lemma 2.**

*Suppose*$ u * \u2208 C \u2113 *$

*and*$ v \u2208 C \u2113$

*. Then*

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

**Theorem 2.**

*With respect to the natural evaluation duality pairing,*$ ( \sigma \u2113 | V \u2113 ) \u2020 = I \u2113 *$

*. Moreover,*

*Hence,*$ I \u2113 *$

*approximates*$ \sigma \u2113 \u2020$

*in the following sense:*