# Time-Dependent Lorentz-Mie-Debye Formulation for Electromagnetic Scattering from Dielectric Spheres.

(Invited Paper)

1. INTRODUCTION

In the last century, technologies that rely on exploiting solutions to Maxwell's equation [1] have been revolutionary in terms of their impact on society. The wealth and breadth of the impact of these technologies, which have affected almost every aspect of modern society, can be attributed to the succinct framework of electromagnetic phenomena via Maxwell's equations and methods (both theoretical and experimental) developed glean insight into the physics of electromagnetics. Some of the earliest solutions to these equations were developed for canonical geometries; a comprehensive catalog of canonical solutions can be found in [2]. A textbook example of a solution to a canonical problem is scattering by a spherical object in the frequency domain [3] and has been known for more than a century. It is based on the existence of analytic eigenfunction of Helmholtz wave equation in spherical coordinate system, and was formulated based on separation of variables. Interestingly, the history of analytical solutions for a spherical system goes back much further, perhaps, the most common being those of Lorentz [4] and Debye [5], both of which are almost contemporary to Mie's work. However, as is rather common, the solution was "reinvented" several times. An excellent review article in 1965 by Logan [6] went through most (if not all) of contributions related with scattering from spheres and uncovered the "lost" contributions of Clebsh [7], Nicholson [8], Bromwich [9], Proudman [10], Doodson [10], Kennedy [10], and White [11]. Widespread lack of awareness of these other work has largely resulted in calling these solutions either Mie series or Mie approach or Lorentz-Debye-Mie theory. A more detailed discussion on this topic can be found in [6] and [12]. It comes as no surprise that applications of this theory have had a widespread impact in a variety of frequency regimes [12,13].

Since the development of Lorentz-Mie-Debye formulation theory for a single sphere of homogeneous medium, research has focused on different applications ranging from scattering for different excitations [14] to scattering from spheres with different constitutive properties (dielectric [15], uniaxial [16], chiral [17]) to spheroidal shapes [18] to scattering from clusters [19]. While the classical Lorentz-Mie-Debye theory was largely relegated to studying scattering of plane wave incident on an object, more recent work has focused on the Generalized Lorentz-Mie theory (GLMT) [14, 20] that extends the approach to analyzing sphere-arbitrary wave interactions by expressing the incident field in terms of vector spherical wave functions. The GLMT has been used to study beam scattering by spheres and spheroids [14,18,21]. Although Mie series approach has been extended to spheroidal particle with spheroidal wave functions, it's difficult (if not impossible) to find analytical eigenfunctions with convergent closed form for scatterers of other shapes (noncanonical or even generalized ellipsoid). The null-field method based transition matrix (T-Matrix) method proposed in 1965 by Waterman [22-25] addresses this problem. The research on this approach as applied to multiple scattering is extensive and has found application in a variety of fields; for a very partial listing, see [26-32]. It should be noted that all the discussions are limited to three dimensional case and Lorentz-Mie- Debye formulation and its equivalent methods in two dimensional case also exist.

Thus far, the historical data presented have largely been restricted to solutions to these equations in the Fourier domain, i.e., wherein the time dependence is assumed to be of the form exp[j[omega]t]. In rather stark contrast, time-dependent Mie solution or Generalized Lorentz-Mie theory or T-matrix method has not been studied as much. Typically, when transient response is desired, one computes these from an inverse Fourier transform of the frequency domain solutions. The rationale for doing so is rather transparent; analytic expressions for transient solutions do not exist.

