Unconditionally stable FDTD algorithm for 3-D electromagnetic simulation of nonlinear media

Mohammad Moradi; Seyyed-Mohammad Pourangha; Vahid Nayyeri; Mohammad Soleimani; Omar M. Ramahi

doi:10.1364/OE.27.015018

1. Introduction

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 [1], four-wave mixing (FWM) [2], amplifiers [3], frequency converters [4], modulators[5], 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 [26] and OptiFDTD [27] 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 $Δ_{m i n} / (\sqrt{N_{d i m}} c_{0})$ , where $Δ_{m i n}$ is the minimum spatial resolution in the computational domain, $N_{d i m} =$ 1, 2, or 3 is the number of space dimensions, and c₀ is the light speed in free space. This time step is known as the Courant-Friedrichs-Lewy (CFL) stabilitycriterion [28]. For FDTD simulations in linear media, a spatial resolution of $λ_{m i n} / 20$ - $λ_{m i n} / 10$ , where $λ_{m i n}$ is the minimum wavelength of interest, is commonly used [7]. However, to achieve accurate results in nonlinear media, it has been shown that a much finer spatial resolution of less than $λ / 100$ is required [10]. In most nonlinear FDTD simulations reported in the literature, a spatial resolution between $λ / 300$ and $λ / 100$ 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 [39] 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 [40], 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 [29], 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 [43] 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:

\nabla \times E = - μ_{0} \frac{\partial H}{\partial t},

\nabla \times H = \frac{\partial}{\partial t} [ε_{0} ε_{r} E + P_{T}],

where E and H are the electric and magnetic fields and

P_{T}

denotes the total polarization vector which is the summation of the linear (

P_{L}

) and nonlinear (P_NL) polarization vectors, i.e.

P_{T} = P_{L} + P_{NL}

. Eq. (2) can be rewritten in terms of polarization current as:

\nabla \times H = ε_{0} ε_{r} \frac{\partial E}{\partial t} + J_{T},

where

J_{T} = J_{L} + J_{NL} = \frac{\partial P_{T}}{\partial t} = \frac{\partial P_{L}}{\partial t} + \frac{\partial P_{NL}}{\partial t}

is the total polarization current vector and

J_{L}

and J_NL are the linear and nonlinear polarization current vectors respectively.

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,

J_{L} = \sum_{p = 1}^{m} J_{Lp},

where each J_Lp denotes the polarization current due to a single pole of the Sellmeier equation whose expression in the frequency domain is given by [21],

{\tilde{J}}_{Lp} = ε_{0} β_{p} ω_{p}^{2} (\frac{j ω}{ω_{p}^{2} - ω^{2}}) \tilde{E},

and β_p and ω_p are the strength and frequency of the pth resonance, respectively [44].

Transferring (6) into the time domain, one can obtain,

ω_{p}^{2} J_{Lp} + \frac{\partial^{2} J_{Lp}}{\partial t^{2}} = ε_{0} β_{p} ω_{p}^{2} \frac{\partial E}{\partial t} .

Discretizing (7) and centering at time-step $n + \frac{1}{2}$ give the following updating equation for J_Lp,

J_{Lp}^{n + \frac{1}{2}} = α_{p} J_{Lp}^{n - \frac{1}{2}} - J_{Lp}^{n - \frac{3}{2}} + γ_{p} (E^{n} - E^{n - 1})

where

\begin{array}{l} α_{p} = 2 - {(ω_{p} Δ t)}^{2}, \\ γ_{p} = ε_{0} β_{p} ω_{p}^{2} Δ t . \end{array}

Then, by substituting (8) in (5), an updating equation for $J_{L}$ is obtained as,

J_{L}^{n + \frac{1}{2}} = \sum_{p = 1}^{m} {α_{p} J_{Lp}^{n - \frac{1}{2}} - J_{Lp}^{n - \frac{3}{2}} + γ_{p} (E^{n} - E^{n - 1})} .

Besides, considering a third-order nonlinearity, the nonlinear polarization is given by,

