Simulation of the nonlinear Kerr and Raman effect with a parallel local time-stepping DGTD solver

Tiancheng Zhang; Yan Peng; Zhou Dai; Huaguang Bao; Zelong Xiao; Xuewen Chen; Dazhi Ding

doi:10.1364/OE.478225

1. Introduction

Maxwell’s equations for nonlinear media are widely used to describe the nonlinear optical phenomena, in which the electric polarization is determined by the electric field [1]. The numerical implementations of nonlinear Maxwell’s equations are of utmost important to obtain the optical responses in the field of optical communication, data processing and other scientific and industrial applications [2–4]. As an example, theoretical analysis has demonstrated that four-wave mixing interactions in specific configurations pave a new way for advanced active photonic systems are targeted for applications in quantum communications and information processing [5]. Furthermore, doping a nonlinear medium in the photonic crystal can generate solitons to perform high-reliability optical communication [6].

To develop the numerical methods for nonlinear optical analysis, the finite difference time-domain (FDTD) method was developed to solve the nonlinear Maxwell’s equations in the time domain for the first time [7]. Traditional implicit FDTD method for nonlinear systems demands to achieve convergence over the entire numerical space at each time step, which limits the modeling of large-scale optical problems [8]. Thus, the FDTD method based on an explicit methodology was proposed to avoid this [8], and was extended for the analysis of second- and third-order optical nonlinearity [9]. However, the explicit leap-frog time stepping scheme for nonlinear systems must satisfy the stability condition, which results in a much smaller time step increment to greatly reduce the simulation efficiency [10]. Recently, an unconditionally stable FDTD algorithm has been introduced for nonlinear optical analysis [11]. With the alternating-direction implicit (ADI) time-splitting scheme, a sparse matrix equation caused by the coupling of the field components is formed and needs to be solved at each time step, which is a challenge to seek an efficient calculation [12]. Crucially, both implicit and explicit FDTD methods suffer from the staircase approximation and numerical dispersion leading to the low-order accuracy, which requires extremely fine mesh elements for spatial discretization [13,14]. To overcome the aforementioned limitations of the FDTD method, the finite element time-domain (FETD) method demonstrates the enormous potential for the flexible capacities of geometric modeling and time integration [14]. However, due to the high computational complexity of the FETD method, the computing scales are greatly limited. Fortunately, by introducing the discontinuous surfaces between the adjacent elements into the FETD method, the discontinuous Galerkin time-domain (DGTD) method can transform the global matrix into a series of weakly coupled submatrices, which can thus perform low computational complexity and highly parallelized computing [15,16]. A novel DGTD method combined with parametric variational principle is used to solve Maxwell’s equations with nonlinear media, which presents greater stability and flexibility compared to the FDTD method [17,18]. In fact, this approach is not suitable for practical and complex nonlinear optics, mainly because it treats the nonlinear constitutive relations as a series of linear complementary ones [18]. Thus, the DGTD method demonstrates great potential for the analysis of nonlinear optical phenomena, providing a more flexible, efficient and accurate approach for multi-scale nonlinear optical problems.

Different from the solution of linear Maxwell’s equations in the DGTD method, the system matrices of nonlinear Maxwell’s equations need to be updated at each time step [19]. Two solution methods are alternatives to deal with such nonlinear systems. One is Newton’s method, which leads to a quadratic nonlinear convergence [20]. However, this approach requires an expensive LU factorization of the Jacobian matrix at each Newton iteration step and destroys the intrinsic local features. Compared to Newton’s method, the iteration method approximates the nonlinear variables with known quantities, which can transform the nonlinear equations into linear ones and iteratively revises the equations [21,22]. Hence, the inherently local characteristic of the DGTD method is preserved allowing for massive parallelization. In contrast, the disadvantage of this approach is that it can only achieve linear convergence, which generally leads to much more time-consuming computations if strong nonlinearities are involved [23].

In this paper, an arbitrary high-order derivatives (ADER)-DGTD method with the local time-stepping (LTS) scheme is proposed for the analysis of third order nonlinear optical phenomena, which can greatly improve the simulation efficiency and has previously been demonstrated to be superior in terms of accuracy and efficiency for nonlinear field-circuit co-simulations [24]. This proposed method eliminates the restriction of the global stability conditions on the linear region and enables the nonlinear and linear subdomains to update themselves with independent time step increments. In addition, with the iteration solution scheme, the mass matrix of the nonlinear optical system preserves the block diagonal form. In this way, the proposed method can be combined with the massive parallelization technique to further enhance the computational efficiency. For the above reasons, numerical evaluations indicate that the proposed method provides a low computational cost tool for practical nonlinear optical analysis.

2. Principle and formulation

2.1 ADER-DGTD method for nonlinear electromagnetism

