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

(2.1) | ||||

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

(2.3a) | |||

(2.3b) |

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

(2.4a) | |||

(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:

(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

(2.6) |

where the energy is defined as

(2.7) |

###### Theorem 2.

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

(2.8) |

where the energy is defined as

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

(3.1) |

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

(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

(3.3) |

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

(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

(3.5a) | ||||

(3.5b) |

and the discrete fourth order finite difference operators,

(3.6a) | ||||

(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

(3.7a) | ||||

(3.7b) |

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

(3.8a) | ||||

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

(3.9a) | ||||

(3.9b) | ||||

(3.9c) | ||||

(3.9d) |

### 3.3 Modified equation approach

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

(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

leading to

(3.11) |

Finally the discrete version of (3.11) becomes:

(3.12a) | |||

Similarly, we obtain the following equations for the other variables | |||

(3.12b) | |||

(3.12c) | |||

(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

(3.13a) | ||||

(3.13b) | ||||

(3.13c) | ||||

(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

(4.1) |

where the energy is defined as

(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

(4.3) |

where , and the matrix operators and are defined by

(4.4a) |

with

Comments

There are no comments yet.