P_{NL} (r, t) = ε_{0} \int_{- \infty}^{t} \int_{- \infty}^{t} \int_{- \infty}^{t} χ^{(3)} (t - τ_{1}, t - τ_{2}, t - τ_{3}) E (r, τ_{1}) E (r, τ_{2}) E (r, τ_{3}) d τ_{1} d τ_{2} d τ_{3},

where

χ^{(3)}

is the third-order susceptibility tensor. Using the Born-Oppenheimer approximation, (11) is simplified to [44],

P_{NL} (t) = ε_{0} χ_{0}^{(3)} E (t) \int_{- \infty}^{t} g (t - τ) {| E (τ) |}^{2} d τ

where

χ_{0}^{(3)}

is a constant value representing the third-order nonlinear susceptibility and

g (t)

is the nonlinearity causal response function which can be written in terms of Kerr and Raman nonlinearity responses as,

g (t) = (1 - α) g_{Raman} (t) + α δ (t),

In (13), $δ (t)$ is a Dirac delta function that models instantaneous Kerr nonresonant virtual transitions, $g_{Raman} (t)$ denotes the Raman nonlinearity response and α is a real valued constant in the range of $0 \leq α \leq 1$ representing the relative strength of the Kerr and Raman interactions [21]. Substituting (13) in (12), we get

P_{NL} = P_{R} + P_{K}

where the Raman and Kerr polarizations are given by

P_{R} (t) = ε_{0} (1 - α) [χ_{0}^{(3)} g_{Raman} (t) * {| E (t) |}^{2}] E (t),

P_{K} (t) = ε_{0} α χ_{0}^{(3)} {| E (t) |}^{2} E (t),

respectively. In (15), the star sign stands for a time domain convolution.

For the Raman polarization, an auxiliary variable $(S_{R})$ can be defined as [21, 29]

S_{R} (t) = (1 - α) χ_{0}^{(3)} g_{Raman} (t) * {| E (t) |}^{2} .

Transferring (17) into the frequency domain yields

\tilde{S_{R}} = (1 - α) χ_{0}^{(3)} {\tilde{g}}_{Raman} (ω) {| \tilde{E} |}^{2},

where

{\tilde{g}}_{Raman}

is the Fourier transform of g_Raman, given by [44]

{\tilde{g}}_{Raman} (ω) = \frac{ω_{R}^{2}}{ω_{R}^{2} + 2 j ω δ_{R} - ω^{2}} .

Substituting (19) in (18) and applying the inverse Fourier-transform, the following differential equation for the auxiliary variable S_R is obtained:

(ω_{R}^{2} + 2 δ_{R} \frac{\partial}{\partial t} + \frac{\partial^{2}}{\partial t^{2}}) S_{R} = (1 - α) χ_{0}^{(3)} ω_{R}^{2} {| E |}^{2} .

Applying a central difference scheme on (20) and enforcing the difference equation at time step n result in the following updating equation for S_R,

S_{R}^{n + 1} = α_{R} S_{R}^{n} + β_{R} S_{R}^{n - 1} + γ_{R} {| E^{n} |}^{2},

where

\begin{array}{l} α_{R} = \frac{2 - ω_{R}^{2} Δ t^{2}}{δ_{R} Δ t + 1}, \\ β_{R} = \frac{δ_{R} Δ t - 1}{δ_{R} Δ t + 1}, \\ γ_{R} = \frac{(1 - α) χ_{0}^{(3)} ω_{R}^{2} Δ t^{2}}{δ_{R} Δ t + 1} . \end{array}

Then, as in (4), the Raman polarization current can be defined as,

J_{R} = \frac{\partial P_{R}}{\partial t} = ε_{0} \frac{\partial}{\partial t} (S_{R} E) .

Discretizing (23) and centering at time-step $n + \frac{1}{2}$ gives the following updating equation for the Raman polarization current,

J_{R}^{n + \frac{1}{2}} = \frac{ε_{0}}{Δ t} (S_{R}^{n + 1} E^{n + 1} - S_{R}^{n} E^{n}) .

The Kerr polarization current is defined as,

J_{K} = \frac{\partial P_{K}}{\partial t} .

Descritizing (25) and enforcing at time step $n + \frac{1}{2}$ , we have

J_{K}^{n + \frac{1}{2}} = \frac{P_{K}^{n + 1} - P_{K}^{n}}{Δ t},

