1. Introduction

Research fields like advanced optical imaging and integrated photonics explore and push forward the terrific possibilities offered by Maxwell’s equations. These developments rely more and more on complex geometries, like photonic crystals [1] or metamaterials [2], as well as on complex field configurations [3], so that full-vector numerical methods play a key role in understanding and designing them. In the past, several techniques have been proposed [4], like finite-difference and finite-element methods, integral equation methods and so forth. Among them, the Finite-Difference Time-Domain (FDTD) method, proposed by K.S. Yee [5], has gained much popularity for several reasons: it is rather easy to implement, the algorithm is intuitive, it can solve Maxwell’s equation for systems of arbitrary shape, it works in space and time domain.

In the FDTD method, space and time are discretized in a way that the derivatives in Maxwell’s curl equations can be written as finite central differences. The structure of the curl equations suggests that the field components are defined on different discrete positions in space, as shown in Fig. 1 for the two-dimensional case. The same holds for the time mesh, where the so-called leapfrog scheme is used. After a few simple rearrangements, the curl equations are finally transformed into loops that simulate the propagation of the electromagnetic field in space and time [6]. There are, however, some intrinsic problems that can make the FDTD method inaccurate, such as staircasing, numerical velocity dispersion and absorbing boundary conditions. The two latters will not be directly addressed in this work; the reader can refer to a reference book on FDTD [6]. Staircasing can be easily understood by looking at Fig. 1(a). Because the mesh can assign only discrete values to position, any object embedded in the simulation domain, will be pixelized in a way that a smooth interface becomes like a staircase. Consequently, the scattering properties of the system will change going from the real shape to the meshed shape, especially if the mesh is coarse and the dielectric contrast is large. This issue has limited the application of the FDTD method when the exact shape and permittivity values are important in determining the electromagnetic response.

A possible solution to staircasing is to depart from the simple Cartesian mesh picture and use non-orthogonal grids or curvilinear coordinates that follow exactly the shape of the objects [7, 8, 9, 10, 11, 12]. However, while improving the accuracy, such approach not only considerably increases the complexity of the algorithm, but can also cause numerical artifacts due to a highly irregular grid, like time instability, velocity dispersion and spurious wave reflection [15]. Better is to exploit a Cartesian mesh as much as possible and introduce distorted cells only when it is really necessary. For the special cells, a Contour Path (CP) FDTD algorithm can be obtained directly from Maxwell’s equations in integral form [13, 14]. The CP-FDTD method can still contain cells that potentially generate instability because of non-reciprocal nearest neighbor borrowing steps [15]. Moreover, the introduction of auxiliary field components and update equations slightly increases memory and CPU time. There have been improvements to CP-FDTD that solve the instability issue [15, 16, 17].

Another possible way for reducing the staircasing error is refining the Cartesian mesh in proximity of the interfaces, the so-called subgridding method [18, 19]. However, this scheme implies modifications at the fields-marching level, making the algorithm more complicated to implement, besides other numerical issues like spurious wave reflection.

A different approach, specific to dielectrics, is using Effective Permittivities (EPs) for the partially filled cells, without any distortion of the Cartesian grid. The question is, what value for the permittivity has to be chosen in order to get the best approximation of the dielectric interface? An early attempt in this direction has been made for modeling thin material sheets [20], even though there is still usage of auxiliary terms for field components normal to the interface and the procedure is limited to rectangular objects aligned with the mesh. A few years later, Kaneda et al. [21] proposed a phenomenological formula for the EP applicable to any kind of interface geometry, including curved surfaces. Their expression matches the rigorous result that can be obtained when the field component is perpendicular or parallel to the interface [22, 23]. These EPs improve the accuracy of the FDTD method, while keeping the same stability and simple structure of the original algorithm. However, there is no guarantee that the formula fulfills the proper boundary conditions at a curved interface or simply at a flat interface tilted with respect to the mesh axes. There are several works presenting other kinds of EPs: a volume average [24], a first-neighbor average [25, 26], Maxwell-Garnett, inverted Maxwell-Garnett and Bruggeman formulae [27], and other phenomenological derivations [28]. Some of them have been tested together for the purpose of comparison [29]. These proposals are not fundamentally more accurate than Kaneda’s approach.

