Open access peer-reviewed chapter

# Electro-magnetic Simulation Based on the Integral Form of Maxwell’s Equations

Written By

Naofumi Kitsunezaki

Submitted: 09 July 2018 Reviewed: 05 September 2018 Published: 27 November 2018

DOI: 10.5772/intechopen.81338

From the Edited Volume

## Recent Advances in Integral Equations

Edited by Francisco Bulnes

Chapter metrics overview

View Full Metrics

## Abstract

Algorithms for a computational method of electromagnetics based on the integral form of Maxwell’s equations are introduced. The algorithms are supported by the lowest- and next-to-the-lowest-order approximations of integrals over a cell surface and edge of the equations. The method supported by the lowest-order approximation of the integrals coincides with the original finite-difference time-domain (FDTD) method, a well-known computational method of electromagnetics based on the differential form of Maxwell’s equations. The method supported by the next-to-the-lowest-order approximation can be considered a correction to the FDTD method. Numerical results of an electromagnetic wave propagating in a two-dimensional slab waveguide using the original and the corrected FDTD methods are also shown to compare them with an analytical result. In addition, the results of the corrected FDTD method are also shown to be more accurate and reliable than those of the original FDTD method.

### Keywords

• Maxwell’s equations
• integral form
• finite-difference time-domain method
• the lowest-order approximation
• next-to-the-lowest-order approximation
• computational method

## 1. Introduction

Maxwell’s equations are considered the fundamental equations of an electromagnetic field. They consist of laws of Faraday, Ampére-Maxwell, and Gauss for magnetic and electric flux densities

t B = × E , E1
t D = × H i , E2
B = 0 , E3
D = ρ , E4

where E and H are electric and magnetic fields, respectively, D and B are electric and magnetic flux densities, respectively, i is current density, ρ is charge density, t f is the time derivative of field f , × A is the rotation of vector field A , A is the divergence of vector A , and ρ is the electric charge density. Taking the divergence of both sides of Eq. (2) and using Eq. (4), law of charge conservation is derived:

t ρ + i = 0 . E5

Maxwell’s equations combined with law of Lorentz are the foundation of electronics, optics, and electric circuits used to understand the physical structure dependence of an electromagnetic field distribution, the interaction between the structure and field, and other relevant characteristics. However, situations having analytical solutions of them are rare. Thus, computational method of electromagnetics is important.

For computational methods of electromagnetics, there are two major types, time domain and frequency domain. In a time-domain method, time is discretized. The field distribution of a particular time step is determined by Maxwell’s equations and by the distribution of the previous time step. In a frequency-domain method, the time derivative is replaced by , where i is the imaginary unit and ω is the angular frequency. Thus, Maxwell’s equations are solved. A user chooses a method by considering the analysis object, calculation accuracy, specifications of his or her computer, and other relevant factors.

The finite-difference time-domain (FDTD) method is a time-domain method used to analyze high-frequency electromagnetic phenomena in optical devices, antennae, and similar devices [1]. Its algorithm is based on the laws of Faraday (1), Ampére-Maxwell (2), and charge conservation (5). In the FDTD method, Gauss’s laws (3) and (4) are not considered except for the initial condition. The reason can be easily understood by taking divergence of both sides of Eq. (1) and (2), and combining the charge conservation law (5) yields

t B = 0 , E6
t D ρ = 0 , E7

the time derivatives of Eq. (3) and (4), respectively. This means that Gauss’s laws of electric and magnetic flux densities are always satisfied when they are initially satisfied.

In the next section, an algorithm of the original FDTD method is shown. Next, a corrected algorithm of the FDTD method based on the integral form of Maxwell’s equations is shown [2, 3]. Then, a numerical result of the propagation of electromagnetic waves in a two-dimensional slab waveguide is shown. In the subsequent section, the accuracy of the original and corrected FDTD methods is compared by showing the differences between the computational and analytical methods. The analytical method is shown in the appendix. The last section is devoted to conclusions.

## 2. Algorithm of the FDTD Method

The FDTD method is a computational method for analyzing the space-time dependence of electromagnetic fields by discretizing space-time variables. This method utilizes a dual lattice called a Yee lattice [1].

Figure 1 shows the Yee lattice. In the figure, there are two cubes called cells. The component parallel to the edge of the electric field is at the center of each edge of a yellow cell. The component perpendicular to the surface of the magnetic field is at the center of each surface of the cell. The cyan cell is placed in such a manner that the component parallel to the edge of the magnetic field is at the center of its edge, and the component normal to the surface of the electric field is at the center of its surface.

The yellow cell is used to calculate the magnetic field at time t = t 0 + Δ t using the magnetic field at t = t 0 and the electric field at t = t 0 + 1 2 Δ t by applying Eq. (1), where t 0 is a particular time and Δ t is the time step. Let us now consider the top surface of the cell. At the center of the surface whose coordinates are x 0 y 0 z 0 , Eq. (1) becomes

t B z t x 0 y 0 z 0 = x E y t x 0 y 0 z 0 + y E x t x 0 y 0 z 0 , E8

where the variables x , y , and z of B and E represent the x -, y -, and z -components of the B and E fields, respectively, and x and y represent the partial derivatives in the x - and y -directions, respectively. Replacing the partial derivatives by the central differences yields

t B z t x 0 y 0 z 0 = B z t 0 + Δ t x 0 y 0 z 0 B z t 0 x 0 y 0 z 0 Δ t + O Δ t 2 , E9
x E y t x 0 y 0 z 0 = E y t 0 + 1 2 Δ t x 0 + 1 2 Δ x y 0 z 0 E y t 0 + 1 2 Δ t x 0 1 2 Δ x y 0 z 0 Δ x + O Δ y 2 , E10
y E x t x 0 y 0 z 0 = E x t 0 + 1 2 Δ t x 0 y 0 + 1 2 Δ y z 0 E x t 0 + 1 2 Δ t x 0 y 0 1 2 Δ y z 0 Δ y + O Δ x 2 , E11

