Open access peer-reviewed chapter

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

Written By

Naofumi Kitsunezaki

Submitted: July 9th, 2018 Reviewed: September 5th, 2018 Published: November 27th, 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

tB=×E,E1
tD=×Hi,E2
B=0,E3
D=ρ,E4

where Eand Hare electric and magnetic fields, respectively, Dand Bare electric and magnetic flux densities, respectively, iis current density, ρis charge density, tfis the time derivative of field f, ×Ais the rotation of vector field A, Ais 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 iis 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

tB=0,E6
tDρ=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=t0+Δtusing the magnetic field at t=t0and the electric field at t=t0+12Δtby applying Eq. (1), where t0is a particular time and Δtis the time step. Let us now consider the top surface of the cell. At the center of the surface whose coordinates are x0y0z0, Eq. (1) becomes

tBztx0y0z0=xEytx0y0z0+yExtx0y0z0,E8

where the variables x, y, and zof Band Erepresent the x-, y-, and z-components of the Band Efields, respectively, and xand yrepresent the partial derivatives in the x- and y-directions, respectively. Replacing the partial derivatives by the central differences yields

tBztx0y0z0=Bzt0+Δtx0y0z0Bzt0x0y0z0Δt+OΔt2,E9
xEytx0y0z0=Eyt0+12Δtx0+12Δxy0z0Eyt0+12Δtx012Δxy0z0Δx+OΔy2,E10
yExtx0y0z0=Ext0+12Δtx0y0+12Δyz0Ext0+12Δtx0y012Δyz0Δy+OΔx2,E11

where t=t0+12Δt. Then, Bzt0+Δtx0y0z0are derived from Eqs. (9), (10), and (11) as the following:

Bz(t0+Δt,x0,y0,z0)=Bz(t0,x0,y0,z0)+Δt[Ey(t0+12Δt,x0+12Δx,y0,z0)Δx+Ey(t0+12Δt,x012Δx,y0,z0)Δx+Ex(t0+12Δt,x0,y0+12Δy,z0)ΔyEx(t0+12Δt,x0,y012Δy,z0)Δy+O(Δxy)2]+O(Δt)3.E12

where OΔxΔy2means m+n2,m,n2OΔxmΔyn. Bxand Byat t=t0+Δtcan also be derived similarly. Usually, the Hfield can be derived from the Bfield. For example, in vacuum, air, or a dielectric

H=1μ0B,E13

where μ0is the vacuum permeability with the value 4π×107Vs/Am. In a magnetic material, the relationship between Hand Boften becomes nontrivial, but this is beyond the scope of this book. However, in small Hand Bregions, 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, μ0is 1.

The cyan cell is used to calculate the electric field at t=t0+12Δtusing the electric field at t=t012Δtand the magnetic field at t=tΔtapplying 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 x1y1z1, Eq. (2) becomes

tDytx1y1z1=zHxtx1y1z1xHztx1y1z1iytx1y1z1.E16

Replacing partial derivatives with central differences yields

tDytx1y1z1=Dyt0+Δtx1y1z1Dyt0x1y1z1Δt+OΔt2,E17
xHztx1y1z1=Hzt0+12Δtx1+12Δxy1z1Hzt0+12Δtx112Δxy1z1Δx+OΔx2,E18
zHxtx1y1z1=Hxt0+12Δtx1y1z1+12ΔzHxt0+12Δtx1y1z112ΔzΔz+OΔz2,E19

where t=t0. Then, Dyt0x1y1z1is derived as follows:

Dy(t0+12Δt,x1,y1,z1)=Dy(t012,x1,y1,z1)+Δt[Hx(t0,x1,y1,z1+12Δz)ΔzHx(t0,x1,y1,z112Δz)ΔzHz(t0,x1+12Δx,y1,z1)Δx+Hz(t0,x112Δx,y1,z1)Δxiy(t0,x1,y1,z1)+O(Δxz)2]+O(Δt)3.E20

Dxand Dzat t=t0+Δtcan also be derived similarly. Typically, the Efield can be derived from the Dflux density. For example, in vacuum, air, or magnetic material

E=1ε0D,E21

