{'fraction_non_alphanumeric': 0.03300859127718173, 'fraction_numerical': 0.01311300201422288, 'mean_word_length': 4.362133568437294, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An enhanced finite difference time domain method for two dimensional Maxwell's equations

An efficient finite-difference time-domain (FDTD) algorithm is built to solve the transverse electric 2D Maxwell's equations with inhomogeneous dielectric media where the electric fields are discontinuous across the dielectric interface. The new algorithm is derived based upon the integral version of the Maxwell's equations as well as the relationship between the electric fields across the interface. It is an improvement over the contour-path effective-permittivity algorithm by including some extra terms in the formulas. The scheme is validated in solving the scattering of a dielectric cylinder with exact solution from Mie theory and is then compared with the above contour-path method, the usual staircase and the volume-average method. The numerical results demonstrate that the new algorithm has achieved significant improvement in accuracy over the other methods. Furthermore, the algorithm has a simple structure and can be merged into any existing FDTD software package very easily.


Introduction
The finite-difference time-domain (FDTD) algorithm is one of the most popular numerical methods to solve Maxwell's equations, as proposed by Yee [1]. 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 [2]. 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 O(∆x 2 ) given a grid size of ∆x for cells that are homogenous. However, in cells that contain a medium interface, the local error becomes O (1). 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 O(∆x).
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 [5]. 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) [6]. 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 [12], 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, CE-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.
In fact, such deficiency of CP-EP has also been identified in [13,14] 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 [15]. Further efforts have been reported in [16] 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 [17] 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. Conclusion and some future research consideration will be provided in the end.

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, E(x, y) = E x i + E y j is located in the incident xy-plane while H(x, y) = H z k is along zaxis. Moreover, we restrict the discussion to dielectric and non-magnetic media such that D(x, y) = ǫE(x, y) and B(x, y) = H(x, y) where ǫ stays as a constant in different media and ǫ 0 = µ 0 = 1. 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 D(x, y) and E(x, y) across the interface of these two media. It is well known that H(x, y) is continuous in the whole domain. Meanwhile, across the media interface, the tangential component ) are the corresponding electric fields. Suppose n = (n x , n y ) is the unit norm vector at A from region 1, then the corresponding tangent unit vector is τ = (n y , −n x ). Therefore, we can calculate E τ and D n at point A in different regions as follows. Meanwhile, Based on the fact that E τ1 = E τ2 and D n1 = D n2 , we have Therefore, by using n 2 x + n 2 y = 1, we get If the interface Γ is parallel or orthogonal to x-axis, either n y = 0 or n x = 0 so n x · n y = 0. Then the relation (5) involving E 1 and E 2 will be simplified such that E x1 and E y1 will only depend on E x2 or E y2 , respectively. Same holds for D 1 and D 2 in relation (6).
By using (6), we can express D x2 in region 2 based on its neighbor D x1 in region 1 and its corresponding D y2 in the same region 2. Similar formula can also be derived for D y2 based on D y1 and D x2 .
(5) will be used to derive the new scheme for Faraday's law while (7) and (8) will be used to derive the new scheme for Ampere's law. 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 E n x and E n y as below. To update E x , we choose S in (9) to be the rectangle in yz- Figure 2. Discretization of Ampere's law across Γ plane whose projection in xy-plane is the line segment OC. Suppose the interface Γ intersect OC at point A and OA = f as in Figure 2. Applying time discretization to Ampere's law along S at time t n+ 1 2 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 dx and dy are ignored for the clarification of the formulas. Similarly, Furthermore, By putting (12), (13) and (14) into (11), we obtain which shows improved discretization of Ampere's law for updating E x . 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 H z values in the last term should stay close to the interface. That is, if the normal direction of the interface n satisfies n x · n y > 0 as shown in Figure 2, we use (H z (i−1,j) ) since these two H z points are closer to the interface, while for the other n x · n y ≤ 0 case, we use (H (16) instead since these two H z nodes are closer to the interface than the other pair. This improvement is feasible since H z in continuous across the interface.