Maxwell’s curl equations for nonlinear media in the time domain are

(1)$$\left\{ {\begin{array}{c} {\nabla \times \mathbf{H} = \frac{{\partial \mathbf{D}}}{{\partial t}}}\\ {\nabla \times \mathbf{E} ={-} {\mu_0}\frac{{\partial \mathbf{H}}}{{\partial t}},} \end{array}} \right.$$

where $\mathbf{E}$ and $\mathbf{H}$ are the electric and magnetic fields, and $\mu$ represents the permeability. The electric displacement $\mathbf{D}$ can be expressed as linear and nonlinear functions of [14]

(2)$$\mathbf{D} = \varepsilon {{\mathbf E}} + \sum {{\mathbf{P}_{NL}}} \textrm{ = }\varepsilon {{\mathbf E}} + \sum {{\varepsilon _0}({\chi ^{(n)}}{{{\mathbf E}}^n})} ,$$

in which $\varepsilon$ is the linear permittivity for general and dispersive dielectric media, and the detailed process for dispersive media in the DGTD method can be found in our previous work [25]. $\chi$ is the nonlinear susceptibility. n is the order of nonlinearity.

Take the nonlinear Raman model as an example, the formulation for the nonlinear electromagnetism is described in detail. With the Born-Oppenheimer approximation, the polarization caused by the Raman effect can be expressed as [26]

(3)$${\mathbf{P}_{R\textrm{ama}n}}(t) = {\varepsilon _0}\mathbf{E}(t)[\chi _{Raman}^{(3)}(t) \ast {\mathbf{E}^2}(t)] = {\varepsilon _0}\mathbf{E}(t)\mathbf{S}(t),$$

where ${\ast} $ represents the convolution integral and $\mathbf{S}(t)$ is the auxiliary variable defined as

(4)$$\mathbf{S}(t)\textrm{ = }\chi _{Raman}^{(3)}(t) \ast {\mathbf{E}^2}(t).$$

Thus, the Fourier transform of Eq. (4) is

(5)$$\mathbf{S}(\omega ) = \chi _{Raman}^{(3)}(\omega ) \cdot \mathrm{{\cal F}}[{{\mathbf{E}^2}(t)} ],$$

where $\mathrm{{\cal F}}$ represents the Fourier transform and

(6)$$\chi _{Raman}^{(3)}(\omega ) = \frac{{(1 - a)\chi _0^{(3)}\omega _{Raman}^2}}{{\omega _{Raman}^2 + 2j\omega {\delta _{Raman}} - {\omega ^2}}},$$

with

(7)$$\left\{ {\begin{array}{c} {{\omega_{Raman}} = \sqrt {\frac{{\tau_1^2 + \tau_2^2}}{{\tau_1^2\tau_2^2}}} }\\ {{\delta_{Raman}} = \frac{1}{{{\tau_2}}},} \end{array}} \right.$$

here, $\alpha$, $\chi _0^{}$ and $\tau $ are parameters of the third-order nonlinearity.

Therefore, the auxiliary equation can be rewritten as

(8)$$\omega _{Raman}^2\mathbf{S}(\omega ) + 2j\omega {\delta _{Raman}}\mathbf{S}(\omega ) - {\omega ^2}\mathbf{S}(\omega ) = (1 - a)\chi _0^{(3)}\omega _{Raman}^2\mathrm{{\cal F}}[{{\mathbf{E}^2}(t)} ].$$

Transform the equation to the time domain by the relation $j\omega \leftrightarrow {\partial / {\partial t}}$, and employ the central difference scheme to discretize the time derivative, the iteration format to update $\mathbf{S}(t)$ is given by

(9)$$\mathbf{S}(t) = \frac{{2 - {{(\Delta t)}^2}\omega _{Raman}^2}}{{{\delta _{Raman}}\Delta t + 1}}\mathbf{S}(t - \Delta t) + \frac{{{\delta _{Raman}}\Delta t - 1}}{{{\delta _{Raman}}\Delta t + 1}}\mathbf{S}(t - 2\Delta t) + \frac{{(1 - \alpha )\chi _0^{(3)}\omega _{Raman}^2{{(\Delta t)}^2}}}{{{\delta _{Raman}}\Delta t + 1}}{(\mathbf{E}(t - \Delta t))^2},$$

where $\Delta t$ represents the time step increment. For other nonlinear effects including second-harmonic generation, optical parametric oscillation, intensity-dependent refractive index, or even quantum nonlinear effects, only the auxiliary variable $\mathbf{S}(t)$ needs to be updated and then introduced into Maxwell’s equations to realize the simulation of nonlinear optical phenomena.

With the conventional Newton’s method to solve this nonlinear system, the mass matrix formed by the DGTD method no longer performs the block diagonal characteristic, which affects the calculation efficiency. Hence, to improve the efficiency, the iteration method is used to solve Eq. (1) and Eq. (9) to transform the nonlinear equation into a linear one.

Thus, the semidiscretized form of the DGTD method for nonlinear electromagnetism [24] is

(10)$$\frac{{\partial \mathbf{u}}}{{\partial t}} = {\mathbf{Q}^{ - 1}}{\mathbf{L}_1}\mathbf{u} + {\mathbf{Q}^{ - 1}}{\mathbf{L}_2}{\mathbf{u}^ + },$$

where $\mathbf{u} = {[{e,h} ]^T}$ is the unknown vector. The superscript “+” implies variables from the adjacent element.

(11)$$\mathbf{Q} = \left[ {\begin{array}{cc} {{\mathbf{T}_{ee}}}&0\\ 0&{{\mathbf{T}_{hh}}} \end{array}} \right],$$

(12)$${\mathbf{L}_1} = \left[ {\begin{array}{cc} { - \mathbf{S}{\mathbf{s}_{ee}}}&{{\mathbf{P}_{eh}} + {\mathbf{S}_{eh}}}\\ { - {\mathbf{P}_{he}} - {\mathbf{S}_{he}}}&{ - \mathbf{S}{\mathbf{s}_{hh}}} \end{array}} \right],$$

(13)$${\mathbf{L}_2} = \left[ {\begin{array}{cc} {\mathbf{Ss}_{ee}^ + }&{\mathbf{S}_{eh}^ + }\\ { - \mathbf{S}_{he}^ + }&{\mathbf{Ss}_{hh}^ + } \end{array}} \right].$$

The detailed definitions of the system matrices are given in the Appendix. Different from linear electromagnetism, the mass matrix ${\mathbf{T}_{ee}}$ of nonlinear media varies as

(14)$${[{\mathbf{T}_{ee}}]_{ij}} = (\varepsilon \textrm{ + }{\varepsilon _0}\mathbf{S}(t))\int_V {\mathbf{N}_i^e \cdot \mathbf{N}_j^e} dV,$$

where $\mathbf{N}$ is the hexahedral hierarchical vector basis function [27]. i and j are the number of the subdomains.

Due to the localized interactions of basis functions [15], the global mass matrix generated by the ADER-DGTD method is a block diagonal form, which means that the inverse of the mass matrix of each element can be obtained independently. Therefore, the proposed method can be highly paralleled, which is a significant improvement over other solutions for nonlinear electromagnetism.

Finally, with the recursive derivations of Eq. (10), the variables $\mathbf{u}$ can be obtained with the arbitrary order by expanding in terms of a Taylor series. Thus, the fully discrete ADER-DGTD system with arbitrary order can be formulated as

(15)$$\mathbf{u}(t + \Delta t) = \mathbf{u}(t) + \sum\limits_{i = 1}^{{N_t}} {\frac{{\Delta {t^i}}}{{i!}}} \frac{{{\partial ^i}\mathbf{u}(t)}}{{\partial {t^i}}}\textrm{ = }\mathbf{u}(t) + \sum\limits_{i = 1}^{{N_t}} {\frac{{\Delta {t^i}}}{{i!}}} {({{\mathbf{Q}^{ - 1}}{\mathbf{L}_1}} )^i}\mathbf{u}(t) + \sum\limits_{i = 1}^{{N_t}} {\frac{{\Delta {t^i}}}{{i!}}} {({{\mathbf{Q}^{ - 1}}{\mathbf{L}_2}} )^i}{\mathbf{u}^ + }(t),$$

where ${N_t}$ is the order of ADER scheme.

2.2 Implementation of LTS scheme

The stability criterion of the ADER-DGTD method for a linear system is determined by the characteristic size of the smallest element, which is defined as [28]

(16)$$\Delta {t_l} < \frac{1}{{2M + 1}} \cdot \frac{\rho }{{{c_0}}},$$

where $\Delta {t_l}$ is the time step size of linear electromagnetic region. M is the order of basis function used in the DGTD method. ${c_0}$ is the speed of light and $\rho$ represents the radius of the insphere of the discretized element. Due to the nonlinear effect, the time step increment should be constrained by its nonlinearity and has to be set much smaller to achieve desired accurate results. Thus, the nonlinear optical subdomains should be updated with a small time step increment [13,14,26,29] as $\Delta {t_s} = \Delta {t_l}/N$.

With the traditional DGTD method, the selection of the time step increment must satisfy the global stability condition, which makes the calculation time prohibitively long. Hence, to enhance the efficiency, the LTS scheme is introduced to make the nonlinear and linear subdomains update themselves independently to be viable. The interpolation approach in Eq. (17) is adopted to exchange the variables between these time intervals.

(17)$$\frac{{{\partial ^p}{\mathbf{u}^ + }(t + (n - 1)\Delta {t_s})}}{{\partial {t^p}}} = \frac{{{\partial ^p}{\mathbf{u}^ + }(t)}}{{\partial {t^p}}} + \sum\limits_{j = 1}^{{N_t} - p} {\frac{{{{((n - 1)\Delta {t_s})}^j}}}{{j!}}\frac{{{\partial ^{j + p}}{\mathbf{u}^ + }(t)}}{{\partial {t^{j + p}}}}} .$$

Consequently, the time cost in linear subdomains will be significantly decreased, which greatly improves the analysis efficiency.

The following steps are performed for the ADER-DGTD method with the LTS scheme for nonlinear optical analysis.

1) At time t, calculate (10) and (15) for linear subdomains with $\mathbf{S}(t)\textrm{ = }0$.
2) From $t + n\Delta {t_s}$, $n = 0,1,\ldots ,N - 1$.
- a) Calculate Eq. (9) to obtain the auxiliary variable.
- b) Update the variables in nonlinear subdomains with (10) and (15).
- c) Exchange the variables from linear subdomains to nonlinear ones with (17).
- d) $n = n + 1$.
3) $t = t + \Delta {t_l}$, move to step 1) and 2). Break out of the cycle until the total simulation time is reached.

