For the first time, an unconditionally stable finite-difference time-domain (FDTD) method for 3-D simulation of dispersive nonlinear media is presented. By applying a new adopted alternating-direction implicit (ADI) time-splitting scheme and the auxiliary differential equation (ADE) technique, the time-step in the FDTD simulations can be increased much beyond the Courant-Friedrichs-Lewy (CFL) stability limit. Thus, in comparison to the classical nonlinear FDTD method, the computational time for the proposed approach is decreased significantly while maintaining a reasonable level of accuracy. Numerical examples are presented to demonstrate the validity, stability, accuracy and computational efficiency of the proposed method.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
Electromagnetic (EM) wave interaction with media having frequency-dependent and intensity-dependent polarizations, namely dispersive nonlinear media, has gained strong interest in the fields of optics and photonics due to their applications in optical switches , four-wave mixing (FWM) , amplifiers , frequency converters , modulators, amongst other applications. The finite-difference time-domain (FDTD) [6, 7] which is a simple, explicit and computationally robust method for full-wave EM simulations has been also known as a powerful tool for simulation of EM wave in dispersive nonlinear media [8–11] and has been applied to numerous applications [12–25]. Commercial EM solvers based on this method such as Lumerical  and OptiFDTD  have been popularly utilized in applications using nonlinear optics and photonics.
Although the FDTD method has many attractive features, its main drawback is the constraint on the time step to ensure stability of the algorithm. To guarantee the stability of the method, the time step cannot exceed , where is the minimum spatial resolution in the computational domain, 1, 2, or 3 is the number of space dimensions, and c0 is the light speed in free space. This time step is known as the Courant-Friedrichs-Lewy (CFL) stabilitycriterion . For FDTD simulations in linear media, a spatial resolution of - , where is the minimum wavelength of interest, is commonly used . However, to achieve accurate results in nonlinear media, it has been shown that a much finer spatial resolution of less than is required . In most nonlinear FDTD simulations reported in the literature, a spatial resolution between and was used [19, 21, 29, 30]. Having such a fine spatial grid makes the simulation time prohibitively long and impractical for many design problems.
To eliminate the CFL restriction, implicit unconditionally stable FDTD algorithms such as the alternating-direction implicit (ADI) FDTD [31–33], Crank-Nicolson (CN) FDTD [34–36], locally one dimensional (LOD) FDTD [37, 38], and split-step (SS) FDTD  methods were developed. At every iteration in these approaches, the fields updating procedure is performed by solving some explicit and implicit equations. Of course this procedure involves more steps than in the classical explicit FDTD method; however, the time step can be set beyond the CFL criterion resulting in noticeable overall saving in computational time. The development of unconditionally stable FDTD algorithms for dispersive nonlinear media, to our best knowledge, was investigated in only two works in which a dispersive nonlinear medium having a third-order Kerr-Raman type nonlinearity was considered [29, 40]. In , by applying the ADI scheme and the Z-transform method [41, 42], an unconditionally stable nonlinear FDTD algorithm was derived for one-dimensional (1-D) and two-dimensional (2-D) problems. Although the application of Z-transform method simplifies the algorithm, it scarifies the second-order temporal accuracy of the FDTD method. In other words, since in the Z-transform method, the time-derivatives are approximated using a backward difference scheme, a first-order temporal-accurate algorithm is achievable [41, 42]. In , an unconditionally-stable nonlinear FDTD algorithm based on the application of the ADI scheme and the auxiliary differential equation (ADE) technique was proposed for 1-D problems. Unlike the Z-transform method, the ADE technique reported in  preserves the second-ordered temporal accuracy of the FDTD method. To our best knowledge, an unconditionally stable nonlinear FDTD method that is capable of solving 3-D electromagnetic problems has not been developed. In fact, applying the original ADI technique to a 3-D nonlinear problem, the nonlinear polarization current leads to the coupling of the electric field components. As a result, updating the field components requires solving three systems of strongly coupled nonlinear equations whose solution is computationally very heavy and inefficient.
In the present work, we introduce a new unconditionally stable FDTD algorithm for 3-D simulation of dispersive nonlinear media with a multiple-pole linear Lorentz polarization and a third-order nonlinearity combined of Kerr and Raman types. Our approach is based on the application of a new adopted ADI technique preserving the second-order temporal accuracy of the method. In the proposed method, the curl operators are time-split as in the original ADI method; however, a new time-splitting approach for the nonlinear polarization current is introduced. The resultant implicit updating equations can be effectively solved using the tridiagonal matrix algorithm. Finally the validity, stability, accuracy and computational efficiency of the proposed method is demonstrated through numerical examples.
2. Method development
2.1. Polarization currents
The Maxwell curl equations in a nonlinear medium can be written as:Eq. (2) can be rewritten in terms of polarization current as:
In materials with multiple-pole linear Lorentz polarization which models the chromatic dispersion and contains a contribution from multiple resonances, the linear polarization current can be written as,21], 44].
Transferring (6) into the time domain, one can obtain,
Discretizing (7) and centering at time-step give the following updating equation for JLp,
Besides, considering a third-order nonlinearity, the nonlinear polarization is given by,(11) is simplified to ,
In (13), is a Dirac delta function that models instantaneous Kerr nonresonant virtual transitions, denotes the Raman nonlinearity response and α is a real valued constant in the range of representing the relative strength of the Kerr and Raman interactions . Substituting (13) in (12), we get(15), the star sign stands for a time domain convolution.
Transferring (17) into the frequency domain yields44]
Applying a central difference scheme on (20) and enforcing the difference equation at time step n result in the following updating equation for SR,
Then, as in (4), the Raman polarization current can be defined as,
Discretizing (23) and centering at time-step gives the following updating equation for the Raman polarization current,
The Kerr polarization current is defined as,
Descritizing (25) and enforcing at time step , we have(16), the Kerr polarization at time step n is written by,
2.2. A novel ADI-FDTD scheme for third-order nonlinear media(30), JTx, JTy and JTz are the total polarization current along x, y and z directions, respectively, which can be expressed as follows: (8),(24) and (28) in (32), we have
By applying the ADI time-splitting scheme  (in which the advance from time step n to time step is performed in two sub-iterations: first, and then ) on (29) and (30) and a new proposed time-splitting scheme on the nonlinear polarization current (JNL in (30)), the following updating equations are derived for the first sub-iteration (i.e.
Recalling from (33), and can be defined by(36) and (40), respectively) in (37), we have the following updating equations for the electric field components at time step ,
Notice that in the left hand side of (41) we have where a is a constant and is the second order central finite-difference operator along , y and z directions, for instance,
Therefore, updating the electric field components at through (41) requires solving a tridiagonal system of linear equations which can be achieved using the tridiagonal matrix algorithm (Thomas algorithm) .(42)) are nonlinear and implicit (notice that for solving each implicit equation a solution of a system of equations is required), we have three systems of nonlinear equations strongly coupled to each other whose solution is computationally very heavy. To remove this problem, we use the following approximation having a second order temporal accuracy
Therefore, when solving (42), is a known constant since the field components at earlier time steps (i.e., , n) are known. Having as a known constant significantly simplifies solving the updating equations of Ex, Ey and Ez (i.e., (42)) since:
- The updating equations of Ex, Ey and Ez (i.e., (42)) are not longer nonlinear.
- The coupling between the updating equations for Ex, Ey and Ez is removed.
Finally, the field updating procedure in each iteration of the algorithm is summarized in the following seven steps (sub-iterations):
3. Numerical example and discussion
Two numerical examples are provided to illustrate the accuracy, efficiency, and applicability of the proposed method. As a first example, third harmonic generation in plane wave transmission through a nonlinear slab is simulated. Then, as a second example, a nonlinear silicon waveguide is investigated. In both examples, the stability of the proposed method when the time step is much beyond the CFL limit is demonstrated. It should be noted that analytically proving the unconditional stability using the von Neumann Method  is very challenging due to the nonlinearity of the medium [29. 40]. The results obtained from the proposed method are validated by providing a comparison with those obtained from the explicit nonlinear FDTD method addressed in .
3.1. First example
As a first example, plane wave transmission through a highly-nonlinear infinite slab (having a very large value of nonlinear susceptibility, ) was considered. The characteristic parameters of the slab were the same as in ,which are , , fs, fs, and
The slab was illuminated by an x-polorized plane wave having two carrier frequencies modulated with a Gaussian waveform,Fig. 1, the slab was placed in a TEM waveguide made of the perfect electric conductor (PEC) and the perfect magnetic conductor (PMC) walls  in a way that the incident electric field is perpendicular to the PEC walls and parallel to the PMC walls. To this end, as shown in Fig. 1, the computational domain was discretized into cells with a size of nm and terminated by PEC, PMC and Mur’s first-order absorbing boundary condition (ABC) in the x, y, and z directions, respectively .
For comparison, the problem was simulated using the explicit FDTD method reported in  (as a reference) and using the proposed unconditionally-stable FDTD method. The CFL condition limits the time step for the explicit FDTD method to fs. In the proposed unconditionally-stable FDTD method, the time step can be set beyond the CFL limit at the expense of sacrificing a little accuracy. In our simulations, the time step was set to , , and . The simulations were run for a large number of time steps (1,400,000 ). No instability was observed. A comparison between the run-times of the simulations are provided in TABLE 1. As is seen in the table, by increasing the time step to 4, 8, and 16 times of the CFL limit, the run time is respectively reduced to 52%, 22% and 13% of that of the explicit FDTD method.
The transmitted Ex was recorded which is shown in Fig. 2 (for the explicit FDTD method and the proposed method with ). A very good agreement between the results is obvious. Then by applying the discrete Fourier transform (DFT) to the time domain recorded data, the power intensity of the transmitted wave was calculated. Figure 3 shows the normalized power intensity obtained using the explicit FDTD method and the proposed method with , , and . In addition to the peak power intensity at the carrier frequencies, THz and THz, Fig. 3 shows two peaks at 200 THz and 500 THz corresponding to the third harmonic frequencies of and due to the third order nonlinearity of the medium which enables the application of a nonlinear slab as an optical mixer. Very good agreement between the results of the proposed unconditionally-stable FDTD method (for , , and ) and the explicit FDTD method is observed in Fig. 3; however, by setting , they diverge. Notice that although in Fig. 2, the results of the proposed method with is shown, the proposed method with smaller time step (i.e., and ) is in a better agreement with the explicit FDTD method.
3.2. Second example
As a second example, we simulated a nonlinear silicon waveguide. The characteristic parameters of the silicon waveguide are considered the same as in , i.e., , , , fs, and fs. As shown in Fig. 4, the nonlinear silicon medium is surrounded by dielectric sections having a relative permittivity of and the free space (the dielectric section is utilized to excite the waveguide). The computational domain consists of 2000 × 140 cells (along x and y directions) where the spatial resolution is set to nm and is terminated by the ABC.
First, to demonstrate the accuracy of the proposed method, a cosine-modulated Gaussian source having a waveform of19] (as a reference) with fs and the proposed method with , and the values of Ex were recorded at grid (1750, 70) (see Fig. 4). The recorded Ex is shown in Fig. 5 indicating a strong agreements between the results of the proposed and the reference methods.
Next, to show the field distribution over the silicon waveguide, the magnetic field was excited by a hard source at the far-left side of the grid, at x = 10 (see Fig. 4) having an initial transverse profile ofFig. 6 depicts the magnitude of the Hz over the computational domain at the steady state. Figure 6 clearly shows a guiding wave in the waveguide, however, due to the nonlinear characteristics of the silicon waveguide, some expansion and contraction of the filed magnitude is observed.
3.3. Third example
Finally, as a third example, we use the proposed method to simulate scattering of spatial optical solitons by subwavelength air holes. We assume a realistic glass background material characterized by a three-pole Sellmeier linear dispersion, an instantaneous Kerr nonlinearity, and a dispersive Raman nonlinearity. The characteristic parameters of the material were given from [21, 49] as the following. The strengths and resonant frequencies of the linear dispersions from Sellmeier’s equation were assumed to be , , , , , . And the Kerr and Raman nonlinear characteristic parameters were considered as , , , fs, and fs.
A computational domain consisting of 1500 × 600 cells (along x and y directions, respectively) where the spatial resolution is set to nm was considered. The computational domain was terminated by a first-order ABC. The magnetic field was excited by a hard source at the far-left side of the grid, at x = 10 having an initial transverse profile of given in (53), where A/m, THz ( nm) and nm. An air hole with a size of 250 nm × 250 nm was assumed in a way that its center is located at grid and . The time step size was set to 10 times beyond the CFL stability criteria and the simulations was run for a large number of time steps until a steady state was achieved.
Fig. 7 depicts the magnitude of the electric field () over the computational domain at the steady state. The figure shows propagation and scattering of spatial optical solitons which is in consistent with the results presented in .
For the first time, a 3-D unconditionally-stable FDTD algorithm based on applying the ADI scheme with eliminated CFL restriction for numerical electromagnetic analysis of nonlinear optical media was introduced. First the ADI scheme was applied on the Maxwell curl equations, then the polarization current was split into two terms in a way that they could be suitably updated in two sub-iterations of the ADI time-splitting scheme. The updating equations resulted in three systems of nonlinear equations strongly coupled to each other whose solution is computationally very heavy. Then using a second-order accurate approximation, not only the coupling between the updating equations was removed but also tridiagonal systems of linear equations were attained; hence, efficient update of each electric field component could be achieved by simply solving a tridiagonal system of linear equations. The stability of the proposed method was tested for different values of time steps appreciably beyond the CFL stability criteria in a highly nonlinear medium, without observing any instability. The accuracy and computational efficiency of the proposed method were investigated through numerical examples. It was shown that by increasing the time step up to 8-10 times beyond the CFL stability condition, the proposed algorithm becomes much faster than the classical FDTD method while keeping the numerical errors within a reasonable bound.
2. H. Fukuda, K. Yamada, T. Shoji, M. Takahashi, T. Tsuchizawa, T. Watanabe, J. Takahashi, and S. Itabashi, “Four-wave mixing in silicon wire waveguides,” Opt. Express 13, 4629–4637 (2005). [CrossRef] [PubMed]
4. V. Raghunathan, R. Claps, D. Dimitropoulos, and B. Jalali, “Parametric Raman wavelength conversion in scaled silicon waveguides,” J. Light. Technol. 23, 2094 (2005). [CrossRef]
5. C. Manolatou and M. Lipson, “All-optical silicon modulators based on carrier injection by two-photon absorption,” J. Light. Technol. 24, 1433 (2006). [CrossRef]
6. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Transactions on Antennas Propag. 14, 302–307 (1966). [CrossRef]
7. A. Taflove and S. C. Hagness, Computational Electrodynamics: the Finite-Difference Time-Domain Method(Artech House, 2005).
8. P. M. Goorjian and A. Taflove, “Direct time integration of Maxwell’s equations in nonlinear dispersive media for propagation and scattering of femtosecond electromagnetic solitons,” Opt. Lett. 17, 180–182 (1992). [CrossRef]
9. R. W. Ziolkowski and J. B. Judkins, “Full-wave vector Maxwell equation modeling of the self-focusing of ultrashort optical pulses in a nonlinear Kerr medium exhibiting a finite response time,” JOSA B 10, 186–198 (1993). [CrossRef]
10. R. W. Ziolkowski and J. B. Judkins, “Applications of the nonlinear finite difference time domain (NL-FDTD) method to pulse propagation in nonlinear media: Self-focusing and linear-nonlinear interfaces,” Radio Sci. 28,901–911 (1993). [CrossRef]
11. R. W. Ziolkowski and J. B. Judkins, “Nonlinear finite-difference time-domain modeling of linear and nonlinear corrugated waveguides,” JOSA B 11, 1565–1575 (1994). [CrossRef]
12. D. M. Sullivan, “Nonlinear FDTD formulations using Z-transforms,” IEEE Transactions on Microw. Theory Tech. 43, 676–682 (1995). [CrossRef]
13. C. V. Hile and W. L. Kath, “Numerical solutions of Maxwell’s equations for nonlinear-optical pulse propagation,” JOSA B 13, 1135–1145 (1996). [CrossRef]
14. R. M. Joseph and A. Taflove, “FDTD Maxwell’s equations models for nonlinear electrodynamics and optics,” IEEE Transactions on Antennas Propag. 45, 364–374 (1997). [CrossRef]
15. P. M. Goorjian and Y. Silberberg, “Numerical simulations of light bullets using the full-vector time-dependent nonlinear Maxwell equations,” JOSA B 14, 3253–3260 (1997). [CrossRef]
16. A. S. Nagra and R. A. York, “FDTD analysis of wave propagation in nonlinear absorbing and gain media,” IEEE Transactions on Antennas Propag. 46, 334–340 (1998). [CrossRef]
17. D. Sullivan, J. Liu, and M. Kuzyk, “Three-dimensional optical pulse simulation using the FDTD method,” IEEE Transactions on Microw. Theory Tech. 48, 1127–1133 (2000). [CrossRef]
18. S. Nakamura, Y. Koyamada, N. Yoshida, N. Karasawa, H. Sone, M. Ohtani, Y. Mizuta, R. Morita, H. Shigekawa, and M. Yamashita, “Finite-difference time-domain calculation with all parameters of Sellmeier’s fitting equation for 12-fs laser pulse propagation in a silica fiber,” IEEE Photonics Technol. Lett. 14, 480–482 (2002). [CrossRef]
19. M. Fujii, M. Tahara, I. Sakagami, W. Freude, and P. Russer, “High-order FDTD and auxiliary differential equation formulation of optical pulse propagation in 2-D Kerr and Raman nonlinear dispersive media,” IEEE J. Quantum Electron. 40, 175–182 (2004). [CrossRef]
20. S. Nakamura, N. Takasawa, and Y. Koyamada, “Comparison between finite-difference time-domain calculation with all parameters of Sellmeier’s fitting equation and experimental results for slightly chirped 12-fs laser pulse propagation in a silica fiber,” J. Light. Technol. 23, 855 (2005). [CrossRef]
22. C. M. Dissanayake, M. Premaratne, I. D. Rukhlenko, and G. P. Agrawal, “FDTD modeling of anisotropic nonlinear optical phenomena in silicon waveguides,” Opt. Express 18, 21427–21448 (2010). [CrossRef] [PubMed]
23. C. M. Reinke, A. Jafarpour, B. Momeni, M. Soltani, S. Khorasani, A. Adibi, Y. Xu, and R. K. Lee, “Nonlinear finite-difference time-domain method for the simulation of anisotropic, , and optical effects,” J. Light. Technol. 24, 624–634 (2006). [CrossRef]
24. M. A. Alsunaidi, H. M. Al-Mudhaffar, and H. M. Masoudi, “Vectorial FDTD technique for the analysis of optical second-harmonic generation,” IEEE Photonics Technol. Lett. 21, 310–312 (2009). [CrossRef]
25. B. Salski, T. Karpisz, and R. Buczynski, “Electromagnetic modeling of third-order nonlinearities in photonic crystal fibers using a vector two-dimensional FDTD algorithm,” J. Light. Technol. 33, 2905–2912 (2015). [CrossRef]
26. I. Lumerical Solutions Inc., “https://www.lumerical.com,” (2018).
27. I. Optiwave Systems Inc., “https://optiwave.com/optifdtd-overview,” (2018).
28. R. Courant, K. Friedrichs, and H. Lewy, “Über die partiellen differenzengleichungen der mathematischen physik,” Math. Annalen 100, 32–74 (1928). [CrossRef]
29. D. Li and C. D. Sarris, “Time-domain modeling of nonlinear optical structures with extended stability FDTD schemes,” J. Light. Technol. 29, 1003–1010 (2011). [CrossRef]
30. I. S. Maksymov, A. A. Sukhorukov, A. V. Lavrinenko, and Y. S. Kivshar, “Comparative study of FDTD-adopted numerical algorithms for Kerr nonlinearities,” IEEE Antennas Wirel. Propag. Lett. 10, 143–146 (2011). [CrossRef]
31. F. Zheng, Z. Chen, and J. Zhang, “A finite-difference time-domain method without the Courant stability conditions,” IEEE Microw. Guid. Wave Lett. 9, 441–443 (1999). [CrossRef]
32. T. Namiki, “A new FDTD algorithm based on alternating-direction implicit method,” IEEE Transactions on Microw. Theory Tech. 47, 2003–2007 (1999). [CrossRef]
33. F. Zhen, Z. Chen, and J. Zhang, “Toward the development of a three-dimensional unconditionally stable finite-difference time-domain method,” IEEE Transactions on Microw. Theory Tech. 48, 1550–1558 (2000). [CrossRef]
34. G. Sun and C. W. Trueman, “Efficient implementations of the Crank-Nicolson scheme for the finite-difference time-domain method,” IEEE Transactions on Microw. Theory Tech. 54, 2275–2284 (2006). [CrossRef]
35. S.-M. Sadrpour, V. Nayyeri, and M. Soleimani, “A new 2d unconditionally stable finite-difference time-domain algorithm based on the crank-nicolson scheme,” in 2016 IEEE International Conference on Computational Electromagnetics (ICCEM), (IEEE, 2016), pp. 55–57.
36. S.-M. Sadrpour, V. Nayyeri, M. Soleimani, and O. M. Ramahi, “A new efficient unconditionally stable finite-difference time-domain solution of the wave equation,” IEEE Transactions on Antennas Propag. 65, 3114–3121 (2017). [CrossRef]
37. J. Shibayama, M. Muraki, J. Yamauchi, and H. Nakano, “Efficient implicit FDTD algorithm based on locally one-dimensional scheme,” Electron. Lett. 41, 1046–1047 (2005). [CrossRef]
38. E. L. Tan, “Unconditionally stable LOD–FDTD method for 3-D Maxwell’s equations,” IEEE Microw. Wirel. Components Lett. 17, 85–87 (2007). [CrossRef]
39. J. Lee and B. Fornberg, “A split step approach for the 3-D Maxwell’s equations,” J. Comput. Appl. Math. 158, 485–505 (2003). [CrossRef]
40. E. P. Kosmidou and T. D. Tsiboukis, “An unconditionally stable ADI-FDTD algorithm for nonlinear materials,” Opt. Quantum Electron 32, 931–946 (2003). [CrossRef]
41. D. M. Sullivan, “Z-transform theory and the FDTD method,” IEEE Transactions on Antennas Propag. 44, 28–34 (1996). [CrossRef]
42. V. Nayyeri, M. Soleimani, J. R. Mohassel, and M. Dehmollaian, “FDTD modeling of dispersive bianisotropic media using Z-transform method,” IEEE Transactions on Antennas Propag. 59, 2268–2279 (2011). [CrossRef]
43. R. M. Joseph, S. C. Hagness, and A. Taflove, “Direct time integration of Maxwell’s equations in linear dispersive media with absorption for scattering and propagation of femtosecond electromagnetic pulses,” Opt. Lett. 16, 1412–1414 (1991). [CrossRef] [PubMed]
44. G. P. Agrawal, Nonlinear Fiber Optics (Academic Press, 2001).
45. S. D. Conte and C. De Boor, Elementary Numerical Analysis: an Algorithmic Approach (SIAM, 2017), vol. 78, pp. 153–156.
46. J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, vol. 88 (Siam, 2004).
47. O. Corporation, OptiFDTD Technical Background and Tutorials (Optiwave Systems Inc., ON, Canada, 2005).
48. N. K. Hon, R. Soref, and B. Jalali, “The third-order nonlinear optical coefficients of Si, Ge, and in the midwave and longwave infrared,” J. Appl. Phys. 110, 9 (2011). [CrossRef]
49. J. H. Greene and A. Taflove, “Scattering of spatial optical solitons by subwavelength air holes,” IEEE Microw. Wirel. Components Lett. 17, 760–762 (2007). [CrossRef]