Explicit Time Marching Schemes for Solving the Magnetic Field Volume Integral Equation

A method for constructing explicit marching-on-in-time (MOT) schemes to solve the time-domain magnetic field volume integral equation (TD-MFVIE) on inhomogeneous dielectric scatterers is proposed. The TD-MFVIE is cast in the form of an ordinary differential equation (ODE) and the unknown magnetic field is expanded using curl conforming spatial basis functions. Inserting this expansion into the TD-MFVIE and spatially testing the resulting equation yield an ODE system with a Gram matrix. This system is integrated in time for the unknown time-dependent expansion coefficients using a linear multistep method. The Gram matrix is sparse and well conditioned for Galerkin testing and consists of only four diagonal blocks for point testing. The resulting explicit MOT schemes, which call for the solution of this matrix system at every time step, are more efficient than their implicit counterparts, which call for inversion of a fuller matrix system at lower frequencies. Numerical results compare the efficiency, accuracy, and stability of the explicit MOT schemes and their implicit counterparts for low-frequency excitations. The results show that the explicit MOT scheme with point testing is significantly faster than the other three solvers without sacrificing from accuracy.


I. INTRODUCTION
A NALYSIS of electromagnetic scattering from inhomogeneous dielectric objects finds applications in numerous areas ranging from medical diagnostics to geophysical surveys. Simulation tools developed for these applications often rely on finite difference time-domain (FDTD) techniques [1]- [4], frequency and time-domain finite element methods (FEMs) [5]- [8], time-domain discontinuous Galerkin (TD-DG) schemes [9]- [13], or frequency and time-domain volume integral equation (VIE) solvers [14]- [31]. VIE solvers are often preferred over differential equation solvers (such as FDTD, FEM, and TD-DG), for open-region scattering problems, since they require only the scatterer to be discretized and implicitly enforce the radiation condition without the need for (approximate) absorbing boundary conditions to terminate the computation domain [32]. Time-domain VIE (TD-VIE) solvers are preferred over their frequency domain counterparts (FD-VIE solvers) for broadband scattering problems [23]- [31] and/or when the scatterer's permittivity is a (nonlinear) function of the fields [30], [31].
VIEs on dielectric scatterers are constructed by enforcing the fundamental field relation (the total field is equal to the summation of the scattered and incident fields) in the (volumetric) support of the scatterer. The scattered field is represented in terms of the equivalent (unknown) total field/flux induced inside the scatterer. If the electric field is used in the fundamental field relation, then the resulting VIE is termed the electric field VIE (EFVIE) [23]- [31] and the unknown is either the electric flux or the electric field induced inside the scatterer. On the other hand, if the magnetic field is used in the fundamental field relation and the unknown is the total magnetic field, then the magnetic field VIE (MFVIE) is obtained [14].
It is well known that the MFVIE has better convergence characteristics than the EFVIE [16]- [22]. It has also been shown that Galerkin discretization of the EFVIE using Schaubert-Wilton-Glisson (SWG) functions [33] yields a more accurate solution than the MFVIE Galerkin-discretized using the lowest order curl-conforming Nedelec functions [34], [35]. However, this accuracy bottleneck has been circumvented for all practical purposes by using recently developed fully linear curl-conforming (FLC) basis functions to expand the unknown magnetic field [17], [35], [36]. In this case, one can also use point testing, which significantly reduces the computational cost (in comparison to Galerkin testing) without sacrificing from accuracy.
Recent research on FD-VIE solvers has mostly focused on spatial discretization techniques and their effects on the accuracy and conditioning of the resulting matrix systems [14]- [22]. On the other hand, research on TD-VIE solvers has been geared toward developing accurate, efficient, and stable marching-on-in-time (MOT) schemes. TD-VIEs are usually solved using implicit MOT schemes [23]- [26] that call for solution of a matrix system (termed MOT matrix system here) at every time step. These schemes are This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ not subject to a Courant-Friedrichs-Lewy (CFL) constraint; their time step size is determined only by the maximum frequency of the excitation. For high-frequency excitations, i.e., when the product of the speed of light and the time step size is comparable to the discretization length, the MOT matrix system is sparse and it is solved efficiently using an iterative method. However, for low-frequency excitations, the MOT matrix system becomes fuller and it cannot be solved efficiently using an iterative method.
Depending on the spatial and temporal discretization schemes and the time step size, the MOT scheme can also be explicit. Even though classical explicit MOT schemes do not call for a matrix solution at every time step, they suffer from stability issues [27], [28], which might be remedied using a small time step size at the cost of increased computation time (i.e., they are subject to a CFL constraint).
This article describes a method for constructing explicit MOT schemes, which do not suffer from these shortcomings, to efficiently and accurately solve the TD-MFVIE. The proposed method casts the TD-MFVIE in the form of an ordinary differential equation (ODE) that relates the unknown magnetic field induced inside the scatterer to its temporal derivative [37]. The magnetic field is expanded using the FLC basis functions [35], [36]; inserting this expansion in the TD-MFVIE and spatially testing the resulting equation yields a time-dependent ODE system. A predictor-corrector algorithm, PE(CE) m , is used to integrate this system in time for the unknown coefficients of the expansion. To facilitate the computation of the retarded-time integrals, which express the scattered magnetic field in terms of the unknown magnetic field induced inside the scatterer, at discrete time steps as required by the PE(CE) m , the piecewise Lagrange polynomial interpolation functions [38]- [40] are used. The resulting time marching algorithm calls for the solution of a system with a (spatial) Gram matrix at the evaluation (E) step. When Galerkin testing is used, the Gram matrix is sparse and well conditioned, and the solution is obtained using an iterative solver. When point testing is used, the Gram matrix consists of four diagonal sub-matrices. Its inverse (which also consists of four diagonal sub-matrices) is computed and stored before the time marching starts. Consequently, the matrix solution required at the evaluation step is obtained with a simple multiplication of the right-hand side with the inverse of the Gram matrix. The resulting MOT schemes are expected to be more efficient than their implicit counterparts, which call for the inversion of a matrix system that gets fuller as the time step size gets larger with decreasing frequency. Indeed, the numerical results demonstrate that the explicit MOT schemes use the same time step sizes as the implicit MOT schemes without sacrificing from stability, and they are more efficient under low-frequency excitations. Especially, the explicit MOT scheme with point testing is significantly faster than the other three solvers without sacrificing from accuracy.
The rest of this article is organized as follows: Section II provides the details of the formulation underlying the explicit and implicit MOT schemes with Galerkin and point testing and derives expressions for their computational complexity estimates. Section III compares the efficiency, stability, and accuracy of the explicit MOT schemes and their implicit counterparts for low-frequency excitations via numerical experiments and demonstrates that the explicit scheme with point testing is significantly faster than the other three without sacrificing from accuracy. In Section IV, conclusions and future research directions are drawn.