3. Numerical results and discussions

In this section, the accuracy and efficiency of the proposed method are confirmed by several numerical examples.

3.1 Efficiency comparison

To demonstrate the advantages of the proposed method, the wave propagation characteristic in 1D nonlinear Kerr media [18] shown in Fig. 1 is analyzed. The thickness of the Kerr slab is 20 cm, and the nonlinear permittivity is defined as ${\varepsilon _r}(\mathbf{E}) = 2 + 0.15{|\mathbf{E} |^2}$. The periodic boundary conditions (PBC) [27] are adopted in the x and y directions to terminate the computational domain for infinite structure and the Mur’s first-order absorbing boundary condition (ABC) is placed in the z-direction. An electromagnetic pulse illuminates the target with the form of the first derivative of the Blackman-Harris window function, whose central frequency is 0.4 GHz. The time-domain electric field is observed at 1 m away from the slab. The total computational domain for the DGTD method has 192 elements, which contains 40 and 152 subdomains for nonlinear media and linear free space with a mesh size of 2 cm and 5 cm, respectively. For comparison, the FDTD method is used with a mesh size of 1 cm to discretize 840 elements.

Fig. 1. Geometry of a nonlinear infinite slab.

Download Full Size | PDF

To satisfy the stability condition [27], the electromagnetic field in whole computational domain is updated with the time step increment of 2.26ps for the traditional DGTD method with ADER scheme and the FDTD method. While with the LTS scheme, the time step increment does not need to be limited by the global stability condition for the linear subdomains. Thus, the linear and nonlinear subdomains update the unknowns with 2.26ps and 11.3ps, respectively. The ratio of the time step size between the linear and nonlinear subdomains reaches $N = 5$. According to the above settings, the comparisons of the receiver’s transient electric field between the proposed method, the traditional DGTD method and the FDTD method are exhibited in Fig. 2. The absolute error defined as the absolute vaule of the differences between the proposed method and the FDTD method is exhibited in Fig. 3(a). A great agreement has achieved between these methods, which demonstrates the accuracy of the proposed method. Moreover, the convergence rates of the proposed method and the FDTD method are estimated with different space discretization sizes. The relative error is evaluated with the following absolute loss function:

(18)$${Relative Error = }\frac{{\int_T {|{{\mathbf{E}_{FDTD/DGTD}}(\mathbf{r},t) - \mathbf{E}_{DGTD}^\alpha (\mathbf{r},t)} |dt} }}{{\int_T {|{\mathbf{E}_{DGTD}^\alpha (\mathbf{r},t)} |dt} }}$$

where T is the total simulated observation time. ${\mathbf{E}_{FDTD/DGTD}}(\mathbf{r},t)$ represents the transient electric field calculated by the FDTD method and the proposed method at a probe $\mathbf{r}$. The superscript $\alpha$ is the selected discretization step to obtain the convergent solution. In this numerical simulation, $\alpha = 0.2\textrm{cm}$. Intuitively, the accuracy of the LTS-DGTD method is much higher than that of the FDTD method. This is benefit from the high-order basis function [30] and the ADER time-stepping scheme [24], which lead to the high accuracy of the DGTD method in both spatial and temporal domains [31]. Furthermore, demonstrated in Fig. 3(b), the traditional FDTD method and DGTD method exhibit nearly first-order and second-order convergence respectively, which is consistent with the reported conclusions [32–34]. Thus, the number of elements is greatly reduced, which is exhibited in Table 1. Therefore, even without the LTS scheme, the computational efficiency of the DGTD method is similar to that of the FDTD method. More importantly, the computational efficiency can be significantly improved by reducing the simulation time in the linear region with the proposed method, which is confirmed to be a much more efficient solution for nonlinear electromagnetism analysis.