where t = t 0 + 1 2 Δ t . Then, B z t 0 + Δ t x 0 y 0 z 0 are derived from Eqs. (9), (10), and (11) as the following:

B z ( t 0 + Δ t , x 0 , y 0 , z 0 ) = B z ( t 0 , x 0 , y 0 , z 0 ) + Δ t [ E y ( t 0 + 1 2 Δ t , x 0 + 1 2 Δ x , y 0 , z 0 ) Δ x + E y ( t 0 + 1 2 Δ t , x 0 1 2 Δ x , y 0 , z 0 ) Δ x + E x ( t 0 + 1 2 Δ t , x 0 , y 0 + 1 2 Δ y , z 0 ) Δ y E x ( t 0 + 1 2 Δ t , x 0 , y 0 1 2 Δ y , z 0 ) Δ y + O ( Δ x y ) 2 ] + O ( Δ t ) 3 . E12

where O Δ x Δ y 2 means m + n 2 , m , n 2 O Δ x m Δ y n . B x and B y at t = t 0 + Δ t can also be derived similarly. Usually, the H field can be derived from the B field. For example, in vacuum, air, or a dielectric

H = 1 μ 0 B , E13

where μ 0 is the vacuum permeability with the value 4 π × 10 7 Vs / Am . In a magnetic material, the relationship between H and B often becomes nontrivial, but this is beyond the scope of this book. However, in small H and B regions, it can be approximated by the following equation:

H = 1 μ B , E14

where μ depends on the material. Often, a value

μ = μ μ 0 , E15

called the relative permeability is used. In an optical wavelength region, μ 0 is 1.

The cyan cell is used to calculate the electric field at t = t 0 + 1 2 Δ t using the electric field at t = t 0 1 2 Δ t and the magnetic field at t = t Δ t applying Eq. (2) representing Ampére-Maxwell’s law. Let us consider the right-hand surface of the cell. At the center of the surface whose coordinates are x 1 y 1 z 1 , Eq. (2) becomes

t D y t x 1 y 1 z 1 = z H x t x 1 y 1 z 1 x H z t x 1 y 1 z 1 i y t x 1 y 1 z 1 . E16

Replacing partial derivatives with central differences yields

t D y t x 1 y 1 z 1 = D y t 0 + Δ t x 1 y 1 z 1 D y t 0 x 1 y 1 z 1 Δ t + O Δ t 2 , E17
x H z t x 1 y 1 z 1 = H z t 0 + 1 2 Δ t x 1 + 1 2 Δ x y 1 z 1 H z t 0 + 1 2 Δ t x 1 1 2 Δ x y 1 z 1 Δ x + O Δ x 2 , E18
z H x t x 1 y 1 z 1 = H x t 0 + 1 2 Δ t x 1 y 1 z 1 + 1 2 Δ z H x t 0 + 1 2 Δ t x 1 y 1 z 1 1 2 Δ z Δ z + O Δ z 2 , E19

where t = t 0 . Then, D y t 0 x 1 y 1 z 1 is derived as follows:

D y ( t 0 + 1 2 Δ t , x 1 , y 1 , z 1 ) = D y ( t 0 1 2 , x 1 , y 1 , z 1 ) + Δ t [ H x ( t 0 , x 1 , y 1 , z 1 + 1 2 Δ z ) Δ z H x ( t 0 , x 1 , y 1 , z 1 1 2 Δ z ) Δ z H z ( t 0 , x 1 + 1 2 Δ x , y 1 , z 1 ) Δ x + H z ( t 0 , x 1 1 2 Δ x , y 1 , z 1 ) Δ x i y ( t 0 , x 1 , y 1 , z 1 ) + O ( Δ x z ) 2 ] + O ( Δ t ) 3 . E20

D x and D z at t = t 0 + Δ t can also be derived similarly. Typically, the E field can be derived from the D flux density. For example, in vacuum, air, or magnetic material

E = 1 ε 0 D , E21

where ε 0 is the vacuum permittivity, whose value is 8.85418782 × 10 12 As / Vm . In a dielectric material, the relationship between E and D often becomes nontrivial, but this is beyond the scope of this book. However, in small E and D regions, it can be approximated by

E = 1 ε D , E22

where ε depends on the material. Often a value

ε = ε ε 0 , E23

called the relative permittivity and a value

n = ε , E24

called the index, is used.

Figure 2 shows the algorithm of the FDTD method for the case in which Eqs. (14) and (22) are satisfied. Initially, distributions of the E and H fields are given which in turn satisfy Eqs. (3) and (4). When the E field distribution at t = t 0 Δ t / 2 and the H field distribution at t = t 0 are known, the E field at t = t 0 + Δ t / 2 is calculated using Eqs. (20) and (22), given the E field at t = t 0 Δ t / 2 and the H field at t = t 0 . The H field at t = t 0 + Δ t is calculated using Eqs. (12) and (14), given the H field at t = t 0 , having determined the E field at t = t 0 + Δ t / 2 . If the time t is less than t fin , then the time becomes t + Δ t and the flow repeats. If the time t exceeds t fin , the algorithm terminates.

## 3. Integral form of Maxwell’s equation and a correction to the FDTD method

The FDTD method is a numerical method for solving Maxwell’s equations using a computer. Any computer can have a finite number of degrees of freedom because it has a finite memory size. In contrast, an electromagnetic field in continuum space has infinitely many degrees of freedom, because the field exists at every point of the space–time continuum. Therefore, Maxwell’s equations must be suitably approximated for us to be able to calculate them using a computer. The algorithm shown in the previous section appears to be suitable for this purpose, because only finitely many degrees of freedom are used to calculate an electromagnetic field distribution if the calculation area is compact.

