Precise-Integration Time-Domain Formulation for Optical Periodic Media.

A numerical formulation based on the precise-integration time-domain (PITD) method for simulating periodic media is extended for overcoming the Courant-Friedrich-Levy (CFL) limit on the time-step size in a finite-difference time-domain (FDTD) simulation. In this new method, the periodic boundary conditions are implemented, permitting the simulation of a wide range of periodic optical media, i.e., gratings, or thin-film filters. Furthermore, the complete tensorial derivation for the permittivity also allows simulating anisotropic periodic media. Numerical results demonstrate that PITD is reliable and even considering anisotropic media can be competitive compared to traditional FDTD solutions. Furthermore, the maximum allowable time-step size has been demonstrated to be much larger than that of the CFL limit of the FDTD method, being a valuable tool in cases in which the steady-state requires a large number of time-steps.

1. Introduction

Simulation of the electromagnetic wave distribution through periodic media is a popular scenario in diffractive optics and photonics in general. Many optical devices are based on periodicity on at least one dimension (one-dimensionally periodic structures), and they can be used for different applications, i.e., waveguides, splitters, filters, etc. The analysis of this kind of element sometimes can be analysed by analytical closed expressions. However, in some applications, the fine detail of the structure is needed, thus implying a numerical formulation based on finite-element methods (FEM), FDTD or rigorous-coupled-wave theory (RCWT) amongst other numerical formulations. One of the most frequently used is FDTD due to its versatility and accuracy for electromagnetic analysis. However, the FDTD method has two main characteristics that limit the accuracy, stability, and indirectly the applicability of the method. These two factors are the Courant-Friedrich-Levy (CFL) condition and the dispersion error. Limiting the dispersion error implies considering small spatial resolutions compared to the wavelength, thus setting up very fine meshes. In these cases, the CFL condition forces small time-step sizes dramatically, thus resulting in demanding simulations in terms of running time and memory resources.

For improving the performance of FDTD, recently, a new approach for FDTD simulations has arisen. The PITD permits to break through the CFL limit in terms of time-step size. The basic idea follows the same paradigm as the traditional FDTD formulation, where the spatial derivatives are computed through the central finite-difference scheme. This step provides a set of ordinary differential equations (ODEs) that are solved by the precise integration (PI) [1,2]. In recent years, researchers have focused their efforts on taking advantage of this paradigm in different areas using the precise-integration time-domain (PITD) for solving Maxwell’s equations in free and lossy space [3,4,5,6,7,8]. There have been different contributions focused on improving PITD method: e.g., the extension of an unconditionally stable PITD for the numerical solutions of Maxwell’s equations to the circular cylindrical coordinate system was reported by Zhao et al. [9], and the works of Kang et al. related to the application of the PITD method to model the wave propagation in magnetised plasma based on the auxiliary differential equation (ADE) [10,11]. In the last years, some problems related to memory requirements on PITD have been addressed from different perspectives. Zhu et al. arranged the transverse electric/transverse magnetic (TE/TM) components in a matrix, constructing a set of Riccati matrix differential equations (RDEs). This proposal was able to reduce memory requirement compared to the conventional PITD for the analysis of homogenous media [12]. Shao et al. [13] implemented the so-called region-splitting (RS) technique for memory-saving realisations of the perfectly matching layers (PMLs) in the PITD method. However, even considering all these new approaches, the problems related to the memory requirements remain and, in some cases, limit the PITD method’s applicability. It is worth noting the recent contribution of Zhu et al. [14] in which a PITD method with a thresholding scheme is shown in order to reduce the computation cost of the matrix exponential involved on PITD.