Fig. 2. Comparison of transmission electric field between the proposed method and the FDTD method.

Download Full Size | PDF

Fig. 3. (a) Absolute error between the proposed method and FDTD method, and (b) estimated convergence rates of these two method.

Download Full Size | PDF

Table 1. Simulation times with different methods

View Table | View all tables in this article

3.2 Third harmonics in nonlinear media

The optical propagation characteristic of a nonlinear infinite slab for third-harmonic generation is simulated to demonstrate the accuracy and efficiency of the proposed method. The blue region with the thickness of 2.4$\mu m$ represents the nonlinear media with the electromagnetic parameter described by the linear Lorentz model and nonlinear Kerr and Raman models, which is the same as in Refs. [11] and [26]. A modulated Gaussian plane wave with the frequency from 150 THz to 550 THz is illuminating the slab along the opposite direction of the z-axis with the electric field linearly polarized along the x-direction, which is defined as

(19)$$\mathbf{E}\textrm{ = }\exp [\frac{{ - {{(t - {t_0})}^2}}}{{{\tau ^2}}}] \times (\cos (2\pi {f_1}t) + \cos (2\pi {f_2}t)),$$

where ${f_1}\textrm{ = }300\textrm{THz}$, ${f_2}\textrm{ = 4}00\textrm{THz}$, ${t_0}\textrm{ = }86.66fs$ and $\tau \textrm{ = }28.89fs$.

Due to the high accuracy and non-conformal technique of the DGTD method, the computational domain does not need to be discretized in such fine size as ${\lambda _0}/75 = 8nm$ in Ref. [11], in which FDTD solution is used to generate conformal mesh elements. In this simulation, the mesh sizes of ${\lambda _0}/20 = 30nm$ and ${\lambda _0}/60 = 10nm$ are adopted in the z-direction of the linear free space and nonlinear media, respectively. The time interval of other subdomains is selected by the stability condition, which equals to $\Delta {t_l} = 8.0 \times {10^{ - 2}}fs$. To ensure the accuracy and stability of nonlinear subdomains, the time step size for nonlinear media is much smaller than the linear subdomains, and the ratio between these time step sizes reaches $\Delta {t_l}/\Delta {t_s}\textrm{ = }8$.

The comparison of the transmitted electric field between the proposed method and the traditional DGTD method is demonstrated in Fig. 4. Moreover, with the discrete Fourier transform method, the power intensity of the transmitted electric field is obtained, which is shown in Fig. 5. A great agreement between the proposed method and the reported results can be observed. Third harmonic phenomenon has been generated by this nonlinear media with the frequencies of $(2{f_1} - {f_2})$ and $(2{f_2} - {f_1})$, which can be used as an optical mixer.

Fig. 4. Comparison of transmission electric field.

Download Full Size | PDF

Fig. 5. Comparison of normalized power intensity.

Download Full Size | PDF

Importantly, the efficiency of the proposed method to analyze the nonlinear optical effect can be greatly improved as exhibited in Table 2. With the traditional DGTD method, the global stability conditions impose constraints on the time step size, which makes all the subdomains must calculate with $\Delta {t_s}$. This limitation significantly increases the computational cost of linear region. By introducing the LTS scheme, the time step increment of linear subdomains is no longer affected by nonlinearity. Hence, the computational cost of linear subdomains is considerably decreased, which achieves a 3.06 times CPU gain with serial computation. What’s more, the proposed method can realize the element-wise solution and perform highly parallelized calculation. This is the largest advantage over the implicit schemes, which makes the system unconditionally stable, but breaks the inherently local characteristics [35,36]. From Table 2, by increasing the number of CPU cores, the simulation time is respectively reduced to 51.2%, 26.8% and 13.6% of that of the DGTD method with the LTS scheme.

Table 2. Simulation times with different methods

View Table | View all tables in this article

4. Conclusion

This paper demonstrates an efficient parallelized solution for the modeling of electromagnetic plus propagation in nonlinear media. The proposed method confirms the competitive calculation accuracy, while requiring far less computational time than the traditional FDTD method, which benefits from the following techniques. First, the discretized subdomains are greatly reduced by adopting the high precision DGTD method with non-conformal modeling technique. Furthermore, the inherently local characteristic of the mass matrix is preserved with the iteration method especially for nonlinear media, which can realize a highly parallelized algorithm. More importantly, combined with the LTS scheme, the time interval of linear electromagnetic subdomains is released from the stability criterion of nonlinear effect to confine the additional computational costs. Numerical evaluations demonstrate that the proposed method can handle the intrinsically high computational cost of nonlinear optical analysis, making this approach more suitable for practical applications. Future extensions of this method are to allow for the analysis of more general nonlinear optical effects.

