Enhancing Efficiency of Electromagnetic Simulation in Time Domain with Transformation Optics

Featured Application

The main contribution of this work is to propose an electromagnetic finite-difference time-domain method that applies transformation optics for the generation of arbitrary non-uniform physical grids in a simulation system. This method can be applied to the simulation study of light interaction with sub-wavelength structures with efficiency.

Abstract

With sub-wavelength scaled structures in a large system, the conventional finite-difference time-domain method can consume much computational resources since it includes both the spatial and temporal dimension in the scheme. In order to reduce the computational cost, we combine the novel methodology “transformation optics” in the simulation to map a physical coordinate with designated non-uniform grids to a uniform numerical coordinate. For a demonstration, the transmission spectrum through a sub-wavelength metallic aperture with one-dimensional and two-dimensional coordinate transformation is simulated, and compared with uniform-grid cases. We show that the proposed method is accurate, and the computational cost can be reduced remarkably to at most 5.31%, in comparison with the simulation of the finest uniform grids demonstrated. We are confident that it should be helpful to the simulation study in sub-wavelength optics due to its verified accuracy and efficiency.

1. Introduction

Since the extraordinary enhanced optical transmission through a metallic sub-wavelength hole array was discovered two decades ago [1], plasmonics has become an active research area. The academic interest is in the anomalous optical phenomena when light propagates through sub-wavelength structures perforated in metallic film [2,3]. The potential applications include plasmonic circuits [4], bio-sensing [5], optical sensors [6], solar cells [7], and many other optical devices.

The finite-difference time-domain (FDTD) method [8] has been extensively applied for the study of electromagnetism for decades because of its capability to obtain the dynamical changing of the electromagnetic (EM) fields and currents when an EM wave interacts with an object. However, the cell-resolution requirement is high when the characteristic scale of objects to be investigated is smaller than the wavelength of the incident light [9], and the requirement will significantly consume the computational resources. Non-uniform cell algorithms that have been developed for decades, such as the sub-gridding approach [10,11], might be applied to lower the computational requirement; but, the implementation of the source codes is usually not simple.

The novel methodology “transformation optics” (TO) [12,13] is utilized for the coordinate transformation of light, and developed as a tailoring strategy of light for potential optical devices, such as invisible cloak [14,15,16,17,18,19], electromagnetic rotator [20,21], optical black hole or concentrator [22,23,24], planar focusing antennas [25], multi-layered plasmonic filter and coupler [26], panoramic lens [27], etc. The methodology is also applied to study plasmonics, such as the structured plasmonic cylindrical waveguide in a transformed Cartesian coordinate [28] or the frequency- and time-domain response of plasmonic particles under electron beam excitation [29]. For numerical purposes, a procedure for structured nonorthogonal discretization grids implemented in FDTD with the reconsideration of the EM field updating scheme is proposed [30]; the methodology is employed to generate non-uniform grids around the interfaces and openings of a metallic slit in a finite-difference frequency-domain scheme [31]; an additional cylindrical coordinate for the local mesh refinement is introduced in the Cartesian-coordinate-based FDTD simulation [32].

In this work, we propose a method that applies TO for the generation of non-uniform physical grids in the FDTD simulation. The complexity of implementing the non-uniformity will be transformed to the implementation of the derived spatially-varying numerical anisotropic medium (or the so-called optical transformation medium); therefore, the EM field updating scheme is conventional and does not need any further modification for the implementation. The physical coordinate that we study is Cartesian such that the implementation of the coordinate transformation is intuitive and simple. The mapping function between the physical and numerical coordinates can be designed according to the structure of interest. For demonstration, we simulate a rectangular sub-wavelength aperture in a nearly perfect electric conductor (PEC-like) for the transmission spectrum [33], where the Courant condition stability factor is re-considered. Fine cells are designed to be around the aperture while coarse cells are in other areas with the designed multiple linear and nonlinear mapping functions for the one-dimensional (1D) and the two-dimensional (2D) coordinate transformation.

We will show that simulated transmission spectrums are in good agreement with that from the simulation in uniform grids, and in the meantime, the computational costs are remarkably reduced. In the optimized case from the 2D coordinate transformation, only a remarkable 5.31% of the computational cost is required. With this method, we can even utilize finer grids around the aperture. In this case, we obtain the more accurate spectrum while its cost is still less than that in uniform grids. Therefore, we are confident that the method should be helpful to the FDTD simulation study in sub-wavelength optics due to its verified accuracy and efficiency.