To the best of our knowledge, the PITD method has not been explicitly applied to periodic optical media and neither to full anisotropic materials where the tensorial nature of the electrical permittivity is considered. The numerical analysis of periodical media through finite-difference approaches has been previously addressed in the literature. In diffractive optics, usually one-dimensional or two-dimensional structures are considered. Furthermore, there are some applications in which the optical device can have at least one dimension larger in terms of the input wavelength, e. g. holographic volume gratings [15,16]. Here, periodical boundary conditions help to minimise the simulation area, thus increasing the accuracy and performance of the numerical method. It is worth noting the outstanding contributions of Roden et al. [17] introducing a reliable formulation based on the split-field technique for simulating periodic structures at oblique incidence. The extension of SF-FDTD for simulating anisotropic media in 2D and 3D schemes was presented by Oh et al. [18,19], and Miskiewicz et al. [20,21], respectively. The authors have contributed to extending the SF-FDTD for the analysis of nonlinear optical media [22,23] and with the computational acceleration by means of graphic processing units (GPUs) [24]. There are other powerful approaches for the analysis of anisotropic optical media, i.e., finite element method (FEM) [25] or RCWA [26,27].

This work includes details related to the periodic boundary conditions (PBC) and anisotropic media implementation of PITD. In order to validate the approach, a set of well-known optical experiences are simulated, e.g., Young’s double-slit experiment, thin-film filters, dielectric binary diffraction gratings, and twisted-nematic liquid crystal (TNLC) regarding anisotropic optical media. The results show that the PITD method is accurate compared to split-field finite-difference time-domain (SF-FDTD) simulations. Furthermore, some analyses are performed related to the computational performance and also, the capability of using time-step resolution larger than the one established by the CFL condition is corroborated.

2. Formulation

As in the conventional FDTD method, Maxwell’s curl equations for inhomogenous and anisotropic media are expressed as

The electric charge density is defined as . The magnetic field and the magnetic flux density , as well as the electric field and the electric flux density , are connected via the following equations:

The magnetic permeability and the electric permittivity are in the most general case time-dependent and imaginary second-order tensors. The electric current density on Equation (2) is connected to Ohm’s law via: where is the electric conductivity. In order to avoid big differences in the magnitude of the electric and magnetic field, the following change of variable is performed: where and are the relative permittivity, and relative permeability, respectively. For simplifying notation, from now . Taking into account Equations (5) and (6), the new normalized electric field Equation (8) can be introduced into Equations (1)–(4) obtaining the following set of equations: being , and c the speed of waves in vacuum.

If the material has lineal dielectric properties, only three dielectric constants (, and ) and three geometric angles (, , and ) are necessary to specify the full tensor description of : where is the transformation matrix fully defined in Equation (3) in [18] and related to the Euler angles (, , and ). Defining and , the Equations (9) and (10) can be redefined as the following set of equations:

We limit ourselves to the analysis of 2D structures, but it is worth noting that applying this scheme to 3D problems would be done straightforwardly. If the simulation plane is in the plane, Equations (12) and (13) can be extended developing both curl terms in the right side and neglecting the spatial derivatives along z-axis. where , and and are the components of the tensor , and , with , and z, respectively. For solving Maxwell’s equations in FDTD simulations the different electromagnetic field components are staggered in a computational grid space known as Yee’s cell [28,29]. The discrete expressions for Equations (14)–(16) can be written as:

The unit cell has , and sizes. Note that employing a discrete grid implies that some values are not available directly, and hence they must be interpolated from neighbouring grids [18,30]. These additional approximations still hold the second order of accuracy. Note that each component is located in an auxiliary grid-point inside the Yee cell, i.e., is defined in . Hence, the spatial first-order partial differentiation on the right side of Equation (22) must be computed on needing elements outside the grid cell in the grid boundary region. Considering normal incidence, the periodic boundary condition (PBC) scheme is where is the periodicity of the system. The same relationships can be applied for periodical structures along the y-axis. PBC is applied along the periodical axis for one-dimensional periodic structures, whereas perfectly matched layers (PML) are applied along the propagation direction. The implementation of PML for the PITD method is fully detailed in [13].

The discrete expressions for Equations (17)–(19) can be also written as:

The above ODEs can be rewritten as a matrix form: where is a one-column vector containing both all of the electromagnetic field components. The matrix contains the information related to the spatial-step and the medium parameters, and is a column vector with the excitations [3,10]. Here, harmonic plane waves are considered for illuminating the PITD grid. More specifically, a source line is placed parallel to the input plane of the device. This scheme is the same that the one used in standard FDTD simulations, e.g., [18]. The full derivation of the excitation source for polarized plane waves can be found in [31]. PBC is applied in Equations (14)–(19) over the unknown terms next to the boundary (out of the PITD grid). These terms are needed due to the spatial interpolation or the spatial derivative. It is worth mentioning that this paradigm is quite different to the one used in the traditional FDTD scheme, where Equations (23) and (24) can be directly applied after the updating of each electromagnetic field thanks to the leapfrog algorithm [28,29]. Here, this procedure can not be applied since all the components are computed at once when Equation (28) is solved.

Taking into account the theory of ODEs, the analytical solution of Equation (28) can be written as follows. and the discrete from of this equation is where , the time step, and is the exponential matrix which can calculated by using the power rule: where n is a preselected arbitrary integer, such as [7,14]. An approximation of the term between brackets in Equation (31) is given as where is the identity matrix and is computed by the iterative process fully defined in Equations (5)–(7) in [14]. The interested reader in PI technique can find more specific information about the implementation in [7,10]. It is convenient to introduce here the Courant-Friedrichs-Lewy (CFL) condition for standard FDTD formulations to link it with the time step previously defined. The maximum time step in Yee’s FDTD is constrained by the CFL condition, which in 2D and square grid cell reads

The CFL limit becomes particularly restrictive in the presence of small geometrical features, which impose a small cell size and, consequently, a small-time step, leading to long simulations. Here, a smaller value for CFL condition is established as a reference for ensuring stability, so we define . PITD permits overcoming this CFL condition, and the syntax for addressing this feature is defined through the scaling factor as .

Returning to the PITD formulation, once is known from Equations (31) and (32), the right-side of Equation (30) is approximated using the Gauss integration technique

Once matrix is fully defined, a thresholding scheme is applied for the matrix exponential terms in . This procedure is fully detailed in [14] where sparse matrices are used for storing the exponential matrices, and the following thresholding scheme is considered: where the threshold value is and the maximum absolute value of all elements in is In practice, the selection of threshold value is not necessary if is experimentally determined previously through a set of preliminary simulation tests. can be set as a smaller value to achieve higher accuracy. It is worth noting that implies no threshold scheme, hence considering values between and should provide good results in almost all analysis. It should be mentioned that Zhu et al. [14] reported the application of the implicitly restarted ARnoldi method (IRAM) method for dynamically determining the value of . However, this approach has not been considered here due to the computation time that implies on each simulation run. Moreover, the convergence of the IRAM method is sometimes compromised, which is not an optimal solution considering the setup here.

3. Results

The diffraction pattern produced by one single slit or two slits is a very popular experiment in Physics since it shows the wave behaviour of light and is an experience that reinforces the Huygens principle and the relationship between diffraction and interference in Physics. The diffraction efficiency can be easily obtained through analytical closed expressions. However, for simulating this experiment numerically, some add-ons must be included. First, it is worth noting that the diffraction efficiency pattern is measured in the Fraunhofer region. Simulating areas far from the slits in PITD would require a considerable grid cell arrangement being impossible to simulate. Therefore, the near-to-far field transformation [28,29,32] has been implemented for obtaining the irradiance pattern in far-field. This far-field distribution is computed from the near-field electromagnetic field close to the slits computed by the PITD simulation. The analysis of one single slit or even two slits is not strictly a periodic problem. Hence, in this first set of analyses, PBC is disabled. This setup serves as an initial validation of the implementation and benchmark for simulating a typical scenario in diffractive optics where the grid form factor is far from being square and large in one dimension compared to the input wavelength.