Note that Eqs. (12) and (20) are exact on a Yee lattice only after taking a zero cell size limit. This appears to cause no problem, but, there is an example in elementary particle physics showing that the discretized continuum theory is different from the original continuum theory [4]. In this example, a fermion in discretized quantum field theory generates nonphysical fermionic degrees of freedom. This problem is called fermion doubling. In essence, this phenomenon is caused by replacing differentials with differences as in Eqs. (9)(10) and (17)(19). An algorithm in which differentials are not replaced with differences must be considered in order to avoid such problems.

As a result of Stokes’ theorem, Faraday’s and Ampére-Maxwell’s laws in Eqs. (1) and (2) can be written in an integral form as

d dt S d a B = S d s E , E25
d dt S d a D = S d s H S d a i , E26

where S is a compact and connected surface, S is the boundary curve of the surface, d a is a surface element normal to the surface, and d s is a line element parallel to the curve. Moreover, integrating both sides of Eq. (25) over t from t = t 0 to t 0 + Δ t and those of Eq. (26) over t from t = t 0 Δ t / 2 to t 0 + Δ t / 2 yields

S d a B t 0 + Δ t x y z = S d a B t 0 x y z t 0 t 0 + Δ t dt S d s E t x y z , E27
S d a D ( t 0 + Δ t / 2 , x , y , z ) = S d a D ( t 0 Δ t / 2 , x , y , z ) + t 0 Δ t / 2 t 0 + Δ t / 2 d t [ S d s H ( t , x , y , z ) S d a i ( t , x , y , z ) ] . E28

Note that no derivative is used in Eqs. (27) and (28), with the result that problems such as fermion doubling cannot occur. Our problem is how to approximate the integrals in Eqs. (27) and (28).

In general, when f is analytical in a region including ξ 0 Δ ξ / 2 ξ ξ 0 + Δ ξ / 2 , the following relationship is satisfied:

ξ 0 Δ ξ / 2 ξ 0 + Δ ξ / 2 dξf ξ = ξ 0 Δ ξ / 2 ξ 0 + Δ ξ / 2 n = 0 1 n ! d n f ξ 0 d ξ n ξ ξ 0 n , E29

and when g is analytical in a region including ξ 0 Δ ξ / 2 ξ ξ 0 + Δ ξ / 2 and η 0 Δ η / 2 η η 0 + Δ η / 2 , the following relationship is satisfied:

ξ 0 Δ ξ / 2 ξ 0 + Δ ξ / 2 η 0 Δ η / 2 η 0 + Δ η / 2 dηg ξ η = ξ 0 Δ ξ / 2 ξ 0 + Δ ξ / 2 η 0 Δ η / 2 η 0 + Δ η / 2 m , n = 0 1 m ! n ! m + n g ξ 0 η 0 ξ m η n ξ ξ 0 m η η 0 n . E30

The lowest-order approximations of Eqs. (29) and (30) are, respectively,

ξ 0 Δ ξ / 2 ξ 0 + Δ ξ / 2 dξf ξ = f ξ 0 Δ ξ + O Δ ξ 3 , E31
ξ 0 Δ ξ / 2 ξ 0 + Δ ξ / 2 η 0 Δ η / 2 η 0 + Δ η / 2 dηg ξ η = g ξ 0 η 0 Δ ξ Δ η + O Δ ξ Δ η 4 . E32

An algorithm for the FDTD method supported by the lowest-order approximation of the integral form of Maxwell’s equations is derived by applying Eqs. (31) and (32) to Eqs. (27) and (28). When S is the top surface of the yellow cell in Figure 1 , Eq. (27) is approximated as

B z ( t 0 + Δ t , x 0 , y 0 , z 0 ) Δ x Δ y = B z ( t 0 , x 0 , y 0 , z 0 ) Δ x Δ y [ E y ( t + Δ t / 2 , x 0 + Δ x / 2 , y 0 , z 0 ) Δ y E x ( t + Δ t / 2 , x 0 , y 0 + Δ y / 2 , z 0 ) Δ x E y ( t + Δ t / 2 , x 0 Δ x / 2 , y 0 , z 0 ) Δ y +   E x ( t + Δ t / 2 , x 0 , y 0 Δ y / 2 , z 0 ) Δ x ] Δ t + O ( Δ t ) 3 , E33

where x 0 y 0 z 0 is the center coordinates of the surface. Comparison of Eq. (12) with Eq. (33) reveals that they are essentially the same.

When S is the right surface of the cyan cell in Figure 1 , Eq. (28) can be approximated as

D y ( t 0 + Δ t / 2 , x 1 , y 1 , z 1 ) Δ x Δ z = D y ( t 0 Δ t / 2 , x 1 , y 1 , z 1 ) Δ x Δ z + [ H x ( t 0 , x 1 , y 1 , z 1 + Δ z / 2 ) Δ x H z ( t 0 , x 1 + Δ x / 2 , y 1 , z 1 ) Δ z H x ( t 0 , x 1 , y 1 , z 1 Δ z / 2 ) Δ x + H z ( t 0 , x 1 Δ x / 2 , y 1 , z 1 ) Δ z i y ( t 0 , x 1 , y 1 , z 1 ) Δ x Δ z ] Δ t + O ( Δ t ) 3 , E34

where x 1 y 1 z 1 is the center coordinates of the surface. Comparison of Eq. (20) with Eq. (34) reveals that they are essentially the same. Therefore, the original FDTD method, which is based on the differential form of Maxwell’s equations, is the same as the FDTD method supported by the lowest-order approximation of the integral form of those equations.

