# A Fully Fourth Order Accurate Energy Stable FDTD Method for Maxwell's Equations in Metamaterials

We present a novel fully fourth order in time and space finite difference method for the time domain Maxwell's equations in metamaterials. We consider a Drude metamaterial model for the material response to incident electromagnetic fields. We consider the second order formulation of the system of partial differential equations that govern the evolution in time of electric and magnetic fields along with the evolution of the polarization and magnetization current densities. Our discretization employs fourth order staggering in space of different field components and the modified equation approach to obtain fourth order accuracy in time. Using the energy method, we derive energy relations for the continuous models, and design numerical schemes that preserve a discrete analogue of the energy relation. Numerical simulations are provided in one and two dimensional settings to illustrate fourth order convergence as well as compare with second order schemes.

## Authors

• 1 publication
• 2 publications
• 4 publications
11/29/2021

### Exponential integrators for second-order in time partial differential equations

Two types of second-order in time partial differential equations (PDEs),...
02/16/2021

### Thermodynamically Consistent Algorithms for Models of Diblock Copolymer Solutions Interacting with Electric and Magnetic Fields

We derive thermodynamically consistent models for diblock copolymer solu...
04/21/2021

### A convergent numerical scheme for a model of liquid crystal dynamics subjected to an electric field

We present a convergent and constraint-preserving numerical discretizati...
04/29/2021

### A Feynman-Kac based numerical method for the exit time probability of a class of transport problems

The exit time probability, which gives the likelihood that an initial co...
08/21/2021

### High-order accurate schemes for Maxwell's equations with nonlinear active media and material interfaces

We describe a fourth-order accurate finite-difference time-domain scheme...
09/30/2019

### Positive and free energy satisfying schemes for diffusion with interaction potentials

In this paper, we design and analyze second order positive and free ener...
09/04/2021

### The anatomy of Boris type solvers and large time-step plasma simulations

This work gives a Lie operator derivation of various Boris solvers via a...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

In this paper, we present the construction and analysis of fully fourth order in space and time finite difference methods based on staggered finite differences in space and the modified equation approach in time [7] for electromagnetic metamaterial models. Metamaterials have revealed a great opportunity for controlling light in nanophotonic devices, in particular by controlling surface plasmons that appear at the interfaces between the materials, with applications for antennas and optical cloaking (e.g. [4, 16, 6, 11, 14, 9, 1]). Surface plasmons are sub-wavelength and highly sensitive to the interface’s geometry. Therefore there is a need for accurate evaluations of the electromagnetic near fields that capture all the multiple scales inherent to those problems.

There exist several papers on finite difference time domain (FDTD) schemes for Maxwell’s equations in metamaterials (e.g. [17, 12, 10, 13]

). High-order schemes based on the first-order formulation of Maxwell’s equations can be derived, but no associated discrete energy estimates have been proved beyond second order in time.

In this paper, we consider models based on the second order formulation of Maxwell’s equations for electromagnetic wave propagation in materials that are characterized as Drude metamaterials [10]. By converting temporal derivatives to spatial derivatives via the modified equation approach, we construct full fourth-order FDTD schemes that mimic, at the discrete level, the energy estimate of the continuous models. Our construction is motivated by high order finite difference methods for the second order wave equation in the time domain as constructed in [15, 7, 8]. High-order numerical methods using second-order formulations have been recently used for dispersive media [2]. However, to our knowledge, discrete energy estimates have not been proved except for the scalar wave equation [7].

The paper is organized as follows: Section 2 introduces the Maxwell-Drude model and its properties, Section 3 presents the construction of the fourth-order FDTD schemes, Section 4 establishes stability of the semi and fully-discrete schemes in one spatial dimension, section 5 presents numerical results for one and two dimensions, and section 6 gives our concluding remarks.

## 2 Setting

### 2.1 Drude Metamaterial Model: First Order Formulation

We consider Maxwell’s equations in a Drude metamaterial [10], given in the form

 ε0∂tE+J=curlH, (2.1) μ0∂tH+K=−curlE, ∂tJ=ε0ω2peE, ∂tK=μ0ω2pmH,

where and are the electric and magnetic fields, respectively, and are the polarization and magnetization current densities, respectively. The parameters and are the vacuum electric permittivity and magnetic permeability, respectively. Finally, the parameters and are the electric and magnetic plasma frequencies, respectively. System (2.1) has to be closed by adding initial and boundary conditions.

System (2.1

) is linear and thus admits harmonic plane wave solutions. Taking Fourier transforms in time, i.e. assuming dependence of solutions on