Figure 1a shows the scheme of a single slit of length b. is the angle formed between the normal of the slit and the path between the centre of the aperture and the observation point in the same plane.

Figure 1

Schematic diagram of the media considered. (a) Scheme of a single slit of length b. (b) Two thin-film filters of 3 and 5 layers where high and low refractive indices are stacked. (c) Diffraction phase grating [33]. (d) Representation of the spatial variation of and in a TN-LC [18].

Figure 2a shows the diffraction efficiency as a function of the parameter , where b is the length of the aperture, k is the wavenumber. where is the input light amplitude. For two slits, the Equation (36) can be reformulated with an additional term related to the separation of the slits a as: where . The results for two slits are shown in Figure 2b. The PITD setup is summarized in Table 1, where and are the spatial and time step resolutions, respectively. The parameter is the factor applied to to overcome the CFL condition. In this specific case, the relationship between and fits the CFL condition. The time window analysed is , and the grid density (number of cells by wavelength) is defined as . The parameters and define the size of the simulation grid in cells for the y and x axis, respectively. Since the analysis of the diffraction pattern of slits is not a periodic problem, PBC is disabled only here, and PML in all boundaries of the simulation grid are considered. The number of additional cells considered for the absorbing boundary conditions is represented by . The parameters n and are related to the PI formulation.

Figure 2

Simulation of the diffraction pattern of one and two slits of width b by means of PITD. The setup parameters are summarised in Table 1 (a) = b = 80. (b) and .

Table 1

Simulation parameters for Figure 2. for analysis in Figure 2a and for the results shown in Figure 2b.

Δ	α^′	Δt=α^′Δt_0	t_s	N_λ	n_y	n_PML	n	δ
10 nm	1	4.71×10^−8 ns	20 ns	15	150	10	20	10^−5

For modelling the metal wall that defines the aperture a considerably high value for S/F is considered in Equations (17)–(19). For the rest of parameters, with and z, since non-magnetic media is considered, and are defined as unity and with due to be the specific case in which homogeneous and isotropic media is considered.

Figure 3 shows the results of the simulation of one single slit but using a time-step resolution larger than the one established by the CFL condition. More precisely, goes from one (CFL condition) up to six. The sequence Figure 3a–f shows the normalised irradiance as a function of the parameter . It can be seen that the results are accurate in all cases. However, some discrepancies can be identified close to the second-order lobes for . Hence, using larger values for would not be recommendable due to the potential error obtained in that analysis and the drawbacks related to computational costs.

Figure 3

Simulation of the diffraction pattern of one slit of width cells by means of PITD. The setup parameters are summarised in Table 2. Each graph is related to a different parameter: (a) = 1. (b) = 2. (c) = 3. (d) = 4. (e) = 5. (f) = 6.

Table 2 summarises the parameters considered for the PITD simulation in Figure 3, whereas Table 3 shows some parameters related to the computational performance of this set of simulations. is the number of time-steps in order to simulate the 60 ns. As rises, the number of time-steps is reduced. However, considering a constant threshold scheme implies that the size of and also grows since increasing implies larger values in these matrices; hence a larger number of values in this matrix are non-zero and are stored. The running time on Table 3 shows that even reducing the number of time-steps, the rise of the computation time of the matrices implies a significant increase in the total running time. This increase is mainly due to the overhead in the stage of setting up the matrices and and in a lower degree, the impact that these structures induce in the time per iteration that also grows when rises. Simulations performed in server with 2 Intel Xeon E5-2630-2.5 GHz, 15 MB Cache, 32 GB RAM DDR3 1333 MHz, SSD Samsung 840 Pro MZ-7PD512 6Gbps, and Mainboard Asus Z9PE-D8 WS and in MATLAB 2019b. Considering the threshold scheme shown in [14], where a variable parameter is computed is not worth the effort since the calculus of the eigenvalues needed for IRAM introduces a non-negligible penalty in terms of time simulation for the calculus of matrices and . Besides using IRAM, it has been experimentally proven that even considering to scale for each scaling simulation, the overhead introduced by the growth of the memory previously detailed is not compensated. That is why a constant parameter that has been previously determined is a more convenient strategy considering the setup here proposed. For comparing the results in Table 3 a FDTD simulation is performed with an equivalent setup of Table 1. Therefore, the traditional FDTD simulation takes 12.29 s for running a complete simulation with the finest time step for the same grid size (and the same number of time steps). The memory required for storing the simulation is 4.23 MB, and the time per iteration ratio is 9.7. As it can be affirmed, the traditional FDTD scheme is still more competitive than PITD, even the ability of PITD of overcoming the CFL condition considering this setup.

