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

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.


Introduction
Electromagnetic and optical devices often rely on materials that exhibit strong Kerr nonlinearity [1][2][3] to induce various interesting physical phenomena, such as higher-order harmonic generation [4], self-focusing [5], self-phase modulation [6], four-wave mixing [7], and three-state switching [8], in their response to electromagnetic excitation.The dielectric permittivity of these nonlinear materials is mathematically modeled as a power series of the magnitude of the electric field weighted with susceptibility coefficients [1][2][3].
Due to this complicated dependence of the dielectric permittivity on the electric field, design of the nonlinear electromagnetic and optical devices, even when they have simplistic geometries, has to be carried out using an electromagnetic simulator.These simulators often rely on time domain techniques since frequency domain methods operate under the assumption of time-harmonic excitation and can not be used in the presence of strong nonlinearities [1].
Majority of the time domain solvers developed to simulate electromagnetic field interactions on materials with Kerr nonlinearity are based on the finite-difference time-domain (FDTD) method [9][10][11][12][13][14][15][16][17][18][19].This can be explained by the fact that FDTD methods are relatively straightforward to be implemented and they permit easy incorporation of the nonlinear permittivity function through the use of an auxiliary equation [9,11,20].This approach calls for iterative updates between the Maxwell equations and this auxiliary equation, which is often implemented using a Newton-type nonlinear solver or by simple explicit updates of the variables in the equations [9,20].
Time marching schemes that rely on similar approaches have also been adopted into time domain finite element method (TD-FEM) [21][22][23][24][25][26].Unlike FDTD, TD-FEM is not restricted to uniform grids for spatial discretization, and therefore can be used in the simulation of devices with complicated geometries without loss of accuracy and/or efficiency.That said, both FDTD and TD-FEM suffer from several well-known drawbacks of the differential equation solvers: They require the computation domain to be truncated using absorbing boundary conditions or perfectly matched layers, their accuracy is limited by numerical phase dispersion, and their time step size is often restricted by the Courant-Friedrichs-Lewy (CFL) condition [1,9].Time domain volume integral equation (TD-VIE) solvers [27][28][29][30][31][32][33][34][35][36][37][38][39] do not suffer from these drawbacks of FDTD and TD-FEM.This is because they rely on a formulation where the scattered electromagnetic field is represented as a spatio-temporal convolution between the Green function of the background medium and current/field induced in the geometry.However, this convolution operation (in discretized form) is also the reason why computational cost and memory requirements of classical marching-on-in-time (MOT)-based TD-VIE solvers are high [31][32][33][34][35].Furthermore, inaccurate discretization/computation of the retarded-time integrals relevant to this convolution cause late-time instabilities in the solution [36].The former issue has been addressed by the development of plane wave time domain (PWTD) algorithm [31][32][33] and fast Fourier transform (FFT)-based schemes [34,35].Late-time instability issue has been mostly alleviated by using highly-accurate interpolation functions in the temporal discretization [36,[40][41][42][43].
These developments have significantly increased the range of TD-VIE solvers' applicability and enabled their use in transient analysis of electromagnetic scattering from electrically large [33][34][35], lossy [31], dispersive [32], and/or high-contrast [36] scatterers.On the other hand, development of TD-VIE solvers for scatterers with nonlinear permittivity has been limited to only two-dimensional (2D) problems [44].This can be explained by the fact that this 2D method uses an implicit MOT scheme and therefore requires a Newton-like nonlinear solver at every time step, which limits its computational efficiency and applicability to three-dimensional (3D) problems.
In this work, an explicit MOT-based time domain electric field volume integral equation (TD-EFVIE) solver is proposed for transient analysis of electromagnetic scattering from 3D nonlinear dielectric objects with Kerr nonlinearity.The proposed method represents the scattered electric field in the form of a (volumetric) convolution between the time domain Green function of the background medium and two unknowns, namely the total electric field and the total electric flux induced in the scatterer.Then, the time derivative of the fundamental field relation, i.e., the sum of the incident and the scattered fields is equal to the total field, is enforced in the scatterer.This yields TD-EFVIE.The nonlinear constitutive relationship between the two unknowns is used as an auxiliary equation that complements TD-EFVIE.
To numerically solve TD-EFVIE, the scatterer is discretized into a mesh of tetrahedrons.The electric flux and the electric field are spatially expanded using "full" and "half" Schaubert-Wilton-Glisson (SWG) basis functions [31,32,37,45] defined on these tetrahedrons.
Inserting these expansions into TD-EFVIE and testing the resulting equation using half SWG functions yield a system of ordinary differential equations (ODEs) in time-dependent expansion coefficients of the SWG basis functions.A predictor-corrector scheme, more specifically a P E(CE) m scheme, is used to integrate this system of ODEs in time for the unknown electric field expansion coefficients [38,39,43,[46][47][48][49].Here, P E and CE refer to the predictor and corrector stages, and m is the number of corrector updates.Similarly, expansions of the electric field and the electric flux are inserted into the nonlinear constitutive relation and its "inverse" obtained using the Padé approximant [50][51][52].The resulting equations are tested using full SWG functions at discrete time steps.This yields matrix systems that relate the electric flux expansion coefficients to those of the electric field.These matrix equations are used at PE and CE stages of the P E(CE) m scheme during the time integration for the explicit updates of the expansion coefficients.This approach does not call for Newton-like nonlinear solvers as done in the implicit MOT solvers.Even though it requires linear matrix system solutions at every time step, the Gram matrices associated with these systems are always sparse and well-conditioned (independent of the time step size).Therefore, these systems are efficiently solved using an iterative solver.
The accuracy and the applicability of the resulting explicit MOT-based TD-EFVIE solver are demonstrated using several numerical examples.These results clearly show that the proposed method is more accurate than FDTD that is traditionally used for analyzing electromagnetic scattering from nonlinear objects.Note that preliminary versions of the method proposed in this work have been described in [53,54] as conference contributions.
The remainder of this paper is organized as follows.Sections 2.1 and 2.2 present the formulation of TD-EFVIE and the nonlinear constitutive relation and their discretization.Section 2.3 describes the P E(CE) m scheme that is used for the solution of the discretized TD-EFVIE (in the form of ODE systems) and the discretized nonlinear constitutive relation.This is followed by Section 2.4 where several comments about the proposed explicit MOTbased TD-EFVIE solver are provided.Section 3 presents numerical results that demonstrate the accuracy, the stability, and the applicability of the proposed method.Finally, Section 4 summarizes this work and outlines future research directions.

TD-EFVIE
Let V denote the volumetric support of a scatterer that resides in an unbounded background medium with permittivity ε 0 and permeability µ 0 .Electric field E inc (r, t) is incident on the scatterer.It is assumed that E inc (r, t) is vanishingly small for ∀r ∈ V and t ≤ 0, and is essentially band-limited to frequency f max .In response to this excitation, the equivalent volumetric electric current J(r, t) is induced in V , and J(r, t) generates the scattered electric field E sca (r, t).One can express E sca (r, t) in terms of J(r, t) using [1] where ∂ t denotes the derivative with respect to time, R = |r − r | is the distance between the source point r and the observation point r, and c 0 = 1/ √ ε 0 µ 0 is the speed of light in the background medium.In V , J(r, t) is expressed in terms of the total electric field E(r, t) and the electric flux D(r, t) using Here, E(r, t) satisfies the fundamental field relation: Inserting ( 1) and (2) into the time derivative of (3) for r ∈ V yields the time derivative form of TD-EFVIE [in unknowns E(r, t) and D(r, t)] [ [29][30][31][32][34][35][36][37] as Here, L[X](r, t) is the volume integral operator defined as Note that an additional time derivative is applied to (3) to obtain (4) because this equation is in the form of an ODE in time and a P E(CE) m scheme be used to integrate it in time for the unknown E(r, t) [38,39,43,[46][47][48][49].
To discretize (4), first V is divided into a mesh of tetrahedrons, and D(r, t) and E(r, t) are discretized using SWG functions that are defined on triangular patches of this mesh [45].
D(r, t) is approximated using where f D n (r) represents the SWG basis function set, {I D (t)} n = I D n (t) are the unknown timedependent expansion coefficients, and N D is the total number of patches in the tetrahedral mesh.f D n (r) associated with triangular patch S n is defined as [45] Here, V ± n are the tetrahedrons on the two sides of S n , |S n | represents the area of S n , |V ± n | are the volumes of V ± n , and r ± n are the "free" nodes of V ± n , i.e., r Note that if S n is located on the surface of V , there is only one tetrahedron attached to it and f D n (r) is set to the "half" SWG function f + n (r) defined in this single tetrahedron represented by V + n in the first row of (7).Note that this combination of full and half SWG functions used in basis set f D n (r) ensures that the normal component of D(r, t) across any two tetrahedrons in V is continuous and the normal component of D(r, t) on the surface of V is properly accounted for.
On the other hand, E(r, t) should be approximated using a basis set that allows its normal component to be discontinuous across any two tetrahedrons in V .Therefore, E(r, t) is expanded using half SWG basis functions [31,32]: where f E n (r) = f + n (r) are the half SWG basis functions associated with triangular patches (there is one half SWG function for each patch that is located on the surface of V and two for each "internal" patch), I E (t) n = I E n (t) are the unknown time-dependent expansion coefficients, and , where N B is the number of patches located on the surface of V .
Inserting ( 6) and ( 8) into (4) and testing the resulting equation with f E m (r), m = 1, . . ., N E yield the spatially-discretized time-dependent TD-EFVIE as Here, the elements of the Gram matrix G EE , the tested incident field vector V inc (t), and the tested scattered field vectors V sca,E (t) and V sca,D (t) are given by The semi-discretized system in ( 9) is integrated in time using the P E(CE) m scheme described in Section 2.3.This calls for sampling (9) in time with time step size ∆t.To compute the samples of V sca,E (t) and V sca,D (t), the retarded-time integrals and have to be evaluated at discrete times t = j∆t, which consequently means and their temporal derivatives) should be interpolated from the samples of I E (t) and I D (t), respectively.This requires expansion of I E (t) and I D (t) as Here, , where T (t) is the temporal interpolation function, and N t is the number of time steps.Inserting ( 14) and ( 15) into ( 9) and point-testing the resulting equation in time, i.e., sampling it at t = j∆t, j = 1, . . ., N t yield the fully discretized TD-EFVIE as where , and the elements of the matrices Z EE j−i and Z ED j−i are given by One can see from the definition of the SWG function in (7) that the full SWG functions can be expressed as linear combinations of the half SWG functions multiplied by the appropriate signs.In other words, basis set f D n (r) can be constructed by linearly combining functions in basis set f E n (r).This means that Z ED j−i can be expressed in terms of Z EE j−i using where P is the sparse matrix of linear mapping from the basis set f E n (r) to the basis set f D n (r), and its non-zero entries are either 1 or −1.Inserting ( 19) into ( 16) yields the final form of the fully discretized TD-EFVIE as The matrix entries {Z EE j−i } mn in ( 17) can be explicitly written as The order of the singularity in the second double integral in ( 21) is reduced by using the chain rule and the divergence theorem [30,36,45,55], which yields Note that the derivation of the surface integral expression in (22) uses the fact that the normal component of f E m (r) is equal to 1 on S m and 0 on the other three surfaces of V + m .

Nonlinear Constitutive Relation
In addition to TD-EFVIE in ( 4), E(r, t) and D(r, t) for r ∈ V are related to each other via the nonlinear constitutive relation [3,9]: Here, ε(E) is the electric-field dependent permittivity and is expressed as where χ (1) and χ (3) are the linear and the third-order nonlinear coefficients associated with the Kerr nonlinearity, respectively.
The constitutive relation in (23) complements TD-EFVIE in (4).Inserting ( 6) and ( 8) into (23) and testing the resulting equation with Here, the elements of the Gram matrices G DD and G DE j are given by where . Using P, the sparse matrix of linear mapping from the basis set f E n (r) to the basis set f D n (r), one can express G DD in terms of G EE using It is assumed that ε(E) is a piece-wise constant function inside the scatterer, with a constant value in each tetrahedron.This constant value is computed at the center of each tetrahedron.
Let r c n represent the center of V ± n .Inserting (8) into (24) and evaluating the resulting expression at r = r c n and t = j∆t yield where the index l runs over the indices of the half basis functions defined in V + n (there are four of them).Let S E j represent a diagonal matrix with entries Then, G DE j can be expressed in terms of G EE as Inserting ( 28) and ( 31) into (25) yields the final form of the discretized constitutive relation as The P E(CE) m scheme described in Section 2.3 requires I E j to be updated from I D j (during evaluation steps).This calls for the "inversion" of the nonlinear constitutive relation in (23).This is done using the Padé approximant [50][51][52]: Inserting ( 6) and ( 8) into (33) and testing the resulting equation with f E m (r), m = 1, . . ., N E at t = j∆t, j = 1, . . ., N t yield Here, the elements of the Gram matrix G ED j are given by Just like ε(E), it is assumed that ε(D) in ( 33) is a piece-wise constant function inside the scatterer, with a constant value in each tetrahedron.Inserting (6) into the expression of ε(D) [see (33)] and evaluating the resulting expression at r = r c n and t = j∆t yield where the index l runs over the indices of the basis functions that have V + n or V − n as support.
Let S D j represent a diagonal matrix with entries Then, G ED j can be expressed in terms of G EE using Inserting (38) into (34) and eliminating G EE from both sides of the resulting equation yield the final form of the discretized Padé approximant as:

PE(CE)m Scheme
The fully-discretized TD-EFVIE (20) relates unknowns I E j and I D j to the time derivative of the unknown İE j , and is integrated in time using a P E(CE) m scheme to yield the unknown [39,[46][47][48].This scheme uses the discretized constitutive relation (32) and the discretized Padé approximant (39) to update I E j and I D j .The steps of the P E(CE) m are provided as follows.
Step 0: Compute V fix j , the part of the right-hand side of (20) that does not change within the time step j:
Step 7.7: Compute İE,(m) j by solving Step 7.8: Check convergence where PECE is the convergence threshold, and x represents the L 2 -norm of vector x.
End loop over m.
End loop over j.
In the P E(CE) m scheme described above, p and c are the predictor and corrector coefficient vectors of length 2k and 2k +1, respectively [43,49].SOR in (48) Here, ITS is the convergence threshold, b is the right-hand side vector, and I (n) j and I (n−1) j are the solutions at iterations n and n − 1, respectively.

Comments
Several comments about the formulation and the discretization schemes described in Sections 2.1 and 2.2, and the P E(CE) m scheme described in Section 2.3 are in order: 1) The second-order nonlinear term (with coefficient χ (2) ) is not considered in the expression of ε(E) given by ( 24) because it is assumed that the scatterer is centrosymmetric [3].
But the method developed in this work could still be used for permittivity functions with second-and higher-order terms.
2) The most common choices for the temporal interpolation function T i (t) used in (14) and ( 15) are the Lagrange polynomials [56][57][58] and the band limited approximate prolate spheroidal wave (APSW) functions [36,[40][41][42].Both options can be used by the MOT solver described in this work.However, when APSW functions are used, the resulting P E(CE) m scheme is no longer casual, i.e., "future" samples of I E i and I D i are required to compute the summation on the right-hand side of (40).In this case, the causality of the time marching is restored using the extrapolation scheme described in [36].This extrapolation scheme is specifically tailored for the accurate and stable solution of TD-EFVIE.
3) Similarly, several options exist to generate the predictor and corrector coefficients, p and c.They can be obtained using polynomial interpolation, which leads to well-known linear multistep methods such as Adams-Moulton, Adams-Bashforth, and backward difference schemes [49].One can also obtain them numerically under the assumption that oscillating and/or decaying exponential functions can be used to approximate the time-dependence of the solution [43].
4) When MOT-based TD-EFVIE solvers are used to analyze electromagnetic scattering from linear objects, ∆t is selected using ∆t = 1/ (2γf max ) where γ is the oversampling factor.
Ideally, γ can be set to 1.0 due to the Nyquist sampling criterion, but, in general, it is set to a value between 2.5 and 15.0 (depending on the desired level of accuracy).On the other hand, for nonlinear scatterers, there is no explicit guidance criterion to select ∆t.For the numerical examples presented in Section 3, γ is increased to sufficiently resolve the frequency of the higher harmonics.6) Within the framework of the proposed explicit MOT solver, a more accurate expansion of E(r, t) can be used.For example, the half SWG functions can be replaced by the fully linear curl-conforming basis functions as done in [39].These functions are defined on the edges of the tetrahedrons and automatically enforce the tangential continuity of E(r, t), leading to a more accurate solution.Switching to this type of basis functions would not change how the nonlinearity is accounted for and the time integration is carried out using a P E(CE) m scheme during time marching.Having said that, using SWG and half SWG functions to expand D(r, t) and E(r, t), respectively, reduces the computational cost of the solver.This is because full SWG functions can be constructed by linearly combining half SWG functions leading to the sparse mappings between the matrices and the expansion coefficients as shown in ( 19), ( 28), (31), and (38).