A. TD-MFVIE
Let V represent the volumetric support of a linear, non-dispersive, non-magnetic, isotropic, and possibly inhomogeneous dielectric scatterer with permittivity ε(r) and permeability μ 0 . The scatterer resides in an unbounded and homogenous medium with permittivity ε 0 and permeability μ 0 . An incident magnetic field H inc (r, t), which is essentially band limited to f max and vanishingly small ∀r ∈ V and t ≤ 0, excites the scatterer. Upon excitation, an equivalent electric current J(r, t) is induced inside V , which in return generates a scattered magnetic field H sca (r, t). H sca (r, t) is expressed in terms of retarded-time magnetic vector potential A(r, t) as Here R = |r − r | is the distance between source point r and observation point r, and c 0 = 1/ √ ε 0 μ 0 is the speed of light in the background medium. J(r, t) is expressed in terms of the total magnetic field H(r, t) as where κ(r) = 1 − ε 0 /ε(r) is the contrast. Substituting (1) and (2) in the temporal derivative of H(r, t) = H inc (r, t) + H sca (r, t) yields the TD-MFVIE

B. Spatial Basis Functions and Temporal Interpolation
To numerically solve the TD-MFVIE (4), V is divided into a mesh of tetrahedrons. Assume that this mesh has N edges. H(r, t) is approximated in terms of the FLC basis functions [35], [36], each of which is defined along one of these edges as Note that this expansion follows the description in [36], where the FLC basis functions are separated into solenoidal and irrotational edge basis functions. In (4), f 1 n (r) and f 2 n (r) are the first-order irrotational edge basis functions [36] and the lowest mixed-order solenoidal edge basis functions [34], and {H 1 (t)} n and {H 2 (t)} n are their unknown time-dependent coefficients, respectively. f s n (r), s ∈ {1, 2} are expressed as where "+" and "−" signs should be selected for s = 1 and s = 2, respectively, S n = ∪ Q n q=1 S q n is the combined support of all Q n tetrahedrons sharing edge n, d 1 n and d 2 n represent the two nodes of this edge, and λ d n (r) and d ∈ {d 1 n , d 2 n } are the barycentric coordinate functions that change linearly from 1 at d to 0 at the face opposite to d. Note that one can easily show ∇ × f 1 n (r) = 0 and ∇ × f 2 n (r) = 0. To facilitate the discretization and computation of the retarded-time integral in the right-hand side of (4), {H s (t)} n , s ∈ {1, 2} are approximated using (shifted) Lagrange interpolation functions as Here N t is the number of time steps, t is the time step size, T (t) is a piecewise polynomial Lagrange interpolation function [38]- [40], and H s i is the sample of H s (t) at t = i t, i.e., H s i = H s (i t).