Table 2

Simulation parameters for Figure 3.

Δ	α^′	Δt=α^′Δt_0	t_s	N_λ	n_y	n_x	n_PML	n	δ
10 nm	2	4.71×10^−8 ns	60 ns	60	100	600	30	20	10^−5

Table 3

Running time and memory resources of simulations shown in Figure 3.

α^′	N_t	Size of T(MB)	Running TimeT and T_i (s)	Running TimePITD (s)	Time/IterPITD (s)
1	1273	761.54	123.04	417.99	0.23
2	637	1594.44	236.59	519.13	0.44
3	425	2682.50	365.78	682.88	0.75
4	319	3960.6	584.55	915.78	1.04
5	255	5352.6	1090.52	1536.29	1.75
6	213	6889.0	1693.62	2151.43	2.15

Figure 1b shows the scheme of the thin film filters (TFF) considered. Specifically, high-reflective coatings (HRCs) are is a basic type of TFF and is composed of a stack of alternate high- and low-index films, all one-quarter wavelength (considering as the design wavelength for the structure) thick as it has been illustrated in Figure 1b. Light reflected within the high-index layers does not suffer any phase shift while a change of 180 in the low-index layers is produced [34]. It is straightforward to see that the light (with ) produced by reflection at successive boundaries throughout the assembly reappear at the front surface, all in phase so that they recombine constructively.

Figure 4 shows the results obtained from PITD simulations and the theoretical curves by means of the matrix method fully detailed in [34]. The results show a good accuracy between the two methods considered. The PITD parameters setup are summarised in Table 4. The physical parameters for the refractive indices for the high and low layers arr = 2.3, and = 1.38, respectively. The refractive index of the substrate is = 1.52, and the excitation wavelength is = 230 nm. It is worth mentioning that results show in Figure 4 are computed using taking advantage of the PITD capabilities.

Figure 4

Simulation of the reflectance for two thin-film filters based on three layers and five stacks of layers. The setup parameters are summarised in Table 4.

Table 4

Simulation parameters for Figure 4.

Δ	α^′	Δt=α^′Δt_0	t_s	N_λ	n_y	n_x	n_PML	n	δ
1.15 nm	2	5.42×10^−9 ns	10^−8 ns	200	500	20	20	20	10^−5

Figure 5 shows the results of the PITD method applied to binary diffraction gratings and TNLC devices. More specifically, Figure 5a shows the diffraction efficiency of the 0th- and 1st-orders as a function of the normalised pillar length. The scheme of a binary phase grating is shown in Figure 1c where the parameters are the same ones used in the analysis in [33]: , and . The length h has been varied for obtaining the results shown in Figure 5a. PITD simulation is compared to SF-FDTD results [33], showing that PITD provides close values to SF-FDTD. Some minor deviations are in the curves that can be produced by the spatial average that must be carried on in time-domain simulations and differences between spatial and time resolutions in both methods. Table 5 summarises the PITD setup for simulating the results in Figure 5a.

Figure 5

