Nonstandard Finite Difference Time Domain Methodology to Simulate Light Propagation in Nonlinear Materials

James B. Cole

University of Tsukuba, Japan
cole@cavelab.cs.tsukuba.ac.jp

Submitted On: July 17, 2023; Accepted On: April 24, 2024

ABSTRACT

We extend the nonstandard (NS) finite difference time domain (FDTD) methodology, originally developed to solve Maxwell’s equations in linear materials, to nonlinear ones. We validate it by computing harmonics generation in a nonlinear dielectric and comparing with theory. The methodology also applies to the quantum electrodynamics that describes the interaction of charged particles with electromagnetic fields, and also to the Ginzburg-Landau model of superconductivity.

Index Terms: Finite difference time domain (FDTD), nonlinear optics, nonlinear susceptibility, nonstandard FDTD, quantum electrodynamics, superconductivity.

I. INTRODUCTION

The conventional or standard (S) finite difference time domain (FDTD) methodology [1] is widely used for linear electromagnetic calculations, but its accuracy is low relative to the computational cost. At wavelength $\lambda$ for space-step size $h$, its error scales as ${(h\mathrm{/}\lambda )}^4$ and in three dimensions its computational also cost scales as ${(\lambda \mathrm{/}h)}^4$. We [2] have introduced what is called a nonstandard (NS) FDTD methodology [3] for which the error scales as ${(h\mathrm{/}\lambda )}^8$, but computational cost still scales as ${(\lambda \mathrm{/}h)}^4$.

Nonlinear problems are generally difficult to solve analytically, but numerical methods also often fail to yield good solutions unless the discretization steps are small, and even then numerical instability can arise. A classic example is the logistic equation, the discrete form of which is a well-known example of deterministic chaos. In this paper we extend the NS-FDTD methodology to solve Maxwell’s equations for nonlinearmaterials.

Although not the topic of this paper, the methodology is also useful to solve certain problems in quantum mechanics. For example, the Hamiltonian of a changed particle (of mass $m,$charge $q$) in an electromagnetic field [4] (vector potential A, scalar potential $\phi$) is nonlinear in the form:

$$H=\frac{1}{2m}{\left(-i\mathrm{\hslash }\nabla -q\mathbf{A}\right)}^2+q\phi .$$

(1)

II. NONLINEAR OPTICS

In a nonlinear dielectric the electric displacement is $\mathbf{D}={\epsilon }_0\left(\mathbf{E}+{\chi }^{(1)}\mathbf{E}+{\chi }^{(2)}{\mathbf{E}}^2+\mathrm{\cdots }\right)$, where ${\epsilon }_0$ is the vacuum electric permeability, and $\mathbf{E}$ the electric field. We use units in which ${\epsilon }_0=1$ and the vacuum magnetic permeability is ${\mu }_0=1\mathrm{/}{c}^2$ ($c=$vacuum velocity of light). Retaining only second order nonlinearity, and defining $\epsilon =1+{\chi }^{(1)},$ then $\mathbf{D}=\epsilon \mathbf{E}+{\chi }^{(2)}{\mathbf{E}}^2$. Taking the magnetic susceptibility to be ${\mu }_0$ everywhere, Maxwell’s equations become:

$\mu {}_0\partial {}_t\mathbf{H}$	$=-\nabla \times \mathbf{E},$
${\partial }_t\mathbf{D}$	$=\nabla \times \mathbf{H},$	(2)

and the index of refraction is $n=\sqrt{\epsilon }$. In a linear material $\nabla \bullet \mathbf{D}=0$ implies $\nabla \bullet \mathbf{E}=0,$ but this is not true in a nonlinear one. If, however, ${\chi }^{(2)}$ is small, $\nabla \bullet \mathbf{E}=0$ is a good approximation [5], and Maxwells equations reduce to a nonlinear wave equation of the form:

$$\left({\partial }_t^2-\frac{c^2}{\epsilon }{\nabla }^2\right)\mathbf{E}(\mathit{\text{x}},t)=-\frac{{\chi }^{(2)}}{\epsilon }{\partial }_t^2\mathbf{E}{(\mathit{\text{x}},t)}^2.$$

(3)

For simplicity, we first develop the finite difference model of (3) in one dimension, where it reduces to:

$$\left({\partial }_t^2-\frac{c^2}{n^2}{\partial }_x^2\right)E(x,t)=-\frac{{\chi }^{(2)}}{n^2}{\partial }_t^{2}E{(x,t)}^2.$$

(4)

Before proceeding further, we first introduce the nonstandard methodology for the linear wave equation.

