Comprehensive Numerical Analysis of Finite Difference Time Domain Methods for Improving Optical Waveguide Sensor Accuracy

The optical planar waveguide is a main element of photonic devices and performs all the optical signal operations, which can include guiding, coupling, switching, splitting, multiplexing and de-multiplexing. The planar waveguide components, electro-optical components, transmitters, receivers, and driving electronics can be integrated into a single chip using planar technology, similar to microelectronics. Although the operation of optical waveguide sensors and components has been researched and is perfectly understood, their particular performance relies on many parameters, including the geometry, wavelength, initial field distribution, material data, and electro-optic driving conditions. These parameters must be optimized before fabricating a device. For a large-scale optoelectronic circuit, accurate modeling is a predominant issue because of the numerous resources required to fabricate a chip. Optical planar waveguide design relies on simulating the propagation of light signals, waveguide modes, mode coupling, losses and gains. The 1st part of the entry data is used to configure numerical calculations. Ideally, entry systems hide or limit the details of the numerical calculations. The 2nd part of the entry data defines the waveguide device in terms of its geometry, fabrication parameters, and material constants. Using an effective simulator with a project layout that can also handle fabrication parameters is the best way to process waveguide data. However, because waveguide modeling often uses sophisticated numerical algorithms, you must be familiar with some aspects of the underlying numerics. Planar waveguides are the building blocks of photonic circuits. The width of a waveguide is defined perpendicular to the path of the waveguide centerline.

The numerical analysis of the optical waveguide sensor’s accuracy (the analytic models of the planar concentrator with cylindrical lenses) has been already derived by Karp. Accordingly, the geometrical concentration factor can be defined as the ratio of the thickness h to the length l. The analytic planar concentrator’s model can be applied to a single-axis planar concentrator [11]. The concentration parameter of the cylindrical lens is the main difference, which will be very small in that case. All the parameters were defined for various positions P and waveguide coupling angles Φ, which was integrated over all the propagation angles in the optical planar waveguide. In addition, the light coupling efficiency was defined for the waveguide coupling angles Φ along the waveguide at point P. As illustrated in Figure 2 and Equations (1)–(4), the commonly used Kretschmann’s prism beam configuration was applied to the metal film. Kretschmann et al. first demonstrated the attenuated total reflectance (ATR) technique. In principal, the SPW exciting process of an OWS is very similar to the Kretschmann ATR coupler process. The wave rays excite the SPW at the external metal’s interface when a guided mode and the SPW are phase matched. An optical waveguide guides the light waves that enter the region through the use of a thin overlayer of metal; the light evanescently penetrates the metal layer. This behavior can be depicted using the following equation:

The reflection coefficient of the coupling prism and the material absorption is calculated using:

The complete concentrator’s optical efficiency was derived by integrating over all the coupled angle points along the waveguide at point P, as in Equation (3):

Figure 2 shows the parameters used to calculate the efficiency of the optical system. It shows that the light concentrates in the optical planar waveguide.

Ύ is the metallized coupling prism’s reflection coefficient, r is the half pitch of the cylindrical lens, α is the attenuation coefficient of the waveguide material, h is the waveguide thickness, l is the waveguide length, and P is the distance to an exit surface (a specific position inside the waveguide). A cylindrical lens array is placed in the waveguide, and a coupling prism is placed at the central point of the lens for coupling sunlight into the waveguide. Accordingly, the light rays travel along the waveguide and exit on the waveguide’s side; the neglected integration limit was equal to 0, whereas the upper limit was Φ_max, as illustrated in Equation (4) by using the same parameters as in the previous equation:

In the last equation, the variable (f/θ) is the ratio between the focal length and the diameter of the cylindrical lens, the constant n is the waveguide refractive index, and the angle θ is the half angular extension of the sun (± 0.26°). A substantial issue of Equation 1 requires careful note. This issue is the reliance of a ray’s coupling angle Ф on the optical efficiency of the optical waveguide. The ray’s coupling angle Ф impacts the value of $\left(\raisebox{1ex}{$2H$}\!\left/ \!\raisebox{-1ex}{$\tan \Phi $}\right.\right)$, which represents the distance between double backside interactions. Accordingly, larger-angle coupled rays have a much greater impact along the backside surface, which will greatly increase the probability that a ray will escape the waveguide by striking the 2nd coupling prism during its propagation [12]. Usually a dark line will appear in the band when surface plasmon resonance occurs at specific spread angles [13].

Given the behavior of the rays as a function of the input wave vector, this structure will generate diffraction orders and act as a diffraction grating, and terms m= {…,−2,−1,0,1,2,…} will be generated in the reflected light. The optical solution technique based on the real time detection (RTD) method delivers unique results across diverse categories of related research. The resonance causes the light ray to be absorbed. The surface plasmon (SP) was set to resonate with the light when the light beam entered the metal film at a resonance angle. The tunability of the incoming wave vector to match the wave vector is an important prerequisite for surface plasmon (SP) excitation. Currently, kx can be set to any value between the line pairs labeled (L1) and (L2). One such line, labeled (L3), is shown. Lines (L1) and (L2) are the dispersion relations for “normal” light in the 2 mediums (ε₁) and (ε₃), and they both depend on the incidence angle α in the experimental setup shown in the inset. By varying the incidence angle α, any line (L3) between the two lines (L1) and (L2) can be achieved.

A view of the top of the surface plasmon resonance (SPR) sensor and the identical cross-section of the silicon-on-insulator (SOI) rib waveguide are illustrated in Figure 3. It contains four components: the V-shaped bend waveguide, input/output waveguides, a double deeply etched facet coating with a metal layer (a surface plasmon resonance interface) and double air trenches. The third component is used to connect the first component with the second component; therefore, double bends in the silicon-on-insulator (SOI) rib waveguides are formed. Thus, the second component is in a straight line and parallel to the third component.

d is the distance between the 3rd component and the centerline of the 2nd component, and the half-angle of the 1st component is θ. L and W represent the width and length of the fourth component, respectively. The surface plasmon resonance array’s sensor could be composed to implement the wavelength interrogation of the sensor’s output signal using the features of silicon micro-fabrication technology and the silicon-on-insulator rib waveguide. Accordingly, bulky and expensive equipment, such as spectrometers, are not necessary; a spectrometer optimizes the surface plasmon resonance array’s sensor for the accurate integration and miniaturization of the complete sensing system [14]. The surface plasmon resonance is completely excited by the guiding mode at the third component (the bend waveguide interface). Accordingly, a bimetallic configuration was employed. The two-dimensional surface plasmon resonance sensor can be simulated and numerically modeled using the 2D-FDTD method with perfectly matched layers (PMLs). The effective index method (EIM) is the fundamental technique for that purpose. The silicon-on-insulator (SOI) rib waveguide structure was numerically modeled as a two-dimensional structure. It is well known that FDTD modeling results in high accuracy and requires large amounts of computer memory and computational time. Basically, the surface plasmon resonance excitation mechanism of a surface plasmon resonance sensor based on a silicon-on-insulator rib waveguide is similar to that of the Kretschmann prism geometry (classical surface plasmon resonance model); the improvement in the resonance condition is illustrated by the following equations: where ${\lambda }_0$ is the optical wavelength in free space, ε_m is the dielectric constant of the metal solution, ε_d is the dielectric constant of the analyte solution, and n_eff is the effective refractive index of the guiding mode of the silicon-on-insulator rib waveguide. Finally, the surface plasmon resonance sensor (SPR) based on a silicon-on-insulator (SOI) rib waveguide is categorized as a special model of Kretschmann that can include a p-polarized incident beam evolved from the TE-polarized guiding mode of the silicon-on-insulator (SOI) rib waveguide and a virtual prism whose refractive index is equal to the effective refractive index of the guiding mode.