C. Explicit MOT Scheme
Inserting (4) in (4) and testing the resulting equation with functions t 1 m (r) and t 2 m (r), m = 1, . . . , N, yield an ODE matrix system of dimension 2N × 2N, which relates unknown vectors H s (t) to their temporal derivativesḢ Here G ps , p, s ∈ {1, 2} are N × N blocks of the Gram matrix G. Their elements are given by where P p m is the support of t p m (r), p ∈ {1, 2}. Two sets of choices are considered for t 1 m (r) and t 2 m (r), which result in Galerkin and point testing, respectively. The specific choice of the testing scheme changes the sparseness structure of G and consequently affects the efficiency and accuracy of the time marching scheme (see Section II-E).
In (7),Ḣ inc, p (t) andḢ sca, p (t), p ∈ {1, 2} are vectors of spatially tested incident and scattered magnetic fields, respectively. Their entries are given by In (10), κ(r) is assumed to be constant in S q n , i.e., κ q n = κ(r q n ), where r q n is the center of the tetrahedron S q n . Also note that since ∇ × f 1 n (r) = 0, the only contribution toḢ sca, p (t) comes from ∇ × f 2 n (r). The samples of the unknown coefficient vectors H s j = H s ( j t), s ∈ {1, 2} are obtained by integrating the system of ODEs in (7) in time using a PE(CE) m -type linear k-step scheme. This approach calls for sampling (7) in time , one has to account for the time retardation in (10); this is done by using temporal interpolation on samples of H 2 (t). Inserting (6) with s = 2 in (10) and evaluating the resulting expression at j t yielḋ where the elements of the MOT matrices M p j −i are given by Substituting (12) into (11) yields Note that to have a more compact notation, (14) is rewritten as The matrix system (15) is integrated in time using a PE(CE) mtype linear k-step method (similar to the one used in [37] to solve the time-domain magnetic field surface integral equation). Therefore, it requires the values of . . , k − 1 are known, the steps of the resulting explicit MOT scheme are detailed below. At each time step, j = k, . . . , N t .
Step 1: The components of the right-hand side of (15), which are not updated within the time step j , are computedḢ Note thatḢ exp j does not include the contributions from H j , i.e., the matrix-vector product M exp 0 H j .
Step 2: Here p is a vector of dimension 2k, which stores the predictor coefficients.
Step 3: Evaluation (E) step. First compute the right-hand side using the predicted H j Then computeḢ j by solving Step 4: SetḢ Step 4.1: Here c is a vector of dimension 2k+1, which stores the corrector coefficients.
Step 4.2: Evaluation (E) step. First compute the right-hand side using the corrected Then computeḢ (m) j by solving Step 5: Once convergence is reached, solutions are stored to be used at the next time step: j . The predictor and corrector coefficients, p and c used in the above scheme can be obtained by polynomial interpolation between time samples (resulting in well-known schemes such as Adam-Moulton, Adam-Bashfort, or backward difference methods [41]) or numerically under the assumption that the solution can be represented in terms of decaying and oscillating exponentials [42]. In this article, p and c obtained through polynomial interpolation are preferred since k associated with these coefficients is much smaller resulting in a more timeand memory-efficient scheme.
At the beginning of time marching, it is assumed that This assumption does not introduce any significant error since H inc (r, t) is vanishingly small ∀r ∈ V and t ≤ 0. For other types of excitations, the Euler method or spectral-deferred correctiontype methods can be used to initialize H i andḢ i , i = 0, . . . , k − 1 [43], [44].
The method used for solving (19) and (22) is selected based on the sparsity structure of G, which depends on the type of spatial testing used as detailed in Section II-E.

