1. Introduction

Various geometries of interest involve curved interfaces between dielectrics and metals. These include periodic split ring resonators (SRR) and plasmonic nanorods, which have been employed in a wide range of applications related to optical metamaterials [1–5]. Another important geometry is the plasmonic ring resonator which has been recently used to control the transmission of surface plasmon polaritons (SPPs) [6, 7]. These motivate the analysis of the wave propagation properties of these structures with full-wave numerical techniques such as the Finite-Difference Time-Domain (FDTD) method. The FDTD method is a popular numerical tool for photonics on account of its ability for handling complex geometries and dispersive media, and its broadband characteristics. However, the method suffers from a serious accuracy degradation due to pronounced staircasing errors when approximating curved interfaces by its conventional rectangular grid [8, 9]. Notably, the staircasing errors become more severe in plasmonic structures, as surface plasmons are strongly localized at metal-dielectric interfaces, where staircasing is applied.

Previous research on these problems is mainly based on two different approaches: modifying effective permittivities to account for the interface cells [10–13] and unstructured mesh techniques [14, 15]. A scalar effective permittivity technique was adopted in [10, 11] and extended to dispersive media by a z-transform and an auxiliary differential equation (ADE) method, respectively. However, the use of scalar effective parameters actually violates the field boundary conditions at the interface, and the resultant fourth-order time derivatives in [10] need significant memory and execution time. Another effective permittivity theory of interest is the subpixel smoothing technique [16], which employed a tensor effective permittivity. This technique was extended to sloped dispersive interfaces in [12] by splitting the electric fields into three auxiliary variables, instead of one, which greatly increases memory and computational requirements. By employing a coordinate rotation of the effective permittivity tensor, [13] reduced the three variables back to two. Moreover, [14] introduced a triangular mesh FDTD scheme in 2-D, which is provably stable and rapidly convergent. However, the time step limit is severely constrained by the smallest element of the triangular mesh. Along similar lines, [15] derived a linearized non-local dispersion model using the Discontinuous Galerkin Time-domain (DGTD) method. The advantage of these techniques is that the unstructured triangular mesh naturally achieves better conformity to curved interface structures than the rectangular FDTD mesh.

In this paper, we present a different approach that strikes a good balance between accuracy and efficiency, by keeping the structured mesh of FDTD and implicitly enforcing field boundary conditions at the interface Yee cells. This is achieved by formulating a tensor ADE-FDTD technique, extending the tensor FDTD method of [17, 18] from simple lossy dielectrics to dispersive media of Lorentz type dispersion. Through further derivation, the proposed method can also handle thin subcellular curved sheets in a Yee cell. Both [19, 20] independently introduced effective permittivities for dielectric thin sheet inclusions by enforcing boundary conditions at material interfaces. However, when using a structured FDTD mesh, these methods only work with planar interfaces. In this paper, we introduce an effective tensor permittivity to account for the boundary conditions at the two interfaces of thin sheets, thus making it possible to embed arbitrarily inclined and curved thin sheets in a Yee cell.

The rest of the paper is structured as follows: in section II, the effective tensor permittivity for sloped interfaces and thin sheets is derived. In section III, the tensor ADE-FDTD update scheme for Lorentz media is derived. In section IV, periodic silver nanorods are simulated to demonstrate the proposed method’s effectiveness at handling sloped dispersive interfaces. As well, a PEC cavity with dielectric thin sheets is modeled to validate our thin sheet model. Finally, a thin ring of a Lorentz medium is simulated to show the efficiency of the tensor ADE-FDTD for dispersive, curved thin sheets. Section V summarizes our contributions.

2. Effective tensor permittivity for sloped interfaces and thin sheets