Appendix

The definitions of the matrices in Eqs. (12)-(14) are exhibited as follows.

{[{\mathbf{T}_{ee}}]_{ij}} = (\varepsilon + \mathbf{S}(t))\int_V {\mathbf{N}_i^e \cdot \mathbf{N}_j^edV}

{[{\mathbf{T}_{hh}}]_{ij}} = \textrm{ - }\mu \int_V {\mathbf{N}_i^h \cdot \mathbf{N}_j^hdV}

{[{\mathbf{P}_{eh}}]_{ij}} = \int_V {\nabla \times \mathbf{N}_i^e \cdot \mathbf{N}_j^hdV}

{[{\mathbf{P}_{he}}]_{ij}} = \int_V {\nabla \times \mathbf{N}_i^h \cdot \mathbf{N}_j^edV}

{[{\mathbf{S}_{eh}}]_{ij}} = \frac{{{Z^ + }}}{{Z + {Z^ + }}}\int_S {\mathbf{N}_i^e \cdot (\hat{n} \times \mathbf{N}_j^h)dS}

{[{\mathbf{S}_{he}}]_{ij}} = \frac{{{Y^ + }}}{{Y + {Y^ + }}}\int_S {\mathbf{N}_i^h \cdot (\hat{n} \times \mathbf{N}_j^e)dS}

{[\mathbf{S}_{eh}^ + ]_{ij}} = \frac{{{Z^ + }}}{{Z + {Z^ + }}}\int_S {\mathbf{N}_i^e \cdot (\hat{n} \times \mathbf{N}_j^{h + })dS}

{[\mathbf{S}_{he}^ + ]_{ij}} = \frac{{{Y^ + }}}{{Y + {Y^ + }}}\int_S {\mathbf{N}_i^h \cdot (\hat{n} \times \mathbf{N}_j^{e + })dS}

{[\mathbf{S}{\mathbf{s}_{ee}}]_{ij}} = \frac{1}{{Z + {Z^ + }}}\int_S {(\hat{n} \times \mathbf{N}_i^e) \cdot (\hat{n} \times \mathbf{N}_i^e)dS}

{[\mathbf{S}{\mathbf{s}_{hh}}]_{ij}} = \frac{1}{{Y + {Y^ + }}}\int_S {(\hat{n} \times \mathbf{N}_i^h) \cdot (\hat{n} \times \mathbf{N}_i^h)dS}

{[\mathbf{Ss}_{ee}^ + ]_{ij}} = \frac{1}{{Z + {Z^ + }}}\int_S {(\hat{n} \times \mathbf{N}_i^e) \cdot (\hat{n} \times \mathbf{N}_i^{e + })dS}

{[\mathbf{Ss}_{hh}^ + ]_{ij}} = \frac{1}{{Y + {Y^ + }}}\int_S {(\hat{n} \times \mathbf{N}_i^h) \cdot (\hat{n} \times \mathbf{N}_i^{h + })dS}

where $Z = \frac{1}{Y}$ represents the wave impedance of the element and superscript “+” implies variables from the adjacent subdomain.

Funding

National Natural Science Foundation of China (62001231, 62025109, 62201257, 62235006); Primary Research and Developement Plan of Jiangsu Province (BE2022070, BE2022070-1); Natural Science Foundation of Jiangsu Province (BK20200467); National Key Laboratory on Electromagnetic Environmental Effects and Electro-optical Engineering (61422062103); Jiangsu Funding Program for Excellent Postdoctoral Talent.

Disclosures

The authors declare no conflicts of interest.

Data availability

No data were generated or analyzed in the presented research.

References

1. R. W. Boyd, Nonlinear Optics (Academic, 1992).

2. B. Kuyken, T. Ideguchi, S. Holzner, M. Yan, T. W. Hansch, J. V. Campenhout, P. Verheyen, S. Coen, F. Leo, R. Baets, G. Roelkens, and N. Picque, “An octave-spanning mid-infrared frequency comb generated in silicon nanophotonic wire waveguide,” Nat. Commun. 6(1), 6310 (2015). [CrossRef]

3. L. Zhang, S. F. Tao, Z. H. Fan, and R. S. Chen, “Efficient method for evaluation of second-harmonic generation by surface integral equation,” Opt. Express 25(23), 28010–28021 (2017). [CrossRef]