D. Implicit MOT Scheme
Inserting (4) and (6) in (4) and testing the resulting equation with functions t 1 m (r) and t 2 m (r), m = 1, . . . , N, yield a linear system of equations Here H j andḢ inc j are same as those in (14), and M imp l , The implicit MOT scheme operates as briefly described next. For j = 1, H 1 is found by solving (23)  Unlike the explicit MOT scheme, the method used for solving (23) at every time step of the implicit MOT scheme does not depend on the sparsity structure of G. The solution of this matrix equation is always obtained using an iterative solver.

E. Spatial Testing Functions
Two different approaches are used to spatially test the TD-MFVIE: Point and Galerkin testing schemes [36].
The inverse of G can be expressed as G −1 is stored using O(N) memory before the time marching starts. Using the pre-computed G −1 , the solution of (19) and (22) is obtained only in O(N) operations. This makes the explicit MOT scheme with point testing significantly faster than the other explicit and implicit MOT schemes under low-frequency excitations (for large t) as shown by the computational complexity analysis carried out in Section II-G (as also demonstrated by the numerical results presented in Section III).
2) Galerkin Testing: For Galerkin testing, (8), one obtains a summation of integrals, each of which has a second-order polynomial integrand defined over a tetrahedron. These integrals are evaluated exactly using a Gaussian quadrature rule specifically designed for tetrahedrons [45], [46]. Analytical expressions can be derived for these integrals but evaluating those would be more computationally expensive than using a quadrature rule.
When Galerkin testing is used, G is sparse and well conditioned regardless of t. Therefore the solution of (19) and (22) is obtained very efficiently using an iterative scheme. The resulting MOT scheme with Galerkin testing is faster than its implicit counterpart under low-frequency excitations (for large t). This is because M 1 0 and M 2 0 become fuller (see Section II-F) as t increases and the computational cost of solving (23) increases. Implicit and explicit schemes have similar computational costs under high-frequency excitations (for small t). These are shown by the computational complexity analysis carried out in Section II-G (also by the numerical results presented in Section III).

F. Comments
Several comments about the formulation and implementation of the explicit and implicit MOT schemes are in order. 1) In (4) H(r, t).
The effect of the additional temporal derivative has been studied in the context of the MOT solution of the time-domain electric field surface integral equation [47]. When the derivative of a time-domain integral equation is solved by an MOT scheme, a dc component is often observed in the solution even if the zero initial condition is enforced at the beginning of time marching. The reason for this dc component is attributed to numerical errors, especially those arising from the solution of the MOT matrix system as discussed in [47]. The numerical results presented in Section III-B demonstrate that a similar discussion applies in the case of the MOT schemes developed in this article to solve (4). The results show that the dc component is stable with a small amplitude, and this amplitude can be reduced by increasing the accuracy of the matrix solution and/or correction updates. 2) Temporal interpolation function T (t) is discretely causal: T (t) = 0 for t ≤ − t. This means that, during time marching (as executed by both the explicit and implicit MOT schemes), are not required to compute H j . Additionally, where t max is the order of the polynomial interpolation. Note that T (t) can be a noncausal interpolation function.
For example, one can use the approximate prolate spheroidal wave functions (APSWFs) [48], [49] since they can interpolate band-limited functions with exponentially increasing accuracy and they have continuous derivatives everywhere along their support. A non causal T (t) means that M p j −i = 0, i > j and consequently the time marching requires H i , i > j to be able to compute H j . The causality of the time marching can be restored using various extrapolation schemes that estimate H i , i > j using H i , i ≤ j [26], [49].
3) The MOT matrices M [see (24)] higher than that ofḢ exp j [see (16)]. However, since t max D max /(c 0 t) for small t (under high-frequency excitations) and G 11 and G 21 are sparser than M 1 j −i and M 2 j −i for large t (under low-frequency excitations), this difference can be ignored in the computational complexity analysis.

