Enhancing Efficiency of Electromagnetic Simulation in Time Domain with Transformation Optics

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.


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.

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]: 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: and det is the determination of a matrix.The electric and magnetic fields in both coordinates have the following relationship: 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): 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.
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: and det is the determination of a matrix.The electric and magnetic fields in both coordinates have the following relationship: 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 Lx, as shown in Figure 1a.The points {x−n, …, xi, …, xn} are employed to define the domain of each grid, where i is an integer from −n to n.The system length Lx and the points {x−n, …, xi, …, xn} 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 Lu, as shown in Figure 1b.In this time, the length Lu 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, …, ui, …, un} can be used to represent the domain of each grid, where ui = −Lu/2 + (i + n)Δu.Assume that f(x) is known.Then, for a physical grid centered at xi+1/2 = (xi+1 + xi)/2 and the corresponding numerical grid centered at ui+1/2 = (ui+1 + ui)/2, the central difference method can be applied for their partial derivative in Equation (2): where Δxi = (xi+1 − xi) 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 Δxi and Δu, or, simply speaking, the slope of f(x), at x = xi+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.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 Lx = Lz = 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 {X1, …, Xj, …, Xm} are defined as the end points of the regions starting from x = 0, where j is an integer varying from 1 to m, and Xm = Lx/2.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.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 Lx = Lz = 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 {X1, …, Xj, …, Xm} are defined as the end points of the regions starting from x = 0, where j is an integer varying from 1 to m, and Xm = Lx/2.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 1 , . . ., 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: where ∆X j represents the cell size in the region j, and: In this case, X 0 = U 0 = 0.The mapping function should be odd.For x = 0, f (x) = 0; for x < 0, we let 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: 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: where (U m−1 + L u /2)/∆u − n < |i| < n.

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 = λ −1 , 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: Appl.Sci.2018, 8, 1133 6 of 13 where 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: 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 17 rad/s) and γ D = 0.In the FDTD simulation, the EM field updating scheme uses the auxiliary differential equation [8]: 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: where  [36], the altered time-step limitation is given as: 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.: 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: 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.

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.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.Normalized transmittance spectrums from: the cases of uniform grids when Δ = 10 nm (green curve), 5 nm (cyan curve), and 2.5 nm (blue curve), respectively, the converged profile (purple curve) from the linear curve fitting (see text), theoretical prediction (black curve) from [33], the cases of 1D non-uniform grids with Δx1 = 2.5 nm when the mapping functions are multiple-linear (long dashed brown curve) and linear-nonlinear (long and short dashed golden curve), respectively, and the cases of 2D non-uniform grids with Δx1 = Δz1 = 2.5 nm (dashed red curve) and 0.625 nm ( long dashed black curve), respectively, when both of the mapping functions are multiple-linear.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.

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 Δxj = 2.5 × 2 j−1 nm for j = 1 to m.The end points of the divided regions {X1, X2, X3, X4} 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) −1 Δ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 Normalized transmittance spectrums from: the cases of uniform grids when ∆ = 10 nm (green curve), 5 nm (cyan curve), and 2.5 nm (blue curve), respectively, the converged profile (purple curve) from the linear curve fitting (see text), theoretical prediction (black curve) from [33], the cases of 1D non-uniform grids with ∆x 1 = 2.5 nm when the mapping functions are multiple-linear (long dashed brown curve) and linear-nonlinear (long and short dashed golden curve), respectively, and the cases of 2D non-uniform grids with ∆x 1 = ∆z 1 = 2.5 nm (dashed red curve) and 0.625 nm ( long dashed black curve), respectively, when both of the mapping functions are multiple-linear.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.

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 1 , X 2 , X 3 , X 4 } 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) −1 ∆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.
Appl.Sci.2018, 8, x FOR PEER REVIEW 9 of 13 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: where nc is the number of parallel-computing cores and ts is the total hours of simulation.The nc, ts, and C for the uniform and non-uniform cases are listed in Table 2.The nc 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, Δx1 and Δz1 varies from 0.625 nm to 10 nm, and Δxj = Δzj = 2.5 × 2 j−1 for j = 2, 3, and 4. The end points of the regions {X1, X2, X3, X4} are as the same as that used in the previous demonstration while {Z1, Z2, Z3, Z4} are 240 nm, 280 nm, 320 nm, and 4000, respectively.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: 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.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 1 and ∆z 1 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 1 , X 2 , X 3 , X 4 } are as the same as that used in the previous demonstration while {Z 1 , Z 2 , Z 3 , Z 4 } are 240 nm, 280 nm, 320 nm, and 4000, respectively.
The spectrum yielded when ∆x 1 = ∆z 1 = 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 1 = ∆z 1 = 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.The spectrum yielded when Δx1 = Δz1 = 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 nc = 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 Δx1 = Δz1 = 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 nc = 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].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.

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.