Comparison between PITD and SF-FDTD results of two periodic media. (a) Diffraction efficiency of the zeroth and first-order of a binary diffraction grating with period , fill factor of 50% and refractive index of [33]. SF-FDTD curve is reproduced with permission from Francés et al., In Proceedings of SPIE Optical Modelling and Design II, Vol. 8429, 2012. (b) Normalized Stokes parameters of a 90 twisted-nematic liquid crystal cell between with parameters: = 0.2, = 1.5 [18].

Table 5

Simulation parameters for Figure 5a.

Δ	α^′	Δt=α^′Δt_0	t_s	N_λ	n_y	n_x	n_PML	n	δ
52.7 nm	1	87.9×10^−9 ns	125·10^−9 s	12	120	24	10	20	10^−5

To demonstrate the potential of the PITD implementation here proposed an anisotropic device is simulated. More precisely, a TN-LC layer is considered. The scheme of this device is shown in Figure 1d. A TN-LC layer is based on LC layer of thickness d that has anisotropy along the y-axis (propagation direction) defined in Equation (38). The TNLC here considered is fully defined in section 3.1 in [18], being , , and . Therefore the matrix is defined as follows: with , being h the thickness of the TNLC layer. Interested reader into the details of the binary phase grating and the TNLC device here considered can find more specific information into [18,33], respectively.

Figure 5b shows the Stokes parameters obtained in the output plane of the TN-LC layer. In both cases (Figure 5a,b), the input light is linearly polarised along the x-axis with a wavelength = 633 nm.

Figure 5b shows the calculated polarization state of outgoing light from the TNLC cell in terms of the normalized Stokes parameters. Table 6 summarises the PITD setup for simulating the results in Figure 5b. The PITD results are compared with those obtained through SF-FDTD simulations and show a good agreement. It is worth noting that the results are consistent with those shown in Figure 3a in [18]. The differences between PITD here are related to the time-step resolution since the computation of Stokes parameter implies the estimation of the phase difference between electric field components. The resolution of the phase difference is closely related to the time resolution. In SF-FDTD, the Stokes parameter was directly related to the inner electric field components since the phase term is implicitly included in the split-field formulation. In spite of this PITD shows a good agreement in this case.

Table 6

Simulation parameters for Figure 5b.

Δ	α^′	Δt=α^′Δt_0	t_s	N_λ	n_y	n_x	n_PML	n	δ
25 nm	1	41.8×10^−9 ns	125·10^−9 s	20	260	12	20	20	10^−5

4. Conclusions

In this paper, the implementation of the PITD method for the analysis of periodic optical media is shown. More precisely, the formulation can simulate anisotropic periodic optical media. The scheme is fully defined, and some canonical examples are simulated and compared with the SF-FDTD method. PITD for optical media is validated through the diffraction efficiency computation of a single and two slits. Here, PITD is tested with time-step resolutions larger than the one established by the CFL condition. The stability of the method and accuracy is good enough for values of the time-step larger than six times than the largest time-step that fits with the CFL condition. A binary phase grating is also analysed regarding periodic optical media, obtaining good results compared to the SF-FDTD curves. Anisotropic periodic media is also covered through the analysis of a 90° TNLC layer. Here, the output polarisation state is analysed and compared to the ones obtained in the literature, showing that PITD is a good approach for simulating this type of devices. It is worth noting that simulating optical media where at least one or two dimensions are larger than the input wavelength can be afforded in the current state without any problems in terms of computational resources.

PITD has shown good accuracy for the analysis of periodic optical media. However, the method still has some drawnbacks related to the formulation. Even using sparse matrices, the memory requirements for time-steps larger than the one established by the CFL condition imply larger values for the matrices, and an experimental tuning of the threshold scheme must be carried out to minimise the overload in terms of computation time and memory requirements. Nevertheless, the authors consider that PITD can become a reference for the analysis of electromagnetic and optical materials where a fine mesh is needed, and hence larger time steps can be applied to reach the steady-state of the structure faster. However, more research and improvements must be developed to outcome the disadvantages that currently has this formulation.