The adopted tensor effective permittivity model was proposed in [21] for handling simple, dielectric interfaces. Then the technique was improved by [17], which adjusted the tensor effective permittivity to field polarization and was extended to simple, lossy media interfaces in [18], via a z-transform method. In this paper, we adopt the tensor model of [18] and consider the 2-D case of Fig. 1(a), where fields (E_x, E_y, H_z) are sought for. The FDTD update scheme will be derived based on the integral form of Maxwell’s equations. The FDTD update equation for $D_y\left(i,j+1/2\right)$, for instance, can be derived by formulating a loop (Fig. 1(b)), which has length of $\mathrm{\Delta }l$ along the $z-$direction, as in the $z-$direction both fields and modeled geometries are considered to be invariant. Thus, the flux density $D_y\left(i,j+1/2\right)$ can be updated as the following expression:

Note that $D_y\left(i,j+1/2\right)$ and $E_y\left(i,j+1/2\right)$ are at the same grid point. If the medium inside the face is homogeneous, the determination of $E_y\left(i,j+1/2\right)$ from $D_y\left(i,j+1/2\right)$ only requires the medium permittivity. However, if the medium is inhomogeneous, boundary conditions must be enforced. In the following, we will first derive the effective tensor permittivity for sloped interfaces and then extend it to sloped thin sheets.

 

Fig. 1 (a) Yee cell of 2-D TE mode (E_x, E_y, H_z). (b) Update of flux density in y-direction from a line integral of magnetic field around a loop centered at (iΔx, (j+1/2)Δy).

Download Full Size | PDF

Considering two Yee cells centered at the grid point $\left(i\mathrm{\Delta }x,\left(j+1/2\right)\mathrm{\Delta }y\right)$, the background dielectric medium ε₁ is interfaced with a dispersive medium ${\epsilon }_2\left(\omega \right)$, as in Fig. 2(a). In the following, we discuss how to find the tensor effective permittivity for the update of the electric field node $E_y\left(i,j+1/2\right)$.

In the middle face centered at $\left(i\mathrm{\Delta }x,\left(j+1/2\right)\mathrm{\Delta }y\right)$ and between the two magnetic nodes, we consider the electric field ${\mathbf{E}}^{\left(1\right)}$ in the background medium and the electric field ${\mathbf{E}}^{\left(2\right)}$ in the dispersive medium. The relationship between ${\mathbf{E}}^{\left(1\right)}$ and ${\mathbf{E}}^{\left(2\right)}$ depends on the interface normal $\widehat{\mathbf{n}}$ and must satisfy the following boundary conditions:

 

Fig. 2 (a) Sloped interface across two Yee cells centered at (iΔx, (j + 1/2)Δy). (b) Sloped thin sheet across two Yee cells centered at (iΔx, (j + 1/2)Δy).

Download Full Size | PDF

The matrix $\left[A\right]$ is the same for the two pairs of fields as the interface normal is the same and introduces no difference to the matrix $\left[A\right]$ in Eq. (4). As a result, we can derive that ${\mathbf{E}}^{\left(1\right)}={\mathbf{E}}^{\left(3\right)}$. By replacing ${\mathbf{E}}^{\left(3\right)}$ with ${\mathbf{E}}^{\left(1\right)}$, the electric field relation of Eq. (4) can be directly adopted in the thin sheet model. Notably, sloped interfaces and thin sheets share the same relationship between the electric field in the dispersive medium and that in the background medium. Once we obtain the relationship between the electric fields ${\mathbf{E}}^{\left(1\right)}$ and ${\mathbf{E}}^{\left(2\right)}$, we can derive the effective tensor permittivity. The electric flux through a face of a Yee cell can be expressed in terms of either D or E as:

We can apply the above equation to $D_y\left(i,j+1/2\right)$ to a face similar to Fig. 1(b). Figures. 3(a) and 3(b) show the corresponding faces for the sloped interface and the thin sheet. The part of the face $\mathrm{\Delta }x\mathrm{\Delta }l$ in the dispersive medium is depicted as gray, and their widths are indicated as g for both cases, which are also marked in Fig. 2. The flux density $D_y\left(i,j+1/2\right)$ can then be expressed as follows:

Hence:

 

Fig. 3 Calculation of the tensor effective permittivity through surface integration at the grid point (iΔx, (j + 1/2)Δy) (a) Faces for the sloped interface. (b) Faces for the thin sheet.