Numerical Results
This section presents several numerical examples that demonstrate the accuracy, stability, and applicability of the proposed explicit MOT solver in characterizing electromagnetic field interactions on scatterers with Kerr nonlinearity.In all the examples considered here, the scatterer resides in free space with permittivity ε 0 and permeability µ 0 .In all simulations, the excitation is a plane wave with electric field where E 0 is the amplitude of the electric field, and p = x and k = ẑ are the unit vectors that represent the directions of the electric field and the propagation of the plane wave, respectively.In (56), P (t) is a band-limited pulse that describes the time dependence of the excitation.

Linear Sphere
In the first example, electromagnetic scattering from a "linear" sphere is analyzed using the proposed method.The sphere is centered at the origin and has a radius of length 1.0 m.The coefficients of the permittivity function of the sphere are χ (1) = 2.0 and χ (3) = 0.The time dependence of the plane wave excitation in ( 56) is a modulated Gaussian pulse expressed as Here, f 0 , t p , and σ are the modulation frequency, time delay, and duration of the pulse, respectively.Let f bw represent the effective bandwidth, then choosing σ = 3/ (2πf bw ) ensures that 99.997% of P (t)'s power is within the frequency band [f min , f max ], where f min = f 0 − f bw and f max = f 0 + f bw [35].
For this example, f 0 = 5.0 MHz, f bw = 2.5 MHz, σ = 0.1910 µs, and t p = 8σ.E(r, t) and D(r, t) induced inside the sphere are discretized using N E = 3 844 and N D = 2 114 spatial basis functions, respectively.The MOT scheme is executed for N t = 1 000 time steps with ∆t = 6.0 ns.The SOR coefficient in ( 48) is selected as α = 0.3.