G. Computational Complexity
In this section, computational complexity of the explicit MOT scheme described in Section II-C is analyzed in detail and compared to that of its implicit counterpart described in Section II-D. Let the computational costs of explicit schemes with point and Galerkin testing and the implicit scheme be represented by C exp PT N t + C, C exp GT N t + C, and C imp N t + C, respectively. Note that the implicit MOT scheme can also be implemented using Galerkin or point testing, but the expressions for these implementations' computational complexity would be the same. That is why C imp does not distinguish between these two implementations.
Here C represents the total cost of computingḢ fixed j for all time steps j = 1, . . . , N t and is dominated by the cost of computingḢ exp j . As explained in Section II-F, the cost of computingḢ exp j (by the explicit MOT schemes) is same as that ofḢ imp j (by the implicit MOT schemes). Therefore, C is assumed same for all schemes. Note that this computation could be accelerated using the time-domain adaptive integral method [40], [50]- [54] or (multilevel) plane wave time-domain algorithm [24], [25], [28], [39], [55], [56].
The differences between the explicit schemes and their implicit counterparts are the other operations executed at a given time step. The computational costs of these operations are represented by C The implicit MOT scheme always uses an iterative method to solve (24), which results in iter is the number of iterations and F iter is the number of matrix-vector multiplications required at a given iteration (it is assumed that the explicit and implicit schemes use the same iterative solver).
In the complexity estimates above, k depends on the order/type of the PE(CE) m ; therefore, it is considered as a user-defined input. Also, N G iter is always small since G is well conditioned and sparse regardless of t. Assuming m max is the same for explicit schemes with point and Galerkin testing, the former scheme is faster since 4 N G iter F iter 2δ. Numerical results in Section III show that the value of m max averaged over all time steps is similar for both the schemes.
Under high-frequency excitations when c 0 t is comparable to the spatial discretization length, γ N and a direct comparison of C exp PT and C exp GT to C imp becomes challenging since it is difficult to accurately estimate which contributions discussed above are dominant.
Under low-frequency excitations when c 0 t is comparable to or larger than the size of the scatterer, γ ∼ N, which means that M 1 0 and M 2 0 become fuller. Consequently, ). This means that the explicit schemes are faster than their implicit counterparts as long as m max < N imp iter F iter . Numerical results presented in Section III show that this condition is indeed satisfied.
Note that one can use the low-frequency extension of the time-domain adaptive integral method [57] to accelerate the matrix-vector multiplication M imp 0 H j required by the iterative method to solve (23). However, the same extension can also be used to accelerate the computation of the matrix-vector multiplication M exp 0 H j required by the explicit scheme during the predictor updates. Therefore, the conclusions drawn above are still applicable even when acceleration methods are used.

III. NUMERICAL RESULTS
This section presents numerical examples to demonstrate the advantages of the proposed explicit MOT schemes. In all examples, the scatterer is illuminated by a plane wave traveling in theẑ direction with aŷ-directed magnetic field where H 0 = √ 0 /μ 0 A/m is the amplitude and G(t) = cos[2π f 0 (t − t p )]e −(t −t p ) 2 /(2σ 2 ) is the modulated Gaussian pulse. Here σ = 3/(2π f bw ) is the duration, f bw is the effective bandwidth, f 0 is the center frequency, and f max = f 0 + f bw is the maximum frequency of the pulse. It is assumed that the scatterer resides in free space. In all examples, the order of the piecewise polynomial Lagrange interpolation function T (t) t max = 4 [40], and the volume integrals present in the entries of the matrices in (8) and (13), and the vector in (9) are computed using the third-order Gauss-Legendre quadrature rule [45], [46].
The accuracy, efficiency, and stability of the four MOT schemes are compared. The implicit scheme with point testing, the explicit scheme with point testing, the implicit scheme with Galerkin testing, and the explicit scheme with Galerkin testing. For the sake of briefness, in the rest of this section, these schemes are referred to using the notations [MOT] exp GT are terminated when the following stopping criterion is satisfied: Here I u l represents the solution vector at the lth time step and the uth TFQMR iteration or at the lth time step and the uth correction update, and χ is the convergence threshold. χ = 10 −6 unless specified otherwise. The PE(CE) m scheme uses the fourth-order Adam-Bashworth and backward difference coefficients at the prediction and correction steps, respectively [41].
After the time-domain simulations are completed, the solutions are Fourier transformed and divided by the Fourier transform of G(t) to yield the time harmonic magnetic field, H(r, f ). Time harmonic electric fieldẼ(r, f ) and the time harmonic current densityJ(r, f ) are computed by taking the curl of (4). The radar cross section (RCS) is computed using is used. Here type ∈ {imp, exp}, test ∈ {GT, PT}, ref ∈ {Mie, FD}, θ = 0.5 • , and φ = 0 • . In (29), σ Mie (θ, φ, f ) and σ FD (θ, φ, f ) refer to the RCS obtained from the Mie series solution orJ(r, f ) computed by an FD-EFVIE solver. The FD-EFVIE solver uses the same mesh as the MOT schemes but discretizes the electric flux density using SWG basis and testing functions [33]. The entries of the resulting method of moments (MoM) matrix are computed using the third-order Gauss-Legendre quadrature rule [45], [46]. The MoM system is solved using the TFQMR method. The iterations are truncated when ZĨ u −ZĨ u−1 < 10 −6 Ṽ inc , whereĨ u ,Z, andṼ inc represent the solution at the uth iteration, the MoM matrix, and the right-hand side vector, respectively.