where, according to (16), the Kerr polarization at time step n is written by,

P_{K}^{n} = α ε_{0} χ_{0}^{(3)} {| E^{n} |}^{2} E^{n} .

Substituting (27) in (26), the following updating equation for $J_{K}$ is obtained:

J_{K}^{n + \frac{1}{2}} = \frac{α ε_{0} χ_{0}^{(3)}}{Δ t} ({| E^{n + 1} |}^{2} E^{n + 1} - {| E^{n} |}^{2} E^{n}) .

2.2. A novel ADI-FDTD scheme for third-order nonlinear media

Discretizing the Maxwell’s curl equations (1) and (2) in time and space and enforcing them at time step $n + \frac{1}{2}$ , we have:

\begin{array}{l} H_{x}^{n + 1} - H_{x}^{n} = \frac{Δ t}{μ} D_{z} E_{y}^{n + \frac{1}{2}} - \frac{Δ t}{μ} D_{y} E_{z}^{n + \frac{1}{2}}, \\ H_{y}^{n + 1} - H_{y}^{n} = \frac{Δ t}{μ} D_{x} E_{z}^{n + \frac{1}{2}} - \frac{Δ t}{μ} D_{z} E_{x}^{n + \frac{1}{2}}, \\ H_{z}^{n + 1} - H_{z}^{n} = \frac{Δ t}{μ} D_{y} E_{x}^{n + \frac{1}{2}} - \frac{Δ t}{μ} D_{x} E_{y}^{n + \frac{1}{2}}, \end{array}

\begin{array}{l} E_{x}^{n + 1} - E_{x}^{n} = \frac{Δ t}{ε_{0} ε_{\infty}} D_{y} H_{z}^{n + \frac{1}{2}} - \frac{Δ t}{ε_{0} ε_{\infty}} D_{z} H_{y}^{n + \frac{1}{2}} - \frac{Δ t}{ε_{0} ε_{\infty}} J_{T}^{n + \frac{1}{2}}, \\ E_{y}^{n + 1} - E_{y}^{n} = \frac{Δ t}{ε_{0} ε_{\infty}} D_{z} H_{x}^{n + \frac{1}{2}} - \frac{Δ t}{ε_{0} ε_{\infty}} D_{x} H_{z}^{n + \frac{1}{2}} - \frac{Δ t}{ε_{0} ε_{\infty}} J_{T}^{n + \frac{1}{2}}, \\ E_{z}^{n + 1} - E_{z}^{n} = \frac{Δ t}{ε_{0} ε_{\infty}} D_{x} H_{y}^{n + \frac{1}{2}} - \frac{Δ t}{ε_{0} ε_{\infty}} D_{y} H_{x}^{n + \frac{1}{2}} - \frac{Δ t}{ε_{0} ε_{\infty}} J_{T}^{n + \frac{1}{2}}, \end{array}

where

D_{x}

,

D_{y}

and

D_{z}

are the first-order central finite-difference operators along the x, y and z directions, for example,

D_{x} f_{i, j, k} = \frac{1}{Δ x} (f_{i + \frac{1}{2}, j, k} - f_{i - \frac{1}{2}, j, k})

and

f_{i, j, k}

denotes

f_{(i Δ x, j Δ y, k Δ z)}

. In (30), J_Tx, J_Ty and J_Tz are the total polarization current along x, y and z directions, respectively, which can be expressed as follows:

J_{T ρ}^{n + \frac{1}{2}} = J_{N L ρ}^{n + \frac{1}{2}} + J_{L ρ}^{n + \frac{1}{2}} = J_{R ρ}^{n + \frac{1}{2}} + J_{K ρ}^{n + \frac{1}{2}} + J_{L ρ}^{n + \frac{1}{2}},

where

ρ = x

, y, and z. Substituting (8),(24) and (28) in (32), we have

J_{T ρ}^{n + \frac{1}{2}} = \frac{ε_{0}}{Δ t} S_{R}^{n + 1} E_{ρ}^{n + 1} - \frac{ε_{0}}{Δ t} S_{R}^{n} E_{ρ}^{n} + \frac{α ε_{0} χ_{0}^{(3)}}{Δ t} {| E^{n + 1} |}^{2} E_{ρ}^{n + 1} - \frac{α ε_{0} χ_{0}^{(3)}}{Δ t} {| E^{n} |}^{2} E_{ρ}^{n} + J_{L ρ}^{n + \frac{1}{2}} .