Hz Hz
Ey  (9) can be taken care analogously as shown in Figure 3, where the interface Γ intersects the line segment OC, projection of S at xy-plane, at A. Let OA = d. By repeating the same procedure as above while based on formula (8), we can also obtain the discretization of Ampere's law for E y as below: where (H S to be the rectangle centered at (i, j). If the medium interface Γ does not cut through interior of S, the standard Yee scheme should apply without modification. So we only consider the case where Γ intersects with S. 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 t n via finite difference scheme yields since H is always perpendicular to the incident xy-plane and so is always continuous across the interface, then by using midpoint rule, However, the situation for E is much more complicated. In Figure 4, the line integral of E n along the right and top sides can be handled as Standard Yee scheme since they are completely inside Ω 2 . But the line integrals along the bottom and left sides need to be handled separately.
∆x (21) Suppose the portion OA of the bottom side in Ω 1 has length d, and the portion OB of the left side in Ω 1 has length f 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 Ω 1 and Ω 2 interchanged, we have Putting it into (22) yields where E n is the value of E n y evaluated at point (i, j − 1 2 ) in Ω 1 . However, only E x instead of E y is calculated at the mesh point (i, j − 1 2 ). Therefore, we will take the average value of E n y at the neighboring pair of mesh points (i − 1 2 , j − 1) and (i + 1 2 , j) if they are both located within Ω 1 , or take the average value of E n y at the other neighboring pair mesh points (i − 1 2 , j) and (i + 1 2 , j − 1) if they are both located within Ω 1 . If none of the pairs are located within Ω 1 , 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 E y node (i − 1 2 , j) is located within Ω 2 instead of Ω 1 , so the new n should be negative of the old n used above. Therefore, applying integration midpoint rule, Taylor expansion and formula (5) where (n x , n y ) replaced by (−n x , −n y ) and ǫ 1 , ǫ 2 interchanged, we have will be taken care similarly as E n y (i,j− 1

)
from the first term. By putting (24), (25) into (21), we get Finally, by putting (19), (20) and (26) into (18) and dividing both sides by ∆x∆y, we obtain the new discretization of Faraday's law: where the last two terms are expressed by If we set ǫ 1 = ǫ 2 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 (close to second order) 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 [14]. Therefore, just as other methods, BC-EP performs well for straight or curved interfaces instead of interfaces with corners.

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 [18] so as to measure the accuracy of each method. Figure 5. The numerical test setup using a dielectric cylinder, giving enough space for the total-field/scattered-field and integration line before the PML The dielectric cylinder is simulated in an area that is 10 times the radius of the cylinder plus some extra space for PML condition. The cylinder is centered at (( Nx 2 + 0.5)∆x, ( Ny 2 + 0.5)∆y) where N x and N y are the number of the cells along the x and y directions, ∆x and ∆y 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 r = 400nm and another case for a smaller cylinder with r = 150nm. The total-field/scattered-field line is three times the length of the cylinder radius r away from the center, and the integration line for calculating SCS is four times r away from the center, while the PML starts at five times r away from the center.
The source wave is a planar wave in the E y direction with a Gaussian envelope e −(t−t 0 ) 2 2 · 1 2π(c/λ 0 −c/λ 1 ) 2 · cos(2πct), where t 0 is six times ∆t, c is speed of light in a vacuum, and [λ 0 , λ 1 ] is the testing range of wavelengths. In our simulation, we calculate SCS for 601 wavelengths ranging from λ 0 = 400nm to λ 1 = 1000nm 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 ∆x = ∆y to simplify the simulation. The time step is set to ∆t = S∆x c where S = 0.98 √ 3 to ensure stability [2]. Also for stability concerns, ∆x were chosen to divide the smallest tested wave length (λ 0 = 400nm) 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 mthods for a circular cylinder with radius r = 400nm and permittivity ǫ = 3. The background media ia always set as air with ǫ = 1. The mesh size is originally set to ∆x = ∆y = 10nm. Each method has been run to calculate the total SCS with different mesh sizes. As seen in Figure 6, the relative errors between the numerical values from all methods and the true solution are quite small. But we can still see that BC-EP has the smallest error while V-EP has the largest error. Meanwhile, for each mesh size ∆x, 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 N λ 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 ǫ = 6, 10. The numerical results are reported in Figure 8 and Figure 9 for for ǫ = 6, and in Figure 10 and Figure 11 for ǫ = 10. 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.
To further investigate the order of convergence of all the methods, we estimate the convergence order by computing the slopes from the log figures for all the given permittivity. Table 1 shows BC-EP converges significantly faster than all the other methods. It approaches second-order in most cases. CE-EP is the second in terms of convergence order and outperforms both staircasing and V-EP as demonstrated in [12] already.
To complete the final evaluation of the performance for all the methods, the whole test suite is then repeated for another circular cylinder with smaller radius    Figure 11. The average relative error of a cylinder in log scale, ǫ = 10 and radius r = 400nm. N λ = log( 400nm ∆x ). r = 150nm. All the results are very similar and for brevity, we only report the case of ǫ = 6, as seen in Figure 12 and Figure 13. It can be seen that BC-EP has a higher convergence order than others. To be more specific, BC-EP, CP-EP, Staircasing and V-EP have the order of convergence 1.4005, 0.92869, 0.92802 and 1.0059, respectively, based on Figure 13.  Therefore, BC-EP is more accurate and converges faster than the other three popular FDTD methods.
Secondly, we investigate the stability of BC-EP since the stability of other methods have been demonstrated in the literature. The algorithm has been run for 200, 000 iterations under FDTD main loop which is long enough for the electric and magnetic waves completely leave the region as seen in Figure 14. This process has been repeated for ǫ = 3, 10, 30 and for mesh size ∆x = 10nm, 7.3nm, 4.6nm, and 2.8nm. At every 2, 000 iterations, a calculation of the SCS for all 601 wavelengths is conducted. Then the average relative error between the calculated SCS at that iteration and the exact SCS over those wavelengths is calculated. We observe that the result are very similar for different ǫ and different ∆x in the sense that the relative errors drop down to zero very quickly. This shows that there is no electromagnetic power left in the region and the fields return back to zero as excepted, meaning no late time instability even at very high contrast.

Conclusions and Future Research
We have built a new BC-EP method for the challenging two-dimensional transverse electric Maxwell's equations involving multiple dielectric media. Based on a relation for the electric fields across the dielectric media interface and the integral version of the Maxwell's equations, we derived this method which represents Maxwell's equations accurately while the previous numerical methods skipped some terms from the boundary conditions for the purpose of algorithm simplification. We then apply BC-EP together with three other popular FDTD methods, BC-EP, staircase and V-EP, to solve a numerical example of a dielectric cylinder with given analytic solution. The numerical results clearly demonstrated that it always achieves the highest convergence (close to second order) compared with other methods. The stability of BC-EP has also been verified numerically. Due to its simplicity, BC-EP can be merged into other FDTD software packages easily.
Our next research consideration is to extend BC-EP to solve two-dimensional Maxwell's equations involving multiple dispersive media. Numerical performance comparison of this method with other reported dispersive methods ( [19,20,21]) will be provided.