Nonlinear Cube
In the second example, electromagnetic scattering from a nonlinear cube is analyzed using the proposed solver.The cube is centered at the origin and has an edge of length 0.1 m.
The coefficients of the permittivity function of the cube are χ  From this perspective, one can argue that the accuracy of the proposed solver is actually higher than MEEP since the late-time solution obtained by the proposed solver reaches a lower level [39,47,48,60].This is further investigated by studying the Fourier transform of the solutions as described next.The following error is used as a measure of the convergence: where E sca,ref 2n−1,k is the scattered electric far-field computed by the simulation with the densest mesh.Fig. 2(d) plots err versus l av and shows that the solution obtained by the proposed solver converges with increasing mesh density.The time dependence of the plane wave excitation, P (t) given by ( 60).(c) x-component of E(r, t) computed at the feeding end [r 0 = (0, 0, −2.6 µm)] and the trailing end [r 0 = (0, 0, 2.6 µm)] in the first simulation, where both layers are linear.(d) x-component of E(r, t) at the trailing end [r 0 = (0, 0, 2.6 µ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.

Four-wave Mixing
In this example, four-wave mixing frequency conversion of electromagnetic fields [15] is analyzed using the proposed method.The scatterer is a sphere that is centered at the origin and has a radius of length 1.0 m.Two simulations are carried out to clearly demonstrate that the nonlinearity results in four-wave mixing.In the first simulation, the sphere is linear and the coefficients of its permittivity function are χ (1) = 1.5 and χ (3) = 0.In the second simulation, the sphere is nonlinear and the coefficients of its permittivity function are χ (1) = 1.5 and χ (3) = 0.075.In both simulations, the time dependence of the plane wave excitation in (56) is a sum of two modulated Gaussian pulses and is expressed as where These peaks are observed in the electromagnetic response because of the four-wave mixing frequency conversion generated as a result of the χ (3) -term in the permittivity function (i.e., Kerr nonlinearity) [15].
[r 0 = (0, 0, 2.6 µ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.The figure shows that with the introduction of the nonlinearity, the electric field at the trailing end is enhanced.
This can be explained by the fact that the stop band for the linear Bragg grating can be partially closed by introducing a negative Kerr nonlinearity [12].

Conclusion
An explicit MOT-based TD-EFVIE solver is developed to analyze 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.Discretizing TD-EFVIE using SWG functions yields a system of ODEs in time-dependent SWG expansion coefficients.This system is integrated in time using a P E(CE) m scheme to obtain the expansion coefficients of the electric field.Similarly, the nonlinear constitutive relation and its inverse obtained using the Padé approximant are discretized using SWG functions.The resulting matrices are used to carry out explicit updates of electric field and electric flux expansion coefficients at the predictor (PE) and the corrector (CE) stages.This approach produces an explicit MOT scheme that does not call for any Newton-like nonlinear solver but only requires solution of sparse and well-conditioned Gram matrix systems at every step.These solutions are done very efficiently using an iterative solver.
The accuracy and the applicability of the explicit MOT-based TD-EFVIE solver are demonstrated using several numerical examples.These results clearly show that the proposed method is more accurate than FDTD that is traditionally used for analyzing electromagnetic scattering from nonlinear objects.
Extension of the proposed scheme to analyze nonlinear problems involving high-contrast

5 )
SOR in(48) [Step 7.2  ] balances between the stability and the convergence of the corrector updates (CE) m , and it does not affect the accuracy of the solution (assuming convergence).Reducing the SOR coefficient α increases the number of corrector updates, which in return increases the computation time.On the other hand, increasing α might result in unstable corrector updates leading to an unstable solution.Additionally, for stronger nonlinearities, one must reduce α to maintain the stability, which again comes with increased computation time.

Fig. 1 (
Fig.1(a) plots the x-component of E(r, t) computed by the proposed solver at the center of the sphere (the sphere is centered at the origin).The figure shows that the proposed solver provides stable results for the entire simulation duration.After the time domain simulation is completed, the Fourier-transformed solution is used to compute the radar cross section (RCS) σ MOT (θ) at f = 5.0 MHz on the φ = 0 plane (θ ∈ [0 • , 180 • ]).Let σ Mie (θ) represent the RCS computed at the same frequency and on the same plane using the Mie series solution[1].Fig. 1(b) plots σ Mie (θ) and |σ Mie (θ) − σ MOT (θ)| versus θ.The figure clearly shows that the result obtained using the proposed solver is accurate.

Figure 2 :Figure 3 :Fig. 2 (
Figure 2: Scattering from a nonlinear cube.(a) x-component of D(r 0 , t) computed by the proposed solver (MOT) and the FDTD-based solver MEEP at the center of the cube [r 0 = (0, 0, 0)].(b) Zoomed version of (a) in the time range [32.5, 41.5] ns.(c) Fourier transform of the x-component of D(r 0 , t) computed by the proposed solver (MOT) and the FDTD-based solver MEEP.(d) Convergence in err defined by(58) with increasing mesh density (decreasing average edge length l av ).

Fig. 2 (
Fig.2(c) compares the Fourier transform of the solutions in Fig.2(a) in the frequency range f ∈ [0, 30.0]GHz.The figure clearly shows that several higher-order harmonics are generated due to the nonlinearity of the dielectric permittivity.Fourier transformed solutions match well up to frequencies where the numerical error in the MEEP solution becomes high enough to significantly effect the solution.Indeed, the mismatch between the two Fourier transformed solutions around the fifth harmonic is on the order of the difference in the levels of the two solutions in the late time as shown in Fig.2(a).Next, it is demonstrated that the solution obtained by the proposed solver converges with increasing mesh density.Four different meshes with average edge length l av = {1.896,1.716, 1.634, 1.553} cm are considered.This results in N E = {5 860, 7 596, 9 108, 10 784} and N D = {3 110, 4 002, 4 812, 5 668} for four simulations, respectively.All other parameters are kept the same.After each time domain simulation, Fourier transformed solutions are used to compute the scattered electric far-field E sca 2n−1,k for φ = 0 and θ = k∆θ, ∆θ = 1.0 • , k = 1, 2, . . ., 180 at the center frequencies of the first ten harmonics (2n − 1)f 0 , n = 1, 2 . . ., 10.

Figure 4 :
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 (60).(c) x-component of E(r, t) computed at the feeding end [r 0 = (0, 0, −2.6 µm)] and the trailing end [r 0 = (0, 0, 2.6 µm)] in the first simulation, where both layers are linear.(d) x-component of E(r, t) at the trailing end [r 0 = (0, 0, 2.6 µ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.

Fig. 3 (
Fig. 3(c) compares the x-component of D(r, t) computed by the proposed solver at the center of the sphere in the first (χ (1) = 1.5, χ (3) = 0) and the second (χ (1) = 1.5, χ (3) = 0.075) simulations.The figure shows that both solutions are stable.Fig. 3(d) compares the Fourier transform of these solutions in the frequency range f ∈ [0, 20] MHz.The two peaks observed at frequencies f = f 1 and f = f 2 in both solutions match the peaks in the Fourier transform of the excitation pulse shown in Fig. 3(b) (f 1 and f 2 are the modulation frequencies of the two Gaussian pulses added in (59)).However, the Fourier transform of the solution in the second simulation has two extra peaks at frequencies f = 2f 1 − f 2 and f = 2f 2 − f 1 .