Download Full Size | PDF

Note that in Eq. (9), the flux density D_y is only dependent on the electric field in medium 1. However, the electric field in the $x-$ direction also contributes to D_y, which means that the electric field in the $y-$ direction cannot be directly derived by D_y. Here, we introduce another flux component $D_x\left(i,j+1/2\right)$, at the same grid point as $D_y\left(i,j+1/2\right)$. By repeating the above procedure via Eqs. (6)–(9) for $D_x\left(i,j+1/2\right)$, the relationship between the flux density D at the grid point $\left(i\mathrm{\Delta }x,\left(j+1/2\right)\mathrm{\Delta }y\right)$ and the electric field ${\mathbf{E}}^{\left(1\right)}$ is formulated as:

3. ADE-FDTD update scheme

In the following, we will first consider the medium in region 2 as a simple dielectric medium and give the normal update scheme. Then, we will consider the medium in region 2 as a Lorentz dispersive medium, and finally give the ADE-FDTD update scheme.

3.1. Electric field update

In this section, we focus on the update scheme for the electric field using the derived tensor effective permittivity in Eq. (10). We present the update of the electric field node $E_y\left(i,j+1/2\right)$. Let us first consider the medium in region 2 as a simple dielectric medium. In every FDTD time step, we can readily update D having the magnetic field at time step n+1/2. Then, D can be used to derive ${\mathbf{E}}^{\left(1\right)}$ from:

To derive ${\mathbf{E}}^{\left(1\right)}$ from the above equation, we need all components of D at the interface grid node. $D_y\left(i,j+1/2\right)$ is found by the normal update equation of Eq. (1) and $D_x\left(i,j+1/2\right)$ can be found by interpolating the adjacent D_x components:

Once ${\mathbf{E}}^{\left(1\right)}$ is derived from the flux density using ${[\epsilon ]}^{-1}$, we can use Eq. (4) to find ${\mathbf{E}}^{\left(2\right)}$. Then, we need to check whether the electric field node $E_y\left(i,j+1/2\right)$ is in the dielectric host or in the dispersive medium. If it is within the dielectric, $E_y\left(i,j+1/2\right)$ is updated by the y-component of ${\mathbf{E}}^{\left(1\right)}$; otherwise $E_y\left(i,j+1/2\right)$ is updated by the y-component of ${\mathbf{E}}^{\left(2\right)}$. Note that the y-components of both ${\mathbf{E}}^{\left(1\right)}$ and ${\mathbf{E}}^{\left(2\right)}$ need to be stored and will be used in the update of the magnetic field.

Then, we consider the ADE-FDTD update scheme for cases when the region 2 is filled with a dispersive medium. The following Lorentz dispersion model is adopted to describe the permittivity of the dispersive medium:

Then, Eq. (10) is cast in the following form:

Mapping Eq. (14) to the time-domain and discretizing it with second-order finite-differences yields update equations of the form:

The coefficients are given by:

Combining Eqs. (15) and (16), we derive explicit update equations for the electric field in the background medium:

Then we need to derive the electric field in the dispersive medium via Eq. (4). However, ${\mathbf{E}}^{\left(2\right)}$ cannot be directly derived, due to its frequency-dependent relation to ${\mathbf{E}}^{\left(1\right)}$. Here we introduce another auxiliary vector $\mathbf{Q}=Q_x{\mathbf{a}}_{\mathbf{x}+Q_y{\mathbf{a}}_{\mathbf{y}}}$ as:

Then, Eq. (4) becomes:

Mapping the auxiliary vector to the time-domain and discretizing it with second-order finite-differences yields update equations of the following form (where the coefficients $c_{1,2,3}$ and $e_{1,2,3}$ are the same as in Eq. (17)):

Combining Eqs. (20) and (21), the following explicit update equations for the electric field in the dispersive medium are derived:

Finally, we check the location of the electric node $E_y\left(i,j+1/2\right)$ to decide whether it is updated by $E_y^{\left(1\right)}$ or by $E_y^{\left(2\right)}$. We also need to store both $E_y^{\left(1\right)}$ and $E_y^{\left(2\right)}$ for the update of the magnetic field.

 