A. Accuracy of the FLC and Nedelec Basis Functions
In this example, the scatterer is a sphere with radius 2 m and permittivity ε 0 . Since κ(r) = 0, H sca (r, t) = 0 and H(r, t) =   (31) are used. Here r k represent the centers of the tetrahedrons and N v is their number, {Ẽ inc (r, f ),H inc (r, f )} are the time harmonic incident electric and magnetic fields, and f = f 0 . Fig. 1(a) and (b) plot err H and err E versus λ 0 /l av , respectively, for the simulations with only f 2 n (r) executed for three different t. Fig. 1(a) and (b) show that err H decreases with increasing mesh density and decreasing t, while err E remains high even for the densest mesh and the smallest t. Fig. 1(c) and (d) do the same comparison for simulations with f 1 n (r) ∪ f 2 n (r). Both err H and err E decrease with increasing mesh density and decreasing t. Fig. 1(c) and (d) clearly show that using only f 2 n (r) results in an inaccurate representation of E(r, t), while using f 1 n (r)∪f 2 n (r) renders E(r, t) as accurate as H(r, t). In other words, f 2 n (r) accurately represent the solution, but the curl of the resulting solution is not accurate. When κ(r) = 0, the curl of the solution is needed to computė H sca, p j , p ∈ {1, 2} [see (10)]. This means that using only f 2 n (r) makes the MOT solution inaccurate (and consequently unstable). The results and the discussion presented in this section clearly justify why f 1 n (r) ∪ f 2 n (r) are used by the MOT schemes developed in this article. Fig. 1(c) and (d) also show that err H and err E decrease with a rate roughly between (l av ) −1 (for large t) and (l av ) −2 (for small t) for λ 0 /32 > l av > λ 0 /45 and with a rate around (l av ) −1 for l av < λ 0 /45. This decrease in the convergence rate can be explained by the fact that for small l av , the dominant error comes from the temporal discretization. This is demonstrated in Fig. 1(c) and (d); for smaller t, the "flattening" of err H and err E curves starts at smaller l av . In other words, for small l av , the accuracy can further be increased by reducing t.

B. Late Time Stability
For this example, scattering from a sphere with radius 1 m and permittivity 10ε 0 is analyzed for different values of the convergence threshold χ. The sphere is discretized using 5350 tetrahedrons resulting in N = 13 494 unknowns. The excitation parameters f 0 = 10 MHz and f bw = 5 MHz. The average, minimum, and maximum edge length of the mesh are l av = λ min /33.28, l min = λ min /62.0, and l max = λ min /19.76, respectively. Here λ min = c 0 /( √ 10 f max ) is the wavelength at f max inside the scatterer. All four MOT schemes are executed for N t = 1200 with t = 6.667 ns (0.1/ f max ) and three different χ ∈ {10 −6 , 10 −8 , 10 −10 }. Fig. 2(a)-(c) plot H(r, t) computed by these schemes at point r = (0.51, −0.64, 0.12) m for χ = 10 −6 , χ = 10 −8 , and χ = 10 −10 , respectively, and show that all four schemes provide stable results with a very small dc component and that the amplitude of the dc component can be further reduced by decreasing χ.      Table I's first group of rows, verify this result. Note that in Table I, the fourth column T fill is the time required to compute all relevant matrices, the fifth column T MOT refers to C exp PT N t + C, C exp GT N t + C, or C imp N t + C (see Section II-G) depending on the scheme used, and the sixth column T tot is the sum of the previous two.
For the same set of simulations with t = 6.667 ns, Fig. 3 Mie (θ, φ, f ), all of which are computed for 0 • < θ < 180 • and φ = 0 • at f = f 0 , and shows that all four MOT schemes practically have the same accuracy. Additionally, the last column of the first group of rows in Table I provides err RCS computed using (29) for this set of simulations and confirms that all four MOT schemes have the same level of accuracy. Table I's second and third groups of rows compares the efficiency and accuracy of the MOT schemes for the sets of simulations with t = 10 ns and t = 13.333 ns. Same conclusions can be drawn: All four MOT schemes have the same level of accuracy, and [MOT] exp PT is significantly faster than the other three. Table I also shows that the accuracy of all schemes increases with decreasing t.
Next, the permittivity of the sphere is set to 50ε 0 . The sphere is discretized using 11697 tetrahedrons, resulting in N = 28970 unknowns. The average, minimum, and maximum edge length of the mesh are l av = λ min /19.51, l min = λ min /39.34, and l max = λ min /11.55, respectively. Here is the wavelength at f max inside the scatterer.

