The finite-difference time-domain (FDTD) algorithm is one of the most popular numerical methods to solve Maxwell’s equations, as proposed by Yee . It is a time marching forward method which easily makes a visual representation of the fields. It also scales very well since the number of computation required is proportional to the size of the model and the method requires no large-scale linear algebra computation . The FDTD method requires a completely structured grid also known as the Yee grid. However, Many real world problems have complicated geometries and numerical inaccuracy happens when the modeled objects do not fit the grid very well.
In this paper, we focus on two-dimensional transverse electric (TE) Maxwell’s equation within multiple non-magnetic media where the electric fields are not continuous across the media interface so the numerical schemes need to be designed to take care of such discontinuity. The other case is the transverse magnetic equation where the electric field is always continuous across the media interface and the original algorithm will work as in the homogeneous medium.
The original FDTD algorithm proposed in Yee’s paper assigns each of the field component an electric permittivity and a magnetic permeability solely based on the material properties at its location. This is commonly known as the staircasing or pixelation method. In general, the staircasing method has an error that scales with given a grid size of for cells that are homogenous. However, in cells that contain a medium interface, the local error becomes . Even though the number of cells that contain an interface is often a small fraction of the total cells, these local errors can cause a rough global error of .
Subsequent research in modern optics and electromagnetism dealing with the Maxwell’s equations have focused on more and more complex geometries with different materials so that the complicated dielectric tensor can cause serious pixelation issues. In turn, there have been many improvements in FDTD algorithms to overcome these pixelation issues. One way to overcome the pixelation problem is by fitting the mesh to the device, thereby generating the so-called non-orthogonal FDTD[3, 4]. The second method is to refine the Cartesian mesh around the interfaces and is thus called subgridding method . Even though these two methods converge faster than the standard FDTD method, their complexity are larger than original FDTD method and they also have time stability issues.
In order to keep a completely structured grid and maintain all the benefits of FDTD, another method redefines Maxwell’s integral equations around the interface of the media. This is known as the effective-permittivity (EP) method. This method uses sub-pixel smoothing techniques to change the permittivities of the field components around the interface to produce better results by taking account of many different factors, such as the interface conditions, permittivities of adjacent field components and so on. The revised effective permittivities smooth out the pixel error. The main benefit of these methods is that there is a negligible increase of the numerical load since all the revised permittivities have been calculated before the main FDTD algorithm loop.
The first method of sub-pixel smoothing is the volume average effective permittivity (V-EP) . This method assigns effective permittivities for all the field components in a cell that contains media interface. The effective permittivities are calculated by taking a weighted average based on how much percent of the volume each medium occupies. This method is simple for implementation and also stabilizes the error fluctuation that occurs in the staircasing method. However, it does not decrease the error rapidly and its error can be worse than staircasing method sometimes due to ignorance of the interface orientation.
The second method proposes dynamic formulas for the EP based on the orientation of the media interface and their expressions are accurate when the interface is perpendicular or parallel to the mesh axes [7, 8, 9]. These EPs improve the accuracy of the FDTD method, while keeping the simple structure of the original algorithm. However, their performance will deteriorate for a curved or flat interface not perpendicular or parallel to the mesh axes.
The third method uses the reflection coefficients to derive effective permittivities and shows the method can also achieve second-order convergence with several special slanted angles between the interface and the Yee grid rather than just orthogonal or parallel to the Yee grid [10, 11]. However, its derivation is quite complicated and has no ways to be extended to arbitrary interfaces.
In a landmark paper , Mohammadi et al. proposed the so-called contour-path effective permittivities (CP-EP) method. The CP-EP method first addresses the orthogonal and parallel cases as shown above, then extends the idea to handle with more general geometry. Unlike the previous methods, it could handle arbitrary boundary with any orientation of the boundary. Compared with traditional staircasing methods, CP-EP has almost no additional runtime cost since determination of the effective permittivities has been done before the main FDTD loop.
However, we will demonstrate that CP-EP does not incorporate one important term from the interface conditions into the FDTD algorithm. Then we will develop a revised EP algorithm where the missing term will be re-considered in order to produce a more accurate and stable algorithm. Since this method incorporates the boundary conditions very well, it will be denoted as BC-EP.
where effective dielectric permittivity tensors are constructed to take care of both isotropic and anisotropic materials so as to achieve second-order convergence under the open source software MEEP. However, such proposed method which satisfies the interface conditions for electromagnetic fields has been shown to have late-time instabilities, and many possibilities to average the effective dielectric tensor are explored to avoid late time instabilities. Further efforts have been reported in  to construct a new second-order scheme by taking the average of eight triplets. Although the new method is highly accurate, its effective dielectric permittivity tensor can still be asymmetric thus unstable for certain conditions. Therefore, the last improvement of the effective dielectric tensor has been given in  to make it symmetric and stable. The numerical results show that this scheme gives the best result in general and the error in practice still lies in between first and second-order in most cases.
In the following, we will introduce the new BC-EP method by adding a few terms to the established CP-EP method for the cells around the interface so as to keep the numerical load increase as small as possible and improve the accuracy significantly. Then the numerical tests verify that BC-EP has a much better performance than CP-EP, V-EP and staircasing methods while still maintaining numerical stability. BC-EP has a very simple structure and can be merged into any existing FDTD packages easily, compared with other established FDTD software [14, 17] which can be run only within their framework. We provide C++ source code for BC-EP together other three methods as supporting information for public access. Conclusion and some future research consideration will be provided in the end.
2. New Algorithm Design
In this section, we will propose a new effective permittivity method to solve Maxwell’s equations in a domain composed of two different non-magnetic media with different electric permittivities. To simply the algorithm description, we will focus on the most important transverse electric (TE) scenario. That is, is located in the incident xy-plane while is along z-axis. Moreover, we restrict the discussion to dielectric and non-magnetic media such that and where stays as a constant in different media and
. Based on the integral version of the Maxwell’s equations, we will derive the new algorithm for any orientation of the interface. This is a significant improvement over all current methods which are only accurate under specific angles between the interface and Cartesian coordinate system, such as parallel or orthogonal to Yee-axis.
As a preliminary step, let us set up the relation of and across the interface of these two media. It is well known that is continuous in the whole domain. Meanwhile, across the media interface, the tangential component of the electric field and the normal component of the electric flux are continuous as well, as shown in Figure 1 where
is the outward unit norm vector from region 1.
Choose . and are the electric fluxes at A from different media, and and are the corresponding electric fields.
Suppose is the unit norm vector at from region 1, then the corresponding tangent unit vector is Therefore, we can calculate and at point in different regions as follows.
Based on the fact that and , we have
Therefore, by using , we get
Furthermore, by using at , we obtain
If the interface is parallel or orthogonal to -axis, either or so . Then the relation (5) involving and will be simplified such that and will only depend on or , respectively. Same holds for and in relation (6).
By using (6), we can express in region 2 based on its neighbor in region 1 and its corresponding in the same region 2. Similar formula can also be derived for based on and .
Now let us express Ampere’s law in integral form
and Faraday’s law integral form
Firstly, let us discretize Ampere’s law (9) around the interface to update and as below.
To update , we choose in (9) to be the rectangle in -plane whose projection in -plane is the line segment . Suppose the interface intersect at point and as in Figure 2. Applying time discretization to Ampere’s law along at time yields
Based on formula (7), we can compute the integrals of left hand side as follows. From now on, all the higher-order terms of and are ignored for the clarification of the formulas.
which, based on in , becomes
which shows improved discretization of Ampere’s law for updating
It should also be noted that the last term on the right hand side of (16) have been ignored by CP-EP for algorithm simplification but it indeed provides necessary corrections for the interface conditions across multiple media. Furthermore, the numerical tests have also verified that in order to make the algorithm for Ampere’s Law more accurate, these two values in the last term should stay close to the interface. That is, if the normal direction of the interface satisfies as shown in Figure 2, we use since these two points are closer to the interface, while for the other case, we use in (16) instead since these two nodes are closer to the interface than the other pair. This improvement is feasible since in continuous across the interface.
Similarly, update of at point by Ampere’s law (9) can be taken care analogously as shown in Figure 3, where the interface intersects the line segment , projection of at -plane, at . Let . By repeating the same procedure as above while based on formula (8), we can also obtain the discretization of Ampere’s law for as below:
where in the last term will be replaced by when holds.
Secondly, let’s discretize Faraday’s law (10) around the interface .
We choose to be the rectangle centered at . If the medium interface does not cut through interior of , the standard Yee scheme should apply without modification. So we only consider the case where intersects with . It’s obvious that will only intersect with at most two sides, e.g. the bottom and left side as in Figure 4. All other cases can be handled in the same way. Similarly as for Ampere’s law, all the higher-order terms of dx and dy are ignored for the clarification of the formulas. Discretization of the left hand side term of (10) at the time via finite difference scheme yields
since is always perpendicular to the incident -plane and so is always continuous across the interface, then by using midpoint rule,
However, the situation for is much more complicated. In Figure 4, the line integral of along the right and top sides can be handled as Standard Yee scheme since they are completely inside . But the line integrals along the bottom and left sides need to be handled separately.
Suppose the portion OA of the bottom side in has length , and the portion OB of the left side in has length as in Figure 4, then the first two terms on the right hand side of (21) can be handled as below. For the first term on the right hand side of (21), we have
By using formula (5) with and interchanged, we have
Putting it into (22) yields
where is the value of evaluated at point in . However, only instead of is calculated at the mesh point . Therefore, we will take the average value of at the neighboring pair of mesh points and if they are both located within , or take the average value of at the other neighboring pair mesh points and if they are both located within . If none of the pairs are located within , this term will be set to zero as in CP-EP. Therefore, BC-EP will try to make corrections for the missing terms in CP-EP if possible.
Similarly, we can calculate the second term on the right hand side of (21). Since the node is located within instead of , so the new should be negative of the old used above. Therefore, applying integration midpoint rule, Taylor expansion and formula (5) where replaced by and interchanged, we have
where the last two terms are expressed by
If we set in the new algorithm (16), (17) for Amepre’s law and (27) for Faraday’s law, we can retrieve the standard FDTD scheme in homogeneous medium. Therefore, BC-EP is a simple extension of the original schemes in order to take care of the media interface. Meanwhile, the extra terms added in the formulas will greatly improve the accuracy compared with CP-EP. Furthermore, BC-EP has a very simple structure and can be merged into any FDTD software by just revising the effective permittivities and adding some extra terms when necessary.
It should be mentioned that due to complicated expressions of those ignored higher-order terms in the newly derived formulas, it is difficult to demonstrate the convergence order of BC-EP theoretically. However, subsequent numerical results will show that BC-EP can achieve the highest convergence compared with other methods. Another point is that it is still an open challenge to accurately handle objects with sharp corners, where the resulting field singularities are known to degrade the accuracy of all numerical schemes . Therefore, just as other methods, BC-EP performs well for straight or curved interfaces instead of interfaces with corners.
3. Numerical Results
In this section, we will demonstrate convergence order and stability of BC-EP numerically. To accomplish this task, we implement the algorithm to solve the Maxwell’s equations on a dielectric cylinder, together with the staircasing, V-EP, and CP-EP methods. All the implementations are set up for the 2D FDTD TE case. The total scattering cross sections (SCS) are calculated by all four methods and then compared with the well-known analytic solution by Mie Theory  so as to measure the accuracy of each method.
The dielectric cylinder is simulated in an area that is times the radius of the cylinder plus some extra space for PML condition. The cylinder is centered at where and are the number of the cells along the and directions, and are the horizontal and vertical mesh sizes of each cell as seen in Figure 5. We run the test for two cases: one case for a larger cylinder with radius and another case for a smaller cylinder with . The total-field/scattered-field line is three times the length of the cylinder radius away from the center, and the integration line for calculating SCS is four times away from the center, while the PML starts at five times away from the center.
The source wave is a planar wave in the direction with a Gaussian envelope , where is six times , is speed of light in a vacuum, and is the testing range of wavelengths. In our simulation, we calculate SCS for wavelengths ranging from to
of visual light spectrum with equal distance. Fast Fourier transform (FFT) is then applied to transform the electric and magnetic fields from spatial domain into light frequency domain so as to calculate SCS over all related wavelengths.
We use square unit cells with to simplify the simulation. The time step is set to where to ensure stability . Also for stability concerns, were chosen to divide the smallest tested wave length () by at least twenty times. The total SCS is calculated by using the Poynting Vector along with the integration line as seen in Figure 5. Meanwhile, Mie theory is used to calculate the true value of total SCS for each wavelength. The number of iterations for the main FDTD loop is set large enough in order to give enough time for the electric and magnetic waves to leave the simulated area for stable SCS calculation.
Firstly, we investigate the accuracy and the convergence order of BC-EP and compare it with other similar methods. We test all the methods for a circular cylinder with radius and permittivity . The background media ia always set as air with . The mesh size is originally set to . Each method has been run to calculate the total SCS with different mesh sizes. As seen in (A) of Figure 6, the numerical values from all methods match with the true solution very well. But we can still see that BC-EP has the smallest relative error while V-EP has the largest relative error on most wavelengths shown in (B) of Figure 6. To make a reasonable comparison, for each mesh size , the average relative error of SCS is calculated among all the wavelengths for each method. In order to observe the error convergence order clearly, the average relative errors versus the mesh sizes are converted into a log scale where denotes the log of the number of mesh points for the radius of the cylinder, as seen in Figure 7.
This whole process is then repeated for other larger permittivity values . The numerical results for the average relative errors versus the mesh sizes are reported in Figure 8 for , and in Figure 9 for . Other values have similar outcomes.
It can be seen that BC-EP gives the least amount of error for most cases, followed by CP-EP. Staircasing results are quite erratic and its error doesn’t go down in a uniform manner due to random nature of some mesh sizes conforming with the cylinder better than others. V-EP has the worst error but gives a uniform error decrease with smaller mesh sizes, and thus gives more consistent results over staircasing.