By applying the ADI time-splitting scheme [33] (in which the advance from time step n to time step $n + 1$ is performed in two sub-iterations: first, $n \to n + \frac{1}{2}$ and then $n + \frac{1}{2} \to n + 1$ ) on (29) and (30) and a new proposed time-splitting scheme on the nonlinear polarization current (J_NL in (30)), the following updating equations are derived for the first sub-iteration (i.e.

\begin{array}{l} H_{x}^{n + \frac{1}{2}} - H_{x}^{n} = \frac{Δ t}{2 μ} D_{z} E_{y}^{n + \frac{1}{2}} - \frac{Δ t}{2 μ} D_{y} E_{z}^{n}, \\ H_{y}^{n + \frac{1}{2}} - H_{y}^{n} = \frac{Δ t}{2 μ} D_{x} E_{z}^{n + \frac{1}{2}} - \frac{Δ t}{2 μ} D_{z} E_{x}^{n}, \\ H_{z}^{n + \frac{1}{2}} - H_{z}^{n} = \frac{Δ t}{2 μ} D_{y} E_{x}^{n + \frac{1}{2}} - \frac{Δ t}{2 μ} D_{x} E_{y}^{n}, \end{array}

\begin{array}{l} E_{x}^{n + \frac{1}{2}} - E_{x}^{n} = \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{y} H_{z}^{n + \frac{1}{2}} - \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{z} H_{y}^{n} - \frac{Δ t}{ε_{0} ε_{\infty}} J_{1 x}^{n}, \\ E_{y}^{n + \frac{1}{2}} - E_{y}^{n} = \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{z} H_{x}^{n + \frac{1}{2}} - \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{x} H_{z}^{n} - \frac{Δ t}{ε_{0} ε_{\infty}} J_{1 y}^{n}, \\ E_{z}^{n + \frac{1}{2}} - E_{z}^{n} = \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{x} H_{y}^{n + \frac{1}{2}} - \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{y} H_{x}^{n} - \frac{Δ t}{ε_{0} ε_{\infty}} J_{1 z}^{n}, \end{array}

and the bellow equations for the second sub-iteration (i.e.

n + \frac{1}{2} \to n + 1

),

\begin{array}{l} H_{x}^{n + 1} - H_{x}^{n + \frac{1}{2}} = \frac{Δ t}{2 μ} D_{z} E_{y}^{n + \frac{1}{2}} - \frac{Δ t}{2 μ} D_{y} E_{z}^{n + 1}, \\ H_{y}^{n + 1} - H_{y}^{n + \frac{1}{2}} = \frac{Δ t}{2 μ} D_{x} E_{z}^{n + \frac{1}{2}} - \frac{Δ t}{2 μ} D_{z} E_{x}^{n + 1}, \\ H_{z}^{n + 1} - H_{z}^{n + \frac{1}{2}} = \frac{Δ t}{2 μ} D_{y} E_{x}^{n + \frac{1}{2}} - \frac{Δ t}{2 μ} D_{x} E_{y}^{n + 1}, \end{array}

\begin{array}{l} E_{x}^{n + 1} - E_{x}^{n + \frac{1}{2}} = \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{y} H_{z}^{n + \frac{1}{2}} - \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{z} H_{y}^{n + 1} - \frac{Δ t}{ε_{0} ε_{\infty}} J_{2 x}^{n + 1}, \\ E_{y}^{n + 1} - E_{y}^{n + \frac{1}{2}} = \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{z} H_{x}^{n + \frac{1}{2}} - \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{x} H_{z}^{n + 1} - \frac{Δ t}{ε_{0} ε_{\infty}} J_{2 y}^{n + 1}, \\ E_{z}^{n + 1} - E_{z}^{n + \frac{1}{2}} = \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{x} H_{y}^{n + \frac{1}{2}} - \frac{Δ t}{2 ε_{0} ε_{\infty}} D_{y} H_{x}^{n + 1} - \frac{Δ t}{ε_{0} ε_{\infty}} J_{2 z}^{n + 1} . \end{array}

Notice that the summation of (34) (and (35)) and (36) (and (37)) results in (29) (and (30)), thus

J_{T ρ}^{n + \frac{1}{2}} = J_{1 ρ}^{n} + J_{2 ρ}^{n + 1}, (ρ = x, y, z) .

Recalling $J_{T ρ}^{n + \frac{1}{2}}$ from (33), $J_{1 ρ}^{n}$ and $J_{2 ρ}^{n + 1}$ can be defined by

J_{1 ρ}^{n} = - \frac{ε_{0}}{Δ t} S_{R}^{n} E_{ρ}^{n} - \frac{α ε_{0} χ_{0}^{(3)}}{Δ t} {| E^{n} |}^{2} E_{ρ}^{n} + J_{L ρ}^{n + \frac{1}{2}},

J_{2 ρ}^{n + 1} = \frac{ε_{0}}{Δ t} S_{R}^{n + 1} E_{ρ}^{n + 1} + \frac{α ε_{0} χ_{0}^{(3)}}{Δ t} {| E^{n + 1} |}^{2} E_{ρ}^{n + 1} .

In (35), by substituting $H_{ρ}^{n + \frac{1}{2}}$ and $J_{1 ρ}^{n}$ ( $ρ = x, y, z$ ) from (34) and (39), respectively, the updating equations for the electric field components at time step $n + \frac{1}{2}$ are obtained as,

\begin{array}{l} (1 - a D_{y}^{2}) E_{x}^{n + \frac{1}{2}} = A^{n} E_{x}^{n} - a D_{x} D_{y} E_{y}^{n} + b D_{y} H_{z}^{n} - b D_{z} H_{y}^{n} - 2 b J_{L x}^{n + \frac{1}{2}}, \\ (1 - a D_{z}^{2}) E_{y}^{n + \frac{1}{2}} = A^{n} E_{y}^{n} - a D_{y} D_{z} E_{z}^{n} + b D_{z} H_{x}^{n} - b D_{x} H_{z}^{n} - 2 b J_{L y}^{n + \frac{1}{2}}, \\ (1 - a D_{x}^{2}) E_{z}^{n + \frac{1}{2}} = A^{n} E_{z}^{n} - a D_{z} D_{x} E_{x}^{n} + b D_{x} H_{y}^{n} - b D_{y} H_{x}^{n} - 2 b J_{L z}^{n + \frac{1}{2}}, \end{array}

and in a similar way, by substituting

H_{ρ}^{n + 1}

and

J_{1 ρ}^{n}

(given in (36) and (40), respectively) in (37), we have the following updating equations for the electric field components at time step

n + 1

,

\begin{array}{l} (A^{n + 1} - a D_{z}^{2}) E_{x}^{n + 1} = E_{x}^{^{n + \frac{1}{2}}} - a D_{x} D_{z} E_{z}^{^{n + \frac{1}{2}}} + b D_{y} H_{z}^{^{n + \frac{1}{2}}} - b D_{z} H_{y}^{^{n + \frac{1}{2}}}, \\ (A^{n + 1} - a D_{x}^{2}) E_{y}^{n + 1} = E_{y}^{^{n + \frac{1}{2}}} - a D_{y} D_{x} E_{x}^{^{n + \frac{1}{2}}} + b D_{z} H_{x}^{^{n + \frac{1}{2}}} - b D_{x} H_{z}^{^{n + \frac{1}{2}}}, \\ (A^{n + 1} - a D_{y}^{2}) E_{z}^{n + 1} = E_{z}^{^{n + \frac{1}{2}}} - a D_{z} D_{y} E_{y}^{^{n + \frac{1}{2}}} + b D_{x} H_{y}^{^{n + \frac{1}{2}}} - b D_{y} H_{x}^{^{n + \frac{1}{2}}}, \end{array}

where

A^{n + 1} = 1 + c S_{R}^{n + 1} + d {| E^{n + 1} |}^{2}

and

a = \frac{Δ t^{2}}{4 μ ε_{0} ε_{\infty}}, b = \frac{Δ t}{2 ε_{0} ε_{\infty}}, c = \frac{1}{ε_{\infty}}, d = \frac{α χ_{0}^{(3)}}{ε_{\infty}} .

Notice that in the left hand side of (41) we have $(1 - a D_{ρ}^{2})$ where a is a constant and $D_{ρ}^{2}$ is the second order central finite-difference operator along $ρ = x$ , y and z directions, for instance,

D_{x}^{2} f_{i, j, k} = \frac{1}{{(Δ x)}^{2}} (f_{i + 1, j, k} - 2 f_{i, j, k} + f_{i - 1, j, k}) .

Therefore, updating the electric field components at $n + \frac{1}{2}$ through (41) requires solving a tridiagonal system of linear equations which can be achieved using the tridiagonal matrix algorithm (Thomas algorithm) [45].

Recalling the definition of $A^{n + 1}$ in (43), in the left hand side of (42), there is a term of ${| E^{n + 1} |}^{2} E_{ρ}^{n + 1}$ ( $ρ = x, y, z$ ), where,

{| E^{n + 1} |}^{2} = {(E_{x}^{n + 1})}^{2} + {(E_{y}^{n + 1})}^{2} + {(E_{z}^{n + 1})}^{2},

thus, the updating equations for

E_{ρ}^{n + 1}

are strongly coupled to each other. Considering that these updating equations (the three equations given in given in (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

E^{n + \frac{1}{2}} = \frac{E^{n + 1} + E^{n}}{2} .

According to (47) , in (43), $E^{n + 1}$ can be substituted by

E^{n + 1} = 2 E^{n + \frac{1}{2}} - E^{n}

which results in

A^{n + 1} = 1 + c S_{R}^{n + 1} + d {| 2 E^{n + \frac{1}{2}} - E^{n} |}^{2} .

Therefore, when solving (42), $A^{n + 1}$ is a known constant since the field components at earlier time steps (i.e., $n + \frac{1}{2}$ , n) are known. Having $A^{n + 1}$ as a known constant significantly simplifies solving the updating equations of E_x, E_y and E_z (i.e., (42)) since:

The updating equations of E_x, E_y and E_z (i.e., (42)) are not longer nonlinear.
The coupling between the updating equations for E_x, E_y and E_z is removed.

In such a way, the tridiagonal matrix algorithm [45] can be applied to efficiently solve (42).

Finally, the field updating procedure in each iteration of the algorithm is summarized in the following seven steps (sub-iterations):

Explicitly update $J_{L x}^{n + \frac{1}{2}}$ , $J_{L y}^{n + \frac{1}{2}}$ , $J_{L z}^{n + \frac{1}{2}}$ using (10).
Implicitly update $E_{x}^{n + \frac{1}{2}}$ , $E_{y}^{n + \frac{1}{2}}$ , $E_{z}^{n + \frac{1}{2}}$ using (41).
Explicitly update $H_{x}^{n + \frac{1}{2}}$ , $H_{y}^{n + \frac{1}{2}}$ , $H_{z}^{n + \frac{1}{2}}$ using (34).
Explicitly update $S_{R}^{n + 1}$ using (21).
Explicitly update $A^{n + 1}$ using (49).
Implicitly update $E_{x}^{n + 1}$ , $E_{y}^{n + 1}$ , $E_{z}^{n + 1}$ using (42).
Explicitly update $H_{x}^{n + 1}$ , $H_{y}^{n + 1}$ , $H_{z}^{n + 1}$ using (36).

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 [46] 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 [19].

Fig. 1 A nonlinear slab in a TEM waveguide made of PEC and PMC walls.

Download Full Size | PDF

Fig. 2 Recorded values of the transmitted E_x through the 3-D nonlinear slab.

Download Full Size | PDF

Fig. 3 Normalized power intensity of the transmitted E_x through the 3-D nonlinear slab.

Download Full Size | PDF

Table 1. Run Times of Solving the 3-D Nonlinear Slab

View Table

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, $χ_{0}^{(3)}$ ) was considered. The characteristic parameters of the slab were the same as in [19],which are $ε_{r} = 2.25 α = 0.7$ , $χ_{0}^{(3)} = 7 \times 10^{- 2} m^{2} / V^{2}$ , $τ_{1} = 12.2$ fs, $τ_{2} = 32$ fs, and

ω_{R}^{2} = \frac{τ_{1}^{2} + τ_{2}^{2}}{τ_{1}^{2} τ_{2}^{2}}, δ_{R} = \frac{1}{τ_{2}} .

The slab was illuminated by an x-polorized plane wave having two carrier frequencies modulated with a Gaussian waveform,

E_{x} = exp [\frac{- {(t - t_{0})}^{2}}{τ^{2}}] \times (cos (2 π f_{1} t) + cos (2 π f_{2} t))

where

f_{1} = 300

THz,

f_{2} = 400

THz,

t_{0} = 86.66

fs and

τ = 28.89

fs. To simulate plane wave transmission, as shown in Fig. 1, the slab was placed in a TEM waveguide made of the perfect electric conductor (PEC) and the perfect magnetic conductor (PMC) walls [47] 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

20 \times 20 \times 1300

cells with a size of

Δ x = Δ y = Δ z = 8

nm

\approx λ_{0} (f = 400 THz) / 75

and terminated by PEC, PMC and Mur’s first-order absorbing boundary condition (ABC) in the x, y, and z directions, respectively [47].

For comparison, the problem was simulated using the explicit FDTD method reported in [19] (as a reference) and using the proposed unconditionally-stable FDTD method. The CFL condition limits the time step for the explicit FDTD method to $Δ t = Δ t_{CFL} = Δ x / \sqrt{3} c_{0} = 1.5 \times 10^{- 2}$ 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 $Δ t = 2 Δ t_{CFL}$ , $4 Δ t_{CFL}$ , $8 Δ t_{CFL}$ and $16 Δ t_{CFL}$ . The simulations were run for a large number of time steps (1,400,000 $Δ t_{CFL}$ ). 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 E_x was recorded which is shown in Fig. 2 (for the explicit FDTD method and the proposed method with $Δ t = 8 Δ t_{CFL}$ ). 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 $Δ t = 2 Δ t_{CFL}$ , $4 Δ t_{CFL}$ , $8 Δ t_{CFL}$ and $16 Δ t_{CFL}$ . In addition to the peak power intensity at the carrier frequencies, $f_{1} = 300$ THz and $f_{2} = 400$ THz, Fig. 3 shows two peaks at 200 THz and 500 THz corresponding to the third harmonic frequencies of $2 f_{1} - f_{2}$ and $2 f_{2} - f_{1}$ 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 $Δ t = 2 Δ t_{CFL}$ , $4 Δ t_{CFL}$ , and $8 Δ t_{CFL}$ ) and the explicit FDTD method is observed in Fig. 3; however, by setting $Δ t = 16 Δ t_{CFL}$ , they diverge. Notice that although in Fig. 2, the results of the proposed method with $Δ t = 8 Δ t_{CFL}$ is shown, the proposed method with smaller time step (i.e., $Δ t = 2 Δ t_{CFL}$ and $4 Δ t_{CFL}$ ) 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 [48], i.e., $α = 0.7$ , $χ_{0}^{(3)} = 2 \times 10^{- 19} m^{2} / V^{2}$ , $ε_{r} = 13.66$ , $τ_{1} = 12.2$ fs, and $τ_{2} = 32$ fs. As shown in Fig. 4, the nonlinear silicon medium is surrounded by dielectric sections having a relative permittivity of $ε_{r} = 13.66$ 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 $Δ x = Δ y = 8$ nm and is terminated by the ABC.

First, to demonstrate the accuracy of the proposed method, a cosine-modulated Gaussian source having a waveform of

H_{z} = H_{0} exp [\frac{- {(t - t_{0})}^{2}}{τ^{2}}] cos (2 π f_{3} t)

is applied at grid (450,70), where

H_{0} = 10^{8}

A/m,

f_{3} = 375

THz,

τ = 9.63

fs and

t_{0} = 3 τ

. The simulations were performed using the explicit FDTD method [19] (as a reference) with

Δ t = Δ t_{CFL} = Δ x / \sqrt{2} c_{0} = 1.8 \times 10^{- 2}

fs and the proposed method with

Δ t = 10 Δ t_{CFL}

, and the values of E_x were recorded at grid (1750, 70) (see Fig. 4). The recorded E_x is shown in Fig. 5 indicating a strong agreements between the results of the proposed and the reference methods.

Fig. 4 A nonlinear silicon waveguide configuration.

Download Full Size | PDF

Fig. 5 Recorded values of E_x at grid (1750,70).

Download Full Size | PDF

Fig. 6 Power intensity distribution over the silicon waveguide.

Download Full Size | PDF

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 $Δ x$ (see Fig. 4) having an initial transverse profile of

H_{z} = H_{0} cos (2 π f_{0} t) sech (y / w)

where

H_{0} = 10^{8}

A/m,

f_{0} = 375

THz (

λ_{0} = 800

nm) and

w = 50

nm. The simulations were run for a large number of time steps until a steady state was achieved. Fig. 6 depicts the magnitude of the H_z 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.

Fig. 7 Scattering of spatial optical solitons by subwavelength air hole with a size of 150 nm × 250 nm whose center is located at x = 5 μm and y = 0 μm. The figure indicates |E| in the computational domain.

Download Full Size | PDF

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 $β_{1} = 0.69617$ , $β_{2} = 0.40794$ , $β_{3} = 0.89748$ , $γ_{1} = 2.7537 \times 10^{16} rad / s$ , $γ_{2} = 1.6205 \times 10^{16} rad / s$ , $γ_{3} = 1.90342 \times 10^{14} rad / s$ . And the Kerr and Raman nonlinear characteristic parameters were considered as $α = 0.7$ , $χ_{0}^{(3)} = 1.89 \times 10^{- 22} m^{2} / V^{2}$ , $ε_{r} = 13.66$ , $τ_{1} = 12.2$ fs, and $τ_{2} = 32$ fs.

A computational domain consisting of 1500 × 600 cells (along x and y directions, respectively) where the spatial resolution is set to $Δ x = Δ y = 10$ 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 $Δ x$ having an initial transverse profile of given in (53), where $H_{0} = 2 \times 10^{8}$ A/m, $f_{0} = 692$ THz ( $λ_{0} = 433$ nm) and $w = 215$ nm. An air hole with a size of 250 nm × 250 nm was assumed in a way that its center is located at grid $x = 5 μ m$ and $y = 0 μ m$ . 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 ( $| E |$ ) 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 [49].

4. Conclusion

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.

References

1. L. Yin, J. Zhang, P. M. Fauchet, and G. P. Agrawal, “Optical switching using nonlinear polarization rotation inside silicon waveguides,” Opt. Lett. 34, 476–478 (2009). [CrossRef] [PubMed]

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]

