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

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.


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 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: 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 iω, 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 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.

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 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 xand y-directions, respectively. Replacing the partial derivatives by the central differences yields 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: where O Δx; Δy Þ . 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 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: where μ depends on the material. Often, a value 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 ; Replacing partial derivatives with central differences yields where t ¼ t 0 . Then, D y t 0 ; x 1 ; y 1 ; z 1 À Á is derived as follows: 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 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 where ε depends on the material. Often a value called the relative permittivity and a value 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.

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 where S is a compact and connected surface, ∂S is the boundary curve of the surface, da is a surface element normal to the surface, and ds 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 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: and when g is analytical in a region including ξ 0 À Δξ=2 ≤ ξ ≤ ξ 0 þ Δξ=2 and η 0 À Δη=2 ≤ η ≤ η 0 þ Δη=2, the following relationship is satisfied: The lowest-order approximations of Eqs. (29) and (30) are, respectively, 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 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 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-lowestorder 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, 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 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 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 for any function f is applied. Applying Eq. (39) to Eq. (37) yields Applying Eq. (39) to Eq. (38) yields Note that points x i AE Δx; y i ; z i À Á and x i ; y i ; z i AE Δ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 where and δ p, q is the Kronecker delta defined as The inverse operator "σ À1 " is defined as Using σ À1 enables Eq. (40) to be rewritten as In addition, Eq. (40) can also be rewritten as 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.

Numerical results
In this section, numerical results of electromagnetic wave transmission in a twodimensional 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 where "anl" indicates that this is an analytical solution. u and w satisfy w ¼ n cl n co 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 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 ð Þ : Brown curves represent which is antisymmetric under the parity transformations x↦ À x as H y t; Àx; y; z ð Þ¼À H y t; x; y; z ð Þ : (60) Figure 5. 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:30m, the core width d as 0:30m, 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 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        FDTD method, the middle column is calculated using the corrected FDTD method, and the right-hand column is the analytical solution wherein the region 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   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.

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.
The discontinuity of n x ð Þ at x ¼ Àd=2 and x ¼ d=2, shown in Figure 16, requires that the boundary condition at x ¼ AEd=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. Defining where P ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Q ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi where Q ¼ ÀP coth Pd 2 : As a result of Eqs. (94) and (95), 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.

Author details
Naofumi Kitsunezaki Amashiro Science, Nagoya, Japan *Address all correspondence to: n.kitsunezaki@amashiro-science.jp © 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.