4. F. Ye, J. Y. Huang, M. S. Arunna Gandhi, and Q. Li, “Nearly self-similar pulse compression of high-repetition-rate pulse trains in tapered silicon waveguides,” J. Lightwave Technol. 39(14), 4717–4724 (2021). [CrossRef]

5. J. W. You, Z. H. Lan, and N. C. Panoiu, “Four-wave mixing of topological edge plasmons in graphene metasurfaces,” Sci. Adv. 6(13), 3910 (2020). [CrossRef]

6. M. Soljačić and J. D. Joannopoulos, “Enhancement of nonlinear effects using photonic crystals,” Nat. Mater. 3(4), 211–219 (2004). [CrossRef]

7. R. M. Joseph and A. Taflove, “FDTD Maxwell’s equations models for nonlinear electrodynamics and optics,” IEEE Trans. Antennas and Propag. 45(3), 364–374 (1997). [CrossRef]

8. C. Varin, R. Emms, G. Bart, T. Fennel, and T. Brabec, “Nonlinear Lorentz Model for Explicit Integration of Optical Nonlinearity in FDTD,” in 2020 International Applied Computational Electromagnetics Society Symposium (ACES), Monterey, USA (2020), pp. 1–2.

9. C. Varin, R. Emms, G. Bart, T. Fennel, and T. Brabec, “Explicit formulation of second and third order optical nonlinearity in the FDTD framework,” Comput. Phys. Commun. 222, 70–83 (2018). [CrossRef]

10. R. W. Ziolkowski and J. B. Judkins, “Applications of the nonlinear finite difference time domain (NL-FDTD) methodto pulse propagation in nonlinear media: Self-focusing and linear-nonlinear interfaces,” Radio Sci. 28(5), 901–911 (1993). [CrossRef]

11. M. Moradi, S. M. Pourangha, V. Nayyeri, M. Soleimani, and O. M. Ramahi, “Unconditionally stable FDTD algorithm for 3-D electromagnetic simulation of nonlinear media,” Opt. Express 27(10), 15018–15031 (2019). [CrossRef]

12. H. Bao and R. Chen, “An efficient domain decomposition parallel scheme for leapfrog ADI-FDTD method,” IEEE Trans. Antennas Propagat. 65(3), 1490–1494 (2017). [CrossRef]

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

14. I. S. Maksymov, A. A. Sukhorukov, A. V. Lavrinenko, and Y. S. Kivshar, “Comparative study of FDTD-adopted numerical algorithms for Kerr nonlinearities,” Antennas Wirel. Propag. Lett. 10, 143–146 (2011). [CrossRef]

15. J. Chen and Q. H. Liu, “Discontinuous Galerkin time-domain methods for multiscale electromagnetic simulations: A review,” Proc. IEEE 101(2), 242–254 (2013). [CrossRef]

16. S. Yan, C. P. Lin, R. R. Arslanbekov, V. I. Kolobov, and J. M. Jin, “A discontinuous Galerkin time-domain method with dynamically adaptive Cartesian mesh for computational electromagnetics,” IEEE Trans. Antennas Propagat. 65(6), 3122–3133 (2017). [CrossRef]

17. S. B. Zeng, J. F. Chen, B. Zhu, and Q. Ren, “A Nonlinear DG-FETD Scheme Based on Parametric Variational Principle,” in 2019 IEEE International Conference on Computational Electromagnetics (ICCEM), Shanghai, China (2019), pp. 1–3.

18. B. Zhu, H. W. Yang, and J. F. Chen, “A novel finite element time domain method for nonlinear Maxwell’s equations based on the parametric quadratic programming method,” Microw. Opt. Technol. Lett. 57(7), 1640–1645 (2015). [CrossRef]

19. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes: The Art of Scientific Computing,3rd ed (Cambridge University Press, 2007).

20. S. Yan, J. M. Jin, C. F. Wang, and J. Kotulski, “Numerical study of a time-domain finite element method for nonlinear magnetic problems in three dimensions,” Prog. Electromagn. Res. 153(11), 69–91 (2015). [CrossRef]

21. T. C. Zhang, C. Liu, H. G. Bao, D. Z. Ding, and R. S. Chen, “Simulation of Electromagnetic Wave Propagation in Plasma Under High-Power Microwave Illumination,” Antennas Wirel. Propag. Lett. 21(1), 144–148 (2022). [CrossRef]

22. F. I. Hantila, G. Preda, and M. Vasiliu, “Polarization method for static fields,” IEEE Trans. Magn. 36(4), 672–675 (2000). [CrossRef]

23. S. Yan and J. M. Jin, “Three-dimensional time-domain finite element simulation of dielectric breakdown based on nonlinear conductivity model,” IEEE Trans. Antennas Propagat. 64(7), 3018–3026 (2016). [CrossRef]