[MOT]
exp GT is executed three times for N t = 600 with t = 6.667 ns (0.1/ f max ), N t = 400 with t = 10 ns (0.15/f max ), and N t = 300 with t = 13.333 ns (0.2/ f max ). Fig. 4(a) plots H(r, t) computed at point r = (0.51, −0.64, 0.12) m for the simulation with t = 6.667 ns. Note that for this problem, the other three schemes do not produce stable results. For the same simulation, Fig. 4

D. Piecewise Slab
In the last example, scattering from a piecewise dielectric slab is analyzed. The slab consists of two equal volumes with permittivities 3ε 0 and 9ε 0 [as shown in the inset of Fig. 5(a)]. The slab is discretized using 7905 tetrahedrons resulting in N = 20 570 unknowns. The excitation parameters f 0 = 10 MHz and f bw = 5 MHz. The average, minimum, and maximum edge length of the mesh are l av = λ min /31.7, l min = λ min /66.6, and l max = λ min /19, respectively. Here is the wavelength at f max inside the right side of the slab. All four schemes are executed three times for N t = 210 with t = 6.667 ns (0.1/ f max ), N t = 140 with t = 10 ns (0.15/ f max ), and N t = 105 with t = 13.333 ns (0.2/ f max ). For all simulations, the sparseness factor of M 1 0 and M 2 0 is γ = N and the sparseness factor of G is δ = 0.0033N. Fig. 5(a) plots H(r, t) computed by all four schemes at point r = (0.23, 0.14, 0.57) m for the set of simulations with t = 6.667 ns. The results agree very well. For the same set of simulations, Fig. 5(b) plots the number of correction imp GT to achieve the convergence criterion in (28) at every time step. The number of iterations required to solve (19) and (22)    For the same set of simulations with t = 6.667 ns,  Table III provides err RCS computed using (29) for this set of simulations and confirms that all four MOT schemes have the same level of accuracy. Table III's second and third groups of rows compares the efficiency and accuracy of the MOT schemes for the sets of simulations with t = 10 ns and t = 13.333 ns. The results show that all four MOT schemes have the same level of accuracy, and [MOT] exp PT is significantly faster than the other three. Also, as expected, the accuracy of all schemes increases with decreasing t.

IV. CONCLUSION
A method for constructing explicit MOT schemes to solve the TD-MFVIE enforced on dielectric scatterers is developed. The TD-MFVIE is first cast in the form of an ODE and the unknown magnetic field is expanded using the FLC basis functions. The expansion is inserted into the TD-MFVIE and the resulting equation is Galerkin or point tested in space. This yields an ODE matrix system, which is integrated in time using a PE(CE) m scheme for the (unknown) expansion coefficients. The resulting MOT scheme calls for the solution of a Gram matrix system at the evaluation (E) steps of every time step. This can be done very efficiently since the Gram matrix is always well conditioned and sparse (for Galerkin testing) or consists of only four diagonal blocks (for point testing). Numerical results demonstrate that the explicit MOT scheme with point testing is significantly faster without sacrificing from accuracy for low-frequency problems.
An extension of the proposed MOT scheme to enable the analysis of electromagnetic scattering from nonlinear dielectric objects is underway.