III. FINITE DIFFERENCE MODELS FOR THE LINEAR WAVE EQUATION

A. Notation and definitions

Define the partial difference operator ($d_t$) by $d_{t}f(x,t)=$ $f(x,t+\mathrm{\Delta }t/2)-f(x,t-\mathrm{\Delta }t/2)$. Then it is easy to show that $d_t^2=d_t{d}_t$ is given by:

$$d_t^{2}f(x,t)=f(x,t+\mathrm{\Delta }t)+f(x,t-\mathrm{\Delta }t)-2f(x,t).$$

(5)

The second derivative is thus approximated by:

$$\mathrm{\Delta }{t}^2{\partial }_t^{2}f(x,t)\cong d_t^{2}f(x,t).$$

(6)

The operator $d_x^2$ is defined analogously to $d_t^2$ and:

$$\mathrm{\Delta }{x}^2{\partial }_x^{2}f(x,t)\cong d_x^{2}f(x,t).$$

(7)

We now construct finite difference models of the homogeneous wave equation:

$$\left({\partial }_t^2-v^2{\partial }_x^2\right)\psi (x,t)=0.$$

(8)

B. Standard finite difference model

Writing $\overline{v}=v\mathrm{\Delta }t\mathrm{/}\mathrm{\Delta }x,$ and replacing the derivatives in (8) with the FD expressions (6) and (7), gives the conventional or standard (S) finite difference (FD) model of (8):

$$\left(d_t^2-{\overline{v}}^2{d}_x^2\right)\psi (x,t)=0.$$

(9)

General solutions of (8) are$f(x\mp vt)$where $f$ is arbitrary. Substituting $f(x\mp vt)$ into the S-FD model above, we find:

$$\left[d_t^2-{\overline{v}}^2{d}_x^2\right]f(x\mp vt)=-\left(\frac{1-{\overline{v}}^4}{12}\right)\mathrm{\Delta }{x}^4{f}^{(4)}(x\mp vt)+\mathrm{\cdots }.$$

(10)

The right side of (10) is the model error. Although it vanishes for $\overline{v}=1$, across multiple media in which $v$ varies, maintaining $\overline{v}=1$ by adjusting $\mathrm{\Delta }x$ or $\mathrm{\Delta }t$ gives rise to large errors on the media boundaries and is thus of limited use in practice.

C. Nonstandard finite difference model

It is, however, possible to construct an exact FD model with respect to harmonic waves, $\phi (x,t)=e^{i\left(kx-\omega t\right)},$where $\omega =vk$. Writing $\overline{\omega }=\omega \mathrm{\Delta }t$, $\overline{k}=k\mathrm{\Delta }x$, and substituting $\phi$ into (9) gives:

$$\left[d_t^2-{\overline{v}}^2{d}_x^2\right]\phi =4\left[-{\sin }^2\left(\overline{\omega }\mathrm{/}2\right)+{\overline{v}}^2{\sin }^2\left(\overline{k}\mathrm{/}2\right)\right]\phi .$$

(11)

The right side of (11) can be made to vanish with the substitution $\overline{v}\to \tilde{v}$, where:

$$\tilde{v}=\sin \left(\overline{\omega }\mathrm{/}2\right)\mathrm{/}\sin \left(\overline{k}\mathrm{/}2\right).$$

(12)

Thus, an exact model of the wave equation with respect to harmonic waves is:

$$\left(d_t^2-{\tilde{v}}^2{d}_x^2\right)\psi (x,t)=0.$$

(13)

This FD model is exact because $\phi$ is a solution of both the wave equation (8) and its model (13). This is an example of what is called an NS model [3]. Expanding $d_t^2\psi (x,t)$ via (5) and solving for $\psi (x,t+\mathrm{\Delta }t)$, we obtain an exact NS FDTD algorithm:

$$\psi (x,t+\mathrm{\Delta }t)=-\psi (x,t-\mathrm{\Delta }t)+\left[2+{\tilde{v}}^2{d}_x^2\right]\psi (x,t).$$

(14)

D. Numerical stability and accuracy

The numerical stability condition [2, 6] for the S- and NS-FDTD algorithms is:

$$\frac{v\mathrm{\Delta }t}{\mathrm{\Delta }x}\le 1.$$

(15)

In the case $v\mathrm{\Delta }t\mathrm{/}\mathrm{\Delta }x=1$, the S- and NS-FD models are equivalent and exact with respect to any waveform. However, whatever the value of $0<v\mathrm{\Delta }t\mathrm{/}\mathrm{\Delta }x\le 1$, for a supposition different frequencies, the shortest period ($T_{\min }$) corresponding to the highest frequency (${\nu }_{\max }$) must satisfy the Nyquist sampling criterion [7]:

$$T_{\min }\mathrm{/}\mathrm{\Delta }t>2,$$

(16)

and the minimum wavelength (${\lambda }_{\min }$) must satisfy:

$${\lambda }_{\min }\mathrm{/}\mathrm{\Delta }x>2.$$

(17)

E. Wave equation with a source

To iterate the FDTD algorithm, two initial fields are needed. They can be generated by turning on sources at time = 0. The wave equation with a source $s(x,t)$ is:

$$\left({\partial }_t^2-v^2{\partial }_x^2\right)\psi (x,t)=s(x,t).$$

(18)

Standard finite difference model

Substituting FD expressions of the derivatives in (18) the S-FD model is:

$$\left(d_t^2-{\overline{v}}^2{d}_x^2\right)\psi (x,t)=\mathrm{\Delta }{t}^{2}s(x,t).$$

(19)

It is interesting to note that while the stability conditions of (15), (16) and (17) still hold, the S-FDTD algorithm for the wave equation with a source is not exact even for $v\mathrm{\Delta }t\mathrm{/}\mathrm{\Delta }x=1$.

Nonstandard finite difference model

To derive the NS-FD model, we examine the analytic solution of (18). Imposing the initial conditions that $\psi$ and its first time derivative vanish for $t\le 0$, ${\psi (x,t)|}_{t\le 0}=$ ${{\partial }_t\psi (x,t)|}_{t\le 0}=0$, the Green’s function that solves (18) is:

$$G(x-x^{\prime },t-t^{\prime })=\frac{1}{2v}\mathrm{\Theta }\left(v[t-t^{\prime }]-\left|x-x^{\prime }\right|\right),$$

(20)

where the step function is defined by $\mathrm{\Theta }(t)=0$ for $t<0$ and $\mathrm{\Theta }(t)=1$ for $t\ge 0$. Using (20) it can be shown that the harmonic point source:

$$s_0(x,t)=2v\omega \delta (x)\mathrm{\Theta }(t)\cos (\omega t),$$

(21)

generates an outgoing unit sine wave:

$${\psi }_0(x,t)=\mathrm{\Theta }\left(\omega t-k\left|x\right|\right)\sin \left(t-\left|x\right|\mathrm{/}v\right).$$

(22)

Modeling $\delta (x)$ as ${\delta }_{x,0}\mathrm{/}\mathrm{\Delta }x$, where ${\delta }_{x,0}=1$ for $x=0$ and ${\delta }_{x,0}=0$ for $x\ne 0$, and putting $s\to s_0$ in (19), the S-FD model becomes:

$$\left(d_t^2-{\overline{v}}^2{d}_x^2\right)\psi (x,t)=2v\omega \frac{\mathrm{\Delta }{t}^2}{\mathrm{\Delta }x}{\delta }_{x,0}\mathrm{\Theta }(t)\cos (\omega t).$$

(23)

We now postulate the NS-FD model to be:

$$\left(d_t^2-{\tilde{v}}^2{\partial }_x^2\right)\psi (x,t)=2v\omega \tilde{A}{\delta }_{x,0}\mathrm{\Theta }(t)\cos (\omega t),$$

(24)

where $\mathrm{\Delta }{t}^2\mathrm{/}\mathrm{\Delta }x\to \tilde{A}$, which is to be determined. Requiring that ${\psi }_0$ be a solution of both (18), with $s=s_0$, and of (24), we find [3]:

$$\tilde{A}=\frac{2}{\omega v}\frac{{\sin }^2(\overline{\omega }\mathrm{/}2)}{\tan (\overline{k}\mathrm{/}2)}.$$

(25)

When a source is abruptly switched on in FDTD calculations it produces extraneous frequency components which give rise to large errors [3]. This is remedied by replacing the step function with a slow switch-on function. A commonly used one is:

$${\mathrm{\Theta }}_{\mathrm{g}}(t)=\mathrm{\Theta }(t)\left[1-e^{-{\beta }^2{t}^2}\right].$$

(26)

Taking $1\mathrm{/}\beta$ to be several wave periods suffices to suppress the errors. The NS source model which generates ${\psi }_0$ in a NS-FDTD calculation is thus:

$${\tilde{s}}_0(x,t)=4\frac{{\sin }^2(\overline{\omega }\mathrm{/}2)}{\tan (\overline{k}\mathrm{/}2)}{\delta }_{x,0}{\mathrm{\Theta }}_{\mathrm{g}}(t)\cos (\omega t).$$