Fig. 4 Contour path integration for the update of magnetic field, E⁽¹⁾ and E⁽²⁾ are both needed.

Download Full Size | PDF

3.2. Magnetic field update

Once all the electric field components have been computed, the magnetic field updates can be derived by contour path integration along the perimeter of a standard cell, as shown in Fig. 4. The lengths inside the dispersive medium on the left and bottom cell edges are marked as l_y and l_x for both models. Then, the line integration of the electric field along the left cell edge is expressed as:

3.3. Implementation of tensor ADE-FDTD for a Lorentz dispersive medium

Before beginning the FDTD time loop, the coefficients of Eq. (17) are determined and stored. Then, at every cell on the curved interface or sloped thin sheet, the matrices ${\overline{\epsilon }}_{\alpha }$ ,${\overline{\epsilon }}_{\beta }$, ${\overline{r}}_{\alpha }$ and ${\overline{r}}_{\beta }$ are stored for the updates of the electric field.

For electric field updates: 1) Compute D_x and D_y at every interface grid node via the standard FDTD update equations and interpolate, at $\mathrm{\Delta }/2$ intervals, to get the D_y and D_x at each other’s grid points. 2) Compute ${\mathbf{E}}^{\left(1\right)}$ at every interface cell according to Eq. (18), then compute ${\mathbf{E}}^{\left(2\right)}$ at every interface cell via Eq. (22). 3) Update the electric field by ${\mathbf{E}}^{\left(1\right)}$ or ${\mathbf{E}}^{\left(2\right)}$ based on the location of the corresponding electric field nodes. 4) Update the auxiliary vectors P and Q using Eq. (16) and Eq. (21).

For magnetic field updates: 1) Calculate the line integrals via Eq. (23). 2) Update the magnetic field using the previously calculated line integrals via Eq. (24).

4. Numerical results

In this section, we demonstrate the proposed tensor ADE-FDTD for sloped interfaces by modelling a one-dimensional array of plasmonic cylinders. Then, we apply the thin sheet version of the method to model a thin Lorentzian ring in a waveguide.

 

Fig. 5 1-D array of periodic silver nonarods and its correspondent simulation structure realized by Bloch’s boundary conditions.

Download Full Size | PDF

4.1. 1-D periodic silver nanorods

We will consider an array of silver nanorods (Fig. 5), an example also treated in [10, 14]. The one dimensional periodic structure is modeled by Bloch’s periodic boundary conditions. In our two dimensional TE computational domain, E_x and H_z are involved in the update equations and the periodic boundary conditions, as tangential to the direction of periodicity (here defined by the y-axis):

 

Fig. 6 Spectrum comparison between staircased FDTD and tensor ADE-FDTD method. The cell size of the former varies from 2.5 nm to 0.625 nm; the cell size for the latter is 2.5 nm.

Download Full Size | PDF

Table 1. Comparison of frequencies (THz) of the first five modes, computed by staircased FDTD, tensor ADE-FDTD and EPs (effective permitivities) [10] for $k_y\alpha =\pi$.

View Table | View all tables in this article

 

Fig. 7 Magnetic field magnitude and electric field vector of the five modes of the nanorods array when k_y α = π, using the tensor ADE-FDTD method with Δ = 2.5 nm. (a) f = 864 THz (b) f=990 THz. (c) f=1036 THz. (d) f=1069 THz. (e) f=1114 THz.

Download Full Size | PDF

 

Fig. 8 Magnetic field magnitude and electric field vector of the five modes of the nanorods array when k_y α = π, using the staircased FDTD method with Δ = 2.5 nm. (a) f = 828 THz (b) f=969 THz. (c) f=1053 THz. (d) f=1081 THz. (e) f=1115 THz.

Download Full Size | PDF