The main problem with the formulation of EPs resides in the vectorial nature of the electromagnetic field. In fact, the same discontinuity can lead to quite different EP values depending on the orientation of the electric field with respect to the interface [22, 23]. Therefore, it is crucial that in the derivation of the EP, not only the geometry, but also the proper boundary conditions are taken into account. Along this line a non-diagonal EP-tensor can be obtained via the homogenization of a partially filled cell [30, 31]. However, its implementation requires the usage of both E and D, implying more storage and CPU time. Moreover, because in the FDTD method the field components are not defined in the same position, a nearest-neighbor average is required for linking E with D. Such average, can wash out the fulfillment of the boundary conditions at the interface. Recently, there have been other original ideas that improve the accuracy of the FDTD method under a rigorous treatment of the electromagnetic field at the dielectric interface, even though they increase the complexity of the algorithm more than the derivation of EPs does [32, 33, 34, 35].

 

Fig. 1. (a) FDTD mesh showing the staircasing effect for a curved interface (blue line); (b) location of the field components for H-modes and integration lines for Ampére law (blue segment) and Faraday law (red segment) for E_x|_i,j-1/2. ∆x and ∆y are the cell dimensions (mesh’s pitch); (i, j) refers to the position of H_z, while E_x and E_y are at (i, j - 1/2) and (i -1/2, j), respectively. Notice that the cells associated to the different field components partially overlap.

Download Full Size | PDF

In this work, using a contour path approach, we obtain EPs starting from Maxwell’s equations in integral form. These permittivities fulfill the boundary conditions for the electromagnetic field at a dielectric interface and, at the same time, they do not require any modification of the FDTD algorithm. Indeed, the main idea is to keep FDTD as simple as possible while reducing the staircasing error. Moreover, we propose another phenomenological formula of EP and compare it with the previous one. In the next section, we discuss in detail the formal derivation of our EPs for a set of partially filled cells. They do not represent the all possible situations, but they are enough to clearly show how one should proceed for any new case. Then, we test our effective-permittivity model against Mie theory [36] and compare it also with staircasing and the Volume-average Effective Permittivity (V-EP) [24]. We consider only the two-dimensional FDTD method for H-modes, those with the magnetic field perpendicular to the mesh, because the formalism is easier to understand and to implement. The case of E-modes, those with the electric field perpendicular to the mesh, simply has the V-EP as correct one, because the electric field is always parallel to the dielectric interface; for this reason, it will not be discussed.

2. The Effective Permittivities

In order to obtain EPs, we start from Ampére and Faraday laws in integral form:

Moreover, we restrict the discussion to dielectric non-magnetic media, so that D = εE and B = H, where the electric permittivity and the magnetic permeability are set to one (ε₀ = μ₀ = 1) for simplicity and e is the relative permittivity. In the two-dimensional FDTD method [6], the xy plane is discretized using space increments ∆x and ∆y along x and y, respectively, and the notation (i, j) represents a rectangular cell, with area ∆_x ∆y, centered at the position (i∆x,j∆y) on the mesh. According to the Yee scheme [5, 6], the non-zero field components for H-modes are positioned as follows: E_x|_i,j-1/2,E_y|_j-1/2,j,D_x|_i,j-1/2,D_y|_i-1/2,j,H_z|_i,j,B_z|_i,j, where the notation F|_i,j is used to name the field component F located in the cell (i,j), see Fig. 1(b). Thus, the electromagnetic field is associated to two orthogonal dual meshes. H_z is at the nodes of one mesh, while E_x and E_y are at the edges of the dual mesh, as it can be noticed from Fig. 1(a) and (b). This captures the topological structure of Maxwell’s equations and provides a geometrical interpretation of the FDTD method [37].

When Eqs. (1) are applied on these meshes, we keep following the Yee scheme for choosing the surface and line integrals: the flux is always computed through a surface normal to the field component, whereas the circulation is computed along the same direction of the field component, as shown in Fig. 1(b) with blue and red lines, respectively. Because the fields are homogeneous along the z direction, Eqs. (1) become:

for Ampère law and

