A Time Domain Volume Integral Equation Solver to Analyze Electromagnetic Scattering from Nonlinear Dielectric Objects

TL;DR

This work addresses transient EM scattering from Kerr nonlinear dielectrics in three dimensions by introducing an explicit marching-on-in-time time-domain electric-field volume integral equation (TD-EFVIE) solver. The method couples the TD-EFVIE with a nonlinear constitutive relation $\mathbf{D}=\varepsilon(\mathbf{E})\mathbf{E}$ where $\varepsilon(\mathbf{E})=\varepsilon_{0}[\chi^{(1)}+\chi^{(3)}|\mathbf{E}|^{2}]$, and inverts it via a Padé approximant to update $\mathbf{E}$ from $\mathbf{D}$, assuming piecewise-constant permittivity inside the scatterer. Time integration is performed with a $PE(CE)^m$ predictor-corrector scheme, using Schaubert-Wilton-Glisson (SWG) basis functions for both $\mathbf{E}$ and $\mathbf{D}$, which yields sparse Gram-matrix systems solved iteratively at each step without Newton-like nonlinear solvers. The approach is validated on several 3D nonlinear scattering problems, showing higher accuracy than conventional FDTD for nonlinear objects and capturing nonlinear phenomena such as higher harmonics and four-wave mixing. These results indicate a practical, efficient tool for transient analysis of nonlinear dielectric scatterers with potential extensions to high-contrast media.

Abstract

A time domain electric field volume integral equation (TD-EFVIE) solver is proposed for analyzing electromagnetic scattering from dielectric objects with Kerr nonlinearity. The nonlinear constitutive relation that relates electric flux and electric field induced in the scatterer is used as an auxiliary equation that complements TD-EFVIE. The ordinary differential equation system that arises from TD-EFVIE's Schaubert-Wilton-Glisson (SWG)-based discretization is integrated in time using a predictor-corrector method for the unknown expansion coefficients of the electric field. Matrix systems that arise from the SWG-based discretization of the nonlinear constitutive relation and its inverse obtained using the Pade approximant are used to carry out explicit updates of the electric field and the electric flux expansion coefficients at the predictor and the corrector stages of the time integration method. The resulting explicit marching-on-in-time (MOT) scheme does not call for any Newton-like nonlinear solver and only requires solution of sparse and well-conditioned Gram matrix systems at every step. Numerical results show that the proposed explicit MOT-based TD-EFVIE solver is more accurate than the finite-difference time-domain method that is traditionally used for analyzing transient electromagnetic scattering from nonlinear objects.

Figures (4)

• Figure 1: Scattering from a linear sphere. (a) $x$-component of $\mathbf{E}(\mathbf{r}, t)$ computed by the proposed solver at the center of the sphere [$\mathbf{r}_0=(0,0,0)$]. (b) $\sigma_{\mathrm{Mie}}(\theta)$ and $|\sigma_{\mathrm{Mie}}(\theta)-\sigma_{\mathrm{MOT}}(\theta)|$, where $\sigma_{\mathrm{MOT}}(\theta)$ and $\sigma_{\mathrm{Mie}}(\theta)$ are the RCS computed at $f=5.0\,\mathrm{MHz}$ on the $\phi=0$ plane using the Fourier-transformed time-domain solution and the Mie series solution, respectively.
• Figure 2: Scattering from a nonlinear cube. (a) $x$-component of $\mathbf{D}(\mathbf{r}_0,t)$ computed by the proposed solver (MOT) and the FDTD-based solver MEEP at the center of the cube [$\mathbf{r}_0=(0,0,0)$]. (b) Zoomed version of (a) in the time range $[32.5, 41.5]\,\mathrm{ns}$. (c) Fourier transform of the $x$-component of $\mathbf{D}(\mathbf{r}_0,t)$ computed by the proposed solver (MOT) and the FDTD-based solver MEEP. (d) Convergence in $err$ defined by \ref{['eq:errorcalc']} with increasing mesh density (decreasing average edge length $l_{\mathrm{av}}$).
• Figure 3: Four-wave mixing frequency conversion. (a) The time dependence of the plane wave excitation, $P(t)$ given by \ref{['eq:4wave-mixing']}, and (b) its Fourier transform. (c) $x$-component of $\mathbf{D}(\mathbf{r}_0,t)$ computed by the proposed solver at the center of the sphere in the first (linear sphere with $\chi^{(1)}=1.5$, $\chi^{(3)}=0$) and the second (nonlinear sphere with $\chi^{(1)}=1.5$, $\chi^{(3)}=0.075$) simulations [$\mathbf{r}_0=(0,0,0)$]. (d) Fourier transform of the $x$-component of $\mathbf{D}(\mathbf{r}_0,t)$ computed in the two simulations.
• Figure 4: Transmission through a Bragg grating. (a) Description of the grating geometry. (b) The time dependence of the plane wave excitation, $P(t)$ given by \ref{['eq:nonlinear_bragg']}. (c) $x$-component of $\mathbf{E}(\mathbf{r}, t)$ computed at the feeding end [$\mathbf{r}_0=(0,0,-2.6\,\mu \mathrm{m})$] and the trailing end [$\mathbf{r}_0=(0,0,2.6\,\mu \mathrm{m})$] in the first simulation, where both layers are linear. (d) $x$-component of $\mathbf{E}(\mathbf{r}, t)$ at the trailing end [$\mathbf{r}_0=(0,0,2.6\,\mu \mathrm{m})$] computed by the proposed solver in the first (both layers are linear) and the second (one layer is linear, the other one is nonlinear) simulations.