Next, an algorithm for the FDTD method supported by the next-to-the-lowest-order approximation of the integral form of Maxwell’s equation is derived. In this case, the next-to-the-lowest-order approximation is applied only in the spatial directions, and the lowest-order approximation is applied in the time direction. Hereafter, the FDTD method supported by the next-to-the-lowest-order approximation of integrals is called the corrected FDTD method.

In general, the next-to-the-lowest-order approximations of Eqs. (29) and (30) are, respectively,

ξ 0 Δ ξ / 2 ξ 0 + Δ ξ / 2 dξf ξ = f ξ 0 Δ ξ + 1 24 d 2 f ξ 0 d ξ 2 Δ ξ 3 + O Δ ξ 5 , E35
ξ 0 Δ ξ / 2 ξ 0 + Δ ξ / 2 η 0 Δ η / 2 η 0 + Δ η / 2 dηg ξ η = g ξ 0 η 0 Δ ξ Δ η + 1 24 2 g ξ η ξ 2 Δ ξ 3 Δ η + 2 g ξ η η 2 Δ ξ Δ η 3 + O Δ ξ Δ η 6 . E36

When S is the top surface of the yellow cell in Figure 1 , Eq. (27) is approximated by applying Eqs. (35) and (36) to yield

B z t 0 + Δ t x 0 y 0 z 0 Δ x Δ y + 1 24 x 2 B z t 0 + Δ t x 0 y 0 z 0 Δ x 3 Δ y + y 2 B z t 0 + Δ t x 0 y 0 z 0 Δ x Δ y 3 = B z t 0 x 0 y 0 z 0 Δ x Δ y + 1 24 x 2 B z t 0 x 0 y 0 z 0 Δ x 3 Δ y + y 2 B z t 0 x 0 y 0 z 0 Δ x Δ y 3 E y t + Δ t / 2 x 0 + Δ x / 2 y 0 z 0 Δ y + 1 24 y 2 E y ( t + Δ t / 2 x 0 + Δ x / 2 y 0 z 0 ) Δ y 3 E x t + Δ t / 2 x 0 y 0 + Δ y / 2 z 0 Δ x 1 24 x 2 E x t + Δ t / 2 x 0 y 0 + Δ y / 2 z 0 Δ x 3 E y t + Δ t / 2 x 0 Δ x / 2 y 0 z 0 Δ y 1 24 y 2 E y t + Δ t / 2 x 0 Δ x / 2 y 0 z 0 Δ y 3 + E x t + Δ t / 2 x 0 y 0 Δ y / 2 z 0 Δ x + 1 24 x 2 E x ( t + Δ t / 2 x 0 y 0 Δ y / 2 z 0 ) Δ x 3 Δ t + O Δ t 3 . E37

When S is the right-hand surface of the cyan cell in Figure 1 , Eq. (28) is approximated by applying Eqs. (35) and (36) to yield

D y t 0 + Δ t / 2 x 1 y 1 z 1 Δ x Δ z + 1 24 x 2 D y t 0 + Δ t / 2 x 1 y 1 z 1 Δ x 3 Δ y + y 2 D y t 0 + Δ t / 2 x 1 y 1 z 1 Δ x Δ y 3 = D y t 0 Δ t / 2 x 1 y 1 z 1 Δ x Δ z + 1 24 x 2 D y t 0 Δ t / 2 x 1 y 1 z 1 Δ x 3 Δ y + y 2 D y t 0 Δ t / 2 x 1 y 1 z 1 Δ x Δ y 3 + H x t 0 x 1 y 1 z 1 + Δ z / 2 Δ x + 1 24 x 2 H x ( t 0 x 1 y 1 z 1 + Δ z / 2 ) Δ x 3 H z t 0 x 1 + Δ x / 2 y 1 z 1 Δ z 1 24 z 2 H z t 0 x 1 + Δ x / 2 y 1 z 1 Δ z 3 H x t 0 x 1 y 1 z 1 Δ z / 2 Δ x 1 24 x 2 H x t 0 x 1 y 1 z 1 Δ z / 2 Δ x 3 + H z t 0 x 1 Δ x / 2 y 1 z 1 Δ z + 1 24 z 2 H z t 0 x 1 Δ x / 2 y 1 z 1 Δ z 3 i y t 0 x 1 y 1 z 1 Δ x Δ z 1 24 x 2 i y ( t 0 x 1 y 1 z 1 ) Δ x 2 Δ z 1 24 z 2 i y ( t 0 x 1 y 1 z 1 ) Δ x Δ z 2 Δ t + O Δ t 3 . E38

There are second derivatives in Eqs. (37) and (38), but they are not calculated in the FDTD method. Therefore, the second derivatives are determined from the calculated electromagnetic field. To determine the second derivatives, the relationship

f ξ 0 + Δ ξ + f ξ 0 Δ ξ = 2 f ξ + d 2 f ξ 0 d ξ 2 Δ ξ 2 + O Δ ξ 4 , E39

for any function f is applied. Applying Eq. (39) to Eq. (37) yields

