High Spatial Order Energy Stable FDTD Methods for Maxwell’s Equations in Nonlinear Optical Media in One Dimension

Abstract

In this paper, we consider electromagnetic (EM) wave propagation in nonlinear optical media in one spatial dimension. We model the EM wave propagation by the time-dependent Maxwell’s equations coupled with a system of nonlinear ordinary differential equations (ODEs) for the response of the medium to the EM waves. The nonlinearity in the ODEs describes the instantaneous electronic Kerr response and the residual Raman molecular vibrational response. The ODEs also include the single resonance linear Lorentz dispersion. For such model, we will design and analyze fully discrete finite difference time domain (FDTD) methods that have arbitrary (even) order in space and second order in time. It is challenging to achieve provable stability for fully discrete methods, and this depends on the choices of temporal discretizations of the nonlinear terms. In Bokil et al. (J Comput Phys 350:420–452, 2017), we proposed novel modifications of second-order leap-frog and trapezoidal temporal schemes in the context of discontinuous Galerkin methods to discretize the nonlinear terms in this Maxwell model. Here, we continue this work by developing similar time discretizations within the framework of FDTD methods. More specifically, we design fully discrete modified leap-frog FDTD methods which are proved to be stable under appropriate CFL conditions. These method can be viewed as an extension of the Yee-FDTD scheme to this nonlinear Maxwell model. We also design fully discrete trapezoidal FDTD methods which are proved to be unconditionally stable. The performance of the fully discrete FDTD methods are demonstrated through numerical experiments involving kink, antikink waves and third harmonic generation in soliton propagation.

Appendices

Appendix A. 2M Order Spatial Discretizations

In this section, we provide a proof for Theorem 4.1 by following the exposition in [2, 6].

Proof