(27)

With $\tilde{A}$ given by (25) the NS-FD model (24) is exact. As expected, in the limits $\mathrm{\Delta }x\to 0$ and $\mathrm{\Delta }t\to 0$, the NS-FD model reduces to the S-FD one.

The NS-FDTD model of the wave equation with a harmonic point source is thus:

$$\begin{array}{ccc}\hfill \psi (x,t+\mathrm{\Delta }t)& \hfill =\hfill & -\psi (x,t-\mathrm{\Delta }t)\hfill \\ & & +\left[2+{\tilde{v}}^2{d}_x^2\right]\psi (x,t)+\tilde{s}(x,t).\hfill \end{array}$$

(28)

When $\tilde{s}={\tilde{s}}_0$ the iteration of (28) generates the outgoing unit sine wave given by (22), where $k$ and $\omega$ are related by:

$$\sin \left(\omega \mathrm{\Delta }t\mathrm{/}2\right)=\tilde{v}\sin \left(k\mathrm{\Delta }x\mathrm{/}2\right).$$

(29)

Initialization and iteration

The FDTD calculation is initialized by taking

$$\psi (x,-\mathrm{\Delta }t)=\psi (x,0)=0,$$

(30)

and switching on the source at $t=0$ generates the incident field.

The boundary conditions at material interfaces are determined by the wave equation itself, viz. continuity of both the field and its first partial derivative with respect to position. Since null fields obviously satisfy these conditions, and because FDTD derives directly from the wave equation, the generated fields automatically satisfy the boundary conditions as they impinge upon material interfaces when the algorithm is iterated.

F. Multi-frequency NS-FDTD

It might seem that NS-FDTD is applicable only to monochromatic waves, but it is also valid for multiple frequencies. For a fixed value of $\tilde{v}$, models (13) and (24) are exact for any angular frequency-wavenumber pair $({\omega }_i,k_i)$ related by:

$$\sin \left({\omega }_i\mathrm{\Delta }t\mathrm{/}2\right)=\tilde{v}\sin \left(k_i\mathrm{\Delta }t\mathrm{/}2\right).$$

(31)

Thus, the NS-FD model is exact with respect to a multi-frequency wave of the form:

$${\phi }_{\mathrm{\Sigma }}(x,t)=\sum _i{a}_i{e}^{i\left(k_{i}x\pm {\omega }_{i}t\right)},$$

(32)

which is produced by multi-frequency source superposition. For example:

$${\tilde{s}}_{\mathrm{\Sigma }}(x,t)=4{\delta }_{x,0}{\mathrm{\Theta }}_{\mathrm{g}}(t)\sum _i\frac{{\sin }^2({\overline{\omega }}_i\mathrm{/}2)}{\tan ({\overline{k}}_i\mathrm{/}2)}\cos ({\omega }_{i}t)$$

(33)

generates a frequency superposition of unit sine waves. Such a superposition is useful for high accuracy computations of reflection or transmission spectra. The maximum frequency is limited by the time step according to the Nyquist sampling criterion given by (16).

G. Nonstandard model of refractive index

Let $k_0$ be the vacuum wave number and ${\omega }_0=c{k}_0$ be the angular frequency of a light wave. In a medium of refractive index $n$ the wavenumber is $k=n{k}_0$. From (12) the algorithmic vacuum velocity of light in the NS-model is:

$$\tilde{c}=\sin \left({\overline{\omega }}_0\mathrm{/}2\right)\mathrm{/}\sin \left({\overline{k}}_0\mathrm{/}2\right).$$

(34)

In a medium of refractive index $n$ the algorithmic velocity of light is:

$$\tilde{v}=\sin \left({\overline{\omega }}_0\mathrm{/}2\right)\mathrm{/}\sin \left(n{\overline{k}}_0\mathrm{/}2\right).$$

(35)

Define the nonstandard index of refraction to be$\tilde{n}=\tilde{c}\mathrm{/}\tilde{v}$:

$$\tilde{n}=\frac{\sin (n{\overline{k}}_0\mathrm{/}2)}{\sin ({\overline{k}}_0\mathrm{/}2)}.$$

(36)

The nonstandard FDTD algorithm (28) is thus:

$$\begin{array}{ccc}\hfill \psi (x,t+\mathrm{\Delta }t)& \hfill =\hfill & -\psi (x,t-\mathrm{\Delta }t)\hfill \\ & & +\left[2+\frac{{\tilde{c}}^2}{{\tilde{n}}^2}{d}_x^2\right]\psi (x,t)+\tilde{s}(x,t).\hfill \end{array}$$