Section 2 introduces the transformation for EM waves between two coordinates to obtain the transformation medium and the design of the mapping functions to simulate the transmittance spectrum of a sub-wavelength aperture with the non-uniform grids. The simulation system setup and the reconsidered Courant condition is outlined in Section 3. Section 4 shows the convergence of the simulated transmittance spectrum with uniform grids. We present the resultant transmittance spectrum with the non-uniform grids and compare them with uniform grids for discussing the accuracy and efficiency enhancement in Section 5. Section 6 is the summary.

2. Coordinate Transformation Applied in Simulation and Mapping Functions

The methodology of transforming EM waves between two coordinates can be implemented by utilizing the optical transformation medium, of which the optical properties can be obtained from the following equation [34]:

$${\epsilon }^{\prime }=\frac{J\epsilon J^T}{\det J},{\mu }^{\prime }=\frac{J\mu J^T}{\det J},$$

(1)

where ε and μ are the relative permittivity and permeability tensors of the real medium in the physical coordinate (x, y, z), while ε' and μ' are the tensors of the optical transformation medium in the numerical coordinate (u, v, w), J is a 3 × 3 Jacobian matrix, which can be expressed as:

$$J=\left[\begin{array}{ccc}\frac{\partial u}{\partial x}& \frac{\partial u}{\partial y}& \frac{\partial u}{\partial z}\\ \frac{\partial v}{\partial x}& \frac{\partial v}{\partial y}& \frac{\partial v}{\partial z}\\ \frac{\partial w}{\partial x}& \frac{\partial w}{\partial y}& \frac{\partial w}{\partial z}\end{array}\right],$$

(2)

and det is the determination of a matrix. The electric and magnetic fields in both coordinates have the following relationship:

$$E^{\prime }={(J^T)}^{-1}E,H^{\prime }={(J^T)}^{-1}H.$$

(3)

Hence, the optical transformation medium (usually anisotropic) depends on the derivative of these coordinates. To utilize the methodology in the finite-difference simulation, we consider 2n grids, with various cell sizes, distributed in the physical coordinate x with the system length L_x, as shown in Figure 1a. The points {x_−n, …, x_i, …, x_n} are employed to define the domain of each grid, where i is an integer from −n to n. The system length L_x and the points {x_−n, …, x_i, …, x_n} are predefined according to the structure of interest. On the other hand, the same number of grids, each with the constant cell size Δu, can be used in the numerical coordinate u with the system length L_u, as shown in Figure 1b. In this time, the length L_u is decided by a designated mapping function, i.e., u = f(x). Tailoring the function depends on the non-uniformity of the physical grids (which we will discuss later). Similarly, the points {u_−n, …, u_i, …, u_n} can be used to represent the domain of each grid, where u_i = −L_u/2 + (i + n)Δu. Assume that f(x) is known. Then, for a physical grid centered at x_i+1/2 = (x_i+1 + x_i)/2 and the corresponding numerical grid centered at u_i+1/2 = (u_i+1 + u_i)/2, the central difference method can be applied for their partial derivative in Equation (2):

$${\frac{\partial u}{\partial x}|}_{x=x_{i+1/2}}={f^{\prime }(x)|}_{x=x_{i+1/2}}\approx \frac{u_{i+1}-u_i}{x_{i+1}-x_i}=\frac{\Delta u}{\Delta x_i},$$

(4)

where Δx_i = (x_i+1 − x_i) is the cell size of the physical grid. Hence, the properties of the optical transformation medium to fill the numerical grid can be derived by the relationship between the cell sizes Δx_i and Δu, or, simply speaking, the slope of f(x), at x = x_i+1/2. The similar relationship applies for the other two directions.

Based on the theoretical work in [33], the simulation with the proposed method to obtain the transmission spectrum through a rectangular sub-wavelength aperture, as shown in Figure 2, is considered for the demonstration. The aperture size is fixed to w = 100 nm, l = 300 nm, and h = 100 nm. The lengths of the metallic film in both directions are L_x = L_z = 8 μm.

To properly resolve the aperture in the physical coordinate (x, y, z), fine grids should be applied in the central area while coarse grids can be applied to the side area; thus, the total cell can be intensively reduced in order to save the computational resource. Then, the mapping function can be designed as that depicted in Figure 3, where the slope in the central area is sharp while that in the side area is smooth.