Generally, this method is an accurate second-order technique for analyzing an optical waveguide sensor (OWS) case. Moreover, another method, the total residual chlorine (TRC) method [7], has already been developed and applied to many optical sensor applications. This method is illustrated in [3] for a second-order Lorentz dielectric. Furthermore, an improved piecewise linear recursive convolution (PLRC) method has been developed and applied to dispersive dielectrics [4]. The improved method can be easily applied to optical waveguide sensors. The capability functions of the time-domain (TD) approach from [3] allowed the analysis of the prerequisite differential equation’s parameters, coefficients, and arithmetic operations, and the formulas presented in [4], which form a solid basis for the PLRC and TRC methods, have much greater accuracy, efficiency, and precision than the conventional recursive convolution method [8]. Similarly, numerical experiments have illustrated that the conventional recursive convolution method has far lower accuracy, efficiency, and precision than the other methods mentioned herein. Therefore, we will consider the PLRC method in this analysis. For the piecewise linear recursive convolution technique, (Eⁿ) usually has piecewise linear functional dependence per time period (Δt). Eⁿ/Δt is constant. The updated equation of the electrical field in the metal can be expressed as: where:

The Young’s direction implicit method is based on several approximations, including substituting ${ℐ}^{n+1/2}$ for $ℐ$(t) in the 2nd Maxwell equation; from Maxwell’s equations: where (v) and (w_p) are the electron-natural collision frequency and photon frequency, respectively. A completely equivalent equation form was used in [2], which led to the following difference equation:

Young’s direction implicit method requires that the electrical field component (Eⁿ) is stored at only a single time. Young’s DI method uses very few multiple components, but it uses more operations/Δt than any of the other methods, as illustrated in Nickisch- Franke’s direction implicit (NFDI) method [2]. This method is slightly different from the original method because the term $\frac{J^{n+1}+J^n}{2}$ was included for J(t):

Equations (10) and (11) are accurate second order approximations [15]. Furthermore, the storage requirements and capacities for the NFDI method are larger than those for the original DI method. The derivation of the RC method is completely different; for double ($ℐ$) time levels, it might be fully stored, but the iteration of the result requires a coupled set of first, second and third order differential equations, as in the Young’s direction implicit case. The new technique was based on the integral relating the time domain, namely, (D) and (E), rather than on the ordinary relations of the differential equation, namely, ($ℐ$) and (E).

The included comparisons are quite similar to the first order dispersive dielectric FDTD method in the Steven direction implicit method, which has already been derived in [16]. The NFDI method, which was called the new direction implicit technique, including the term (${ℐ}^{\mathrm{n}+1}+ {ℐ}^{\mathrm{n}} ) /2$ for $ℐ$(t), can finally be presented:

The principle of this method that has been used for solving parabolic differential equations is utilized. However, a modified alternating direction implicit technique is developed and compared to the typical ADI techniques. The mixed coordinates are altered at each time step rather than each coordinate at every time step. Consequently, in the 3D case, the requirement is simply two changes in the solution identification. Therefore, the Courant-Fredericks-Levy (CFL) stability constraint is completely mitigated by the finite difference time domain methods, resulting in the modified ADI algorithm. The CFL stability conditions do not limit the FDTD time step as the modeling precision of the FDTD methods because of the mixing between the ADI and FDTD techniques. The ADI-FDTD approach is unconditionally stable based on a theoretical proof, and this proposed approach is efficient and effective, as predicted by numerical results. For a similar numerical precision, the proposed ADI-FDTD requires at least four times less numerical iteration than that of a typical FDTD method. The FDTD technique is very significant in electrical engineering. Eventually, the number of iterations and CPU time were decreased [15]. Nevertheless, the CFL stability constraint and numerical dispersion errors restrict the high Q structures and geometrical efficiency. The ADI-FDTD approach was developed to mitigate the Courant-Fredericks-Levy stability constraint. The error minimized ADI-FDTD and enhanced ADI-FDTD techniques were created to enhance the accuracy. Two approaches, PSTD and MRTD, were employed to minimize the numerical dispersion errors. Recently, the LOD-FDTD approach has been discussed as an unconditionally stable method, along with the CN-FDTD approach. The Crank-Nicolson method is computationally expensive but is more exact than the ADI-FDTD approach. Since Yee created the FDTD method, it has been successfully employed for a large number of problems. In addition, the CFL stability condition reduced the computational efficiency of FDTD, including for high Q structures in particular or structures where field variations occur rapidly. Therefore, the unconditionally stable CN-FDTD and ADI-FDTD approaches were proposed to address that restriction. Recently, the popularity of both techniques has increased because of unconditional stability. However, although ADI-FDTD is computationally effective, it has large errors when large time steps are used. This drawback of ADI-FDTD with unconditional stability compromises its utilization. However, high performance can be achieved using the CN-FDTD approach with large time steps, but it requires a larger computational period. Eventually, the benefits of both approaches are needed. In fact, from advance studies of both methods, it can be assumed that ADI-FDTD is an altered form of the CN-FDTD approach, which is known as the splitting error term. Accordingly, antagonistic techniques that can solve the Crank-Nicolson method in an alternating direction implicit platform have been illustrated. These techniques incorporate an iteration loop for each advance of Δt in the FDTD method. The ADI-FDTD algorithm [15] in the “Novel unconditionally stable FDTD method for electromagnetic and microwave modeling” developed by Zheng is numerically stable, even beyond the CFL point of the conventional FDTD. The original ADI technique has been applied to solve the parabolic differential equation and has proven to be stable. The ADI technique applied to the existing ADI-FDTD is unlike the original one. Rather than alternating the coordinate directions, from which ADI gets its name, the modified ADI technique used in ADI-FDTD is applied by alternating the order of the expressions on the right-hand side of the equations. Hence, the stability of the existing approach must be rigorously examined to demonstrate that the existing scheme is stable without any conditions. Preliminary numerical experiments show the stability of the existing ADI-FDTD even for time steps that are hundreds of times the threshold of the conventional FDTD. The conventional FDTD approach was developed by Yee [17]. The CFL stability condition limits the time step of the FDTD method, which is a marching-in-time method. Consequently, in some applications, the FDTD may require a large number of iterations when small time steps are needed, e.g., structures with fine geometries such as small via. The existing ADI-FDTD scheme does not require the CFL stability condition. The CFL condition does not restrict the time step but, rather, the accuracy of the modeling results and the Nyquist sampling limit. Hence, the number of iterations can be significantly decreased compared to the conventional FDTD when simulating the same structure. It will be shown that ADI-FDTD has an advantage in this case. Similar to any numerical method, numerical dispersion is inherent in the FDTD scheme. It is the result of the spatial discretization where the phase velocity of the numerical models in the FDTD is forced to change with the grid size, time step, wave propagation direction, and modal wavelength. These nonphysical results are produced by numerical dispersion, for example, pulse distortion and artificial anisotropy. To estimate the errors caused by numerical dispersion, the numerical dispersion relation of discrete lattices must be studied. For the conventional FDTD scheme, the numerical dispersion has been investigated extensively. In the homogenous waveguide case, the air-filled cavity has a dimension of 9 mm × 6 mm × 15 mm. For the ADI-FDTD and Yee’s FDTD approaches, a network with Δ L = 0.6 mm was used, resulting in 15 × 10 × 25 grid points. Note that that the time step Δt FDTD = 0.8 ps and the larger time step Δt i = Δt FDTD = 2.4 ps were chosen for FDTD and ADI-FDTD, respectively. The errors in both methods are comparable in Table 1, but the suggested ADI-FDTD yielded a 1.55× savings in CPU time.

