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.
- Maxwell’s equations
- integral form
- finite-difference time-domain method
- the lowest-order approximation
- next-to-the-lowest-order approximation
- computational method
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 and are electric and magnetic fields, respectively, and are electric and magnetic flux densities, respectively, is current density, is charge density, is the time derivative of field , is the rotation of vector field , is the divergence of vector , 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 , where 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 . 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
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 .
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 using the magnetic field at and the electric field at by applying Eq. (1), where is a particular time and is the time step. Let us now consider the top surface of the cell. At the center of the surface whose coordinates are , Eq. (1) becomes
where the variables , , and of and represent the -, -, and -components of the and fields, respectively, and and represent the partial derivatives in the - and -directions, respectively. Replacing the partial derivatives by the central differences yields
where means . and at can also be derived similarly. Usually, the field can be derived from the field. For example, in vacuum, air, or a dielectric
where is the vacuum permeability with the value . In a magnetic material, the relationship between and often becomes nontrivial, but this is beyond the scope of this book. However, in small and 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, is 1.
The cyan cell is used to calculate the electric field at using the electric field at and the magnetic field at 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 , Eq. (2) becomes
Replacing partial derivatives with central differences yields
where . Then, is derived as follows:
and at can also be derived similarly. Typically, the field can be derived from the flux density. For example, in vacuum, air, or magnetic material
where is the vacuum permittivity, whose value is . In a dielectric material, the relationship between and often becomes nontrivial, but this is beyond the scope of this book. However, in small and 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 and fields are given which in turn satisfy Eqs. (3) and (4). When the field distribution at and the field distribution at are known, the field at is calculated using Eqs. (20) and (22), given the field at and the field at . The field at is calculated using Eqs. (12) and (14), given the field at , having determined the field at . If the time is less than , then the time becomes and the flow repeats. If the time exceeds , 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 . 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.
where is a compact and connected surface, is the boundary curve of the surface, is a surface element normal to the surface, and is a line element parallel to the curve. Moreover, integrating both sides of Eq. (25) over from to and those of Eq. (26) over from to yields
In general, when is analytical in a region including , the following relationship is satisfied:
and when is analytical in a region including and , the following relationship is satisfied:
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 is the top surface of the yellow cell in Figure 1 , Eq. (27) is approximated as
where 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.
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
Note that points and are at adjacent cells to the cell including . 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.
and is the Kronecker delta defined as
The inverse operator “” is defined as
Using enables Eq. (40) to be rewritten as
In addition, Eq. (40) can also be rewritten as
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 and , respectively. The core region is extended infinitely in the - and -directions and has width in the -direction. The cladding region is the rest of space. An electromagnetic wave propagates in the -direction, and its electromagnetic field is assumed to have no dependence. Because the system has no dependence, it is essentially a two-dimensional system. An analytical solution is known and is derived in the appendix. The solution is
where “anl” indicates that this is an analytical solution. and satisfy
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
defined by Eq. (56) is called the V-parameter, which is determined by the parameters defining the system. and are determined using Figure 5 . In the figure, red curves represent Eq. (55) which is symmetric under the parity transformation as
Brown curves represent
which is antisymmetric under the parity transformations as
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 is called fundamental mode. The number of modes of the system is determined by and increases by one with respect to each .
In the computational methods, the system parameters are set with the wavelength as m, the core width as m, the same with the wavelength, the core index as , and the cladding index as . The lengths of the cell edges and are both , and the time step is 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 and . These parameter values show that the solution is the fundamental mode. A magnetic field is excited at as
Figures 6 – 14 are numerical and analytical results at times at which the values are integer multiples of . In these figures, violet curves represent , 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 , , , and 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
is zero. Moreover, the differences in the results between the FDTD calculations and those of the analytical ones make it clear that 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 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 plane per unit length of the -direction.
Figure 15 shows the 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 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.
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 function, and the corrected method was found to be more accurate and reliable than the original.
The author would like to thank Mr. A. Okabe, who is the collaborator of reference , 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 -dependence in the -direction with angular frequency and propagation constant , the wave number in the -direction, is written as
where 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 are shown, and they are TM modes.
The discontinuity of at and , shown in Figure 16 , requires that the boundary condition at 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, and in Eqs. (76)–(78) are continuous.
Solving Eq. (76) requires considering the following three cases:
In case 1, the solutions of Eq. (76) are
However, this solution does not show electromagnetic wave propagation in the -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 , is a linear combination of
a plane wave whose wavenumber is and , respectively. When , is a linear combination of
a plane wave whose wavenumbers are and , respectively. Therefore, with a suitable linear combination of Eqs. (79) and (80), the solution becomes that a plane wave with wavenumber is incident from the region to be reflected and transmitted by a film of the region with a reflected wave propagated in the region and a transmitted wave propagated in the region. Therefore, this solution is not what we want.
In case 2, the solutions of Eq. (76) are
In case 3, the solutions of Eq. (76) are
is satisfied. When is pure imaginary and is real, the solution reduces to case 2. There is a possibility that is neither real nor purely imaginary. However, such a solution must be attenuated when becomes large. Mathematically, there can be a divergent solution when becomes large, but such a solution cannot conserve energy. Therefore, such a solution is not what we want.