5 6 B z t 0 + Δ t x 0 y 0 z 0 + 1 24 B z t 0 + Δ t x 0 + Δ x y 0 z 0 + B z ( t 0 + Δ t x 0 Δ x y 0 z 0 ) + B z t 0 + Δ t x 0 y 0 + Δ y z 0 + B z ( t 0 + Δ t x 0 y 0 Δ y z 0 ) x Δ y = 5 6 B z t 0 x 0 y 0 z 0 + 1 24 B z t 0 x 0 + Δ x y 0 z 0 + B z ( t 0 x 0 Δ x y 0 z 0 ) + B z ( t 0 x 0 y 0 + Δ y z 0 ) + B z ( t 0 x 0 y 0 Δ y z 0 ) Δ x Δ y 11 12 E y t + Δ t / 2 x 0 + Δ x / 2 y 0 z 0 + 1 24 E y t + Δ t / 2 x 0 + Δ x / 2 y 0 + Δ y z 0 + E y ( t + Δ t / 2 x 0 + Δ x / 2 y 0 Δ y z 0 ) Δ y 11 12 E x t + Δ t / 2 x 0 y 0 + Δ y / 2 z 0 + 1 24 E x t + Δ t / 2 x 0 + Δ x y 0 + Δ y / 2 z 0 + E x ( t + Δ t / 2 x 0 Δ x y 0 + Δ y / 2 z 0 ) Δ x 11 12 E y t + Δ t / 2 x 0 Δ x / 2 y 0 z 0 + 1 24 ( E y ( t + Δ t / 2 x 0 Δ x / 2 y 0 + Δ y z 0 ) + E y ( t + Δ t / 2 x 0 Δ x / 2 y 0 Δ y z 0 ) ) Δ y + 11 12 E x t + Δ t / 2 x 0 y 0 Δ y / 2 z 0 + 1 24 E x t + Δ t / 2 x 0 + Δ x y 0 Δ y / 2 z 0 + E x ( t + Δ t / 2 x 0 Δ x y 0 Δ y / 2 z 0 ) Δ x Δ t + O Δ t 3 . E40

Applying Eq. (39) to Eq. (38) yields

5 6 D y t 0 + Δ t / 2 x 1 y 1 z 1 + 1 24 ( D y ( t 0 + Δ t / 2 x 1 + Δ x y 1 z 1 ) + D y ( t 0 + Δ t / 2 x 1 Δ x y 1 z 1 ) + D y t 0 + Δ t / 2 x 1 y 1 z 1 + Δ z + D y ( t 0 + Δ t / 2 x 1 y 1 z 1 Δ z ) ) Δ x Δ z = 5 6 D y t 0 Δ t / 2 x 1 y 1 z 1 + 1 24 ( D y ( t 0 Δ t / 2 x 1 + Δ x y 1 z 1 ) + D y ( t 0 Δ t / 2 x 1 Δ x y 1 z 1 ) + D y t 0 Δ t / 2 x 1 y 1 z 1 + Δ z + D y ( t 0 Δ t / 2 x 1 y 1 z 1 Δ z ) ) Δ x Δ z + 11 12 H x t 0 x 1 y 1 z 1 + Δ z / 2 + 1 24 ( H x ( t 0 x 1 + Δ x y 1 z 1 + Δ z / 2 ) + H x ( t 0 x 1 Δ x y 1 z 1 + Δ z / 2 ) ) Δ x 11 12 H z t 0 x 1 + Δ x / 2 y 1 z 1 + 1 24 H z t 0 x 1 + Δ x / 2 y 1 z 1 + Δ z + H z ( t 0 x 1 + Δ x / 2 y 1 z 1 Δ z ) Δ z 11 12 H x t 0 x 1 y 1 z 1 Δ z / 2 + 1 24 H x t 0 x 1 + Δ x y 1 z 1 Δ z / 2 + H x ( t 0 x 1 Δ x y 1 z 1 Δ z / 2 ) Δ x + 11 12 H z t 0 x 1 Δ x / 2 y 1 z 1 + 1 24 H z t 0 x 1 Δ x / 2 y 1 z 1 + Δ z + H z ( t 0 x 1 Δ x / 2 y 1 z 1 Δ z ) Δ z 5 6 i y t 0 x 1 y 1 z 1 + 1 24 i y t 0 x 1 + Δ x y 1 z 1 + i y ( t 0 x 1 Δ x y 1 z 1 ) + i y ( t 0 x 1 y 1 z 1 + Δ z ) + i y ( t 0 x 1 y 1 z 1 Δ z ) Δ x Δ z Δ t + O Δ t 3 . E41

Note that points x i ± Δ x y i z i and x i y i z i ± Δ z are at adjacent cells to the cell including x i y i z i . Therefore, all terms in Eqs. (40) and (41) are fields defined on the Yee lattice. However, in contrast to the original FDTD method, the left-hand sides (LHSs) of these equations are a linear combination of fields at five points. Therefore, it is impossible to directly determine the values of fields at new times using these equations.

The LHSs of Eqs. (40) and (41) can be written symbolically as

m , n σ x 0 y 0 x 0 + m Δ x y 0 + n Δ y B z t 0 + Δ t x 0 + m Δ x y 0 + n Δ y z 0 , E42
m , n σ x 1 z 1 x 1 + m Δ x y 1 z 1 + n Δ z D y t + Δ t / 2 x 1 + m Δ x y 1 z 1 + n Δ z , E43

where

σ ξ η ξ + m Δ ξ η + n Δ η = 5 6 δ m , 0 δ n , 0 + 1 24 δ m , 1 δ n , 0 + δ m , 1 δ n , 0 + δ m , 0 δ n , 1 + δ m , 0 δ n , 1 , E44

and δ p , q is the Kronecker delta defined as

δ p , q = 1 p = q 0 p q . E45

The inverse operator “ σ 1 ” is defined as

p , q σ ξ η ξ + p Δ ξ η + q Δ η σ 1 ξ + p Δ ξ η + q Δ η ξ + m Δ ξ η + n Δ η = δ m , 0 δ n , 0 . E46

Using σ 1 enables Eq. (40) to be rewritten as