where ε0is the vacuum permittivity, whose value is 8.85418782×1012As/Vm. In a dielectric material, the relationship between Eand Doften becomes nontrivial, but this is beyond the scope of this book. However, in small Eand Dregions, 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 Eand Hfields are given which in turn satisfy Eqs. (3) and (4). When the Efield distribution at t=t0Δt/2and the Hfield distribution at t=t0are known, the Efield at t=t0+Δt/2is calculated using Eqs. (20) and (22), given the Efield at t=t0Δt/2and the Hfield at t=t0. The Hfield at t=t0+Δtis calculated using Eqs. (12) and (14), given the Hfield at t=t0, having determined the Efield at t=t0+Δt/2. If the time tis less than tfin, then the time becomes t+Δtand the flow repeats. If the time texceeds tfin, 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

ddtSdaB=SdsE,E25

where Sis a compact and connected surface, Sis the boundary curve of the surface, dais a surface element normal to the surface, and dsis a line element parallel to the curve. Moreover, integrating both sides of Eq. (25) over tfrom t=t0to t0+Δtand those of Eq. (26) over tfrom t=t0Δt/2to t0+Δt/2yields

SdaBt0+Δtxyz=SdaBt0xyzt0t0+ΔtdtSdsEtxyz,E27

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 fis analytical in a region including ξ0Δξ/2ξξ0+Δξ/2, the following relationship is satisfied:

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

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

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

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

ξ0Δξ/2ξ0+Δξ/2dξfξ=fξ0Δξ+OΔξ3,E31
ξ0Δξ/2ξ0+Δξ/2η0Δη/2η0+Δη/2dη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 Sis the top surface of the yellow cell in Figure 1 , Eq. (27) is approximated as

Bz(t0+Δt,x0,y0,z0)ΔxΔy=Bz(t0,x0,y0,z0)ΔxΔy[Ey(t+Δt/2,x0+Δx/2,y0,z0)ΔyEx(t+Δt/2,x0,y0+Δy/2,z0)ΔxEy(t+Δt/2,x0Δx/2,y0,z0)Δy+ Ex(t+Δt/2,x0,y0Δy/2,z0)Δx]Δt+O(Δt)3,E33

where x0y0z0is the center coordinates of the surface. Comparison of Eq. (12) with Eq. (33) reveals that they are essentially the same.

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

Dy(t0+Δt/2,x1,y1,z1)ΔxΔz=Dy(t0Δt/2,x1,y1,z1)ΔxΔz+[Hx(t0,x1,y1,z1+Δz/2)ΔxHz(t0,x1+Δx/2,y1,z1)ΔzHx(t0,x1,y1,z1Δz/2)Δx+Hz(t0,x1Δx/2,y1,z1)Δziy(t0,x1,y1,z1)ΔxΔz]Δt+O(Δt)3,E34

where x1y1z1is 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+Δξ/2dξfξ=fξ0Δξ+124d2fξ0dξ2Δξ3+OΔξ5,E35
ξ0Δξ/2ξ0+Δξ/2η0Δη/2η0+Δη/2dηgξη=gξ0η0ΔξΔη+1242gξηξ2Δξ3Δη+2gξηη2ΔξΔη3+OΔξΔη6.E36

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

Bzt0+Δtx0y0z0ΔxΔy+124x2Bzt0+Δtx0y0z0Δx3Δy+y2Bzt0+Δtx0y0z0ΔxΔy3=Bzt0x0y0z0ΔxΔy+124x2Bzt0x0y0z0Δx3Δy+y2Bzt0x0y0z0ΔxΔy3Eyt+Δt/2x0+Δx/2y0z0Δy+124y2Ey(t+Δt/2x0+Δx/2y0z0)Δy3Ext+Δt/2x0y0+Δy/2z0Δx124x2Ext+Δt/2x0y0+Δy/2z0Δx3Eyt+Δt/2x0Δx/2y0z0Δy124y2Eyt+Δt/2x0Δx/2y0z0Δy3+Ext+Δt/2x0y0Δy/2z0Δx+124x2Ex(t+Δt/2x0y0Δy/2z0)Δx3Δt+OΔt3.E37

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