The mapping function could be composed of multiple functions with their neighboring end points connected to divide the system into a number of regions properly. Since the aperture is symmetrical, we discuss the x ≥ 0 domain at first. In this demonstration, the domain is divided into m regions, and the points {X₁, …, X_j, …, X_m} are defined as the end points of the regions starting from x = 0, where j is an integer varying from 1 to m, and X_m = L_x/2.

For simplicity, we assume that each region has a constant cell size. Then, the mapping function can be expressed as the combination of the multiple-linear functions. For x > 0, we obtain:

$$u=f(x)=\{\begin{array}{cc}\frac{\Delta u}{\Delta X_1}x,\hfill & 0<x\le X_1\hfill \\ \frac{\Delta u}{\Delta X_2}x+U_1-\frac{\Delta u}{\Delta X_2}{X}_1,\hfill & X_1<x\le X_2\hfill \\ ⋮& \\ \frac{\Delta u}{\Delta X_j}x+U_{j-1}-\frac{\Delta u}{\Delta X_j}{X}_{j-1},\hfill & X_{j-1}<x\le X_j\hfill \\ ⋮& \end{array}$$

(5)

where ΔX_j represents the cell size in the region j, and:

$$U_j=U_{j-1}+\frac{\Delta u}{\Delta X_j}(X_j-X_{j-1}).$$

(6)

In this case, X₀ = U₀ = 0. The mapping function should be odd. For x = 0, f(x) = 0; for x < 0, we let f(x) = −f(−x).

In another demonstration, we could also use a nonlinear function to generate x-varying cell size in the most outside region X_m−1 < x < X_m since it is considered less important to the accuracy. Since the gradual decrease of the slope is required, an inverse tangent function might be suitable. Thus, in the region, the function can be expressed as:

$$u=\frac{1}{\alpha }{\tan }^{-1}(\frac{x}{\gamma }),$$

(7)

where γ and α can be determined by the end point of the previous region x = X_m−1 and the minimum slope we design for the boundary x = X_m. In this case, U_m is decided mathematically by these conditions. As a result, the x-varying cell size of the grid centered at x = x_i+1/2 in this region is:

$${\Delta x_i|}_{x=x_{i+1/2}}=\frac{\alpha \gamma }{\cos (\alpha u_{i+1/2})}\Delta u,$$

(8)

where (U_m−1 + L_u/2)/Δu − n < |i| < n.

3. Simulation System and Courant Condition

A FDTD parallel computing package, Meep [35], is employed for our simulation. In this package, light speed in vacuum and metric scale are normalized to unity and 1 μm, respectively. Thus, we can define the unit of time as “Meep-time,” a time interval that light in vacuum can travel 1 μm. Accordingly, the relationship between the frequency f and wavelength λ of an EM wave in vacuum is simply as f = λ⁻¹, where the units of λ are in micrometers.

The system height in both physical and numerical coordinates is fixed to 0.2 μm because we are interested in the transmission at the aperture exit. While L_x = L_z = 8 μm, the system lengths L_u in the u and L_w in the w direction are estimated according to the mapping functions.

A broadband incident EM wave source is located above the metallic film and extends in the entire xz plane for estimating the normalized transmittance spectra, and the incident electric field is polarized along the x axis. Its distribution in the wavelength domain is set to be a rectangular function ranged from λ_a = 0.1 μm to λ_b = 1.2 μm; from the inverse Fourier transform, the source in the time domain is given as:

$$A(t)=\frac{1}{\pi (t-t_R/2)}\sin [2\pi f_0(t-t_R/2)]\sin (2\pi \Delta f(t-t_R/2),$$

(9)

where f₀ = (1/λ_a + 1/λ_b) / 2, Δf = (1/λ_a − 1/λ_b), t_R is the running-time, and t_R/2 is the time-delay. The singularity t = t_R/2 in the denominator can be avoided since we use the E_x component as the wave source, and it is generated in every half integer time-step. The Fourier transform will be performed for the simulated electric and magnetic field in the time domain to the wavelength domain. According to our test, t_R should be equal to or larger than 16 Meep-time for the convergence of the transform.

The perfect-matched-layer (PML) boundary condition is employed to the v direction with a thickness of 10 layers. In the u and w directions, however, the periodic boundary condition is applied. In the condition, the coupling due to the periodicity will be too weak to experience an influence because the system size is large enough, On the other hand, if the PML condition is used, then the edge effect at the boundaries due to the discontinuity of the source fields will generate non-physical waves that travel back and forth in the system, and thus cause unwanted resonance.

The film is assumed to be PEC in [33]. However, for simulations with TO, the built-in PEC material is not suitable for the application because its permittivity is fixed to −∞, which does not meet the need of altering the properties of medium. Instead, an artificial dispersive medium simulated with the Drude mode is employed as the PEC-like material. In the Drude mode, the relative permittivity of the medium is:

$$\epsilon ={\epsilon }_{\infty }-\frac{{\omega }_p^2}{{\omega }^2+i{\gamma }_D\omega },$$

(10)

where ε_∞ is the relative permittivity at infinite frequency, ω_p is the plasma frequency, and γ_D is the damping coefficient. To approximate the PEC, we let ε_∞ = 1, ω_p be sufficiently high (e.g., 6π × 10¹⁷ rad/s) and γ_D = 0. In the FDTD simulation, the EM field updating scheme uses the auxiliary differential equation [8]:

$$\frac{\partial J_p}{\partial t}+{\gamma }_D{J}_p={\epsilon }_0{\omega }_p^{2}E$$

(11)

in synchronism with Maxwell’s curl equations to yield a composite self-consistent system, where J_p is the polarization current associated with the wave interaction with the dispersive medium. According to our test, the numerical transmittance spectrums without coordinate transformation show no difference when compared with those using the built-in PEC material. For the coordinate transformation, the dispersive medium will become anisotropic. The updating scheme for this kind of medium has been implemented in Meep. Therefore, the conventional FDTD algorithm does not need any further modification to employ the simulation with non-uniform grids.

The maximum time-step is limited by the cell size in the x, y, and z directions, known as the Courant condition [8], to avoid numerical instability. For the optical transformation medium used in the numerical grids to map to smaller physical grids, the numerical phase velocity of light is greater than that in the physical space because the coordinate transform is invariant of time. Under the circumstance, the Courant condition should be re-considered. Assuming that ε′ and μ′ are obtained from Equation (1) by the mapping functions u = f(x), v = g(y), and w = h(z), then they can be written as the diagonal matrixes:

$${\epsilon }^{\prime }=\left[\begin{array}{ccc}{\epsilon }_u& 0& 0\\ 0& {\epsilon }_v& 0\\ 0& 0& {\epsilon }_w\end{array}\right],{\mu }^{\prime }=\left[\begin{array}{ccc}{\mu }_u& 0& 0\\ 0& {\mu }_v& 0\\ 0& 0& {\mu }_w\end{array}\right],$$

(12)

where ε_u = ε(∂u/∂x)/[(∂v/∂y)(∂w/∂z)], ε_v = ε(∂v/∂y)/[(∂u/∂x)(∂w/∂z)], ε_w = ε(∂w/∂z)/[(∂u/∂x)(∂v/∂y)], μ_u = μ(∂u/∂x)/[(∂v/∂y)(∂w/∂z)], μ_v = μ(∂v/∂y)/[(∂u/∂x)(∂w/∂z)], and μ_w = μ(∂w/∂z)/[(∂u/∂x)(∂v/∂y)]. From [36], the altered time-step limitation is given as:

$$\Delta t\le \frac{\Delta }{c}{\left[\frac{1}{{(\Delta u/\Delta )}^2\sqrt{{\epsilon }_v{\epsilon }_w{\mu }_v{\mu }_w}}+\frac{1}{{(\Delta v/\Delta )}^2\sqrt{{\epsilon }_u{\epsilon }_w{\mu }_u{\mu }_w}}+\frac{1}{{(\Delta w/\Delta )}^2\sqrt{{\epsilon }_u{\epsilon }_v{\mu }_u{\mu }_v}}\right]}^{-1/2},$$

(13)

where Δ is the normalizing cell size. If no coordinate transformation is taken (x = u, y = v, z = w), and the medium is vacuum (ε = μ = 1), this inequality can be reduced to conventional Courant as indicated in [8], i.e.:

$$\Delta t\le \frac{\Delta }{c}{\left[\frac{1}{{(\Delta x/\Delta )}^2}+\frac{1}{{(\Delta y/\Delta )}^2}+\frac{1}{{(\Delta z/\Delta )}^2}\right]}^{-1/2}.$$

(14)

The Courant condition stability factor can be defined as S = cΔt/Δ. If we let the mapping functions be u = 5x, y = v, and w = 5z, and let the cell size be Δu = Δv = Δw = Δ, we obtain the stability factor S = 51^−1/2 (~ 0.14). On the other hand, in the physical coordinate, the corresponding cell size in the x, y and z directions becomes Δx = 0.2Δ, Δy = Δ, and Δz = 0.2Δ, respectively. From Equation (14), the same stability factor S is obtained. Therefore, Equation (13) is valid. Note that in Meep, the time-step Δt is decided by the stability factor S and Δ.

According to the definition in [33], the transmittance spectrum is the amount of the transmitted power flux that passes through the aperture exit and normalized to that passes the same area in the absence of the film. In the simulation, the power flux in the v direction is:

$$P(\omega )=\mathrm{Re}\{-\widehat{v}\cdot {\displaystyle {\int }_s{E}_{\omega }(u,v=-h/2,w)\times H_{\omega }(u,v=-h/2,w)ds}\},$$

(15)

where s is the area of the aperture exit, and E_ω and H_ω is the Fourier transform of the simulated electric and magnetic field vectors; the normalized transmittance is obtained from the measured P(ω) divided by that measured at the same location without the film. Note that the power flux remains the same in both numerical and physical coordinates because: according to Equation (3), the electric and magnetic field in the numerical coordinate should be transformed to that in the physical coordinate while the cell size is also correspondingly changed. Since both the field change and cell size change cancel each other, the power flux is not affected by the coordinate transformation.

4. Convergence of Transmittance Spectrum with Uniform Grids

The accuracy of the simulated spectrum depends on the cell resolution, and it should converge when the cell size approaches to zero. To test the convergence, the simulation in uniform grids with the similar system setup in Section 3 are performed to obtain the transmittance spectrum T(λ); the cell size is set to be constant in all directions, i.e., Δx = Δy = Δz = Δ. In this test, the metal is the built-in PEC. For each uniform case, we use the stability factor S = 0.5, as the default setting in Meep.

Figure 4 shows the results yielded when Δ are 10 nm, 5 nm, and 2.5 nm, respectively. As seen, while the profiles are similar, they converge with Δ. The wavelength at peak transmittance λ_p and the corresponding peak T(λ_p) for each Δ is listed in

Table 1. Wavelength at peak transmittance λ_p and corresponding peak T(λ_p) for the cases of uniform grids, 1D non-uniform grids, and 2D non-uniform grids.

	Uniform Grids			1D Non-Uniform Grids		2D Non-Uniform Grids
	10 nm	5 nm	2.5 nm	Multiple-Linear (2.5 nm)	Linear-Nonlinear (2.5 nm)	Multiple-Linear (2.5 nm)	Multiple-Linear (0.625 nm)
λp (nm)	727	678	657	665	665	665	642
T(λp)	4.14	3.87	3.64	3.83	3.80	3.82	3.68

. The convergence for both λ_p and T(λ_p) is almost linear. Suppose that they have the relationships λ_p = aΔ + b and T(λ_p) = cΔ + d, respectively, where a, b, c, and d are constants. We obtain a = 9.44, b = 632, c = 0.06, and d = 3.51 when the curve fitting is utilized. In other words, if Δ approaches to zero, it is expected that λ_p and T(λ_p) will approximate to 632 nm and 3.51, respectively. To verify, the transmittance spectrum predicted in [33] is reproduced and shown in Figure 4. From the theory, we obtain λ_p = 626 nm and T(λ_p) = 3.49, respectively, which are well agreed with our approximation.

Besides λ_p and T(λ_p), the rest parts of the spectrums in Figure 4 are also seen to converge linearly with the cell size at relative wavelengths, i.e., for the wavelength ratio r, rλ_p and T(rλ_p) have similar linear relationships with Δ to those when r = 1. With curve fitting for each r from 0.6 to 1.3, the approximated transmittance spectrum is shown in Figure 4. The result is close to that from the theoretical prediction. Therefore, it can be considered as the converged profile from the simulation in uniform grids.

5. Results from 1D and 2D Non-Uniform Grids and Their Efficiency

The results from the mapping function x = f(u) to the 1D non-uniform grids are demonstrated, and the numerical cell size in all direction are set to be Δu = Δv = Δw = Δ = 2.5 nm. For the multiple-linear mapping function, we let m = 4, and Δx_j = 2.5 × 2^j−1 nm for j = 1 to m. The end points of the divided regions {X₁, X₂, X₃, X₄} are 80 nm, 120 nm, 160 nm, and 4000 nm, respectively. The mesh of the non-uniform grids around the aperture is shown Figure 5a. The cell size in the x direction increases as x increases. For the linear-nonlinear mapping function, the cell sizes and end points in the first three regions are identical to those used in the multiple-linear function while γ and α used in Equation (8) are 1.85 and 0.79, respectively, when the units of x and u are in micrometers. The mesh of the non-uniform grids is shown Figure 5b. For |x| ≤ 160 nm, the mesh is same as that shown in Figure 5a. For 160 nm < |x| < 200 nm, the slope of Equation (7) in this range is about 0.68. The cell size in the x direction is therefore about Δx = (0.68)⁻¹Δu = 3.68 nm. The slope will decrease with |x| such that the cell size will increase. The film is simulated with the PEC-like material (see Section 3). In these cases, the stability factor is set to be S ≅ 0.45, slightly smaller than the default 0.5, due to the stability requirement of the PEC-like material.

The transmittance spectrums are shown in Figure 4. Both spectrums coincide with each other. The wavelength at peak transmittance λ_p and the corresponding peak T(λ_p) are listed in Table 1. Between these two cases, we obtain that λ_p is the same and T(λ_p) has the less than 1% difference. Moreover, they are in good agreement with those from the uniform 2.5-nm case. Therefore, the results are accurate. We verify that the cell size resolution in the central area is important to the accuracy. It should be noted that the linear-nonlinear case implies the flexibility of the mapping functions for complex configurations.

To quantify the efficiency, we define the computational cost for each simulation:

C = n_c × t_s,

(16)

where n_c is the number of parallel-computing cores and t_s is the total hours of simulation. The n_c, t_s, and C for the uniform and non-uniform cases are listed in

Table 2. Number of parallel-computing cores n_c, total hours of simulation t_s, and computational cost C for the cases of uniform grids, 1D non-uniform grids, and 2D non-uniform grids.

	Uniform Grids			1D Non-Uniform Grids		2D Non-Uniform Grids
	10 nm	5 nm	2.5 nm	Multiple-Linear (2.5 nm)	Linear-Nonlinear (2.5 nm)	Multiple-Linear (2.5 nm)	Multiple-Linear (0.625 nm)
nc	48	48	144	96	144	48	96
ts	0.37	3.39	23.54	8.06	15.62	3.75	21.76
C	18	163	3390	774	2249	180	2089

. The n_c is decided by the total number of cells. While the uniform cases shows the exponential growth of computational cost, the multiple-linear and linear-nonlinear cases only require 22.83% and 66.34% of the cost in the uniform 2.5-nm case. This confirms the efficiency of the non-uniform grid transformation.

The non-uniform grids for the z direction is further implemented while the numerical cell size remains the same. At this time, the mapping functions x = f(u) and z = h(w) are demonstrated only for the multiple-linear case. For both functions, m = 4, Δx₁ and Δz₁ varies from 0.625 nm to 10 nm, and Δx_j = Δz_j = 2.5 × 2^j−1 for j = 2, 3, and 4. The end points of the regions {X₁, X₂, X₃, X₄} are as the same as that used in the previous demonstration while {Z₁, Z₂, Z₃, Z₄} are 240 nm, 280 nm, 320 nm, and 4000, respectively.

The spectrum yielded when Δx₁ = Δz₁ = 2.5 nm is shown in Figure 4; in this case, the stability factor is also S ≅ 0.45 due to the same stability requirement of the PEC-like material. The spectrum coincides with those from the cases of the 1D non-uniform grids. The λ_p and T(λ_p) of the spectrum are listed in Table 1; we find that λ_p is identical to the 1D cases and T(λ_p) has negligible difference. Since the cases of the 2D and 1D non-uniform grids have the same resolution in the central area, this is expected. However, in this 2D case, the computational cost is remarkably reduced to C = 180 while n_c = 48, as listed in Table 2. Therefore, with similar accuracy, the 2D transformation requires only 5.31% of the cost in the uniform 2.5-nm case.

We show the spectrum yielded when Δx₁ = Δz₁ = 0.625 nm in Figure 4; in this case, the stability factor S = 0.17 due to the numerical stability requirement of the transformation medium. As seen, the spectrum moves toward that at convergence previously obtained based on the spectrums of the cases of the uniform grids. Because the case of the 2D non-uniform grids has the higher resolution in the central area, the yielded spectrum expectedly becomes closer to that at convergence. The λ_p and T(λ_p) of the spectrum are listed in Table 1; we find that λ_p = 642 nm and T(λ_p) = 3.68, similar to those of that at convergence, λ_p = 632 nm and T(λ_p) = 3.51. The computational cost is C = 2089 while n_c = 96, as listed in Table 2. Although the number of time-steps increases (due to the lower S) in this case, we still obtain the less computational cost than that of the uniform 2.5-nm case.

To analyze the increasing rate of the computational cost due to cell resolution, the C for the cases of uniform grids and non-uniform grids as functions of cell sizes are shown in Figure 6. Two power functions are used to fit the curves. According to the results, we obtain that the increasing rates are the cell size ratio to the power of −3.55 and −1.86, respectively. For the uniform cases, the power term being close to −4 is expected since the FDTD algorithm deals with both the space (three dimensions) and time (one dimension); all the numbers of cells in the three directions and also the number of time steps increase proportionally with the inverse of cell size. On the other hand, for the non-uniform cases, the cost is increased mainly due to the cell number only in the x and z directions; although the S becomes lower when the cell size is smaller than 2.5 nm, the iteration time in one time step is reduced more than that in the uniform case (because the cell number is reduced). Therefore, we obtain that not only the cost is always lower than that in the uniform case with the same cell resolution, but the increasing rate is also two orders slower.

Real metals, such as gold or silver, are usually considered dispersive mediums and their optical properties can be described with the Drude model, as indicated in Equation (10). However, the plasma frequency is close to the wave frequency of the optical regime and the damping coefficient is not negligible. Since sub-wavelength structures in real metals could exhibit enhanced transmission due to the excitation of surface plasmons [1], a sub-wavelength aperture in a real metal film associated with the properties is often considered for the more extensive purposes [37].

Since the metals are dispersive, the simulation of the EM wave interaction with the metallic structures requires the auxiliary differential equation, as shown in Equation (11), in the EM field updating scheme. With the proposed coordinate-transformation method for the simulation with non-uniform grids, the dispersive medium becomes anisotropic. Although the metal film in our current study is PEC, we have already employed an artificial dispersive medium with sufficiently high plasma frequency and zero damping coefficient to approximate the film (see Section 3), and obtained the accurate transmission spectrums with efficiency. Because the implementation for the real metal films is the same, we believe that this method should also be valid in the study of the light transmission for the case of real metals with accuracy and efficiency.

6. Summary

In summary, we have shown the capability of the FDTD simulation with TO for non-uniform grids to enhance efficiency. The non-uniformity is accomplished by the derived spatially-varying numerical anisotropic medium from the designated mapping function. In the demonstration, we used fine grids around the aperture in a large film of a PEC-like material, and coarse grids at other areas to reduce the computational cost. To obtain converged transmittance spectrum, simulations in uniform grids are performed; the resultant spectrums are used to approximate that when the cell size approaches to zero. In comparison, the approximation is close to that reproduced from the theoretical work in [33]. For the 1D coordinate transformation cases with the multiple-linear and linear-nonlinear functions, the normalized transmittance spectrums agree with that from the uniform-grid simulation when the central cell sizes are the same, but both the cases require only 22.83% and 66.34% of the computational cost, respectively, showing the accuracy and efficiency of the method and implying its flexibility for complex configurations. The efficiency is further improved by performing the 2D coordinate transformations; in this case, the cost is remarkably reduced to 5.31% while the accuracy remains. A finer central cell size of 0.625 nm is utilized for the more accurate spectrum while its cost is still less than that of the uniform case. The increasing rate of the computational cost due to cell resolution is also analyzed; the analysis shows that the increasing rate of the cost is two orders slower than that in the uniform case. According to the results, we are confident that the method should be helpful to the FDTD simulation study in the sub-wavelength optics with accuracy and efficiency.