The frequencies of the first five modes are shown in Table 1, along with those found by applying the isotropic effective permittivities (EPs) method of [10]. The results of the proposed tensor ADE-FDTD method converge well as the mesh is refined, unlike the results of the staircased FDTD which exhibit oscillations with mesh refinement. Then, we compare the resonant patterns of the first five modes using both methods. The accurately resolved modal patterns of the magnetic field along with the electric field vectors are plotted in Fig. 7, using the tensor ADE-FDTD method. Similar graphs are generated using the staircased FDTD method in Fig. 8. The cell size is chosen as 2.5 nm for both methods. The patterns found with the proposed method are clear and symmetrical. In comparison, the staircased FDTD fails to obtain clear modal field patterns, especially for the last three modes, although their corresponding resonant frequencies are close to those of the proposed tensor ADE-FDTD method. Then, we check the convergence rates of the two methods. We compute the root-mean-square (RMS) error of the frequencies of the first two modes, with the following equation:

 

Fig. 9 Convergence comparison between staircased FDTD and tensor ADE-FDTD method (log-log scale). RMS error of the resonant frequency of the first two modes is compared for the two methods; the cell sizes of both methods vary from 5 nm to 0.625 nm.

Download Full Size | PDF

 

Fig. 10 Dispersion diagram of the nanorods array, results in [10] are also compared.

Download Full Size | PDF

Table 2. Comparison of execution time and RMS error in the resonant frequencies of the first two modes, computed by the staircased FDTD and the proposed method using various cell sizes.

View Table | View all tables in this article

4.2. Sloped thin dispersive sheet

In this section, we present two case studies, where the proposed tensor ADE-FDTD technique is applied to sloped thin sheets. Our method is first evaluated in a cavity geometry and then applied to compute the scattering parameters of a thin dispersive ring in a waveguide.

 

Fig. 11 (a) Sloped thin dielectric sheet array loaded cavity, 50 cm × 2.5 cm in size. The thin sheets has a width of w and inclination angle of θ. An actual cell for the tensor method is marked as a square. (b) Spectrum resolved by tensor FDTD (Δ = 5 mm) and FDTD with very fine mesh (Δ = 0.2 mm), with sheet width w = 1 mm and inclination angle θ = 45°.

Download Full Size | PDF

4.2.1. Verification of thin sheet model with simple dielectric media

The accuracy of our proposed method was verified by simulating a 50 cm × 2.5 cm cavity with PEC walls loaded with parallel, sloped thin sheets. The layout of the simulation is shown in Fig. 11(a). The dielectric thin sheets (${ϵ}_r=4$, μ₀) are located in free space (ϵ₀, μ₀). Note that in this case, the permittivity in region 2 is a constant, so the effective tensor permittivity in Eq. (10) is a simple matrix. Hence, the normal update scheme via Eq. (11) can be adopted. The cavity was excited by a Gaussian point source with a bandwidth of 3GHz. To verify the accuracy of the proposed method, we compared the resonant frequency of the structure computed by the tensor FDTD and standard FDTD with a very fine mesh (reference FDTD). In the reference FDTD simulation, the structure was discretized by cells much smaller than the sheet width, with a cell size $\mathrm{\Delta }=0.2$ mm; the same results were obtained with even smaller cells. The cell size of tensor FDTD was chosen as ${\lambda }_{min}/10=v_p/\left(10f_{max}\right)=5$ mm, where v_p is the speed of light inside the dielectric sheets. The width of the thin sheet is $w={\lambda }_{min}/50=1$ mm and the inclination angle is $\theta ={45}^{\circ }$. Fig. 11(b) illustrates the comparison between tensor FDTD and reference FDTD with respect to the spectrum ofthe electric field within the cavity. The results agree well. We also test the accuracy of the method with varying inclination angles of the thin sheets. The sheets with a fixed width of w = 1 mm are tilted from 0° to 45° with a step of ${5.6}^{\circ }$. The errors (compared to the reference FDTD) of three modes (TE${ }_{01}$, TE${ }_{03}$ and TE${ }_{05}$) are shown in Fig. 12(a). The accuracy remains high as θ varies, which proves the effectiveness of our method for arbitrary angles. Finally, we checked the error variation with sheet width. We fixed the sheet angle to $\theta ={45}^{\circ }$ and changed the sheet width from 0.5 mm to 3 mm with a step size of 0.5 mm. Fig. 12(b) illustrates how the error of the three modes varies with sheet width. The results show that the error is very small when the sheet is thin. However, the accuracy decreases as the sheet width increases. This phenomenon is due to the assumption that the electric field $E^{\left(2\right)}$ in region 2 is constant. The variation of the electric field inside the thin sheet can be ignored when the sheet is sufficiently thin. Then, Eq. (5) is satisfied. As the sheet becomes thicker, the electric field differs across the two interfaces, thus Eq. (5) cannot be satisfied. As a result, the accuracy of the method degrades in this case.

 