Dyt0+Δt/2x1y1z1ΔxΔz+124x2Dyt0+Δt/2x1y1z1Δx3Δy+y2Dyt0+Δt/2x1y1z1ΔxΔy3=Dyt0Δt/2x1y1z1ΔxΔz+124x2Dyt0Δt/2x1y1z1Δx3Δy+y2Dyt0Δt/2x1y1z1ΔxΔy3+Hxt0x1y1z1+Δz/2Δx+124x2Hx(t0x1y1z1+Δz/2)Δx3Hzt0x1+Δx/2y1z1Δz124z2Hzt0x1+Δx/2y1z1Δz3Hxt0x1y1z1Δz/2Δx124x2Hxt0x1y1z1Δz/2Δx3+Hzt0x1Δx/2y1z1Δz+124z2Hzt0x1Δx/2y1z1Δz3iyt0x1y1z1ΔxΔz124x2iy(t0x1y1z1)Δx2Δz124z2iy(t0x1y1z1)ΔxΔz2Δt+OΔt3.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Δξ=2fξ+d2fξ0dξ2Δξ2+OΔξ4,E39

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

56Bzt0+Δtx0y0z0+124Bzt0+Δtx0+Δxy0z0+Bz(t0+Δtx0Δxy0z0)+Bzt0+Δtx0y0+Δyz0+Bz(t0+Δtx0y0Δyz0)xΔy=56Bzt0x0y0z0+124Bzt0x0+Δxy0z0+Bz(t0x0Δxy0z0)+Bz(t0x0y0+Δyz0)+Bz(t0x0y0Δyz0)ΔxΔy1112Eyt+Δt/2x0+Δx/2y0z0+124Eyt+Δt/2x0+Δx/2y0+Δyz0+Ey(t+Δt/2x0+Δx/2y0Δyz0)Δy1112Ext+Δt/2x0y0+Δy/2z0+124Ext+Δt/2x0+Δxy0+Δy/2z0+Ex(t+Δt/2x0Δxy0+Δy/2z0)Δx1112Eyt+Δt/2x0Δx/2y0z0+124(Ey(t+Δt/2x0Δx/2y0+Δyz0)+Ey(t+Δt/2x0Δx/2y0Δyz0))Δy+1112Ext+Δt/2x0y0Δy/2z0+124Ext+Δt/2x0+Δxy0Δy/2z0+Ex(t+Δt/2x0Δxy0Δy/2z0)ΔxΔt+OΔt3.E40

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

56Dyt0+Δt/2x1y1z1+124(Dy(t0+Δt/2x1+Δxy1z1)+Dy(t0+Δt/2x1Δxy1z1)+Dyt0+Δt/2x1y1z1+Δz+Dy(t0+Δt/2x1y1z1Δz))ΔxΔz=56Dyt0Δt/2x1y1z1+124(Dy(t0Δt/2x1+Δxy1z1)+Dy(t0Δt/2x1Δxy1z1)+Dyt0Δt/2x1y1z1+Δz+Dy(t0Δt/2x1y1z1Δz))ΔxΔz+1112Hxt0x1y1z1+Δz/2+124(Hx(t0x1+Δxy1z1+Δz/2)+Hx(t0x1Δxy1z1+Δz/2))Δx1112Hzt0x1+Δx/2y1z1+124Hzt0x1+Δx/2y1z1+Δz+Hz(t0x1+Δx/2y1z1Δz)Δz1112Hxt0x1y1z1Δz/2+124Hxt0x1+Δxy1z1Δz/2+Hx(t0x1Δxy1z1Δz/2)Δx+1112Hzt0x1Δx/2y1z1+124Hzt0x1Δx/2y1z1+Δz+Hz(t0x1Δx/2y1z1Δz)Δz56iyt0x1y1z1+124iyt0x1+Δxy1z1+iy(t0x1Δxy1z1)+iy(t0x1y1z1+Δz)+iy(t0x1y1z1Δz)ΔxΔzΔt+OΔt3.E41