B z t 0 + Δ t x 0 y 0 z 0 Δ x Δ y = B z t 0 x 0 y 0 z 0 Δ x Δ y m , n σ 1 x 0 y 0 x 0 + m Δ x y 0 + n Δ y × 11 12 E y ( t 0 + Δ t / 2 x 0 + m + 1 2 Δ x y 0 + n Δ y z 0 ) E y ( t 0 + Δ t / 2 x 0 + m 1 2 Δ x y 0 + n Δ y z 0 ) Δ y E y ( t 0 + Δ t / 2 x 0 + m Δ x y 0 + n + 1 2 Δ y z 0 ) E y ( t 0 + Δ t / 2 x 0 + m Δ x y 0 + n 1 2 Δ y z 0 ) Δ x + 1 24 E y ( t 0 + Δ t / 2 x 0 + m + 1 2 Δ x y 0 + n + 1 Δ y z 0 ) + E y ( t 0 + Δ t / 2 x 0 + m + 1 2 Δ x y 0 + n 1 Δ y z 0 ) E y ( t 0 + Δ t / 2 x 0 + m 1 2 Δ x y 0 + n + 1 Δ y z 0 ) E y ( t 0 + Δ t / 2 x 0 + m 1 2 Δ x y 0 + n 1 Δ y z 0 ) Δ y E x ( t 0 + Δ t / 2 x 0 + m + 1 Δ x y 0 + n + 1 2 Δ y z 0 ) + E x ( t 0 + Δ t / 2 x 0 + m + 1 Δ x y 0 + n 1 2 Δ y z 0 ) E x ( t 0 + Δ t / 2 x 0 + m 1 Δ x y 0 + n + 1 2 Δ y z 0 ) E x ( t 0 + Δ t / 2 x 0 + m 1 Δ x y 0 + n 1 2 Δ y z 0 ) Δ x Δ t + O Δ t 3 . E47

In addition, Eq. (40) can also be rewritten as