Fig. 12 (a) Error (compared with converged FDTD with very fine mesh) of tensor FDTD varies by inclination angle from θ = 0° to θ = 45°, with a fixed sheet width w = 1 mm. (b) Error (compared with converged FDTD with very fine mesh) of tensor FDTD varies by width from w = 0.5 mm to w = 3 mm, with a fixed inclination angle θ = 45°.

Download Full Size | PDF

 

Fig. 13 A thin ring of Lorentz media excited by plane wave in a parallel plate waveguide. Actual cells for the tensor ADE-FDTD method are marked as squares.

Download Full Size | PDF

4.2.2. Thin Lorentz medium ring in a waveguide

Finally, let us consider a thin ring filled with a Lorentz medium, inside a waveguide (Fig. 13). The ring has an inner radius of $r_1=3.00$ nm and an outer radius of $r_2=3.11$ nm. For the Lorentz medium parameters in Eq. (13), we use ${\omega }_t=6\times {10}^{16}$ rad/s, ${\omega }_t=4\times {10}^{16}$ rad/s and $\delta =0.28\times {10}^{16}$ Hz. The ring is placed in a parallel plate waveguide, excited by a plane wave with maximum frequency of $3\times {10}^{16}$ Hz. The scattering parameters of the ring are computed by probing the vertical electric field in front of and behind the ring. With the tensor ADE-FDTD method, we introduce cells much larger than the width of the ring, which are marked with red squares in Fig.13. The S-parameters are computed and compared to the standard, staircased FDTD simulation in fine mesh (Fig. 14). The tensor ADE-FDTD method adopts a cell $\mathrm{\Delta }={\lambda }_{min}/15=v_p/\left(15f_{max}\right)=0.667$ mm, where v_p is the wave velocity in vacuum. The standard FDTD converges well when the cell size is smaller than 0.0222 nm. A good agreement between the proposed method and the reference FDTD is observed throughout the simulated frequency range. Finally, the execution times are compared in Table 3. Our method requires less than 10 seconds to compute the scattering parameters, where as the staircased FDTD requires two orders of magnitude greater execution time for the same level of accuracy.

 

Fig. 14 S-parameters comparison between the tensor ADE-FDTD method and standard FDTD with very fine mesh. The cell size of the former is 0.677 nm; the cell size for the latter varies from 0.037 nm to 0.0111 nm.

Download Full Size | PDF

Table 3. Comparison of execution time for the staircased FDTD using various cell sizes and the proposed method using 0.667 nm cells.

View Table | View all tables in this article

5. Conclusion

A generalized tensor ADE-FDTD method for sloped dispersive interfaces and thin sheets has been presented. The improvement in accuracy for sloped dispersive interfaces is validated by computing the plasmonic modes of a one-dimensional silver nanorods array. Compared with the conventional staircased FDTD, our proposed method shows faster and more stable convergence; plasmonic mode patterns are also well resolved. For the sloped thin dispersive sheet, we have compared the proposed method with an ultra-fine mesh standard FDTD simulation, for the scattering parameters of a thin dispersive ring loaded waveguide. We have shown that our method uses much less time to extract accurate results. The proposed method shows great versatility for arbitrarily shaped interfaces and thin sheets. Hence, it is a useful technique for computational electromagnetics and optics.