3. R. Claps, D. Dimitropoulos, V. Raghunathan, Y. Han, and B. Jalali, “Observation of stimulated Raman amplification in silicon waveguides,” Opt. Express 11, 1731–1739 (2003). [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]

21. J. H. Greene and A. Taflove, “General vector auxiliary differential equation finite-difference time-domain method for nonlinear optics,” Opt. Express 14, 8305–8310 (2006). [CrossRef] [PubMed]

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, $χ^{(2)}$ , and $χ^{(3)}$ 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 $S i_{1 - x} G e_{x}$ 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]

Classical FDTD	This Work
$Δ t_{CFL} = Δ x / \sqrt{3} c_{0}$	$Δ t = 2 Δ t_{CFL}$	$Δ t = 4 Δ t_{CFL}$	$Δ t = 8 Δ t_{CFL}$	$Δ t = 16 Δ t_{CFL}$
2968 (sec)	3001 (sec)	1532 (sec)	667 (sec)	379 (sec)

Unconditionally stable FDTD algorithm for 3-D electromagnetic simulation of nonlinear media

Abstract

1. Introduction

2. Method development

2.1. Polarization currents

2.2. A novel ADI-FDTD scheme for third-order nonlinear media

3. Numerical example and discussion

3.1. First example

3.2. Second example

3.3. Third example

4. Conclusion

References

Cited By

Figures (7)

Tables (1)

Equations (53)

Optics Express