24. T. C. Zhang, H. G. Bao, P. F. Gu, D. Z. Ding, D. H. Werner, and R. S. Chen, “An Arbitrary High Order DGTD Method with Local Time-Stepping for Nonlinear Field-Circuit Co-Simulation,” IEEE Trans. Antennas Propagat. 70(1), 526–535 (2022). [CrossRef]

25. T. Zhang, H. Bao, D. Ding, and R. Chen, “Interior penalty DGTD method for solving wave equation in dispersive media described with GDM model,” IEEE Trans. Antennas Propagat. 69(9), 6105–6110 (2021). [CrossRef]

26. 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(2), 175–182 (2004). [CrossRef]

27. H. G. Bao, L. Kang, S. D. Campbell, and D. H. Werner, “PML implementation in a nonconforming mixed-element DGTD method for periodic structure analysis,” IEEE Trans. Antennas Propagat. 67(11), 6979–6988 (2019). [CrossRef]

28. A. Taube, M. Dumbser, C. D. Munz, and R. Schneider, “A high-order discontinuous Galerkin method with time-accurate local time stepping for the Maxwell equations,” Int. J. Numer. Model. 22(1), 77–103 (2009). [CrossRef]

29. D. Li and C. D. Sarris, “Time-domain modeling of nonlinear optical structures with extended stability FDTD schemes,” J. Lightwave Technol. 29(7), 1003–1010 (2011). [CrossRef]

30. Q. Ren, H. G. Bao, S. D. Campbell, L. J. Prokopeva, A. V. Kildishev, and D. H. Werner, “Continuous-discontinuous Galerkin time domain (CDGTD) method with generalized dispersive material (GDM) model for computational photonics,” Opt. Express 26(22), 29005 (2018). [CrossRef]

31. C. Y. Tian, Y. Shi, and C. H. Chan, “Interior penalty discontinuous Galerkin time-domain method based on wave equation for 3-D electromagnetic modeling,” IEEE Trans. Antennas Propagat. 65(12), 7174–7184 (2017). [CrossRef]

32. M. Fujii, D. Lukashevich, I. Sakagami, and P. Russer, “Convergence of FDTD and Wavelet-Collocation modeling of curved dielectric interface with the effective dielectric constant technique,” IEEE Microw. Wireless Compon. Lett. 13(11), 469–471 (2003). [CrossRef]

33. A. C. Lesina, A. Vaccari, P. Berini, and L. Ramunno, “On the convergence and accuracy of the FDTD method for nanoplasmonics,” Opt. Express 23(8), 10481–10497 (2015). [CrossRef]

34. L. D. Angulo, J. Alvarez, M. F. Pantoja, S. G. Garcia, and A. R. Bretones, “Discontinuous Galerkin time domain methods in computational electrodynamics: State of the art,” Forum Electromagn. Res. Methods Appl. Technol. 10, 1–24 (2015).

35. V. Dolean, H. Fahs, L. Fezoui, and S. Lanteri, “Locally implicit discontinuous Galerkin method for time domain electromagnetics,” J. Comput. Phys. 229(2), 512–526 (2010). [CrossRef]

36. Q. T. Sun, Q. W. Zhan, Q. Ren, and Q. H. Liu, “Wave equation-based implicit subdomain DGTD method for modeling of electrically small problems,” IEEE Trans. Microwave Theory Tech. 65(4), 1111–1119 (2017). [CrossRef]

Method	CPU Time (sec)		Parallel Efficiency (%)
DGTD	10283.74		/
LTS-DGTD	3356.39		100
Parallel LTS-DGTD	2 cores	1720.46	97.5
	4 cores	900.37	93.2
	8 cores	458.48	91.5

Method	CPU Time (sec)		Parallel Efficiency (%)
DGTD	10283.74		/
LTS-DGTD	3356.39		100
Parallel LTS-DGTD	2 cores	1720.46	97.5
	4 cores	900.37	93.2
	8 cores	458.48	91.5

Simulation of the nonlinear Kerr and Raman effect with a parallel local time-stepping DGTD solver

Abstract

1. Introduction

2. Principle and formulation

2.1 ADER-DGTD method for nonlinear electromagnetism

2.2 Implementation of LTS scheme

3. Numerical results and discussions

3.1 Efficiency comparison

3.2 Third harmonics in nonlinear media

4. Conclusion

Appendix

Funding

Disclosures

Data availability

References

Data availability

Cited By

Figures (5)

Tables (2)

Equations (31)

Optics Express

Method	Elements	Time step size (fs)		Computation Time (s)
FDTD	840	2.26		102.422
DGTD	192	2.26		114.599
LTS-DGTD	192	2.26	11.3	28.946