In the inhomogeneous waveguide case, half of the cavity is occupied by a dielectric material with a relative permittivity of ε_r = 64, and the other half is occupied by air. The cavity has a dimension of l m × 2 m × 1.5 m. Again, a uniform mesh with ΔL = 0.1 m was used, which produced 10 × 20 × 15 mesh grid points, and a Δt FDTD = 1.5 × 10⁻¹ s was chosen with the typical FDTD, whereas a higher time step of Δt i = 4 × Δt FDTD = 6 × 10⁻¹⁰ s was selected with the suggested ADI-FDTD. To establish the strength of the ADI-FDTD approach, one homogeneous and one inhomogeneous optical waveguide cavities were determined using the ADI-FDTD and the typical Yee’s FDTD. The cavities were chosen because analytical solutions are readily available for comparison in Table 2.

It is very desirable to create unconditionally stable FDTD approaches that are free from the condition of CFL [18]. The celebrated alternating direction implicit (ADI) method is one of the methods based on FDTD. Recently, several scientists have suggested other unconditionally stable (non-dissipative) approaches such as split-step (SS) [19,20] and 1-D LOD-FDTD approaches [21,22]. LOD–FDTD techniques are used to establish 2D approaches related to the split-step process [19], which is designated as SS1. As normally contributed by scientists, LOD–FDTD produces an additional non-commutative error expression, consequently producing techniques precise to only the 1st stage in time.

For acquiring 2nd stage progressive accuracy, the Strang splitting technique can be implemented following the split-step technique [19], which is designated as SS2. Nevertheless, this technique requires extra arithmetic functions for executing three processes at every time period when calculating the fields, [20] and separate or parallel processes cannot be used in this case. It has been demonstrated that the substitute LOD–FDTD technique with the unconditionally stable property cannot produce the above-mentioned second stage commutative error. The three-dimensional Maxwell’s equation model can be presented as an LOD–FDTD. A numerical explanation will be given to establish its 2nd stage sequential accuracy. Moreover, the current three-dimensional LOD–FDTD approach includes two upgraded procedures that simply support the utilization of parallel and minimized outcome processing. The computational efficiency is higher than those that can be attained using the ADI-FDTD and SS2-FDTD approaches.

The LOD-FDTD approach is an additional example of a 2-D case; it is a two-stage technique, although its formulation is complex. It was observed that the previous approach was comparable to the split-step technique. Open structure problems can be resolved by convolutional perfectly matched layers (CPML), and perfectly matched layers (PML) have been effective for this approach, as shown by Iftikhar Ahmed et al. In the end, this research favors previous FDTD methods rather than the LOD-FDTD method, which includes three phases and is certainly well-defined by the common LOD method in the three spatial coordinate directions. The LOD approach can compute the equations, and the same number of equations can be computed using the typical ADI-FDTD approach, which requires more arithmetic operations. Table 3 compares both methods in term of arithmetic Operation Counts.

Therefore, LOD is more computationally efficient than the ADI-FDTD approach. The results were demonstrated based on the numerical calculations for the comparison of the analytical and typical ADI-FDTD approaches to demonstrate the usefulness of the existing approach. For validation purposes, the existing approach was employed on a 9 × 6 × 15 mm³ cavity. The selected cell sizes were Δx = Δy = Δz = 0.6 mm. The multiple modes had resonant frequencies that were computed and are shown in the following tables, and the analytical solutions and the solutions obtained with the typical ADI-FDTD approaches were compared in Table 4. It can be observed that ADI-FDTD and the current LOD-FDTD had similar relative errors.

In this case, there were eight thousand iterations. Furthermore, with respect to computational cost, the two approaches were dissimilar; the LOD approach required a simulation time of 20 s and 1.336 MB of memory, but the ADI-FDTD approach required 25 s and 1.336 MB of memory. Explicitly, both approaches used the same amount of memory; however, LOD-FDTD was almost 20% faster than the ADI-FDTD approach. The simulation was performed on a Dell Pentium VI computer with 256 MB of RAM, and the program was created in C++ by Ahmed et al. the results of that simulation in terms of time steps, CPU time, and memory size are illustrated in Table 5.