D y t 0 + Δ t / 2 x 1 y 1 z 1 Δ x Δ z = D y t 0 Δ t / 2 x 1 y 1 z 1 + i y ( t 0 x 1 y 1 z 1 ) Δ x Δ z + m , n σ 1 x 1 z 1 x 1 + m Δ x z 1 + n Δ z × 11 12 H x ( t 0 x 1 + m Δ x y 1 z 1 + n + 1 2 Δ z ) H x ( t 0 x 1 + m Δ x y 1 z 1 + n 1 2 Δ z ) Δ x H z ( t 0 x 1 + m + 1 2 Δ x y 1 z 1 + n Δ z ) H z ( t 0 x 1 + m 1 2 Δ x y 1 z 1 + n Δ z ) Δ z + 1 24 H x ( t 0 x 1 + m + 1 Δ x y 1 z 1 + n + 1 2 Δ z ) + H x ( t 0 x 1 + m 1 Δ x y 1 z 1 + n + 1 2 Δ z H x ( t 0 x 1 + m + 1 Δ x y 1 z 1 + n 1 2 Δ z ) H x ( t 0 x 1 + m 1 Δ x y 1 z 1 + n 1 2 Δ z ) Δ x H z ( t 0 x 1 + m + 1 2 Δ x y 1 z 1 + n + 1 Δ z ) + H z ( t 0 x 1 + m + 1 2 Δ x y 1 z 1 + n 1 Δ z H z ( t 0 x 1 + m 1 2 Δ x y 1 z 1 + n + 1 Δ z ) H z ( t 0 x 1 + m 1 2 Δ x y 1 z 1 + n 1 Δ z ) ) Δ z Δ t + O Δ t 3 . E48

Then, the algorithm of the corrected FDTD method, which is supported by the next-to-the-lowest-order approximation, can be obtained by using Eqs. (48) and (47) repeatedly.

## 4. Numerical results

In this section, numerical results of electromagnetic wave transmission in a two-dimensional slab waveguide based on the original and corrected FDTD methods are compared with the analytical result.

Figure 3 shows the slab waveguide used in the computational methods, and Figure 4 shows its calculation domain. This system consists of core and cladding regions whose indices are n co and n cl , respectively. The core region is extended infinitely in the y - and z -directions and has width d in the x -direction. The cladding region is the rest of space. An electromagnetic wave propagates in the z -direction, and its electromagnetic field is assumed to have no y dependence. Because the system has no y dependence, it is essentially a two-dimensional system. An analytical solution is known and is derived in the appendix. The solution is

E x anl t x y z = β h 0 ω n co 2 ε 0 cos 2 ux d sin ωt βz x d 2 β h 0 ω n cl 2 ε 0 cos u e w 2 x d d sin ωt βz x > d 2 , E49
E y anl t x y z = 0 , E50
E z anl t x y z = 2 uh 0 ω n co 2 ε 0 d sin 2 ux d cos ωt βz x d 2 2 wh 0 ω n cl 2 ε 0 d sign x sin u e w 2 x d d cos ωt βz x > d 2 , E51
H x anl t x y z = 0 , E52
H y anl t x y z = h 0 cos 2 ux d sin ωt βz x d 2 h 0 cos u e w 2 x d d sin ωt βz x > d 2 , E53
H z anl t x y z = 0 , E54

where “anl” indicates that this is an analytical solution. u and w satisfy

w = n cl n co 2 u tan u , E55
v 2 = u 2 + w 2 = n co 2 n cl 2 d 2 π 2 λ 2 , E56

where λ is the wavelength in a vacuum, ω is the angular frequency, and β is the propagation constant. The propagation constant is the propagation directional component of wave number vector and calculated as

β = 2 π n co 2 w 2 + n cl 2 u 2 . E57

v defined by Eq. (56) is called the V-parameter, which is determined by the parameters defining the system. u and w are determined using Figure 5 . In the figure, red curves represent Eq. (55) which is symmetric under the parity transformation x x as

H y t x y z = H y t x y z . E58

Brown curves represent

w = n cl 2 n co 2 u cot u , E59

which is antisymmetric under the parity transformations x x as

H y t x y z = H y t x y z . E60

The blue curve shows Eq. (56). At each intersection of curves Eqs. (55) and (56), there is an independent symmetric mode satisfying Eq. (58), and at each of the curves Eq. (59) and (56), there is an independent antisymmetric mode satisfying Eq. (60). The mode with the lowest u is called fundamental mode. The number of modes of the system is determined by V and increases by one with respect to each π / 2 .

In the computational methods, the system parameters are set with the wavelength λ as 0.30 m, the core width d as 0.30 m, the same with the wavelength, the core index n co as 2.0 , and the cladding index n cl as 1.0 . The lengths of the cell edges Δ x and Δ z are both λ / 20 , and the time step Δ t is 10 12 s. With these parameter values, the parameter values of the analytical solution in Eqs. (49)(54) can be derived as

u = 1.50 , E61
w = 5.23 , E62
v = 5.44 , E63
βλ 2 π = 1.94 , E64

where the LHS of Eq. (64) is called the effective index, a value between n cl and n co . These parameter values show that the solution is the fundamental mode. A magnetic field is excited at z = 0 as

H y t x ,0,0 = h 0 cos 2 ux d sin ωt x d 2 h 0 sign x cos u e w 2 x d d sin ωt x > d 2 , E65

with the parameter values in Eqs. (61) and (62).

Figures 6 14 are numerical and analytical results at times at which the ωt values are integer multiples of 2 π . In these figures, violet curves represent H y / h 0 , and green curves represent core and cladding regions. The region in which the value is 0.5 is the core region with index 2.0, and the region in which the value is −0.5 is the cladding region with index 1.0. In the figures, time goes downward. The time values are 1.0 , 5.0 , 10.0 , and 20.0 ns. The left-hand column is calculated using the original FDTD method, the middle column is calculated using the corrected FDTD method, and the right-hand column is the analytical solution wherein the region

ωt < βz , E66

H z is zero. Moreover, the differences in the results between the FDTD calculations and those of the analytical ones make it clear that H y / h 0 at some points exceeds one in the FDTD calculations, even though the values at any point are equal to or less than 1 in the analytical results. However, the differences in the calculation results between the original and corrected FDTD methods are unclear. This indicates that it is impossible to conclude whether the corrected FDTD method is better than the original one using these figures.

To compare the accuracy and reliability of the original and corrected FDTD methods, we use a function err t defined as

err t = p , q H y num t p Δ x 0 q Δ z H y anl t p Δ x 0 q Δ z n p Δ x 0 q Δ z 2 p , q H y anl t p Δ x 0 q Δ z 2 n p Δ x 0 q Δ z 2 , E67

which shows the error between the FDTD and the analytical calculations at each time. In Eq. (67), the denominator of the right-hand side is proportional to the power of the propagating electromagnetic wave passing through the x y plane per unit length of the y -direction.

Figure 15 shows the err functions of the original and corrected FDTD methods defined by Eq. (67). As shown in the figure, almost every time except for less than 0.14 ns, the err function of the corrected FDTD method is less than that of the original. This means that the corrected method is more accurate than that of the original. In addition, when the time is greater than 6 ns, both curves begin to oscillate. The amplitude of the oscillation of the corrected FDTD method is clearly less than that of the original. This indicates that the corrected method is more reliable.

## 5. Conclusion

In this chapter, a higher-order correction to the original FDTD method supported by the next-to-the-lowest-order approximation of the integral form of Maxwell’s equation was shown. The essence of this method is the approximation of integrals over a cell surface and edge using discretized electric and magnetic fields.

The results of numerical calculations of an electromagnetic wave propagating in a two-dimensional slab waveguide using the corrected and original FDTD methods and analysis were also shown. The differences between the corrected and original FDTD methods were compared using the err function, and the corrected method was found to be more accurate and reliable than the original.

## Acknowledgments

The author would like to thank Mr. A. Okabe, who is the collaborator of reference [2], which is based on this chapter. The author would also like to thank Enago (www.enago.com) for the English language review.

## Conflict of interest

The author declares no conflicts of interest associated with this manuscript.

A.1 Analytical results of two-dimensional slab waveguide

The analytic solution of an electromagnetic wave in a slab waveguide shown in Figure 3 is provided in various textbooks regarding optical waveguides and related fields [5, 6, 7, 8]. In this appendix, the analytical results of Eqs. (49)(54) are derived in accordance with these references.

A propagating electromagnetic wave with no y -dependence in the z -direction with angular frequency ω and propagation constant β , the wave number in the z -direction, is written as

E t x y z = e x e i ωt βz , E68
H t x y z = h x e i ωt βz . E69

Maxwell’s equations in dielectrics using Eqs. (13) and (22) become

i μ 0 ω h x x = e y x E70
i μ 0 ω h y x = e x x + x e z x , E71
i μ 0 ω h z x = x e y x , E72
in x 2 ε 0 ω e x x = h y x , E73
in x 2 ε 0 ω e y x = h x x x h z x , E74
in x 2 ε 0 ω e z x = x h y x , E75

where n x is the index distribution shown in Figure 16 . As shown in Eqs. (70)(75), there are two closed equation classes. The first class contains Eqs. (70), (72), and (74), which have components of electric and magnetic fields transverse to the propagation direction and a longitudinal component of the magnetic field in that direction. The second class contains Eqs. (71), (73), and (75), which have components of electric and magnetic fields transversed to the propagation direction and a longitudinal component of the electric field in that direction. Solution to the first class comprise the transverse electric (TE) mode because the electric field has only a component transversed to the propagation direction, and solution to the second class comprises the transverse magnetic (TM) mode because the magnetic field has only a component transversed to that direction. In Section 4, numerical and analytical results of H y are shown, and they are TM modes.

Hereafter, our discussion is limited to the TM mode. Then, Eqs. (71), (73), and (75) are rewritten as

x 2 h y x = 2 πn x λ 2 β 2 h y x , E76
e x x = β n x 2 ε 0 ω h y x , E77
e z x = 1 in x 2 ε 0 ω x h y x . E78

The discontinuity of n x at x = d / 2 and x = d / 2 , shown in Figure 16 , requires that the boundary condition at x = ± d / 2 be considered. Because of the integral form of Maxwell’s equations, the components of electric and magnetic fields parallel to the boundary surface of the indices are continuous, and the components of the electric and magnetic flux densities normal to the surface are continuous. Consequently, H y x and x h y x in Eqs. (76)(78) are continuous.

Solving Eq. (76) requires considering the following three cases:

1. β < n cl

2. n cl β < n co

3. n co β

In case 1, the solutions of Eq. (76) are

h y x = Q cos Px x d / 2 Q cos Pd 2 cos Qd 2 + P sin Pd 2 sin Qd 2 cos Qx P sin Pd 2 cos Qd 2 Q cos Pd 2 sin Qd 2 sin Qx x d / 2 , E79

and

h y x = Q sin Px x d / 2 sign x { Q sin Pd 2 cos Qd 2 P cos Pd 2 sin Qd 2 cos Qx + P cos Pd 2 cos Qd 2 + Q sin Pd 2 sin Qd 2 sin Qx } x d / 2 , E80

where

P = 2 π n co λ 2 β 2 , E81
Q = 2 π n cl λ 2 β 2 . E82

However, this solution does not show electromagnetic wave propagation in the z -direction but shows a reflection and transmission problem of the film when the incident angle is less than the critical angle. Let us discuss this in more detail. When x d / 2 , H y t x y z is a linear combination of

e i ωt Px βz , and e i ωt + Px βz , E83

a plane wave whose wavenumber is P 0 β and P 0 β , respectively. When x > d / 2 , H y t x y . z is a linear combination of

e i ωt Qx βz , and e i ωt + Qx βz , E84

a plane wave whose wavenumbers are Q 0 β and Q 0 β , respectively. Therefore, with a suitable linear combination of Eqs. (79) and (80), the solution becomes that a plane wave with wavenumber P 0 β is incident from the x < d / 2 region to be reflected and transmitted by a film of the x d / 2 region with a reflected wave propagated in the x < d / 2 region and a transmitted wave propagated in the x > d / 2 region. Therefore, this solution is not what we want.

In case 2, the solutions of Eq. (76) are

h y x = cos Px x d / 2 cos Pd 2 e Q 2 x d d , E85

where

P = 2 π n co λ 2 β 2 , E86
Q = β 2 2 π n cl λ 2 , E87
Qd 2 = n cl n c o 2 Pd 2 tan Pd 2 , E88

and

h y x = sin Px x d / 2 sign x sin Pd 2 e Q 2 x d d x > d / 2 , E89

where

Qd 2 = n cl n c o 2 Pd 2 cot Pd 2 . E90

Defining

u = Pd 2 , E91
w = Qd 2 , E92

yields Eqs. (53), (55), (56), and (59). This is the solution we want.

In case 3, the solutions of Eq. (76) are

h y x = cosh Px x d / 2 cosh Pd 2 e Q x d / 2 x > d / 2 , E93

where

P = β 2 2 π n co λ 2 , E94
Q = β 2 2 π n cl λ 2 , E95
Q = P tanh Pd 2 , E96

and

h y x = sinh Px x d / 2 sign x sinh Pd 2 e Q x d / 2 x > d / 2 , E97

where

Q = P coth Pd 2 . E98

As a result of Eqs. (94) and (95),

P 2 + Q 2 = 2 π n co λ 2 2 π n cl λ 2 , E99

is satisfied. When P is pure imaginary and Q is real, the solution reduces to case 2. There is a possibility that Q is neither real nor purely imaginary. However, such a solution must be attenuated when z becomes large. Mathematically, there can be a divergent solution when z becomes large, but such a solution cannot conserve energy. Therefore, such a solution is not what we want.

## References

1. 1. Yee K. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Transactions on Antennas and Propagation. 1966;14:302-307. DOI: 10.1109/TAP.1966.1138693
2. 2. Kitsunezaki N, Okabe A. Higher-order correction to the FDTD method based on the integral form of Maxwell’s equations. Computer Physics Communications. 2014;185:1582-1588. DOI: 10.1016/j.cpc.2014.02.022
3. 3. Kitsunezaki N. Higher-order Correction to the Finite-Difference Time-Domain Method Based on the Integral Form of Maxwell’s Equations. In: BIT’s 4th Annual Global Congress of Knowledge Economy; 19–21 September 2017; Qingdao, Dalian: BIT Group Global Ltd.; 2017. p. 096
4. 4. Nielsen HN, Ninomiya M. A no-go theorem for regularizing chiral fermions. Physics Letters. 1981;B105:219-223. DOI: 10.1016/0370-2693(81)91026-1
5. 5. Marcuse D. Light Transmission Optics. 2nd ed. New York: Van Nostrand Reinhold Company Inc; 1982. 534 p. ISBN: 0-442-26309-0
6. 6. Haus HA. Waves and Fields in Optoelectronics. Prentice Hall: Englewood Criffs; 1983. 402 p. ISBN: 0-139-460-535
7. 7. Buch JA. Fundamentals of Optical Fibers. Hoboken: Wiley; 2004. 332 p. ISBN: 0-471-221-910
8. 8. Okamoto K. Fundamentals of Optical Waveguides. San Diego: Academic Press; 2006. 561 p. ISBN: 0-12-525096-7

Written By

Naofumi Kitsunezaki

Submitted: 09 July 2018 Reviewed: 05 September 2018 Published: 27 November 2018