In general, compared to time-harmonic case, time dependent Lorentz-Mie-Debye formulation or more general time domain multiple theory has received far less attention from both the physics (EM/optics) and the mathematics community. This is not a result of lack of interest; the first piece of work, to the author's knowledge, was as early as 1919 [9]. In this work, general solution in time domain or time-dependent wave function for scalar wave equation was proposed, where its form involves differential operators on a time-dependent function, which could be reduced to Lorentz-Mie- Debye formulation if time-harmonic factor is assumed. Afterwards, time-domain representations for "spherical-polar-vector" components of dipole field were derived [33] and applied for field extrapolation [34]. In 1966, Granzow generalized earlier results to high order and developed the multipole theory in time domain (MTTD) [35]. Granzow's MTTD is also based on the high order differential operator working on arbitrary time-dependent functions to represent each component in the general solutions to wave equations, which makes the implementation complicated. His method has been used in two types of boundary values problems: finding the exterior field given (1) normal component of the field [35] (impulse function or Green's function approach) and (2) tangential components of the field on the sphere [36] (Laplace and inverse Laplace approach). Other research on time dependent multipole expansions include [37] and their applications [38], most of which, instead of working with time-domain versions of wave functions, focused mainly on radiation due to surface or volume sources. None of those methods provides a time-domain Lorentz-Mie-Debye solution for electromagnetic scattering from spheres. In [39], time-domain acoustic scattering from a sphere is studied with the help of time-dependent wave functions (inverse Fourier transform of the frequency domain spherical wave functions). The same idea is also extended to two-sphere case [40]. Considering the fact that mode matching concept is used to extract the time-domain signatures, this type of approach is direct time-domain variant of conventional Lorentz-Mie-Debye approach. However, the method involves explicit convolutions between inverse Fourier transform (IFT) of spherical Hankel and Bessel functions (the IFT of spherical Hankel function only exists in convolution sense according to [39]), which introduces instability especially for high order modes. Recently, for acoustic scattering, Laplace transform based approach [41, 42] and TDIE based approach [43] have been proposed to solve the time-dependent scattering problem. These two methods used the same expansions to represent the unknown physical quantities and didn't use spherical wave functions at all. Both methods have been extended to electromagnetics [44,45]. Compared to [39], the main advantages of TDIE based method lie in the stability and flexibility, because no explicit time domain wave functions is used. Though stable and self-consistent, Laplace based methods involves finding the roots of special functions (recursion rules have to be used). TDIE based method has explicit modal decoupling, where a set of independent Volterra integral equations will be solved, and extension into multilayer or multiple spheres problems is fairly straightforward.

Another approach, not Lorentz-Mie-Debye formulation based, but the so-called Debye series [46] based, has been proposed to study the behaviors of time-domain scattering from spheres [47, 48]. As in its frequency domain counterparts [46, 49], each term in the Mie series is represented by another set of Debye series. After exchange of the summations and introducing some large argument approximations for spherical Bessel functions, the scattering field could be written as a simple form of Debye series, which could be transformed into time-domain using inverse Fourier transform.

The rest of the paper is organized as follows: In Section 2 we develop the requisite formulation, which is then followed in Section 3 by a novel addition theorem that we will exploit. Using these addition theorems, we develop Mie solution in Section 4. Results and discussions are presented in Sections 5 and 6, respectively.

2. TIME DEPENDENT FORMULATION AND FIELD REPRESENTATION

Frequency domain Mie type solutions are typically obtained starting with eigenfunctions of the underlying differential equation. Early attempts to obtain transient solutions relied using inverse Fourier transform of the frequency domain solutions. The approach espoused in this paper is to use an integral equation based approach; note, in the frequency domain, it can be readily proven that both approaches yield identical results. To this end, consider a sphere with non-dispersive constitutive parameters ([[epsilon].sub.-], [[mu].sub.-]) that occupies a region denoted by [OMEGA] whose surface is denoted by [partial derivative][OMEGA]. Let [[OMEGA].sub.+] = [R.sup.3]/[OMEGA], and [[OMEGA].sub.-] = [OMEGA]/[partial derivative][OMEGA]. Surfaces [partial derivative][[OMEGA].sub.[+ or -]] can be thought of as boundaries of [[OMEGA].sub. [+ or -]] that are conformal to, and just outside and inside [partial derivative][OMEGA], respectively. These surfaces are equipped with a unique outward pointing normal that is denoted using [[??].sub.[+ or -]], where the subscript [+ or -] denotes that the currents that are defined either on [partial derivative][[OMEGA].sub.[+ or -]]. Henceforth, we denote quantities that are associated with [[OMEGA].sub.[+ or -]] using the subscript [+ or -]. Assume that an incident field denoted using {[E.sup.i]([bar.t], t), [H.sup.i]([bar.r], t)} that is approximately bandlimited to [f.sub.max] and quiescent for t < 0. The scattered field can be constructed using the equivalence theorem by first defining equivalent electric and magnetic currents [J.sub.[+ or -]]([bar.r], t) and [M.sub.[+ or -]]([bar.r], t). As usual, the currents are defined as [J.sub.[+ or -]]([bar.r], t) = [[??].sub.[+ or -]] x [H.sub.[+ or -]]([bar.r], t) and [M.sub.[+ or -]]([bar.r], t) = [E.sub.[+ or -]]([bar.r], t) x [[??].sub.[+ or -]]. The total field in regions exterior and interior to [partial derivative][OMEGA] can be defined as follows.

[E.sub.+]([bar.r], t) = [E.sup.i]([bar.r], t) + [E.sup.s.sub.+]([bar.r], t) = [E.sup.i]([bar.r], t) - [L.sub.+]([J.sup.+]([bar.r'], t)) + [K.sub.+]([M.sup.+]([bar.r'], t)) (1a)

[H.sub.+]([bar.r], t) = [H.sup.i]([bar.r], t) + [H.sup.s.sub.+]([bar.r], t) = [H.sup.i]([bar.r], t) - 1/[[eta].sup.2.sub.+] [L.sub.+]([J.sup.+]([bar.r'], t)) + [K.sub.+]([M.sup.+]([bar.r'], t)) (1b)

[E.sub.-]([bar.r], t) = [E.sup.s.sub.-]([bar.r], t) = - [L.sub.-]([J.sup.-]([bar.r'], t)) + [K.sub.-]([M.sup.-]([bar.r'], t)) (1c)

[H.sub.-]([bar.r], t) = [H.sup.s.sub.-]([bar.r], t) = -1/[[eta].sup.2.sub.-] [L.sub.-]([M.sup.-]([bar.r'], t)) - [K.sub.-]([J.sup.-]([bar.r'], t)) (1d)

where [[eta].sub.[+ or -]] = [square root of [[mu].sub.[+ or -]]/[[epsilon].sub.[+ or -]]] is the impedance of the medium.

The operators [L.sub.[+ or -]](*) and [K.sub.[+ or -]](*) are spatial-temporal operators associated with dyadic Green's functions, and are defined as follows:

[L.sub.[+ or -]](X([bar.r'], t))([bar.r], t) = [mu][[partial derivative].sub.t] [[integral].sub.S'] [[??].sub.e[+ or -]] ([bar.r], [bar.r'], t) [cross product] X([bar.r'], t)dS' [equivalent to] [mu][[partial derivative].sub.t]([[??].sub.e[+ or -]] [*.sub.st] J (2)

[K.sub.[+ or -]](X([bar.r'], t))([bar.r], t) = [[integral].sub.S'] [[??].sub.m[+ or -]] ([bar.r], [bar.r'], t) [cross product] X([bar.r'], t)dS' [equivalent to] [[??].sub.e[+ or -]] [*.sub.st] M (3)

In the above definitions, [[??].sub.e[+ or -]] and [[??].sub.m[+ or -]], respectively, denote dyadic Green's functions of electric type and magnetic type in the homogeneous exterior or interior space. The subscript e and m are used to denote electric and magnetic type, respectively. These Green's functions can be expressed as

[[??].sub.e[+ or -]]([bar.r], [bar.r'], t) = (I - [c.sup.2.sub.[+ or -]] [[partial derivative].sup.-2.sub.t] [nabla][nabla]) [delta] (t - R/[c.sub.[+ or -]])/r[pi]R (4a)

[[??].sub.m[+ or -]]([bar.r], [bar.r'], t) = [nabla] x [I [delta] (t - R/[c.sub.[+ or -]])/4[pi]R] = [nabla] [delta] (t - R/[c.sub.[+ or -]])/4[pi]R x I (4b)

where I is the idempotent, R = [absolute value of [bar.r] - [bar.r']] and [delta](t - R/[c.sub.[+ or -]]/4[pi]R is the retarded kernel for scalar potential with [c.sub.[+ or -]] the speed of light in the exterior or interior domain.

One can obtain a set of time-dependent integral equations (TDIE) by following the application of the boundary condition, viz., tangential continuity of fields. The continuity for tangential components is preserved by defining [J.sub.+]([bar.r], t) = -[J.sub.-]([bar.r], t) = J([bar.r], t) and [M.sub.+]([bar.r], t) = - [M.sub.-]([bar,r], t) = M([bar.r], t) in the limit [[OMEGA].sub.[+ or -]] [right arrow] [OMEGA]; note in this limit, [[??].sub.+] = -[[??].sub.-] [right arrow] [??]. Hence, no net artificial currents exists on the interface. In this work, the dependence of currents on position and time will be suppressed unless specified otherwise.

Let [??] x defined trace operator. Then, a linear combination of tangential trace of (1a) and tangential trace of (1c) gives the generalized electric field integral equation (EFIE);

[(1 + [alpha])I + [??] x ([K.sub.+] - [alpha][K.sub.-])] M - [[??] x ([L.sub.+] - [alpha][L.sub.-])] J = -[??] x [E.sup.i] (5)

Likewise, a linear combination of tangential trace of (1b) and tangential trace of (1d) would give a generalized magnetic field integral equation (MFIE);

[[??] x (1/[[eta].sup.2.sub.+] [L.sub.+] - [beta]/[[eta].sup.2.sub.-] [L.sub.- ])] M + [(1 + [beta])I + [??] x ([K.sub.+] - [beta][K.sub.-])] J = [??] x [H.sup.i] (6)

In the above equations, [alpha] and [beta] are appropriately chosen constants [50].

In a similar manner, if trace operator were defied as [??] x [??] x, one gets another set of TDIEs; viz.,

[(1 + [alpha])[??] x I + [??] x [??] x ([K.sub.+] - [alpha][K.sub.-])] M - [[??] x [??] x ([L.sub.+] - [alpha][L.sub.- ])] J = -[??] x [??] x [E.sup.i] (7)

[[??] x [??] x (1/[[eta].sup.2.sub.+] [L.sub.+] - [beta]/[[eta].sup.2.sub.-] [L.sub.-])] M + 1(1 + [beta]) [??] x I + [??] x [??] x ([K.sub.+] - [beta][K.sub.-])] J = [??] x [??] x [H.sup.i] (8)

The above two TDIEs can be regarded as widely used TD-PMCHWT or TD-Muller formulations for particular values for [alpha] and [beta], as in frequency domain [50]. There are other types of formulation that could be used to model the scattering from dielectric objects. All of the associated integral equation formulations might show slight differences in performance and accuracy when discretized. It should be noted that for spherical systems, vector spherical harmonics (VSH) form a natural basis to represent the currents. As shown in [51], using the addition theorem for dyadic Green's functions (the frequency domain counterparts of (4)), the use of VSH leads to a diagonal system. To this end, we seek a representation of the current density as follows:

J = [summation over (nm)] ([J.sup.1.sub.nm](t) [[PSI].sup.m.sub.n]([??]) + [J.sup.2.sub.nm](t) [[PHI].sup.m.sub.n]([??])) (9a)

M = [summation over (nm)] ([M.sup.1.sub.nm](t) [[PSI].sup.m.sub.n]([??]) + [M.sup.2.sub.nm](t) [[PHI].sup.m.sub.n]([??])) (9b)

where [summation over (nm)] is defined as two-fold summation [[infinity].summation over (n=0)][n.summation over (m= -n)]. The vector spherical harmonic functions in (9) are defined as

[[PSI].sup.m.sub.n]([??]) = [[PSI].sup.m.sub.n] ([theta], [phi]) = r/[square root of n(n + 1)] [[nabla].sup.t][Y.sup.m.sub.n] ([theta], [phi]) (10)

[[PHI].sup.m.sub.n]([??]) = [[PHI].sup.m.sub.n] ([theta], [phi]) = r/[square root of n(n + 1)] [[nabla].sup.t][Y.sup.m.sub.n] ([theta], [phi]) (11)

where

[Y.sup.m.sub.n] ([theta], [phi]) = [square root of 2n + 1/4[pi] (n - m)!/(n + m)!] [P.sup.m.sub.n] (cos [theta]) [e.sup.jm[phi]] (12)

Additionally, a set of orthogonal basis to represent the normal component of the field in terms of VSHs can be written as

Y([??]) = Y([theta], [phi]) = [??][Y.sup.m.sub.n]([theta], [phi]) (13)

In the above expressions, the unknown quantities, [Q.sup.i.sub.nm](t) for Q denoting either J/M, and i = 1,2 need to be recovered to provide the necessary transient behavior. Our approach to constructing these will proceed as follows: (i) use representations (9) in the appropriate TDIEs, (ii) use a novel addition theorem for the dyadic Green's function, and (iii) obtain a one-dimensional integral equation (as a function of time) for each of coefficients. The approach can be easily extended to both multiple spheres as well as multilayered spheres.

3. ADDITION THEOREMS FOR TIME DOMAIN GREEN'S FUNCTION

In this section, a set of addition theorems for time domain Green's function will be introduced. The spherical expansion for the time-domain scalar Green's function in free space will form the basis to construct the set of addition theorems for time domain dyadic Green's function. Two types of modified dyadic Green's functions are discussed next. The first relates current sources to trace components of the fields and will be used to solve the IEs. The second type relates the current sources and fields. Henceforth, these representation of the dyadic Green's functions will be referred to as current propagators and radiators, respectively.

3.1. Spherical Expansion for Time-Domain Scalar Green's Function

Like its variant in frequency domain, the addition theorem for the time domain scalar Green's function in free space can be expressed in terms of spherical harmonics [43] as

[delta](t - R/c)/4[pi]R = [summation over (nm)] c/2rr' [P.sub.n] ([r.sup.2] + [r'.sup.2] - [c.sup.2][t.sup.2]/2rr') [P.sub.[alpha][beta](t) [Y.sup.m.sub.n] ([theta]', [phi]') = [summation over (nm)] [c.sub.n] (t) [Y.sup.m.sub.n] ([theta], [phi]) [Y.sup.m*.sub.n] ([theta]', [phi]') (14)

where (r, [theta], [phi]) and (r', [theta]', [phi]') are the locations of source and observation points in spherical coordinate, respectively. Note, this expression is very different from those obtained using an inverse Fourier transform of expressions in [43]. The detailed derivation of (14) and its application in deriving Mie-like solution for transient scattering of acoustic waves from hard and soft spherical objects can be found in [43]. In (14), the time t lies in the range [[[absolute value of r - r']]/c, [absolute value of r + r']/c] [43], which is tantamount to constraining the argument of [P.sub.n](*) to within [-1,1]. This requirement corresponds to the causality of the scattering system. One can consider the Green's function in time domain as the impulse response [??](t) of the scattering system, and it vanishes when t < [absolute value of r - r']/c or t > [absolute value of r + r']/c, where the radial position for input and output signal are r and r', respectively.

Each term in the expansion can be seen as one mode kernel, which has time dependence as a Legendre polynomial with [t.sup.2] in its argument. For example, by assuming that [r.sub.1] = 3, [r.sub.2] = 5 and c = 1, Fig. 1 gives several plots of time-dependent coefficients cn(t) in the expansion of time-domain Green's function.

In addition to the fact that the above expansion is the time-domain counterpart of the conventional addition theorem for time-harmonic case, it also has following properties:

(i) It is symmetric in terms of source and observation position.

(ii) For fixed r and r', space and time dependence can be separated.

(iii) No singularity exists if it is used in integral equations.

(iv) The time-domain dependent coefficient [c.sub.n](t) in the expansion corresponds to the Green's function for the orthogonal spatial mode [Y.sup.m.sub.n].

(v) No explicit singular time-domain outgoing wave function is involved.

(vi) Operators associated with time or radial (such or normal derivatives) on the time domain Green's function can be done on the time-dependent coefficient [c.sub.n](t).

Property (i) is unique in time domain. Properties (ii) to (v) are essential to derive analytic or semi-analytic Mie type of solutions from integral equations of first kind (which is associated with single layer potential operator). Property (vi) can ensure the existence of addition theorem for [partial derivative]/[[partial derivative].sub.r] G(t, r, r') or [partial derivative]/[[partial derivative].sub.r] [partial derivative]/[[partial derivative].sub.r'] G(t, r.r') in system with spherical symmetry, and thus it can lead to analytic or semi- analytic Mie type of solutions from integral equations involving those kernels with higher order singularity.

The above addition theorem provides the basis of deriving the spherical expansion for time-domain dyadic Green's function. It should be noted that the (14) is not new; it was first presented in [40] when the authors tried using this expression to derive analytical expression for the coefficients.

3.2. Spherical Expansion for Tangential Trace of Dyadic Green's Function

As in scalar case, the dyadic Green's function in time domain can be also expanded in the form of addition theorem [45]. As shown in [45], when working with surface integral equations it is easier to work with a modified version of dyadic Green's function specifically its tangential trace form [G.sub.e0]. A quick review of the derivation is given below.

The tangential trace forms of frequency-domain dyadic Green's functions of electric and magnetic types can be written as

[mathematical expression not reproducible] (15)

and

[mathematical expression not reproducible] (16)

In the above expansions, it is assumed r > r', and [N.sup.(p).sub.nm] (k, [bar.r]) and [M.sup.(p).sub.nm](k, [bar.r]) are vector spherical wave functions (VSWF) of degree n and order m, respectively, [z.sup.(1).sub.n](kr') is the spherical Bessel funcition and [z.sup.(4).sub.n] (kr) denotes the spherical Hankel function of second type with k being the wave number in the medium. The above derivations used the relations between vector spherical harmonics and vector spherical wave functions [45].

Inverse Fourier transforms of (15) and (16) lead to time-domain additions theorems for tangential form of dyadic Green's functions that can be written as

[mathematical expression not reproducible] (17)

and

[mathematical expression not reproducible] (18)

where the following auxiliary function [K.sup.(0).sub.n] is defined as

[K.sup.(0).sub.n] (r,r',t) [equivalent to] [F.sup.-1] [-jk[z.sup.(4).sub.n] (kr) [z.sup.(1)*.sub.n] (kr')] = c/2rr' [P.sub.n] ([r.sup.2] + [r'.sup.2] - [c.sup.2][t.sup.2]/2rr')[P.sub.[alpha][beta]](t) (19)

It should be noted that one can derive the above addition theorems for the modified time-domain dyadic Green's function using the approach given in [43], though a little more complicated than the one given here based on the result of scalar case. Spherical expansions for time-domain modified dyadic Green's function share almost all of the same properties associated with the time-dependent addition theorem of scalar case.

3.3. A Simplified Current-Field Radiator in Spherical Coordinate System

The addition theorems presented above are sufficient to solve the TDIEs and obtain the current coefficients for each harmonic. Furthermore, it is possible to use these expressions to find trace fields on any spherical surface (say for scattering from multiple spheres or spheres from with multiple layers). What is also necessary is to compute radiated fields wherein one may need all components of both the electric and magnetic fields. To derive the field propagators, we start with

[mathematical expression not reproducible] (20)

and

[mathematical expression not reproducible] (21)

Transforming the above two relations results in the two corresponding field propagators:

[mathematical expression not reproducible] (22)

and

[mathematical expression not reproducible] (23)

Thus, given the equivalent current density on the interface, one can construct the radiated electric or magnetic field through (1a) or (1b), respectively. As defined in Section 2, the current densities are expanded in terms of vector spherical harmonics. Consequently, considering the orthogonality properties of VSH, the spatial integral can be evaluated analytically, thus reducing the four dimensional spatial-temporal convolution between source currents and the time-domain dyadic Green's function into a weighted sum of a set of temporal convolutions. One can evaluate the electric field observed at the spherical surface with radius r due to the currents on the dielectric scatterer (with radius r'). Without losing generality, assuming that [[epsilon].sub.+] = [[epsilon].sub.0], the scattered fields can be written as

[mathematical expression not reproducible] (24)

where the property [[??].sub.e/m0] * X = -[[??].sup.ft.sub.e/m0] * X is used. Substituting for the spherical expansion for both currents densities and the modified dyadic Green's functions, we get

[mathematical expression not reproducible] (25)

with equivalent time domain kernels defined as

[K.sup.1.sub.n](r, r', t) = -[mu] c[[partial derivative].sup.-1.sub.t]/rr' [[partial derivative].sub.r][[partial derivative].sub.r'] [rr' [K.sup.(0).sub.n] (r, r', t)] (26a)

[K.sup.2.sub.n](r, r', t) = -[mu]/c [partial derivative]/[partial derivative]t [K.sup.(0).sub.n] (r, r', t)] (26b)

[K.sup.3.sub.n](r, r', t) = -[[partial derivative].sub.r']/r' [r' [K.sup.(0).sub.n] (r, r', t)] (26c)

[K.sup.4.sub.n](r, r', t) = -[[partial derivative].sub.r]/r [r [K.sup.(0).sub.n] (r, r', t)] (26d)

[K.sup.5.sub.n](r, r', t) = [square root of n(n + 1)] [mu] c[[partial derivative].sup.-1.sub.t]/rr' [[partial derivative].sub.r'] [r' K.sup.(0).sub.n] (r, r', t)] (26e)

[K.sup.6.sub.n](r, r', t) = [square root of n(n + 1)] [K.sup.0.sub.n] (r, r', t)/r (26f)

It follows, that one can use a similar procedure to map the fields on [partial derivative][[OMEGA].sub.+] to fields on the boundary of any spherical surface. Thus, we have expressions for two propagators. The first can be used to solve TDIEs (as examined next) and the second can be used to post-process data.

4. TIME-DEPENDENT LORENTZ-MIE-DEBYE SERIES SOLUTION

In this section, TDIEs governing the time-dependent scattering from dielectric spheres are to be solved. The spatial expansion mentioned in Section 2 and the addition theorems in Section 3 will be used to derive the reduced semi-discrete system. Afterwards, the reduced system are to be solved using an hp-discontinuous Galerkin approach.

4.1. Reduced System

As in usual discretization of integral equations of time-harmonic case, after substituting the representations of currents defined in (9) into (5) and (7) and using expansion for the tangential trace of dyadic Green's function, one can test the TDIEs with [[PSI].sup.*.sub.nm] and [[PHI].sup.*.sub.nm], respectively, using the Galerkin testing scheme. The following reduced system can be obtained for each mode set [[PSI].sub.nm] and [[PHI].sub.nm]:

[(1 + [alpha])[delta](t) + ([K.sup.4.sub.+nm] - [alpha][K.sup.4.sub.+nm])] [*.sub.t] [M.sup.1.sub.nm] - [[K.sup.2.sub.+nm] - [alpha][K.sup.2.sub.+nm]] [*.sub.t] [J.sup.2.sub.nm] = [f.sub.1] (27a)

[1.5mm] [(1 + [alpha])[delta](t) + ([K.sup.3.sub.+nm] - [alpha][K.sup.3.sub.+nm])] [*.sub.t] [M.sup.2.sub.nm] - [[K.sup.1.sub.+nm] - [alpha][K.sup.1.sub.+nm]] [*.sub.t] [J.sup.1.sub.nm] = [f.sub.2] (27b)

[1.5mm] [[K.sup.2.sub.+nm]/[[eta].sup.2.sub.+] - [beta][K.sup.2.sub._nm]/[[eta].sup.2.sub.+]] [*.sub.t] [M.sup.2.sub.nm] + [(1 + [beta])[delta](t) + ([K.sup.4.sub.+nm] - [beta][K.sup.4.sub.+nm])] [*.sub.t] [J.sup.1.sub.nm] = [f.sup.3] (27c)

[1.5mm] [[K.sup.1.sub.+nm]/[[eta].sup.2.sub.+] - [beta][K.sup.1.sub._nm]/[[eta].sup.2.sub.+]] [*.sub.t] [M.sup.1.sub.nm] + [(1 + [beta])[delta](t) + ([K.sup.3.sub.+nm] - [beta][K.sup.3.sub.+nm])] [*.sub.t] [J.sup.2.sub.nm] = [f.sup.4] (27d)

where the dependence on time is suppressed, and the right hand side terms can be written as

[f.sup.1] = - [[integral].sub.S'] [[PSI].sup.*.sub.nm] * ([??] x [E.sup.i])dS'; [f.sup.2] = - [[integral].sub.S'] [[PHI].sup.*.sub.nm] * ([??] x [E.sup.i])dS'; [f.sup.3] = - [[integral].sub.S'] [[PSI].sup.*.sub.nm] * ([??] x [H.sup.i])dS'; [f.sup.4] = - [[integral].sub.S'] [[PSI].sup.*.sub.nm] * ([??] x [H.sup.i])dS'. (28)

Compared to the PEC case in [45], the above reduced system is coupled with four unknowns instead of one. After careful observation, (27a) and (27d) are independent from (27b) and (27c); in another words, the above system can be further reduced into two smaller systems. All the reductions of the semidiscrete version of original TDIEs are based on the orthogonality of VSH, and the spatial convolution is taken care of analytically. Again, no radial argument is involved during the reduction. It is apparent from the above equations that a direct Mie type analytical solution to these coefficients is difficult (if not impossible).

If the tangential component of excitation can be represented using

[??] x [E.sup.i] (t, [bar.r]) = [summation over (nm)] ([E.sup.1.sub.nm](t)[[PSI].sub.nm] + [E.sup.2.sub.nm](t)[[PHI].sub.nm]) and [??] x [H.sup.i] (t, [bar.r]) = [summation over (nm)] ([H.sup.1.sub.nm](t) [[PSI].sub.nm] + [H.sup.2.sub.nm](t) [[PHI].sub.nm]]),

the integrals in (28) can be done analytically. The summations need to be truncated so that finite number of reduced equations are to be solved. Fortunately, if the incident field is approximately band-limited in both space and time, then the expansion in the excitation and field representation be truncated, and one can develop error bounds that can be correlated to truncation. Typically, the highest degree of the modes chosen is proportional to the electric size of the sphere being studied, very similar to the time-harmonic case. For transient systems, this amounts to the largest electrical size of the object. For pulses, the coefficients [E.sup.1.sub.nm] and [E.sup.2.sub.nm] are obtained using the properties of [[PSI].sub.nm] and [[PHI].sub.nm].

4.2. Marching-on-in-Time Solution

The current coefficients, [J.sup.i.sub.nm], [M.sup.i.sub.nm] for i = 1,2 are obtained using temporal deconvolution, i.e., solve the reduced system in (27a)-(27d). It is done first by representing the time- dependent coefficients using temporal basis functions, such as shifted Lagrange polynomials, piecewise Legendre polynomials or approximate prolate function. Then, using either Galerkin testing or collocation, one can get a discrete convolution system. The procedure is similar to the one given in [43, 45]. Minor differences include that (1) the kernels involved are linear combinations of kernels corresponding to the exterior and interior domains and (2) the spatial system dimension is only two by two for homogeneous dielectric problem instead of one by one in perfectly conducting case.

Since the system is greatly decoupled (hence the size for the deconvolution problem is greatly decreased), marching-on-in-time or time-stepping approach [52] is easy to be applied here. While details are not given in this paper, readers are referred to [43,45] for a more complete description. It has to be noted that the stable MOT solutions for four different kernels have been obtained for conducting sphere case [45].

5. NUMERICAL EXAMPLES

In this section, a simple numerical example is given to test the formulation derived earlier for electromagnetic scattering from a dielectric sphere of unit radius. The incidence field is a modulated Gaussian pulse [mathematical expression not reproducible], where [mathematical expression not reproducible] denotes the direction of propagation, [sigma] = 3/(2[pi] B) with B denoting the bandwidth, and [t.sub.p] = 20 [sigma]. Center frequency and bandwidth are set as 0.5 GHz and 0.3 GHz respectively in the example. Fig. 2(a) plots the time-dependent coefficients ([J.sup.1](t) and [J.sup.2](t) respectively) for the modes [[PSI].sub.41] and [[PHI].sub.41] in the electric current density, and Fig. 2(b) plots the time-dependent coefficients ([M.sup.1](t) and [M.sup.2](t) respectively) for the same two modes in the magnetic current density.

The current density values are transformed into Fourier domain to compare with frequency-domain Lorentz-Mie-Debye series. Fig. 3(a) plots the spectrum of the transient responses (magnitude of frequency domain coefficients) for the modes [[PSI].sub.41] and [[PHI].sub.41] in the electric current density, along with the results from frequency domain Lorentz-Mie-Debye series. Similarly, Fig. 3(b) plots the comparison in the magnetic current density data. From the comparisons, the agreement between timed-dependent Lorentz-Mie-Debye series based solution and its Fourier domain counterpart is evident.

6. DISCUSSIONS AND CONCLUSIONS

In this contribution, we construct a time-dependent Lorentz-Mie-Debye formulation for dielectric spherical objects based on the new set of addition theorems for time domain dyadic Green's function. With the new current and field propagators, the presented approach can help build self-consistent time-domain spherical multipole theory, which can be widely applied in both the microwave and optics regime.

Furthermore, it can be easily verified that the above formulation can be trivially extended to solve scattering from multiple spheres as well as multilayered spheres. For both the cases, at each interface there would be a set of IEs like (5) and (6). Each interface will have an unknown function associated with it to represent the time signature of the currents on it. For concentrically layered case, the spatial coupling will be in block diagonal form. When analyzing multiple spheres, the spatial system would be fully coupled except the block corresponding to each interface. In this case, fast method like plane-wave time-domain [53] can be used to reduce the memory and time complexity.

Another example of applications of time-dependent Mie formulation includes an exact nonreflection global boundary condition [54, 55] for partial differential equation based simulation. That is straightforward in time-harmonic case. The open domain could be truncated within a spherical artificial domain, on which boundary integral is employed to satisfy the Sommerfeld radiation boundary condition. This observation allows the possible applications in many more complicated problems with non-canonical shapes, inhomogeneous materials or nonlinear behaviors.

Finally, an interesting discussion item is the stability of the time marching system. This has been numerically demonstrated in [45]. However, ongoing work is focused on studying late time stability for this system and extending it to analyzing stabilities of TDIEs as applied to electromagnetics--a long standing problem where significant literature exists [56-62].

ACKNOWLEDGMENT

The authors wish to acknowledge the HPCC facility at Michigan State University and support from NSF CMMI-1250261, and support from the Department of Defense CREATE-RF program. B. Shanker is grateful for support from the Department of Defense High Performance Computing Modernization Program Office (HPCMPO) and its Computational Research and Engineering for Acquisition Tools and Environments (CREATE) Program.

REFERENCES

[1.] Maxwell, J. C., "A dynamical theory of the electromagnetic field," Philosophical Transactions of the Royal Society of London, Vol. 155, 459-512, 1865.

[2.] Uslenghi, P. L. E., "Scattering by an impedance sphere coated with a chiral layer," Electromagnetics, Vol. 10, 201-211, Jan. 1990.

[3.] Mie, G., "Beitrage zur Optik truber Medien, speziell kolloidaler Metallosungen," Annalen der Physik, Vol. 330, No. 3, 377-445, 1908.

[4.] Lorenz, L., "Lysbevaegelsen i og uden for en af plane Lysbolger belyst Kugle," Det Kongelige Danske Videnskabernes Selskabs Skrifter, Vol. 6, No. 6, 1-62, 1890.

[5.] Debye, P., "Der Lichtdruck auf Kugeln von beliebigem material," Annalen der Physik, Vol. 30, No. 1, 57-136, 1909.

[6.] Logan, N., "Survey of some early studies of the scattering of plane waves by a sphere," Proceedings of the IEEE, Vol. 53, 773-785, Aug. 1965.

[7.] Clebsch, A., "Ueber die reflexion an einer Kugelflache," Journal fur die Reine Und Angewandte Mathematik, Vol. 61, 195-262, 1863.

[8.] Nicholson, J., "XVII. On the diffraction of short waves by a rigid sphere," Philosophical Magazine Series 6, Vol. 11, 193-205, Feb. 1906.

[9.] Bromwich, T., "X. Electromagnetic waves," Philosophical Magazine Series 6, Vol. 38, 143-164, Jul. 1919.

[10.] Proudman, J., A. T. Doodson, and G. Kennedy, "Numerical results of the theory of the diffraction of a plane electromagnetic wave by a perfectly conducting sphere," Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, Vol. 217, 279- 314, Jan. 1918.

[11.] White, F. P., "The diffraction of plane electromagnetic waves by a perfectly reflecting sphere," Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, Vol. 100, 505-525, Feb. 1922.

[12.] Wriedt, T., "Mie theory: A review," The Mie Theory, W. Hergert and T. Wriedt, eds., Vol. 169, 53-71, Springer Series in Optical Sciences, Springer Berlin Heidelberg, 2012.

[13.] Born, M. and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffaction of Light, Pergamon Press, 1959.

[14.] Lock, J. A. and G. Gouesbet, "Generalized lorenzmie theory and applications," Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 110, No. 11, 800-807, Light Scattering: Mie and More Commemorating 100 years of Mie's 1908 publication, 2009.

[15.] Aden, A. L. and M. Kerker, "Scattering of electromagnetic waves from two concentric spheres," Journal of Applied Physics, Vol. 22, No. 10, 1951.

[16.] Geng, Y.-L., X.-B. Wu, L.-W. Li, and B.-R. Guan, "Mie scattering by a uniaxial anisotropic sphere," Phys. Rev. E, Vol. 70, 056609, Nov. 2004.

[17.] Bohren, C. F., "Light scattering by an optically active sphere," Chemical Physics Letters, Vol. 29, 458-462, Dec. 1974.

[18.] Xu, F., K. Ren, G. Gouesbet, G. Grehan, and X. Cai, "Generalized lorenz- mie theory for an arbitrarily oriented, located, and shaped beam scattered by a homogeneous spheroid," J. Opt. Soc. Am. A, Vol. 24, 119-131, Jan. 2007.

[19.] Xu, Y.-L. and B. Gustafson, "A generalized multiparticle mie-solution: Further experimental verification," Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 70, No. 4-6, 395-419, Light Scattering by Non-Spherical Particles, 2001.

[20.] Hough, J. and G. Gouesbet, "Generalized LorenzMie theories, the third decade: A perspective," Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 110, No. 14, 1223-1238, 2009.

[21.] Han, Y., G. Grehan, and G. Gouesbet, "Generalized lorenz-mie theory for a spheroidal particle with off-axis gaussian-beam illumination," Appl. Opt., Vol. 42, 6621-6629, Nov. 2003.

[22.] Waterman, P. C., "New formulation of acoustic scattering," The Journal of the Acoustical Society of America, Vol. 45, No. 6, 1417, 1969.

[23.] Strom, S., "T matrix for electromagnetic scattering from an arbitrary number of scatterers with continuously varying electromagnetic properties," Phys. Rev. D, Vol. 10, 2685- 2690, Oct. 1974.

[24.] Strom, S. and W. Zheng, "The null field approach to electromagnetic scattering from composite objects," IEEE Transactions on Antennas and Propagation, Vol. 36, 376-383, Mar. 1988.

[25.] Waterman, P. C., "The T-matrix revisited," J. Opt. Soc. Am. A, Vol. 24, 2257-2267, Aug. 2007.

[26.] Nilsson, A. M., P. Alsholm, A. Karlsson, and S. Andersson-Engels, "T- matrix computations of light scattering by red blood cells," Appl. Opt., Vol. 37, 2735-2748, May 1998.

[27.] Mishchenko, M. I., G. Videen, N. G. Khlebtsov, and T. Wriedt, "Comprehensive T-matrix reference database: A 20122013 update," Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 123, 145-152, 2013.

[28.] Mishchenko, M. I., M. Kahnert, D. W. Mackowski, and T. Wriedt, "Peter Waterman and his scientific legacy," Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 123, 1, 2013.

[29.] Mishchenko, M. I. and P. A. Martin, "Peter Waterman and T-matrix methods," Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 123, 2-7, 2013.

[30.] Mishchenko, M. I., L. Liu, and D. W. Mackowski, "T-matrix modeling of linear depolarization by morphologically complex soot and soot-containing aerosols," Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 123, 135-144, 2013.

[31.] Khlebtsov, N. G., "T-matrix method in plasmonics: An overview," Journal of Quantitative Spectroscopy and Radiative Transfer, Vol. 123, 184-217, 2013.

[32.] Wang, Y. and W. Chew, "A recursive t-matrix approach for the solution of electromagnetic scattering by many spheres," IEEE Transactions on Antennas and Propagation, Vol. 41, 1633-1639, Dec. 1993.

[33.] Abraham, M. and R. Becker, The Classical Theory of Electricity and Magnetism, Blackie & Son Ltd., London, 1937.

[34.] Wicklund, J., "Extrapolation of the electromagnetic field," Tech. Rep., Diamond Ordnance Fuze Laboratories, 1962.

[35.] Granzow, K. D., "Multipole theory in the time domain," Journal of Mathematical Physics, Vol. 7, 634, Dec. 1966.

[36.] Granzow, K. D., "Time-domain treatment of a spherical boundary-value problem," Journal of Applied Physics, Vol. 39, 3435, Nov. 1968.

[37.] Davidon, W. C., "Time-dependent multipole analysis," Journal of Physics A: Mathematical, Nuclear and General, Vol. 6, 1635-1646, Nov. 1973.

[38.] Campbell, W., J. Macek, and T. Morgan, "Relativistic time-dependent multipole analysis for scalar, electromagnetic, and gravitational fields," Phys. Rev. D, Vol. 15, 2156-2164, Apr. 1977.

[39.] Buyukdura, O. M. and S. S. Koc, "Two alternative expressions for the spherical wave expansion of the time domain scalar free-space Green's function and an application: Scattering by a soft sphere," Journal of the Acoustical Society of America, Vol. 101, No. 1, 87-91, 1997.

[40.] Azizoglu, S., S. Koc, and O. Buyukdura, "Spherical wave expansion of the time-domain free- space dyadic Green's function," IEEE Transactions on Antennas and Propagation, Vol. 52, 677-683, Mar. 2004.

[41.] Sauter, S. and A. Veit, "Retarded boundary integral equations on the sphere: Exact and numerical solution," IMA Journal of Numerical Analysis, 2013.

[42.] Greengard, L., T. Hagstrom, and S. Jiang, "The solution of the scalar wave equation in the exterior of a sphere," Journal of Computational Physics, Vol. 274, 191-207, Oct. 2014.

[43.] Li, J., D. Dault, and B. Shanker, "A quasianalytical time domain solution for scattering from a homogeneous sphere," The Journal of the Acoustical Society of America, Vol. 135, 1676-1685, Apr. 2014.

[44.] Greengard, L., T. Hagstrom, and S. Jiang, "Extension of the lorenzmiedebye method for electromagnetic scattering to the time-domain," Journal of Computational Physics, Vol. 299, 98-105, 2015.

[45.] Li, J. and B. Shanker, "Time-dependent debye-mie series solutions for electromagnetic scattering," IEEE Transactions on Antennas and Propagation, Vol. 63, 3644-3653, Aug. 2015.

[46.] Hovenac, E. A. and J. A. Lock, "Assessing the contributions of surface waves and complex rays to far-field mie scattering by use of the debye series," J. Opt. Soc. Am. A, Vol. 9, 781-795, May 1992.

[47.] Lock, J. A. and P. Laven, "Mie scattering in the time domain. Part 1. The role of surface waves," J. Opt. Soc. Am. A, Vol. 28, 1086-1095, Jun. 2011.

[48.] Lock, J. A. and P. Laven, "Mie scattering in the time domain. Part II. The role of diffraction," J. Opt. Soc. Am. A, Vol. 28, 1096-1106, Jun. 2011.

[49.] Li, R., X. Han, H. Jiang, and K. F. Ren, "Debye series for light scattering by a multilayered sphere," Appl. Opt., Vol. 45, 1260-1270, Feb. 2006.

[50.] Li, J., D. Dault, N. Nair, and B. Shanker, "Analysis of scattering from complex dielectric objects using the generalized method of moments," J. Opt. Soc. Am. A, Vol. 31, 2346- 2355, Nov. 2014.

[51.] Hsiao, G. and R. Kleinman, "Mathematical foundations for error estimation in numerical solutions of integral equations in electromagnetics," IEEE Transactions on Antennas and Propagation, Vol. 45, 316-328, Mar. 1997.

[52.] Shanker, B., A. Ergin, K. Aygun, and E. Michielssen, "Analysis of transient electromagnetic scattering from closed surfaces using a combined field integral equation," IEEE Transactions on Antennas and Propagation, Vol. 48, 1064-1074, Jul. 2000.

[53.] Shanker, B., A. Ergin, M. Lu, and E. Michielssen, "Fast analysis of transient electromagnetic scattering phenomena using the multilevel plane wave time domain algorithm," IEEE Transactions on Antennas and Propagation, Vol. 51, 628-641, Mar. 2003.

[54.] Grote, M. J. and J. B. Keller, "Exact nonreflecting boundary conditions for the time dependent wave equation," SIAM Journal on Applied Mathematics, Vol. 55, No. 2, 280-297, 1995.

[55.] Alpert, B., L. Greengard, and T. Hagstrom, "Nonreflecting boundary conditions for the timedependent wave equation," J. Comput. Phys., Vol. 180, 270-296, 2002.

[56.] Tijhuis, A. G., "Toward a stable marching-on-in-time method for two- dimensional transient electromagnetic scattering problems," Radio Science, Vol. 19, No. 5, 1311-1317, 1984.

[57.] Sadigh, A. and E. Arvas, "Treating the instabilities in marching-on-in- time method from a different perspective [electromagnetic scattering]," IEEE Transactions on Antennas and Propagation, Vol. 41, 1695-1702, Dec. 1993.

[58.] Ha-Duong, T., B. Ludwig, and I. Terrasse, "A Galerkin bem for transient acoustic scattering by an absorbing obstacle," International Journal for Numerical Methods in Engineering, Vol. 57, No. 13, 1845-1882, 2003.

[59.] Wang, X., R. Wildman, D. S. Weile, and P. Monk, "A finite difference delay modeling approach to the discretization of the time domain integral equations of electromagnetics," IEEE Transactions on Antennas and Propagation, Vol. 56, 2442-2452, Aug. 2008.

[60.] Pray, A., N. Nair, and B. Shanker, "Stability properties of the time domain electric field integral equation using a separable approximation for the convolution with the retarded potential," IEEE Transactions on Antennas and Propagation, Vol. 60, 3772-3781, Aug. 2012.

[61.] Van't Wout, E., D. van der Heul, H. van der Ven, and C. Vuik, "The influence of the exact evaluation of radiation fields in finite precision arithmetic on the stability of the time domain integral equation method," IEEE Transactions on Antennas and Propagation, Vol. 61, 6064-6074, Dec. 2013.

[62.] Van't Wout, E., D. R. van der Heul, H. van der Ven, and C. Vuik, "Stability analysis of the marching-on-in-time boundary element method for electromagnetics," Journal of Computational and Applied Mathematics, Vol. 294, 358-371, 2016.

Jie Li (1) and Balasubramaniam Shanker (2), *

Received 14 December 2015, Accepted 31 December 2015, Scheduled 5 January 2016

Invited paper for the Commemorative Collection on the 150-Year Anniversary of Maxwell's Equations.

* Corresponding author: Balasubramaniam Shanker (bshanker@msu.edu).

(1) Department of Electrical and Computer Engineering, Michigan State University, USA. (2) Department of Physics and Astronomy, Michigan State University, USA.

Caption: Figure 1. Time-dependent coefficients with degree (a) n = 0, (b) n = 3 and (c) n = 10. ([r.sub.1] = 3, [r.sub.2] = 5 and c = 1.)

Caption: Figure 2. Transient response (coefficients) of modes ^41 and \$41 in (a) electric currents and (b) magnetic currents. "Re" and "Im" denote real and imaginary part, coefficients that are not listed are zero.

Caption: Figure 3. Spectrum of temporal responses of modes [[PSI].sub.41] and [[PHI].sub.41] due to the modulated Gaussian pulse incidence. (a) Electric currents and (b) magnetic currents. The curves from Fourier domain and time-dependent Lorentz-Mie-Debye formulation are denoted by "FD" and "TD" respectively.