for Faraday law. In Eqs. (2) and (3), the subscripts with x or y mean that the field component along the whole integration path is needed, while the subscripts with i, j, i+1/2 and so forth, represent the field component only at the corresponding position on the mesh. Within each integral, if the medium is homogeneous and the mesh is sufficiently fine, the field components can be assumed to be constant, i.e. E_x|_x,j-1/2 = E_x|_i,j-1/2. With little manipulation and with the time discretisation, Eqs. (2) and (3) become the usual Yee algorithm for H-modes in two-dimensions [6].

 

Fig. 2. Partially filled cells: (a) interface parallel to the field component, (b) interface orthogonal to the field component, (c) inclusion without crossing the integration lines, (d) inclusion crossing both integration lines, (e) inclusion crossing only the integration line of Ampére law, (f) inclusion crossing only the integration line of Faraday law. d,f represent the line filling factors, 1 and 2 mean media with ε₁ and ε₂ respectively.

Download Full Size | PDF

In the presence of a dielectric interface that crosses the line integrals, the electric field components cannot be assumed constant any more. The magnetic field is not affected, so that the integral on the left side of Eq. (3) follows the same treatment of the Yee algorithm. Therefore, we focus only on the left sides of Eqs. (2) and on the right side of Eq. (3). The two simplest situations of partially filled cells are sketched in Fig. 2(a) and (b), where a dielectric planar interface is parallel or perpendicular to the field component.

In Fig. 2(a), the interface crosses only the line integral of Eq. (2), while the integral of Eq. (3) can be treated as discussed earlier. The line integral of Eq. (2), which corresponds to the electric flux across the vertical dashed line, can be performed using the continuous tangential electric field, namely E _∥1 = E _∥2 = E _∥, where 1 and 2 refer to the electric field inside materials with electric permittivity equal to ε₁ and ε₂, respectively. E∥ is the value stored in the computer memory for that cell. Inside the cell, the field component in each medium is assumed to be constant. Applying this rule to D_x, for instance, gives:

where d measures the crossing of the dielectric interface in the cell. Therefore, Eq. (4) leads to the original Yee algorithm with the electric permittivity of that cell replaced by an EP ε_∥ = ε₂(d/∆y)+ε₁(1 - d/∆y), as already reported by several works [20,21, 22, 23]. The same result can be obtained for the other field component E_y.

When the dielectric interface is perpendicular to the field component, as shown in Fig. 2(b), only the integral of Eq. (3) needs special treatment. The procedure is similar to the previous case. Now, the physical quantity is the line integral of the electric field along the horizontal dashed line and it can be evaluated using the continuous normal electric displacement, namely D _⊥1 = D _⊥2 = ε ₁ E _⊥, where ε ₁(ε₂) has to be used if the location of the field component is in material 1(2). E _⊥ is the value that should be stored in the computer memory. Considering that the interface cuts only the cells associated to E_x|_{i,j-1 /2} and E_x|_i,j²1/2, the integral in Eq. (3) is performed as:

where f measures the crossing of the dielectric interface in the cell. Eq. (5) can be further simplified if the following change of variable is done for the value that has to be stored in the computer memory: E′_x|_i,j±1/2 = [(f/∆x) + (1 - f/∆x)ε ₂/ε ₁]E_x|_i,j±1/2. With this substitution, Eq. (5) leads to the original Yee algorithm, while Eq. (2) leads to the original Yee algorithm as well, if E′_x|_i,j-1/2 is used in place of E_x|_i,j-1/2 and the electric permittivity of the corresponding cell is replaced by the EP, ε_⊥ = [(f/∆x)/ε₂ + (1 -f/∆x)/ε₁]^-1. The same permittivity is discussed in [21, 22, 23] and more recent works. Notice that the value stored in the computer memory corresponds neither to the electric field before nor to the one after the dielectric interface. A similar approach holds for the E_y component. The EPs derived so far have been shown to preserve second-order accuracy in the Yee algorithm [22, 23].