, one finds that the Drude metamaterial model is characterized by a dispersive permittitivity and permeability given as

 ε(ω)=ε0(1−ω2peω2),μ(ω)=μ0(1−ω2pmω2), (2.2)

following Drude’s law [10]. We call system (2.1) the Maxwell-Drude model.

### 2.2 Drude Metamaterial Model: Second Order Formulation

System (2.1) can be rewritten in second order form in which we get two sets of decoupled equations, with one pair of equations involving only the field variables modeled by the system of equations

 ∂ttE+c2curlcurlE+ω2peE=−c2curlK, (2.3a) ∂ttK+ω2pmK=−ω2pmcurlE, (2.3b)

and the second pair involving the field variables modeled by the system of equations

 ∂ttH+c2curlcurlH+ω2pmH=c2curlJ, (2.4a) ∂ttJ+ω2peJ=ω2pecurlH. (2.4b)

Since the two sub-systems in (2.3) and (2.4) are decoupled they can be solved separately. We note that in 2D and 3D, these systems will be coupled through the divergence conditions. In the 1D case, when the divergence conditions are decoupled from the curl equations, these systems are truly decoupled. In this paper, we do not consider the divergence conditions.

### 2.3 Energy Estimates

To derive energy estimates we need to define the following functional space:

 H#(curl;Ω):={v∈(L2(Ω))3,curlv∈(L2(Ω))3,v periodic on ∂Ω}. (2.5)

From now on denotes the norm in . From the second order formulation one can check that we have the following energy estimates.

###### Theorem 1.

Assume , . Under the assumption of periodic boundary conditions, system (2.3) satisfies the following energy estimate

 dEEKdt=0, (2.6)

where the energy is defined as

 EEK=12(1c2||∂tE||2+1ω2pm ||∂tK||2+ω2pec2||E||2+||curlE+K||2). (2.7)
###### Theorem 2.

Assume , . Under the assumption of periodic boundary conditions, system (2.4) satisfies the energy estimate

 dEHJdt=0, (2.8)

where the energy is defined as

 EHJ=12(1c2||∂tH||2+1ω2pe ||∂tJ||2+ω2pmc2||H||2+||curlH−J||2). (2.9)

The proofs are based on the energy method which employs integration by parts.

## 3 Fourth order FDTD schemes

### 3.1 Staggered High Order Spatial Finite Difference Methods

In this section we describe the construction of second and fourth order approximations to the first and second order derivative operators to get approximation to the vector

operator in 2D and 3D and the scalar curl operator in 1D and 2D. To that aim we need to define approximations of the derivative operators , , . The construction presented in this section follows the exposition in [8, 5].

Consider a domain , with . Given , let be a uniform mesh step size, and define for , 444 denotes the set of integers .. For , define , , . We define several staggered grids on . The primal grid in a direction along one of the axes, , is defined as

 Gp:(Xℓ)ℓ,X∈{x,y,z},l∈⟦0,M⟧. (3.1)

The dual grid in one direction, , on is defined as

 Gd:(Xℓ+1/2)ℓ,X∈{x,y,z},ℓ∈⟦0,M−1⟧. (3.2)

We define staggered grids on as , and , , . For smooth functions and , define , with , , and , with , . On the primal grid we define the discrete space

 V0,h:={vh:=(vℓ),ℓ∈⟦0,M⟧d,hM∑ℓ=0|vℓ|2<∞} + periodic B.C, (3.3)

and on the dual grids , , we define the discrete spaces

 Vm12,h (3.4)

The norms on , and , denoted by , and , respectively, are derived from corresponding scalar products , and . We now define the discrete second order finite difference operators

 F(2)h,m:V0,h→Vm12,h, ∀vh∈V0,h,(F(2)h,mvh)ℓ+12:=vmℓ+1−vℓh, (3.5a) F(2)h,m∗:Vm12,h→V0,h, ∀umh∈Vm12,h,(F(2)h,m∗umh)ℓ:=umℓ+12−umℓ−12h, (3.5b)

and the discrete fourth order finite difference operators,

 F(4)h,m:V0,h→Vm12,h, ∀vh∈V0,h,(F(4)h,mvh)ℓ+12:=98vmℓ+1−vℓh−124vmℓ+2−vmℓ−1h, (3.6a) F(4)h,m∗:Vm12,h→V0,h, ∀umh∈Vm12,h,(F(4)h,m∗umh)ℓ:=98umℓ+12−umℓ−12h−124umℓ+32−umℓ−32h. (3.6b)

