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

by   Puttha Sakkaplangkul, et al.

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.



There are no comments yet.


page 1

page 2

page 3

page 4


Exponential integrators for second-order in time partial differential equations

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

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

We derive thermodynamically consistent models for diblock copolymer solu...

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

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

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

Positive and free energy satisfying schemes for diffusion with interaction potentials

In this paper, we design and analyze second order positive and free ener...

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


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


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


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


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:


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


where the energy is defined as

Theorem 2.

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


where the energy is defined as


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


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


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


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


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


and the discrete fourth order finite difference operators,


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


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


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),


3.3 Modified equation approach

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


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

leading to


Finally the discrete version of (3.11) becomes:

Similarly, we obtain the following equations for the other variables
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


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


where the energy is defined as


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


where , and the matrix operators and are defined by