There are many other possibilities for having partially filled cells, especially when one deals with complex geometries. In the following, we intend to show that the application of the previous simple ideas can always lead to an EP that represents the integrals of Eqs. (2) and (3). For instance, consider the partially filled cells sketched in Fig. 2(c), (d), (e) and (f). The starting point is the question: does the dielectric inclusion cross the line integrals? If not, then the integration is performed in the usual way, otherwise special treatment is required. Case (c) of Fig. 2 does not present any crossing, therefore this cell will have the usual Yee algorithm using ε ₁, because all integrals are performed in the region with material 1. The case (e) is similar to case (a), except the fact that the field component stored in the computer memory resides in material 1 instead of 2. However, the EP that is used is still ε _∥ = ε ₂(d/Δ) + ε ₁ (1 - d/Δ), where Δ means Δx or Δy if one refers to the E_y or the E_x field component, respectively. The case (f) is similar to case (b), except the fact that the field component stored in the computer memory is E_x′|_i,j-1/2 = [(f/Δx)ε ₁/ε ₂ + (1 - f/Δx)E_x|_i,j-1/2, with E_x|_i,j-1,2 in material 1, but the EP that has to be used is always ε _⊥ = [(f/Δx)/ε ₂ + (1 - f/Δx) /ε ₁]^-1, because D_x|_i,j-1/2 = ε ₂ E_x|_i,j-1/2 is replaced by D_x|_i,j-1/2 = ε ₁ E_x|_i,j-1/2.

The case (d), however, is different from the situations discussed so far and requires some more attention. In fact, the dielectric inclusion crosses both integration lines, one with the field component parallel and the other one with the field component perpendicular to the interface. Assuming that we are working with the E_x field component, the integral of Eq. (2) will give an EP ε _∥ = ε ₂ (d/Δy) + ε ₁ (1 - d/Δy). Besides, the integral of Eq. (3) will give the field transformation E_x′|_i,j-1/2 = [(f/Δx)+(1 - f/Δx)ε ₂/ ε ₁]E_x|_i,j-1/2. The combination of these effects results in a new permittivity for E′_x|_i,j-1/2, the value stored in the computer memory, given by

Once again, this kind of partially filled cell can be modeled with the usual Yee algorithm using the EP ε _eff = ε _∥ ε_⊥|ε ₂.

When the dielectric interfaces are neither parallel nor perpendicular to the field component, but tilted by a certain angle, the previously derived permittivities need to be modified. Again, we consider first the situation where only one integration line is crossed, then both. Figure 3(a) shows a tilted interface that affects only the integrals in Eqs. (2). n is the normal to the interface and d is measured exactly at the crossing with the integration line. In order to fulfill the boundary conditions, the field component has to be projected into normal and parallel components to the interface, using n, and re-projected to the usual component because of the scalar products in Eqs. (1). Therefore, if n is the projection of n along the field component, the integrals in Eqs. (2) become, referring to the D_x component,

 

Fig. 3. Partially filled cells with tilted interfaces: (a) tilted interface crossing only the integration line of Ampère law, (b) tilted interface crossing only the integration line of Faraday law, (c) tilted interface crossing both integration lines, (d) curved interface crossing only the integration line of Ampère law, (e) curved interface crossing only the integration line of Faraday law, (f) curved interface crossing both integration lines. n,m represent unit vectors normal to the interface, d, f represent the line filling factors, 1 and 2 mean media with ε ₁ and ε ₂, respectively.

Download Full Size | PDF

The EP ε _∥,n = [ε ₁ n ² + ε ₂(1 - n ²)](d/Δy)+ ε ₁ (1 - d/Δy) reduces to ε _∥ for n = 0. In a similar way, if d > Δy/2, the permittivity is ε _∥,n = ε ₂(d/Δy) + [ε ₂ n ² + ε ₁(1 - n ²)](1 - d/Δy). The treatment for the other field component is completely analogous.

The opposite situation is shown in Fig. 3(b), where only the integral in Eq. (3) crosses the dielectric interface. By application of the boundary conditions and considering for the moment only the E_x component the integral is computed as follows,

After exploiting a field transformation similar to the one performed for Eq. (5), the result is once again an EP ε _⊥,n = {(f/Δx)[n ²/ε ₂ + (1 - n ²)/ε ₁]+(1 - f/Δx)/ε ₁)^-1, which becomes ε _⊥ for n = 1. In a similar way, if f > Δx/2, the permittivity is ε _⊥,n = {(f/Δx)/ε ₂+(1 - f/Δx)[n ²/ε ₁ + (1-n ²)/ε ₂]}^-1.

Now that one knows how to deal with tilted interface, the case when two crossings occurs (see Fig. 3(c)) can be readily solved using the previous result: the integrals of Eqs. (2) and (3) are implemented using an EP ε _eff,n = ε _∥,n ε _⊥,n/ε ₂ if the center of the cell is located in material 2, like in Fig.3(c), otherwise ε _eff,n = ε _∥,n ε _⊥,n/ε ₁. Notice also that the EPs ε _∥,n and ε _⊥,n change if d ≶ Δ Δ/2 and f ≶ Δ/2, where Δ can be Δx or Δy depending on the field component.

We have seen that also for interfaces that are tilted with respect to the Cartesian FDTD mesh, the integrals involving the electric field and the electric displacement in Eqs. (2) and (3) become the usual Yee algorithm by assigning an EP to a partially filled cell related to a field component. Of course, there exist many more cases for tilted interfaces, for examples with a corner similar to Fig. 2(d). Nevertheless, all occurrences can be tackled following the procedure already shown.

It is worth to briefly discuss another situation that is quite common in FDTD simulations: partially filled cells with curved interfaces. Some examples of them are shown in Fig. 3(d), (e) and (f). In these cases, the angle between the normal to the interface and the field component is not constant over the FDTD cell, because the dielectric inclusion has a curvature. Following our approach, the normal that is needed for the calculation of the EP is only at the crossing between the integration line and the dielectric interface. Therefore, the result will be as for a tilted interface: ε _∥,n = [ε ₁ n ² + ε ₂(1 - n ²)](d/Δy)+ ε ₁(1 - d/Δy) for case (d), ε _⊥,n = {(f/∆x)[n ²/ε₂ + (1 - n ²)/ε₁] +(1 - f/∆x)/ε₁}^-1 for case (e) and ε_eff,m,n = ε_∥,mε_⊥,n/ε₂ for case (f). Notice that for the latter there are two normals, n and m, inside the formula for the EP. Other situations with curved dielectric inclusions can be treated in a similar manner.

The method of the Contour Path Effective Permittivity (CP-EP) is expected to be more accurate than the staircase approximation, because much more information on the actual geometry is passed to the FDTD algorithm. On the other hand, because this information can be represented using EPs, the algorithm remains as fast and memory efficient (except a little more memory to store more permittivity values) as the basic Yee algorithm. In fact, the calculation of the permittivities can be done before time-marching the fields, maybe with external software, then passed to existing FDTD code. This preprocessing time can be very small, especially if the geometry is defined by means of analytical functions, like in vector graphics. The method can be readily extended to three-dimensions, keeping in mind that the line integrals of Eqs. (2) become surface integrals.

Since we are also interested in finding a trade-off between accuracy and complexity of coding, we propose a simple phenomenological formula for EP that still achieves a great improvement with respect to staircasing. First, instead of considering the crossing between the dielectric interface and the integration lines (surfaces in three-dimensions), we compute the cell filling ratio s, i.e. the fraction of material 2 with respect to the cell area (volume in three-dimensions). Secondly, the normal n to the interface is approximated by defining a unit vector along the connecting line between the center of the cell and the center of curvature of the interface. The idea is to weight ε _∥ and ε_⊥ using the following formula:

where n is the projection of n along the field component of that cell and ε_∥ and ε_⊥ are computed using the filling ratio s, in place of d and f. Eq. (9) represents a Volume-average Polarized Effective Permittivity (VP-EP) that accounts for the orientation of the field with respect to the interface. The formula in Eq. (9) becomes en or ε_⊥ if n = 0 or n = 1, respectively. The calculation of the VP-EP is much simpler than the CP-EP, but, as a drawback, it is not able to clearly distinguish among all the possible partial fillings that occur and it is not a rigorous derivation from Eqs. (2) and (3). Notice that this formula is similar to the result obtained by Kaneda et al. [21], but it is different from the approach presented in [30, 31] by the fact that the resulting dielectric tensor is diagonal and its elements are computed at the position of the relative field component; as a result, the averaging step during the time-marching introduced by [30] is not required.

3. Numerical Tests and Discussion

The CP-EP and the VP-EP have been implemented in a two-dimensional FDTD code for testing. For comparison, we have also considered the staircase and the V-EP, which corresponds to ε _∥ with filling ratio s. We are particularly interested in assessing the convergence properties of V-EP, because this approach seems to be quite popular in the FDTD community. There are actually several choices for staircasing a certain structure on the FDTD mesh; we have chosen this one: ε = ε ₂ if for a given field component, the center of the cell is inside medium 2, and ε = ε ₁ otherwise. Like for the EPs, this criterion is applied to each field component, independently, meaning that E_x|_i,j-1/2 and E_y|_i,j-1/2 may have different ε, because they are not exactly in the same position. In the literature, we have seen two main categories of tests related to the FDTD method: calculation of scattering cross sections (SCS) [14, 26, 27, 29, 32], calculation of resonant frequencies [16, 17, 21, 24, 28, 34]. We want to compare our method with Mie theory [36] to have both kinds of tests: total SCS and position of the Mie resonances for a dielectric cylinder. Because the Mie coefficients depend on the boundary conditions of the electromagnetic field at the interface between the background medium and the cylinder, we believe that the result obtained with the FDTD method should be quite sensitive on how the dielectric interface is modeled.

 

Fig. 4. Layout of the FDTD calculation: CPML layers (gray), cylindrical scatterer (orange), the H-polarized incident plane wave is excited using the total field / scattered field method, the integration line is for the calculation of the total SCS.

Download Full Size | PDF

The layout of the FDTD calculation is shown in Fig. 4: the center of the cylinder is positioned at (i ₀ - 1/2, j ₀ - 1/2), where (i ₀, j ₀) is referred to H_z. There are more options to locate the center of the cylinder [26], however we have observed that our method is almost insensitive to the choice of its position in the mesh. An E_y-polarized incident plane wave is excited using the total field / scattered field method [6] and the temporal profile of the source has a Gaussian envelope, exp[- (t - t ₀)²/2σ²]cos(2πct/λ ₀), to excite many wavelengths at the same time. t ₀ ensures that the incident field starts with a negligible value at the beginning of the simulation, σ is chosen to cover the desired wavelength range, c is the speed of light and λ₀ is the central wavelength of the pulse. The time step is set to be ∆t = S∆/c, where S = 0.98/√2 is chosen to enforce numerical stability [6] and ∆ is the spatial discretization (we use a square mesh ∆x = ∆y = ∆). The scattered field exits the computational mesh with negligible reflection through Convolutional Perfectly Matched Layers (CPML) absorbing boundary conditions [39], where a layer thickness of 16 grids was used. The fields are transformed into frequency domain using an on-the-fly discrete Fourier transformation. The resulting wavelength step ∆λ, which is related to the number of points in the discrete Fourier transform N _DFT and the wavelength range by ∆λ = (λ_max - λ_min)N _DFT, is different for each numerical test and it is specified in the corresponding figure caption. Finally, the total SCS is obtained by summing the normal component of the Poynting vector along the integration line shown in Fig. 4. The reference data from Mie theory are given by the scattering code for an infinite dielectric cylinder under normal illumination, as presented in [36]. Because we are mainly interested in the optical and near-infrared wavelength range, our tests were performed for wavelengths between 0.4μm and 1.6μm on cylinders with radii r from a few tens of nm up to a few hundreds of nm and relative dielectric constant ε from 3 to 20. The background medium is set to vacuum, ε = 1. In the following, we discuss a representative selection of them.

The first test considers a low-index dielectric cylinder, ε = 3, with radius r = 400nm. Figure 5(a) shows the total cross section computed using CP-EP, V-EP and Mie theory for N_λ = 25. N_λ is the number of divisions for the shortest wavelength inside the cylinder, λ = 400nm/√ε. For N_λ = 25, a fine FDTD mesh, all EPs agree quite well with Mie theory, so that Fig. 5(a) would look the same also for staircase and VP-EP In order to see the small differences, almost not visible by eye, in Fig. 5(b) we plot the relative error err_λ = (SCS_FDTD - SCS_Mie)/SCS_Mie as a function of wavelength for N_λ = 25. Notice that the largest error accumulates where the SCS exhibits minima and maxima. The fact that the error oscillates between negative and positive values indicates that the SCS computed with FDTD is crossing several times the one computed using Mie theory. While staircase, VP-EP and CP-EP have comparable oscillations, the V-EP clearly exhibit a larger error though still quite small. Since the dielectric contrast is small, we do not expect that the EP converges much better than staircase. However, V-EP is actually doing worse, suggesting that a wrong choice of EP can damage the convergence of the FDTD method. Moreover, in the wavelength region above 0.7 μm, the error associated to staircase, VP-EP and CP-EP is further reduced to a value close to 0.1%, while V-EP is around 0.3%.

 

Fig. 5. Accuracy on the total SCS: (a) total SCS for N_λ = 25, (b) relative error on the total SCS for N_λ = 25, (c) average relative error on the total SCS. N_λ is the number of divisions for the shortest wavelength in the cylinder; i. e. N_λ = (400nm/∆)/√ε. The phase-velocity error in (c) is computed for N_λ. Parameters: ε = 3, r = 400nm and ∆λ = 1nm.

Download Full Size | PDF

We also show how this error depends on N_λ. To this purpose, for each N_λ we compute an average error errN_λ = Σ_λ |err_λ|/N, where N is the total number of wavelengths used in the sum (N = N_DFT). As shown in Fig. 5(c), errN_λ exhibits fluctuations for the staircase case. In other words, a finer mesh does not always imply an improved result. In fact, during staircasing it could happen that a particular N_λ is able to match better the geometry than the successive finer division. Such magic-number like behavior is not predictable, especially when the geometry is more complex than a cylinder. On the other hand, the effective permittivities are always improving the result as the mesh gets finer and finer. However, if the EP does not properly preserve the scattering properties of the object, the convergence can be slower, as shown for the V-EP. In our tests we have noticed that also the two-dimensional Maxwell-Garnett EP [40] has given very good convergence results, but not for high-index contrasts.

The next class of tests regards both the calculation of the total SCS and of resonant wavelengths for a cylinder with parameters ε = 12 and r = 150nm. Figure 6(a) displays the total SCS, computed with V-EP and CP-EP for N_λ = 25, compared to the Mie theory result. By eye, one sees that the CP-EP calculation shows better agreement with Mie theory both in the total SCS and in the position of the resonances. The average relative error on the total SCS, errN_λ, is given in Fig. 6(b) for all cases. Again, the V-EP has the worst performances, but the difference with respect to staircase and the other EPs is even more evident than in Fig. 5(c), because the index contrast is higher. The VP-EP and CP-EP are a few percent better than staircase, with less pronounced oscillations, up to N_λ =25. When the mesh is very fine, staircase is as good as CP-EP.

 

Fig. 6. Accuracy on the total SCS and the resonant wavelengths: (a) total SCS for N_λ =25, (b) average relative error on the total SCS, (c) relative error on the wavelength of the resonance λ ₁ = 675.8nm, (d) relative error on the wavelength of the resonance λ ₂ = 5323nm. N_λ is the number of divisions for λ ₂ in the cylinder; i. e. N_λ = (λ ₂/∆)/√ε. The phase-velocity error in (b)-(d) is computed for N_λ. Parameters: ε = 12, r = 150nm and ∆λ = 0.25nm.

Download Full Size | PDF

We now consider the error associated to the resonant wavelength for two peaks, the first and the second one starting from the right of Fig. 6(a). The relative error, computed as err_peak = (λ _FDTD - λ _Mei)λ _Mie, is shown in Fig. 6(c) and (d) for the first and second peaks, respectively. In both Fig. 6(c) and (d), V-EP is the worst result. VP-EP and CP-EP are nearly the same and perform better than staircase, both in term of error and in terms of error fluctuations. For the finer meshes, the error is as small as the error on the discrete Fourier transform, which is of the order of 0.1%. A similar behavior has been found for the third sharp peak of Fig. 6(a) at λ ≃ 440nm.

 

Fig. 7. Accuracy on resonant wavelengths: (a) resonant peak for N_λ = 11 and (b) N_λ = 19, (c) relative error on the wavelength of the resonance λ_o = 679 4nm. N_λ is the number of divisions for λ_o in the cylinder; i. e. N_λ = (λ_o/∆)/√ε. The phase-velocity error in (c) is computed for N_λ. Parameters: ε = 20, r = 120nm and ∆_λ = 0.2nm.

Download Full Size | PDF

As the last example, we present another calculation of resonant wavelengths for a larger dielectric constant, ε = 20, and slightly smaller radius r = 120nm. The peak in the total SCS is shown for two values of N_λ in Fig. 7(a) and (b). Notice the anomalous behavior of staircasing: the calculation seems to be more accurate for the coarser (see Fig. 7(a)) than for the finer mesh (see Fig. 7(b)). This effect is even more evident in Fig. 7(c), where the relative error exhibits significant jumps even for very fine meshes. Once again, V-EP is much less accurate than VP-EP and CP-EP. For the finest mesh, the error is comparable to the wavelength discretization due to the Fourier transform.

Because of the space-time discretization, the actual phase velocity associated to a wave propagating in the FDTD mesh is different than that determined just by the medium’s refractive index v _ph = c/√ε, both in terms of magnitude and anisotropy [6]. In the previous Figs., we have also plotted the error related to the phase velocity only for propagation along the grid axes, which is given by the formula errv _ph = (v _ph - ṽ _ph)/v _ph = 1 - π/{N_λ sin^-1[sin(πS/N_λ)/S]}, where v _ph is the numerical phase velocity and S = v _ph∆t/∆ [6]. When the mesh is coarse (N_λ is small), err_vph can be so large that it exceeds the FDTD error obtained by comparison with Mie theory. We have not directly addressed this issue here, but it would be interesting to know how much the velocity error affects the calculation for scattering problems and test the EPs using low-dispersion FDTD algorithms [38]. We believe that without low-dispersion algorithms it is difficult to say if the small differences between VP-EP and CP-EP are due to the approximation in the EP or to velocity dispersion.

 

Fig. 8. Color maps of the EPs: (a) staircase, (b) V-EP, (c) VP-EP, (d) CP-EP.

Download Full Size | PDF

Regarding resonance wavelengths, the staircase result jumps around the correct value, while the V-EP approaches to it without fluctuations, but with a larger error. That is why, especially in the absence of analytical solutions or other references, one is led to think that V-EP is better than staircasing. On the other hand, VP-EP and CP-EP are shown to be better than staircasing in terms of accuracy and stability of the error. Moreover, they give very much the same error. To gain more insight on these facts, we look at the permittivity values given by staircase, V-EP, VP-EP and CP-EP. A color map for them is shown in Fig. 8 for the E_x component. Obviously, for staircase (a) only two values are possible, corresponding to black or white. Notice that while V-EP, VP-EP and CP-EP have similar values for the cells close to the top and bottom of the cylinder, V-EP departs from VP-EP and CP-EP for cells close to the left and right sides of the cylinder. Indeed, the interface at the top and bottom is almost parallel to the field component, so that V-EP is close to the correct result. On the contrary, at the right and left sides, the interface is nearly perpendicular to the field component and V-EP is larger than the proper EP. This points out the importance of taking into account the polarization of the field for the derivation of EPs. The fact that VP-EP and CP-EP have nearly the same color maps and yield similar convergence properties suggests that one could use VP-EP, which is very easy to implement, rather than CP-EP.

4. Conclusion

We have proposed two schemes for the calculation of effective permittivities to improve the accuracy of the two-dimensional FDTD method for dielectric media. The CP-EP results from the rigorous application of the field boundary conditions at the crossings between the integration lines and the interface and can be applied to any kind of partially filled cell. The VP-EP is a phenomenological formula that tries to find a compromise between accuracy and complexity of coding. It has been shown that VP-EP is as good as CP-EP over a wide range of tests, so that it can be used in place of CP-EP, being much simpler to implement.

Our results were compared with Mie theory and two other FDTD schemes: staircase and V-EP. All methods approach the Mie-theory result for very fine meshes, but the CP-EP and VP-EP exhibit faster convergence. The major problem of staircasing is the instability of the error, so that one does not know if the mesh is fine enough to fulfill a determined benchmark. On the contrary, V-EP has the advantage of exhibiting a stable error, allowing benchmarking, but it introduces a significant bias that worsens the performances reachable with properly chosen EPs. Even though these tests do not represent a mathematical proof of the convergence of FDTD with EPs, as shown for instance in Ref. [33], we consider them quite convincing and of practical importance.

It would be interesting to extend these ideas to dispersive media and, more importantly, to metals at optical wavelengths, where accurate modeling of complex structures is still a challenging problem with FDTD, especially when surface plasmon resonances are involved.