Following the discussion in Sect. 4, assume that $u,v \in C_{\#}^{2m+3}([0,L])$ with $m\in {\mathbb {N}}$ an integer, and $m \ge 1$, where the subscript $_\#$ is used to indicate periodic boundary conditions. Then, if $v_h \in V_{0,h}$ is a restriction of v to the primal grid and $u_h \in V_{\frac{1}{2},h}$ is a restriction of u to the dual grid, we have the following Taylor expansions [8, p.53]

$$\begin{aligned}&\left( \tilde{{\mathcal {D}}}^{(2)}_{h}u_h\right) _\ell = \frac{\partial u}{\partial x}(x_\ell ) +\sum _{i=1}^m \displaystyle \frac{h^{2i}}{(2i+1)!2^{2i}}\displaystyle \left( \frac{\partial ^{2i+1}u}{\partial x^{2i+1}}\right) (x_\ell ) + {\mathcal {O}}\left( h^{2m+2}\right) , \end{aligned}$$

(A.1)

$$\begin{aligned}&\left( {\mathcal {D}}^{(2)}_{h}v_h\right) _{\ell +\frac{1}{2}} = \frac{\partial v}{\partial x}(x_{\ell +\frac{1}{2}}) +\sum _{i=1}^m \displaystyle \frac{h^{2i}}{(2i+1)!2^{2i}}\displaystyle \left( \frac{\partial ^{2i+1}v}{\partial x^{2i+1}}\right) (x_{\ell +\frac{1}{2}}) + {\mathcal {O}}\left( h^{2m+2}\right) . \end{aligned}$$

(A.2)

By replacing h / 2 by $(2p-1)h/2$ in (A.1) and by inserting the different Taylor expansions obtained from (A.1) into (4.8) with $m = M-1$ we obtain ( [8, p. 53])

$$\begin{aligned} \left( {\tilde{{\mathcal {D}}}^{(2M)}_{h}u_h} \right) _{\ell }= & {} \sum _{p=1}^M \lambda _{2p-1}^{2M} \sum _{i=0}^{M-1}\left[ \displaystyle \left( \frac{(2p-1)h}{2}\right) ^{2i}\displaystyle \frac{1}{(2i+1)!} \displaystyle \left( \frac{\partial ^{2i+1}u}{\partial x^{2i+1}}\right) (x_{\ell })+ {\mathcal {O}}(h^{2M})\right] \nonumber \\= & {} \sum _{i=0}^{M-1}\left[ \displaystyle \frac{h^{2i}}{2^{2i}(2i+1)!} \left( \displaystyle \frac{\partial ^{2i+1}u}{\partial x^{2i+1}}\right) (x_\ell ) \sum _{p=1}^M \lambda _{2p-1}^{2M}(2p-1)^{2i}\right] + {\mathcal {O}}(h^{2M}).\nonumber \\ \end{aligned}$$

(A.3)

Requiring $\left( { \tilde{{\mathcal {D}}}^{(2M)}_{h}u } \right) _\ell$ to approximate $\left( \frac{\partial u}{\partial x}\right) (x_\ell )$ with error ${\mathcal {O}}(h^{2M})$ leads to a system of equations in the coefficients $\lambda _i^{2M}$, for $i = 1,2,\ldots , 2M-1$, given by

$$\begin{aligned} \left. \begin{array}{lllllll} &{}\lambda _1^{2M}&{}+\lambda _3^{2M}&{}+\lambda _5^{2M}&{}+\ldots &{}+\lambda _{2M-1}^{2M} &{}= 1\\ &{}\lambda _1^{2M}&{}+3^2\lambda _3^{2M}&{}+5^2\lambda _5^{2M}&{}+\ldots &{}+(2M-1)^2\lambda _{2M-1}^{2M} &{}= 0\\ &{}\lambda _1^{2M}&{}+3^4\lambda _3^{2M}&{}+5^4\lambda _5^{2M}&{}+\ldots &{}+(2M-1)^4\lambda _{2M-1}^{2M} &{}= 0\\ &{}\vdots &{}\vdots &{}\vdots &{}\vdots &{}\vdots &{}\vdots \\ &{}\lambda _1^{2M}&{}+3^{2M-2}\lambda _3^{2M}&{}+5^{2M-2}\lambda _5^{2M}&{}+\ldots &{}+(2M-1)^{2M-2}\lambda _{2M-1}^{2M} &{}= 0 \end{array}\right. . \end{aligned}$$

(A.4)

To solve system (A.4) and derive explicit formulas for the coefficients $\lambda _i^{2M}$, we rewrite system (A.4) in matrix form as

$$\begin{aligned} \left( \begin{array}{ccccc} 1^0 &{} 3^0 &{} 5^0 &{} \ldots &{} (2M-1)^0 \\ 1^2 &{} 3^2 &{} 5^2 &{} \ldots &{} (2M-1)^2 \\ 1^4 &{} 3^4 &{} 5^4 &{} \ldots &{} (2M-1)^4 \\ \vdots \\ 1^{2M-2} &{} 3^{2M-2} &{} 5^{2M-2} &{} \ldots &{} (2M-1)^{2M-2} \end{array}\right) \left( \begin{array}{c} \lambda _1^{2M}\\ \lambda _3^{2M}\\ \lambda _5^{2M}\\ \vdots \\ \lambda _{2M-1}^{2M}\end{array}\right) = \left( \begin{array}{c} 1\\ 0\\ 0\\ \vdots \\ 0 \end{array}\right) . \end{aligned}$$

(A.5)

Let the matrix of system (A.5) be denoted as $W_{2M}$. We define the vector $\lambda ^{2M} = \left( \lambda _1^{2M},\lambda _3^{2M},\lambda _5^{2M}, \ldots , \lambda _{2M-1}^{2M}\right) ^T$. Multiplying the linear system of equations in (A.5) by any vector $U = (U_1,U_2,\ldots ,U_M)^T \in {\mathbb {R}}^M$, we get the equation

$$\begin{aligned} U^TW_{2M}\lambda ^{2M} = U_1. \end{aligned}$$

(A.6)

Let ${\mathcal {P}}:= {\mathcal {P}}_{2M-2}^{\mathrm {even}}({\mathbb {R}})$ be the set of all even polynomials of degree $2M-2$ with real coefficients. Associated to the vector $U \in {\mathbb {R}}^M$, we define $P_U \in {\mathcal {P}}$ as

$$\begin{aligned} P_U(x) = U_1+U_2(2x-1)^2+U_3(2x-1)^4+\ldots +U_M(2x-1)^{2M-2}, \end{aligned}$$

(A.7)

from which we obtain the equation

$$\begin{aligned} U^TW_{2M} = \left( P_U(1),P_U(2),P_U(3),\ldots ,P_U(M)\right) , \end{aligned}$$

(A.8)

which permits rewriting Eq. (A.6) as

$$\begin{aligned} \sum _{j=1}^M P_U(j)\lambda ^{2M}_{2j-1} = P_U\left( \frac{1}{2}\right) . \end{aligned}$$

(A.9)

Satisfying Eq. (A.6) for any $U \in {\mathbb {R}}^M$ is equivalent to having Eq. (A.9) hold for any polynomial $P \in {\mathcal {P}}$.

For each integer $1\le p\le M$, we now consider the polynomials in ${\mathcal {P}}$ defined as

$$\begin{aligned} Q_p(x) = \prod _{1\le r\le M, r\ne p} \left( 1-\displaystyle \frac{(2x-1)^2}{(2r-1)^2}\right) , \end{aligned}$$

(A.10)

The polynomial $Q_p(x)$ vanishes at all $x = 1,2,3 \ldots , M$ except at $x=p$. Using $P = Q_p$ in (A.9) we have

$$\begin{aligned} Q_p(p)\lambda ^{2M}_{2p-1} = Q_p\left( \frac{1}{2}\right) = 1, \end{aligned}$$

which implies that

$$\begin{aligned} \lambda _{2p-1}^{2M} = \displaystyle \frac{1}{Q_p(p)} = \prod _{1\le r\le M, r\ne p}\left( 1-\displaystyle \frac{(2p-1)^2}{(2r-1)^2}\right) ^{-1}. \end{aligned}$$

(A.11)

We use the following identities (given without proof)

$$\begin{aligned} \prod _{1\le r\le M, r\ne p} \left( 1+\displaystyle \frac{(2p-1)}{(2r-1)}\right) ^{-1} = \displaystyle \frac{2^p(p-1)!(2M-1)!!}{(2M+2p-2)!!}, \end{aligned}$$

(A.12)

and

$$\begin{aligned} \prod _{1\le r\le M, r\ne p} \left( 1-\displaystyle \frac{(2p-1)}{(2r-1)}\right) ^{-1} = \displaystyle \frac{(-1)^{p-1}(2M-1)!!}{2^{p-1}(2M-2p)!!(2p-1)(p-1)!}, \end{aligned}$$

(A.13)

where $p \in {\mathbb {Z}}, 1\le p\le M$.

From Eqs. (A.11), (A.12) and (A.13), and some algebraic manipulations, we can obtain the explicit formula (4.10). $\square$

Remark A.1

The result in (4.10) has been obtained, using other techniques, by other authors in the past (see [11, 15], and [12]). In [2], the authors prove several additional properties of the corresponding coefficients for higher order approximations of the 1D Laplace operator. Similar properties for the coefficients $\lambda ^{2M}_{2p-1}$ can be proved. Some of these properties have been proved in [15] and [12].

Appendix B. Symbol of the Operator ${\mathcal {A}}_h$

Following [4], we define $K = \frac{kh}{2}$, and use Chebyshev polynomial identities to get

$$\begin{aligned} \displaystyle \sin ((2m-1)K) = \sum _{j = 1}^m \alpha _j^m \sin ^{2j-1}(K), \end{aligned}$$

(B.1)

where, for $j = 1,2,\ldots ,m$

$$\begin{aligned} \displaystyle \alpha _j^m:= (-1)^{2m-j-1}\frac{2m-1}{m+j-1} \frac{(m+j-1)!}{(m-j)!}\frac{2^{2j-2}}{(2j-1)!}. \end{aligned}$$

(B.2)

Using the identity (B.1) and rearranging terms, we obtain

$$\begin{aligned} \displaystyle \sum _{p=1}^M \frac{\lambda _{2p-1}^{(2M)}}{2p-1} \sin \left( (2p-1) \frac{kh}{2} \right)&= \displaystyle \sum _{p=1}^M \frac{\lambda _{2p-1}^{(2M)}}{2p-1} \sum _{\ell =1}^{p} \alpha _\ell ^p\sin ^{2\ell -1}\left( \frac{kh}{2} \right) \nonumber \\&= \sum _{\ell =1}^{M}\left( \sum _{p = \ell }^{M} \frac{\lambda _{2p-1}^{(2M)}}{2p-1} \alpha _\ell ^p\right) \sin ^{2\ell -1}\left( \frac{kh}{2} \right) . \end{aligned}$$

(B.3)

By definition of $\displaystyle \lambda _{2p-1}^{(2M)}$ and $\displaystyle \alpha _\ell ^p$, we have

$$\begin{aligned} \displaystyle \sum _{p = \ell }^{M} \frac{\lambda _{2p-1}^{(2M)}}{2p-1} \alpha _\ell ^p = \sum _{p = \ell }^{M} \frac{1}{2p-1}\frac{(-1)^{3p-\ell -2}2^{2\ell -1}[(2M-1)!!]^2 (p+\ell -2)!}{(2M+2p-2)!!(2M-2p)!!(p-\ell )!(2\ell -1)!}. \end{aligned}$$

(B.4)

Using the fact that $(2n)!! = 2^n n!$ and changing index from $p = j+ \ell$ to index j, we get

$$\begin{aligned}\displaystyle \sum _{p = \ell }^M \frac{\lambda _{2p-1}^{(2M)}}{2p-1} \alpha _\ell ^p&= \sum _{j = 0}^{M-\ell } \frac{1}{(2j + 2\ell -1)}\frac{(-1)^{3j+2\ell -2}2^{2\ell -1}[(2M-1)!!]^2(j+2 \ell -2)!}{(2M+2j+2\ell -2)!!(2M-2j-2\ell )!!j!(2\ell -1)!} \\&= \sum _{j=0}^{M-\ell } \frac{1}{2j+2\ell -1} \frac{(-1)^j 2^{2\ell -1}[(2M-1)!!]^2(2\ell +j-2)!}{2^{M+j+\ell -1}(M+j+\ell -1)!2^{M-j-\ell }(M-j-\ell )! j!(2\ell -1)!}\\&= \frac{[(2M-1)!!]^2 2^{2\ell }}{2^{2M}(2\ell -1)!} \sum _{j=0}^{M-\ell } \frac{(-1)^j (2\ell +j-2)!}{(2j+2\ell -1)!(M+j+\ell -1)!(M-j-\ell )!j!}. \end{aligned}$$

Using the following recursion formulas, which we state without proof, for $n \in {\mathbb {N}}$,

$$\begin{aligned} \displaystyle \Gamma (n+1) = n \Gamma (n); \qquad \text {with} \qquad \Gamma \left( n + \frac{1}{2} \right) := \frac{(2n-1)!!\sqrt{\pi }}{2^n}, \end{aligned}$$

(B.5)

we then have

$$\begin{aligned} \displaystyle \sum _{p = \ell }^M \frac{\lambda _{2p-1}^{(2M)}}{2p-1} \alpha _\ell ^p = \frac{[(2M-1)!!]^2 2^{2\ell }}{2^{2M}(2\ell -1)!} \frac{[\Gamma (\ell - \frac{1}{2})]^2}{4[\Gamma \left( M + \frac{1}{2}\right) ]^2} = \frac{[(2\ell -3)!!]^2}{(2\ell -1)!}, \end{aligned}$$

(B.6)

and

$$\begin{aligned} \displaystyle \sum _{p=1}^M \frac{\lambda _{2p-1}^{(2M)}}{2p-1} \sin \left( (2p-1) \frac{kh}{2} \right) = \sum _{\ell = 1}^M \frac{[(2\ell -3)!!]^2}{(2\ell -1)!} \sin ^{2\ell - 1}\left( \frac{kh}{2} \right) . \end{aligned}$$

(B.7)