(37)

Before proceeding to the nonlinear NS-model we introduce a simplified and abbreviated notation. Discretizing space time as $x=\chi \mathrm{\Delta }x$, $t=\tau \mathrm{\Delta }t$ and defining $\psi (x,t)={\psi }_{\chi }^{\tau }$,

($\tau$, $\chi$ integers), (37) is compactly rewritten as:

$${\psi }_{\chi }^{\tau +1}=-{\psi }_{\chi }^{\tau -1}+\left[2+\frac{{\tilde{c}}^2}{{\tilde{n}}^2}{d}_x^2\right]{\psi }_{\chi }^{\tau }+{\tilde{s}}_{\chi }^{\tau },$$

(38)

with the definitions $d_x^2{\psi }_{\chi }^{\tau }={\psi }_{\chi +1}^{\tau }+{\psi }_{\chi -1}^{\tau }-2{\psi }_{\chi }^{\tau }$, $d_t^2{\psi }_{\chi }^{\tau }=$ ${\psi }_{\chi }^{\tau +1}+{\psi }_{\chi }^{\tau -1}-2{\psi }_{\chi }^{\tau }$ (compare with equation 5).

IV. FINITE DIFFERENCE MODELS OF THE NONLINEAR WAVE EQUATION

A. Standard finite difference model

First construct the S-FD model of (4). For notational clarity suppress the spatial dependence and denote time as a subscript, thus $E_{\tau }=E(x,\tau \mathrm{\Delta }t)$. Writing $\overline{c}=c\mathrm{\Delta }t\mathrm{/}\mathrm{\Delta }x$, the S-FD model of nonlinear wave equation (4) is:

$$\left(d_t^2-\frac{{\overline{c}}^2}{n^2}{d}_x^2\right)E_{\tau }=-\frac{{\chi }^{(2)}}{n^2}\left[E_{\tau +1}^2+E_{\tau -1}^2-2E_{\tau }^2\right].$$

(39)

The intractability of this model is immediately evident. Since (39) is quadratic in $E_{\tau +1}$, there are two solutions and it is unclear a priori which one to use.

B. Nonstandard finite difference model

In nonstandard models a term of power $m$ (a positive integer), such as $E{(x,t)}^m$ is modeled as [3, 8]:

$$E_{\tau }^m=E_{\tau }^m{E}_{\tau -1}^m\mathrm{\cdots }{E}_{\tau -m+1}^m.$$

(40)

The NS-model of $E{(x,t)}^2$is thus $E_{\tau }^2=E_{\tau }{E}_{\tau -1}$. Modeling $d_t^{2}E{(x,t)}^2$ as:

$$d_t^2{E}_{\tau }^2=E_{\tau +1}{E}_{\tau }+E_{\tau -1}{E}_{\tau -2}-2E_{\tau }{E}_{\tau -1},$$

(41)

postulate the NS-FD model of (4) to be:

$\left(d_t^2-{\displaystyle \frac{{\tilde{c}}^2}{{\tilde{n}}^2}}{d}_x^2\right)E_{\tau }=-{\displaystyle \frac{{\chi }^{(2)}}{{\tilde{n}}^2}}$
$\mathrm{\hspace{1em}\hspace{1em}}\left[E_{\tau +1}{E}_{\tau }+E_{\tau -1}{E}_{\tau -2}-2E_{\tau }{E}_{\tau -1}\right]+{\tilde{s}}_{\tau },$	(42)

where ${\tilde{s}}_{\tau }$ is a nonstandard source term that generates the fields. Solving for $E_{\tau +1}$ we find:

$\mathrm{\hspace{1em}\hspace{1em}}\mathit{\hspace{1em}\hspace{1em}}{E}_{\tau +1}=$
${\displaystyle \frac{2E_{\tau }\left[{\tilde{n}}^2+{\chi }^{(2)}{E}_{\tau -1}\right]-E_{\tau -1}\left[{\tilde{n}}^2+{\chi }^{(2)}{E}_{\tau -2}\right]+{\tilde{c}}^2{d}_x^2{E}_{\tau }}{{\tilde{n}}^2+{\chi }^{(2)}{E}_{\tau }}}+{\tilde{s}}_{\tau }.$	(43)

Whereas the FDTD algorithms for the linear wave equation require two initial fields, the NS-FDTD for nonlinear wave equation requires three. To iterate (4.2) Take $E_{-2}=E_{-1}=$ $E_0=0$, and switch on the source at $t=0$.

C. Computational model