Note that points xi±Δxyiziand xiyizi±Δzare at adjacent cells to the cell including xiyizi. 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σx0y0x0+mΔxy0+nΔyBzt0+Δtx0+mΔxy0+nΔyz0,E42
m,nσx1z1x1+mΔxy1z1+nΔzDyt+Δt/2x1+mΔxy1z1+nΔz,E43

where

σξηξ+mΔξη+nΔη=56δm,0δn,0+124δm,1δn,0+δm,1δn,0+δm,0δn,1+δm,0δn,1,E44

and δp,qis the Kronecker delta defined as

δp,q=1p=q0pq.E45

The inverse operator “σ1” is defined as

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

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

Bzt0+Δtx0y0z0ΔxΔy=Bzt0x0y0z0ΔxΔym,nσ1x0y0x0+mΔxy0+nΔy×1112Ey(t0+Δt/2x0+m+12Δxy0+nΔyz0)Ey(t0+Δt/2x0+m12Δxy0+nΔyz0)ΔyEy(t0+Δt/2x0+mΔxy0+n+12Δyz0)Ey(t0+Δt/2x0+mΔxy0+n12Δyz0)Δx+124Ey(t0+Δt/2x0+m+12Δxy0+n+1Δyz0)+Ey(t0+Δt/2x0+m+12Δxy0+n1Δyz0)Ey(t0+Δt/2x0+m12Δxy0+n+1Δyz0)Ey(t0+Δt/2x0+m12Δxy0+n1Δyz0)ΔyEx(t0+Δt/2x0+m+1Δxy0+n+12Δyz0)+Ex(t0+Δt/2x0+m+1Δxy0+n12Δyz0)Ex(t0+Δt/2x0+m1Δxy0+n+12Δyz0)Ex(t0+Δt/2x0+m1Δxy0+n12Δyz0)ΔxΔt+OΔt3.E47

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