Operators (3.5)-(3.6) are the second-order and fourth-order discrete approximations of the operator , with stepsize , respectively. Finally we can define the second and fourth order discrete approximations of the 3D operator in terms of , and , . For all , , and for , we have

 (curl(p)hvh)ℓ+12:=((F(p)h,yvzh)ℓ+12−(F(p)h,zvyh)ℓ+12,(F(p)h,zvxh)ℓ+12−(F(p)h,xvzh)ℓ+12,(F(p)h,xvyh)ℓ+12−(F(p)h,yvxh)ℓ+12), (3.7a) (curl(p)h,∗uh)ℓ:=((F(p)h,y∗uzh)ℓ−(F(p)h,z∗uyh)ℓ,(F(p)h,z∗uxh)ℓ−(F(p)h,x∗uzh)ℓ,(F(p)h,x∗uyh)ℓ−(F(p)h,y∗uxh)ℓ). (3.7b)

Note that for one needs scalar and vector curl operators. For , their discrete version is defined as

 (curl(p)hvh)ℓ+12 :=((F(p)h,yvh)ℓ+12,−(F(p)h,xvh)ℓ+12), (3.8a) (curl(p)hvh)ℓ+12 :=(F(p)h,xvyh)ℓ+12−(F(p)h,yvxh)ℓ+12, (3.8b)

and similarly for , . Finally for , the curl operator and its dual, up to a sign, reduce to and , , , respectively.

### 3.2 Semi-discretization

We start with the discretization of the system for the pair and then follow that with the discretization for the second sub-system corresponding to the pair . The discretizations are based on the staggered finite difference methods presented above on dual and primal grids.

Using the discretization introduced in section 3.1, we obtain the following semi-discretizations: for defined on a primal grid, and defined on a dual grid, (similarly for defined on a primal grid, and defined on a dual grid),

 ∂ttEh+c2(curl(4)h,∗curl(4)h)Eh+ω2peEh=−c2curl(4)h,∗Kh, (3.9a) ∂ttKh+ω2pmKh=−ω2pmcurl(4)hEh, (3.9b) ∂ttHh+c2(curl(4)h,∗curl(4)h)Hh+ω2pmHh=c2curl(4)h,∗Jh, (3.9c) ∂ttJh+ω2peJh=ω2pecurl(4)hHh. (3.9d)

### 3.3 Modified equation approach

Given , we define , and , . For all , we have

 ∂ttXnℓ =Xn+1ℓ−2Xnℓ+Xn−1ℓΔt2−Δt212∂4tXnℓ+O(Δt4)=δttXnℓ−Δt212∂4tXnℓ+O(Δt4). (3.10)

Based on [7], we use the modified equation approach to convert fourth order time derivatives into derivatives in space. If the spatial derivatives are then discretized by second order finite differences, we will be able to achieve overall fourth order accuracy in space and time. Starting with (2.3a), and finding an expression of in terms of spatial derivatives, assuming we can interchange time and space differentiation and using (2.3) we have

 ∂4tE =−c2curlcurl∂ttE−ω2pe∂ttE−c2curl∂ttK, = −c2curlcurl(−c2curlcurlE−ω2peE−c2curlK)−ω2pe(−c2curlcurlE−ω2peE−c2curlK) − c2curl(−ω2pmK−ω2pmcurlE),

 ∂4tE=c4(curlcurl)2E+c2(2ω2pe+ω2pm)curlcurlE+ω4peE+c4curlcurlcurlK+c2(ω2pe+ω2pm)curlK. (3.11)

Finally the discrete version of (3.11) becomes:

 ∂4tEh=c4(curl(2)h,∗curl(2)h)2Eh+c2(2ω2pe+ω2pm)(curl(2)h,∗curl(2)h)Eh+ω4peEh (3.12a) +c4(curl(2)hcurl(2)h,∗curl(2)h)Kh+c2(ω2pe+ω2pm)curl(2)h,∗Kh. Similarly, we obtain the following equations for the other variables ∂4tKh=c2ω2pm(curl(2)hcurl(2)h,∗curl(2)h)Eh+ω2pm(ω2pe+ω2pm)curl(2)hEh+c2ω2pm(curl(2)hcurl(2)h,∗)Kh+ω4pmKh, (3.12b) ∂4tHh=c4(curl(2)h,∗curl(2)h)2Hh+c2(2ω2pm+ω2pe)(curl(2)h,∗curl(2)h)Hh+ω4pmHh (3.12c) −c4(curl(2)hcurl(2)h,∗curl(2)h)Jh−c2(ω2pm+ω2pe)curl(2)h,∗Jh, ∂4tJh=−c2ω2pe(curl(2)hcurl(2)h,∗curl(2)h)Hh−ω2pe(ω2pm+ω2pe)curl(2)hHh+c2ω2pe(curl(2)hcurl(2)h,∗)Jh+ω4peJh. (3.12d)