We take the computational domain to be $0\le x\le N\mathrm{\Delta }x$. In Fig. 1 the nonlinear dielectric (magenta) is immersed in vacuum (white) and illuminated by a point harmonic source (S) of angular frequency ${\omega }_0$. S is located at $x=p>\mathrm{\Delta }x$, away from the computational boundary, thus:

$${\tilde{s}}_0(x,t)=4\frac{{\sin }^2({\overline{\omega }}_0\mathrm{/}2)}{\tan ({\overline{k}}_0\mathrm{/}2)}{\delta }_{x,p}{\mathrm{\Theta }}_{\mathrm{g}}(t)\cos ({\omega }_{0}t),$$

(44)

where $c={\omega }_0\mathrm{/}{k}_0$. We choose $\mathrm{\Delta }t$ and $\mathrm{\Delta }x$ such that $\overline{c}=$ $c\mathrm{\Delta }t\mathrm{/}\mathrm{\Delta }x=1$, and thus $\tilde{c}=\overline{c}=1$. Thus (only) in vacuum, as noted in Section II D, the S-FDTD and NS-FDTD algorithms are equivalent and exact. With the choices above in the vacuum:

$$E_{\tau +1}=-E_{\tau -1}+\left[2+d_x^2\right]E_{\tau }+{\tilde{s}}_{\tau },$$

(45)

where $s_{\tau }=s_0(x,t)$. Because the boundaries of the computational domain are vacuum and $\overline{c}=1$, the Mur absorbing boundary [9, 10] is exact, and is given by:

$E(0,t+\mathrm{\Delta }t)$	$=E(\mathrm{\Delta }x,t),$	(46)
$E(N\mathrm{\Delta }x,t+\mathrm{\Delta }t)$	$=E(N\mathrm{\Delta }x-\mathrm{\Delta }x,t).$	(47)

Taking the vacuum wavelength of the incident field be ${\lambda }_0=1200\mathrm{n}\mathrm{m}$, let $\mathrm{\Delta }x={\lambda }_0\mathrm{/}64=18.75\mathrm{nm}$, which implies the wave period $T_0={\lambda }_0\mathrm{/}c$. Setting

$$\mathrm{\Delta }t=T_0\mathrm{/}64=6.25\times {10}^{-17}\sec$$

(48)

gives $\overline{c}=\tilde{c}=1$. For this choice, the source simplifies to:

$${\tilde{s}}_0(x,t)=2\sin ({\overline{\omega }}_0){\delta }_{x,0}{\mathrm{\Theta }}_{\mathrm{g}}(t)\cos ({\omega }_{0}t).$$

(49)

Figure 1: Computational set-up to simulate propagation and harmonic generation in a non-linear dielectric illuminated by a point harmonic source S located in vacuum. The time dependence of the transmitted field is recorded at point O outside the dielectric.

V. RESULTS AND COMPARISON WITH SEMI-ANALYTIC CALCULATION

Take $n$= 1.6, ${\chi }^{(2)}=0.05$, the source amplitude to be 1.1, ${\omega }_0\mathrm{/}2\pi =1\mathrm{/}64\mathrm{\Delta }t$, and the source rise time $1\mathrm{/}\beta =$ $4T_0$, where $T_0=2\pi \mathrm{/}{\omega }_0$.

Figure 2: Light from a source (yellow dot) enters the dielctric material (orange line). As the electric field (black curve) traverses the material, harmonics are generated.

Figure 2 is a snapshot of the position dependence of $E$. A time series of the electric field amplitude was collected at an observation point outside the dielectric after the source switched on and the field had completely traversed the dielectric. The data were analyzed with a discrete Fourier transform (DFT). The DFT amplitudes of the harmonic frequencies are shown in Fig. 3.

Figure 3: Harmonics generation in a nonlinear dielectric. Time steps from $1300\le t\mathrm{/}\mathrm{\Delta }t\le 1940$ were analyzed. The frequency unit is $\mathrm{\Delta }{t}^{-1}=$$16\bullet {10}^{15}$Hz.

In Fig. 4 we compare our simulation with a semi-analytic calculation based on the low depletion approximation [1]. The low depletion model assumes that energy is slowly transferred from the fundamental mode to the higher harmonics. This is the usual case when ${\chi }^{(\mathit{\text{2}})}$ is small.

Figure 4: Amplitudes first harmonic (black) and second harmonic (blue) with position and in the nonlinear material (orange). Comparison with theory (red) in the low depletion approximation.

VI. EXTENSION TO TWO AND THREE DIMENSIONS

The linear homogeneous wave equation in three dimensions is:

$$\left({\partial }_t^2-v^2{\nabla }^2\right)\psi (\mathit{\text{x}},t)=0,$$

(50)

where $\mathit{\text{x}}=(x,y,z)$. Defining $d_y^2$ and $d_z^2$ by analogy with $d_x^2$, define the Laplacian difference operator:

$${\mathbf{d}}^2=d_x^2+d_y^2+d_z^2.$$

(51)

Taking $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{\Delta }z=h$:

$${\nabla }^{2}f(x)\cong \frac{{\mathbf{d}}^2}{h^2}f(x).$$

(52)

The S-FD model for the three-dimensional wave equation is thus:

$$\left(d_t^2-{\overline{v}}^2{\mathbf{d}}^2\right)\psi (x,t)=0,$$

(53)

where $\overline{v}=v\mathrm{\Delta }t\mathrm{/}h$.

To construct the NS-model, define:

$${\tilde{\mathbf{d}}}^2={\mathbf{d}}^2+\left(\frac{1}{6}+\frac{{(k_{0}h)}^2}{180}\right)\left[d_x^2{d}_y^2+d_x^2{d}_z^2+d_y^2{d}_z^2\right]$$

(54)

$$+\frac{1}{30}{d}_x^2{d}_y^2{d}_z^2.$$

The motivation and derivation of this definition is found in [3]. The NS-FD model of wave equation (50) becomes

$$\left(d_t^2-{\tilde{v}}^2{\tilde{\mathbf{d}}}^2\right)\psi (x,t)=0,$$

(55)

where $\tilde{v}$ is given by (12) and $\overline{k}=kh$. The NS-source model remains the same but with $\mathrm{\Delta }x\to h$.

Table 1 lists the stability conditions for the S- and NS- algorithms.

Table 1: Stability conditions for the S- and NS- algorithms

	Maximum Stable Value
	Theoretical/Practical
	S-FDTD $\overline{v}$	NS-FDTD $\tilde{v}$
1-D	$1/1$	$1/1$
2-D	${\displaystyle \frac{\sqrt{2}}{2}}\cong 0.70/0.67$	$0.86/0.80$
3-D	${\displaystyle \frac{\sqrt{3}}{3}}\cong 0.57/0.45$	$0.80/0.70$

For S-FDTD, the maximum value of $\overline{v}$ is given, while for the NS-FDTD algorithm the maximum value of $\tilde{v}$ is given. The practical stability limits are somewhat lower due to the termination of the computational boundary. Details of the derivation are found in [2]. The greater stability of NS-FDTD allows the solution of problems using fewer time steps. The advantage over S-FDTD is greatest in three dimensions.

The NS-FD model has been validated against analytic solutions of Mie scattering [11, 12], as depicted in Fig. 5. An infinite plane wave (not shown) is incident from the left and the scattered field intensitycomputed.

Figure 5: Mie scattering of an infinite dielectric cylinder (white). Scattered intensity (shades of red) in the exterior was computed using the NS-FDTD and S-FDTD algorithms and compared with Mie theory. Index of refraction = 1.7, vacuum wavelength = 800 nm, cylinder radius = 600 nm, grid spacing = 100 nm.

In the vacuum (black) the wavelength is ${\lambda }_0=$800 nm and in the dielectric (white) it is ${\lambda }_d={\lambda }_0\mathrm{/}1.7$, thus in the vacuum $h\mathrm{/}{\lambda }_0=1\mathrm{/}8=0.125$, but in the dielectric $h\mathrm{/}{\lambda }_d=1.7\mathrm{/}8=$ $0.2125$, which is just slightly greater that the minimum allowed by the Nyquist criterion for a 2-dimensional (uniform) grid where $h\mathrm{/}\lambda$ must satisfy $h\mathrm{/}\lambda <1\mathrm{/}(2\sqrt{2})$ $\cong 0.3536$. Nonetheless the NS-FDTD error remains low. The theoretical error of the NS-FDTD calculation is ${\epsilon }_{\mathrm{NS}}\cong$${\left(kh\right)}^8\mathrm{/}438840$, while that of S-FDTD is ${\epsilon }_{\mathrm{S}}\cong {\left(kh\right)}^4\mathrm{/}48$.

The three-dimensional NS-FDTD algorithm for the nonlinear dielectric derives from the one-dimensional form (42) with the substitution $d_x^2\to {\tilde{\mathbf{d}}}^2$:

$\left(d_t^2-{\displaystyle \frac{{\tilde{c}}^2}{{\tilde{n}}^2}}{\tilde{\mathbf{d}}}^2\right)E_{\tau }=$	(56)
$\mathrm{\hspace{1em}\hspace{1em}}-{\displaystyle \frac{{\chi }^{(2)}}{{\tilde{n}}^2}}\left[E_{\tau +1}{E}_{\tau }+E_{\tau -1}{E}_{\tau -2}-2E_{\tau }{E}_{\tau -1}\right]+{\tilde{s}}_{\tau },$

where $E_{\tau }$ stands for the $x-$, $y-$, or $z-$component of $\mathbf{E}(\mathit{\text{x}},\tau \mathrm{\Delta }t)$.

VII. SUMMARY AND CONCLUSIONS

We introduced a high precision finite difference time domain algorithm derived from a nonstandard finite difference model to simulate electromagnetic propagation in nonlinear dielectrics. We validated the results of our simulation against an analytic calculation based on the low depletion approximation [5].

This NS-methodology can also be applied to other nonlinear problems, such as quantum electrodynamics in magnetic fields, and to higher order nonlinearities.

We introduced our methodology in one-dimension and extended it to two and three dimensions and have verified its high accuracy and numerical stability [2].

REFERENCES

[1] K. S. Kunz and R. J. Luebbers, The Finite Difference Time Domain Method for Electromagnetics. Boca Raton: CRC Press, 1993.

[2] J. B. Cole and S. Banerjee, Computing the Flow of Light, Bellingham: SPIE Press, 2017.

[3] R. Mickens, Advances in the Applications of Nonstandard Finite Difference Difference Schemes, Singapore: World Scientific, 2005.

[4] E. Merzbacher, Quantum Mechanics, 2nd ed. Hoboken: John Wiley, 1970.

[5] R. W. Boyd, Nonlinear Optics, 3${}^{rd}$ ed. Cambridge, MA: Academic Press, 2008.

[6] J. Strikwerda, Finite Difference Schemes and Partial Differential Equations, 2nd ed. Philadelphia: SIAM, 2004.

[7] J. R. Pierce, An Introduction to Information Theory: Symbols, Signalsand Noise, 2nd ed. New York: Dover Publications, 1980.

[8] A. Kiran Güçoğlu, The Solution of Some Differential Equations by Nonstandard Finite Difference Method, MS Dissertation, İzmir Institute of Technology, Türkiye, 2005.

[9] B. Engquist and A. Majda, “Absorbing boundary conditions for the numerical simulation of waves,” Mathematics of Computation, vol. 31, pp. 629–651, 1977.

[10] G. Mur, “Absorbing boundary conditions for the finite- difference approximation of the time domain electromagnetic field equations,” IEEE Transactions on Electromagnetic Compatibility, vol. 23, pp. 377–382, 1981.

[11] J. A. Stratton, Electromagnetic Theory. New York: McGraw-Hill Book Company, 1941.

[12] P. W. Barber and S. C. Hill, Light Scattering by Particles: Computational Methods, World Scientific, Singapore, 1990.

BIOGRAPHIES

James B. Cole graduated from the University of Maryland, PhD particle physics. After a post-doctoral fellowship (US National Research Council) at the NASA-Goddard Space Flight Center, he was a research physicist at several US National Laboratories, before joining the faculty of University of Tsukuba (Japan).

After returning to the US, he was a senior research fellow of the National Academy of Sciences, and is now a corporate research physicist. He specializes in mathematical models and high precision algorithm development for applications to computational optics and photonics, quantum mechanics, and machine learning. He is one of the pioneers of the methodology of nonstandard finite differences, and has published numerous papers and a book on the subject.

ACES JOURNAL, Vol. 39, No. 3, 183–188
doi: 10.13052/2024.ACES.J.390303
© 2024 River Publishers

ABSTRACT

I. INTRODUCTION

II. NONLINEAR OPTICS

III. FINITE DIFFERENCE MODELS FOR THE LINEAR WAVE EQUATION

A. Notation and definitions

B. Standard finite difference model

C. Nonstandard finite difference model

D. Numerical stability and accuracy

E. Wave equation with a source

F. Multi-frequency NS-FDTD

G. Nonstandard model of refractive index

IV. FINITE DIFFERENCE MODELS OF THE NONLINEAR WAVE EQUATION

A. Standard finite difference model

B. Nonstandard finite difference model

C. Computational model

V. RESULTS AND COMPARISON WITH SEMI-ANALYTIC CALCULATION

VI. EXTENSION TO TWO AND THREE DIMENSIONS

VII. SUMMARY AND CONCLUSIONS

REFERENCES

BIOGRAPHIES