Dyt0+Δt/2x1y1z1ΔxΔz=Dyt0Δt/2x1y1z1+iy(t0x1y1z1)ΔxΔz+m,nσ1x1z1x1+mΔxz1+nΔz×1112Hx(t0x1+mΔxy1z1+n+12Δz)Hx(t0x1+mΔxy1z1+n12Δz)ΔxHz(t0x1+m+12Δxy1z1+nΔz)Hz(t0x1+m12Δxy1z1+nΔz)Δz+124Hx(t0x1+m+1Δxy1z1+n+12Δz)+Hx(t0x1+m1Δxy1z1+n+12ΔzHx(t0x1+m+1Δxy1z1+n12Δz)Hx(t0x1+m1Δxy1z1+n12Δz)ΔxHz(t0x1+m+12Δxy1z1+n+1Δz)+Hz(t0x1+m+12Δxy1z1+n1ΔzHz(t0x1+m12Δxy1z1+n+1Δz)Hz(t0x1+m12Δxy1z1+n1Δz))ΔzΔt+OΔt3.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 ncoand ncl, respectively. The core region is extended infinitely in the y- and z-directions and has width din 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 ydependence. Because the system has no ydependence, it is essentially a two-dimensional system. An analytical solution is known and is derived in the appendix. The solution is

Exanltxyz=βh0ωnco2ε0cos2uxdsinωtβzxd2βh0ωncl2ε0cosuew2xddsinωtβzx>d2,E49
Eyanltxyz=0,E50
Ezanltxyz=2uh0ωnco2ε0dsin2uxdcosωtβzxd22wh0ωncl2ε0dsignxsinuew2xddcosωtβzx>d2,E51
Hxanltxyz=0,E52
Hyanltxyz=h0cos2uxdsinωtβzxd2h0cosuew2xddsinωtβzx>d2,E53
Hzanltxyz=0,E54

where “anl” indicates that this is an analytical solution. uand wsatisfy

w=nclnco2utanu,E55
v2=u2+w2=nco2ncl2d2π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πnco2w2+ncl2u2.E57

vdefined by Eq. (56) is called the V-parameter, which is determined by the parameters defining the system. uand ware determined using Figure 5 . In the figure, red curves represent Eq. (55) which is symmetric under the parity transformation xxas

Hytxyz=Hytxyz.E58

Brown curves represent

w=ncl2nco2ucotu,E59

which is antisymmetric under the parity transformations xxas

Hytxyz=Hytxyz.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 uis called fundamental mode. The number of modes of the system is determined by Vand increases by one with respect to each π/2.

In the computational methods, the system parameters are set with the wavelength λas 0.30m, the core width das 0.30m, the same with the wavelength, the core index ncoas 2.0, and the cladding index nclas 1.0. The lengths of the cell edges Δxand Δzare both λ/20, and the time step Δtis 1012s. 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 ncland nco. These parameter values show that the solution is the fundamental mode. A magnetic field is excited at z=0as

Hytx,0,0=h0cos2uxdsinωtxd2h0signxcosuew2xddsinωtx>d2,E65

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

Figures 6 14 are numerical and analytical results at times at which the ωtvalues are integer multiples of 2π. In these figures, violet curves represent Hy/h0, 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.0ns. 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

Hzis zero. Moreover, the differences in the results between the FDTD calculations and those of the analytical ones make it clear that Hy/h0at 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 errtdefined as

errt=p,qHynumtpΔx0qΔzHyanltpΔx0qΔznpΔx0qΔz2p,qHyanltpΔx0qΔz2npΔx0qΔz2,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 xyplane per unit length of the y-direction.

Figure 15 shows the errfunctions 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 errfunction 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 errfunction, 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

Etxyz=exeiωtβz,E68
Htxyz=hxeiωtβz.E69

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

iμ0ωhxx=eyxE70
iμ0ωhyx=exx+xezx,E71
iμ0ωhzx=xeyx,E72
inx2ε0ωexx=hyx,E73
inx2ε0ωeyx=hxxxhzx,E74
inx2ε0ωezx=xhyx,E75

where nxis 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 Hyare shown, and they are TM modes.

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

x2hyx=2πnxλ2β2hyx,E76
exx=βnx2ε0ωhyx,E77
ezx=1inx2ε0ωxhyx.E78

The discontinuity of nxat x=d/2and x=d/2, shown in Figure 16 , requires that the boundary condition at x=±d/2be 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, Hyxand xhyxin Eqs. (76)(78) are continuous.

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

1. β<ncl

2. nclβ<nco

3. ncoβ

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

hyx=QcosPxxd/2QcosPd2cosQd2+PsinPd2sinQd2cosQxPsinPd2cosQd2QcosPd2sinQd2sinQxxd/2,E79

and

hyx=QsinPxxd/2signx{QsinPd2cosQd2PcosPd2sinQd2cosQx+PcosPd2cosQd2+QsinPd2sinQd2sinQx}xd/2,E80

where

P=2πncoλ2β2,E81
Q=2πnclλ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 xd/2, Hytxyzis a linear combination of

eiωtPxβz,andeiωt+Pxβz,E83

a plane wave whose wavenumber is P0βand P0β, respectively. When x>d/2, Hytxy.zis a linear combination of

eiωtQxβz,andeiωt+Qxβz,E84

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

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

hyx=cosPxxd/2cosPd2eQ2xdd,E85

where

P=2πncoλ2β2,E86
Q=β22πnclλ2,E87
Qd2=nclnco2Pd2tanPd2,E88

and

hyx=sinPxxd/2signxsinPd2eQ2xddx>d/2,E89

where

Qd2=nclnco2Pd2cotPd2.E90

Defining

u=Pd2,E91
w=Qd2,E92

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

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

hyx=coshPxxd/2coshPd2eQxd/2x>d/2,E93

where

P=β22πncoλ2,E94
Q=β22πnclλ2,E95
Q=PtanhPd2,E96

and

hyx=sinhPxxd/2signxsinhPd2eQxd/2x>d/2,E97

where

Q=PcothPd2.E98

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

P2+Q2=2πncoλ22πnclλ2,E99

is satisfied. When Pis pure imaginary and Qis real, the solution reduces to case 2. There is a possibility that Qis neither real nor purely imaginary. However, such a solution must be attenuated when zbecomes large. Mathematically, there can be a divergent solution when zbecomes 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: July 9th, 2018 Reviewed: September 5th, 2018 Published: November 27th, 2018