Figure 1 .
Figure 1.(a) Number of 2n grids, with various cell sizes, are distributed in the x coordination with length Lx, where {x−n, …, xi, …, xn} defines the domain of each grid, and i is an integer from -n to n,; for a grid centered at xi+1/2 = (xi + xi+1)/2, the physical cell size is Δxi = (xi − xi+1).(b) Number of 2n grids, each with the constant cell size Δu, are distributed in the u coordination with length Lu, where the points {u−n, …, ui, …, un} define the domain of each grid.The length Lu are decided by the designated mapping function u = f(x).

Figure 1 .
Figure 1.(a) Number of 2n grids, with various cell sizes, are distributed in the x coordination with length L x , where {x −n , . . ., x i , . . ., x n } defines the domain of each grid, and i is an integer from -n to n; for a grid centered at x i+1/2 = (x i + x i+1 )/2, the physical cell size is ∆x i = (x i+1 − x i ).(b) Number of 2n grids, each with the constant cell size ∆u, are distributed in the u coordination with length L u , where the points {u −n , . . ., u i , . . ., u n } define the domain of each grid.The length L u are decided by the designated mapping function u = f (x).

Figure 2 .
Figure 2. Schematic side view (a) and top view (b) of the configuration for the finite-difference time-domain (FDTD) simulation employed with transformation optics (TO) to generate non-uniform grids.The aperture width w, length l, and height h are fixed to 100, 300, and 100 nm, respectively.The film extends to the entire xz plane of the physical system, and its lengths are Lx = Lz = 8 μm.The film is illuminated by a broadband incident electromagnetic (EM) wave from the top with the electric field polarized along the x axis.

Figure 3 .
Figure 3. Schematic of the mapping function u = f(x).

Figure 2 .
Figure 2. Schematic side view (a) and top view (b) of the configuration for the finite-difference time-domain (FDTD) simulation employed with transformation optics (TO) to generate non-uniform grids.The aperture width w, length l, and height h are fixed to 100, 300, and 100 nm, respectively.The film extends to the entire xz plane of the physical system, and its lengths are L x = L z = 8 µm.The film is illuminated by a broadband incident electromagnetic (EM) wave from the top with the electric field polarized along the x axis.

Figure 2 .
Figure 2. Schematic side view (a) and top view (b) of the configuration for the finite-difference time-domain (FDTD) simulation employed with transformation optics (TO) to generate non-uniform grids.The aperture width w, length l, and height h are fixed to 100, 300, and 100 nm, respectively.The film extends to the entire xz plane of the physical system, and its lengths are Lx = Lz = 8 μm.The film is illuminated by a broadband incident electromagnetic (EM) wave from the top with the electric field polarized along the x axis.

Figure 3 .
Figure 3. Schematic of the mapping function u = f(x).

Figure 3 .
Figure 3. Schematic of the mapping function u = f (x).

Figure 4 .
Figure 4. Normalized transmittance spectrums from: the cases of uniform grids when Δ = 10 nm (green curve), 5 nm (cyan curve), and 2.5 nm (blue curve), respectively, the converged profile (purple curve) from the linear curve fitting (see text), theoretical prediction (black curve) from[33], the cases of 1D non-uniform grids with Δx1 = 2.5 nm when the mapping functions are multiple-linear (long dashed brown curve) and linear-nonlinear (long and short dashed golden curve), respectively, and the cases of 2D non-uniform grids with Δx1 = Δz1 = 2.5 nm (dashed red curve) and 0.625 nm ( long dashed black curve), respectively, when both of the mapping functions are multiple-linear.

Figure 4 .
Figure 4. Normalized transmittance spectrums from: the cases of uniform grids when ∆ = 10 nm (green curve), 5 nm (cyan curve), and 2.5 nm (blue curve), respectively, the converged profile (purple curve) from the linear curve fitting (see text), theoretical prediction (black curve) from[33], the cases of 1D non-uniform grids with ∆x 1 = 2.5 nm when the mapping functions are multiple-linear (long dashed brown curve) and linear-nonlinear (long and short dashed golden curve), respectively, and the cases of 2D non-uniform grids with ∆x 1 = ∆z 1 = 2.5 nm (dashed red curve) and 0.625 nm ( long dashed black curve), respectively, when both of the mapping functions are multiple-linear.

Figure 5 .
Figure 5. Mesh of the 1D non-uniform grids from (a) the multiple-linear mapping function and (b) the linear-nonlinear mapping function.The thick lines represent the boundaries of the aperture.

Figure 5 .
Figure 5. Mesh of the 1D non-uniform grids from (a) the multiple-linear mapping function and (b) the linear-nonlinear mapping function.The thick lines represent the boundaries of the aperture.

Figure 6 .
Figure 6.Computational cost C from the cases of uniform grids vs. cell size (blue dots) and the cases of non-uniform grids vs. central cell size (red dots), respectively, and the curve-fitting power functions for the cases of uniform grids (blue dashed curve) and the cases of non-uniform grids (red dashed curve), respectively.

Figure 6 .
Figure 6.Computational cost C from the cases of uniform grids vs. cell size (blue dots) and the cases of non-uniform grids vs. central cell size (red dots), respectively, and the curve-fitting power functions for the cases of uniform grids (blue dashed curve) and the cases of non-uniform grids (red dashed curve), respectively.

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.

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.

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

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.