###### Remark 1.

For this model, the schemes follow exactly [7]. The reason is the system is composed of wave-like equations, with a diffusion term.

### 3.4 Full discretization

Plugging the modified equations (3.12) into (3.10) for each field, then into the semi-discrete equations (3.9) we obtain the fully discrete scheme

 δttEnℓ+c2(curl(4)h,∗curl(4)h)Enℓ+ω2peEnℓ+c2curl(4)h,∗Knℓ+1/2−Δt212[c4(curl(2)h,∗curl(2)h)2Enℓ+c2(2ω2pe+ω2pm)(curl(2)h,∗curl(2)h)Enℓ+ω4peEnℓ+c4(curl(2)hcurl(2)h,∗curl(2)h)Knℓ+1/2+c2(ω2pe+ω2pm)curl(2)h,∗Knℓ+1/2]=0, (3.13a) δttKnℓ+1/2+ω2pmKnℓ+1/2+ω2pmcurl(4)hEnℓ−Δt212[c2ω2pm(curl(2)hcurl(2)h,∗curl(2)h)Enℓ+ω2pm(ω2pe+ω2pm)curl(2)hEnℓ+c2ω2pm(curl(2)hcurl(2)h,∗)Knℓ+1/2+ω4pmKnℓ+1/2]=0, (3.13b) δttHnℓ+c2(curl(4)h,∗curl(4)h)Hnℓ+ω2pmHnℓ−c2curl(4)h,∗Jnℓ+1/2−Δt212[c4(curl(2)h,∗curl(2)h)2Hnℓ+c2(2ω2pm+ω2pe)(curl(2)h,∗curl(2)h)Hnℓ+ω4pmHnℓ−c4 (curl(2)hcurl(2)h,∗curl(2)h)Jnℓ+1/2−c2(ω2pm+ω2pe)curl(2)h,∗Jnℓ+1/2]=0, (3.13c) δttJnℓ+1/2+ω2peJnℓ+1/2−ω2pecurl(4)hHnℓ−Δt212[−c2ω2pe(curl(2)hcurl(2)h,∗curl(2)h)Hnℓ−ω2pe(ω2pm+ω2pe)curl(2)hHnℓ+c2ω2pe(curl(2)hcurl(2)h,∗)Jnℓ+1/2+ω4peJnℓ+1/2]=0, (3.13d)

where is the three point time stencil defined in (3.10). The discrete system (3.13) is then fourth order in time and space. We will refer to it as (4,4)-scheme.

###### Remark 2.

In the numerical examples we will compare with the so-called -FDTD scheme and -FDTD scheme. The -scheme is obtained from (3.13) by dropping the terms multiplied by , while the -scheme is obtained from (3.13) by dropping the terms multiplied by and replacing the fourth order discrete operators by second order discrete operators.

## 4 Stability (in 1D)

For simplicity, we detail the energy estimates for the semi-discrete and fully-discrete schemes in one dimension, and for the pair . Similar results hold for the other pair, and similarly for 2D and 3D. We now consider generically , and the scalar operator curl.

### 4.1 Semi-discrete energy estimates

The semi-discrete system for (E,K) in one dimension satisfies the following

###### Theorem 3 (Semi-discrete Energy Estimate).

Assume , . Then the semi-discrete in space system satisfies the energy estimate

 dEh,EKdt=0, (4.1)

where the energy is defined as

 EEK= 12(1c2||∂tEh||20,h+1ω2pm||∂tKh||212,h+ω2pec2||Eh||20,h+||Kh+curl(4)hEh||212,h). (4.2)

For simplicity we denote the discrete operator norms
and , without precising : will be fixed by the considered polarization. The proof is based on the fact that the operator is the adjoint of the discrete operator with respect to the scalar product. The rest follows the same steps as for the continuous energy estimate.

### 4.2 Full-discrete energy estimates

The full fourth-order discretization of the one dimensional version of (3.13a)-(3.13b) can be rewritten in vector form

 PδttWnj+A1Wnj−Δt2c212A2Wnj=0, (4.3)

where , and the matrix operators and are defined by

 A2=[a11a12a21a22], (4.4a)

with

 a11:=(curl(2)h,∗curl(2)h)2−(2ω2p