On the Computation of Power in Volume Integral Equation Formulations

# On the Computation of Power in Volume Integral Equation Formulations

## Abstract

We present simple and stable formulas for computing power (including absorbed/radiated, scattered and extinction power) in current-based volume integral equation formulations. The proposed formulas are given in terms of vector-matrix-vector products of quantities found solely in the associated linear system. In addition to their efficiency, the derived expressions can guarantee the positivity of the computed power. We also discuss the application of Poynting’s theorem for the case of sources immersed in dissipative materials. The formulas are validated against results obtained both with analytical and numerical methods for scattering and radiation benchmark cases.

{keywords}

Electromagnetic scattering, method of moments (MoM), numerical analysis, Poynting’s theorem, volume integral equations.

## I Introduction

Volume integral equation (VIE) formulations have been extensively used over the last decades for the numerical solution of electromagnetic (EM) scattering and radiation problems (here is a non-exhaustive list of references [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]). Admittedly, VIE methods do not hold the workhorse status of their early days, since partial differential equation (PDE)-based methods, such as the finite-difference and finite-element methods, reached a certain level of maturity. Nevertheless, some recent insights have spawned renewed interest in the development of more competitive VIE methods [17, 18, 19, 20, 21, 22, 23].

The objective of this paper is to introduce concise and computationally efficient formulas for the power absorbed, scattered, and radiated by bodies modeled with VIE methods. Of course, in the VIE or any other scattering formalism these quantities could be computed simply by direct numerical cubature of the Poynting vector over appropriate bounding surfaces, an approach we term the “Poynting method” (PM); however, in many cases the PM may not be the best way to perform numerical power computations. One reason is that, in methods such as the current-based VIE (JM-VIE) which solve directly for (volumetric) sources rather than fields, computation of the Poynting vector requires an extra post-processing step to compute the fields at each cubature point. A more urgent problem is that, in many cases, the Poynting-vector cubature over the bounding surface may be badly behaved due to large cancellations from different surface regions, requiring large numbers of cubature points to obtain decent accuracy. This difficulty may be mitigated by using a distant bounding surface far removed from the surface of the scatterers, but this strategy is unavailable in many cases of interest, such as geometries involving interleaved bodies or problems in which we need the power absorbed by just one of two nearby objects.

These issues were recently discussed in the context of surface integral equation (SIE) solvers [24, 25]. More specifically, Ref. 24 noted that the power absorbed by a homogeneous body may be computed solely from knowledge of the tangential currents flowing on its surface. Ref. 25 extended this observation by noting that in fact the surface currents suffice to determine not only the absorbed power but also the scattered power, total power (extinction), force (radiation pressure), and torque on a homogeneous body. Indeed, as discussed in [25], all of these quantities may be expressed compactly as vector-matrix-vector product (VMVP) expressions of the form , where is the vector of surface-current basis-function coefficients obtained in the SIE solution to a scattering problem and is a matrix which assumes different forms for the various different quantities we are calculating. In this work we extend the developments of [25] to the VIE domain, deriving VMVP formulas for the absorbed, scattered, radiated, and total power. Our formulas are analogous to those of [25] in that they compute the power directly from the volume-current coefficient vector (the solution to the linear VIE system), thus bypassing the post-processing step of computing scattered fields and Poynting’s vectors. An advantage of our formulas over their SIE counterparts is that, as we prove rigorously in Section IIIE, the numerical predictions that they yield obey the key physical requirement of positivity of the power absorbed, scattered1, and radiated by passive material bodies, even in their discretized form. In contrast, while the SIE power formulas derived in [25] (see also [26]) are positive in the exact Maxwell equations, this positivity relied on a delicate cancellation that can sometimes break down when they are approximated by a discretized basis (such as boundary elements).

As noted above, the VMVP formulas derived previously in the SIE context for classical scattering [25] have also proven useful for numerical modeling of quantum/statistical-mechanical phenomena [26, 27, 28]. In direct analogy, the formulas presented here for classical problems are key building blocks for an efficient VIE-based numerical approach to computational fluctuation physics; in particular, our VMVP formulas for the absorbed power (36) and the radiated power (54) may be extended to matrix-trace formulas for rates of radiative heat transfer and fluorescence. This correspondence will be addressed in future work.

The remainder of this paper is organized as follows. In Section II we set up the JM-VIE formulation and the associated linear system. In addition, we describe the incorporation of dipole sources in JM-VIE solvers. In Section III, we present the main results of this work: the boxed VMVP expressions for the computation of power in VIE methods, and we prove their positivity. Finally, in Section IV we validate our formulas, and we demonstrate some of their useful properties. Table I lists some notation used in this work.

## Ii Volume integral equations

### Ii-a Formulation

We consider the scattering of time-harmonic EM waves by a penetrable object, occupying the bounded domain in 3-D Euclidean space, . The working angular frequency is and the electric properties are defined as

 ϵ=ϵ0,μ=μ0inR3∖Ω;ϵ=ϵr(r)ϵ0,μ=μr(r)μ0inΩ (1)

Here, the vacuum (or free-space) permittivity and permeability are real positive values, while the relative permittivities and read

 ϵr(r)=ϵ′r(r)−iϵ′′r(r)μr(r)=μ′r(r)−iμ′′r(r) (2)

with and , assuming a time factor .

The total time harmonic fields () in the presence of an isotropic inhomogeneous object can be expressed in terms of equivalent polarization and magnetization currents (), as follows (dropping some function arguments where no confusion exists):

 (eh)=(einchinc)+(escahsca) (3)

where the incident fields () are the fields generated by sources in the absence of the scatterer and the scattered fields () are given by2

 (escahsca)=(1ceL−KK1cmL)(jm)=(1ce(N−I)−KK1cm(N−I))Asca(jm) (4)

where . is simply the convolution operator with the Green function connecting currents to fields in vacuum. More explicitly, the associated integro-differential operators are

 Lf ≜(k20+∇∇⋅)V(f) (5a) Kf ≜∇×V(f) (5b) Nf ≜∇×∇×V(f) (5c)

where

 V(f)≜∫Ωg(r−r′)f(r′)d3r′ (6)

is the volume vector potential and is the fundamental Helmholtz solution,

 g(r)=e−ik0|r|4π|r| (7)

with being the wavenumber in free-space. Also, the equivalent current densities are defined in terms of the fields as follows:

 j(r) ≜ceχe(r)e(r) (8a) m(r) ≜cmχm(r)h(r). (8b)

where

 χe≜ϵr−1,χm≜μr−1 (9)

is the electric and magnetic susceptibility, respectively. Finally, the JM-VIE formulation can be derived by combining (3), (4) and (8) [19, 22],

 (10)

where

 A=(ANeAKe−AKmANm)=(Mϵr−MχeNceMχeK−cmMχmKMμr−MχmN) (11)

and

 (12)

are multiplication operators that multiply by the local parameter functions , while is the identity dyadic tensor.

### Ii-B Linear System

Usually, JM-VIE formulations are numerically solved by means of a Galerkin method, where the equivalent volumetric currents are approximated as expansions in some discrete set of vector-valued square-integrable basis functions, e.g. as in [18, 19, 22]:

 j≈∑αxe,αpα,m≈∑αxm,αpα (13)

The linear system arising from the Galerkin “testing”, i.e. with and , reads

 (ANeAKe−AKmANm)(xexm)=(bebm) (14)

where

 A=(ANeAKe−AKmANm)=(MϵrG−MχeNceMχeK−cmMχmKMμrG−MχmN) (15)

with

 Nαβ =⟨pα,Npβ⟩% \tiny V (16a) Kαβ =⟨pα,Kpβ⟩% \tiny V (16b)

and is the Gram matrix, given by

 Gαβ=⟨pα,pβ⟩% \tiny V. (17)

Also, and are the discrete versions of the associated operators. More specifically, matrices are diagonal for isotropic material with the non-zero values being equal to the material properties at the corresponding element, while matrix is the identity matrix with multiplication pre-factors and for the two diagonal sub-blocks, respectively. Note that the equivalent currents have no continuity constraints (at the interface of the elements) and the support of the basis/testing functions is restricted to single elements. Hence, the associated Gram matrix is diagonal, when non-overlapping basis functions are used [19, 22].

In the above we have used the inner product:

 ⟨f,g⟩\tiny V=∫Ω¯¯¯f⋅g. (18)

Finally, the right-hand side in (14) is given by

 b=(bebm)=CMχ(einchinc) (19)

where

 einc,α =⟨pα,einc⟩\tiny V (20a) hinc,α =⟨pα,hinc⟩\tiny V. (20b)

### Ii-C Dipole Excitation

In most radiation problems, we need to deal with elementary excitations, such as electric (or magnetic) Hertzian oscillating dipoles, . Numerically, the finite discretization means that the 3D dirac delta will have support only within a voxel of size , and will be given by

 δ△V(r)=p△V={1△V, if r∈supp(p)0,otherwise (21)

with being the magnitude of the basis functions for the approximation of the polarization/magnetization densities.

The use of dipole sources in our analysis could prove quite problematic, since current sources are inherently modeled in VIE as uniform distributions throughout the volume elements. Motivated by the limit , we introduce the notion of the distributed-dipole (DD) source. Specifically, the current of a DD source immersed in element , along the direction of the basis reads

 σα(x)=pαVα. (22)

Hence, the impressed current vector is given by

 (dedm)=^G−1(pepm) (23)

where

 ^G=(G00G). (24)

The non-zero elements of the vector in (23) depend solely on the location and direction of the dipole sources under consideration.

Finally, the incident fields of (19) read

 (einchinc)=Asca(dedm) (25)

where is the discrete form of the operator in (4):

 Asca=(1ce(N−G)−KK1cm(N−G)). (26)

In the case of dipole sources located outside the scatterer, the computation of the incident fields is quite straightforward, since there are no singularities in the fields (inside the scatterer). More specifically, the impressed sources are propagated by means of the free-space Green function, and the incident fields are “tested” as in (20).

## Iii Power Formulas

The power flowing into or radiated from a material body can be expressed as an integral of the normal (inward or outward-directed, respectively) component of the (total) Poynting vector over its surface [29]. Similar expressions may also be derived for the scattered and extinguished power by integrating the associated components of the Poynting vector. In this section, we present simple formulas for computing the powers directly from the equivalent polarization and magnetization currents, which are the immediate output of the numerical solution of JM-VIE formulations. The derivation is based on the conservation of energy, or Poynting’s theorem, which relates the energy flowing out through the boundary surfaces of the body to the work done by the fields on the currents [29].

In particular, we find it useful to decompose the time average Poynting vector [29]

 ⟨stot⟩≜12Re(etot×¯¯¯htot) (27)

in terms of incident, scattered, and extinguished components:

 ⟨stot⟩=⟨sinc⟩+⟨ssca⟩+⟨sext⟩ (28)

where

 ⟨sinc⟩ ≜12Re(einc×¯¯¯hinc) (29a) ⟨ssca⟩ ≜12Re(esca×¯¯¯hsca) (29b) ⟨sext⟩ ≜12Re(esca×¯¯¯hinc+einc×¯¯¯hsca) (29c)

are given by the corresponding fields. Poynting’s theorem can therefore be written in the following form3 [29]:

 ∇⋅⟨sϕ⟩=−⟨Wϕ⟩ (30)

where

 ⟨Wϕ⟩=⟨Wfϕ⟩+⟨Wbϕ⟩ (31)

with , is the time average of the work done on total currents by the corresponding fields. More specifically, it is the sum of work done on free currents () and bound (polarization and magnetization) currents ():

 ⟨Wf,bϕ⟩≜12Re(¯jf,b⋅eϕ+¯¯¯¯¯mf,b⋅hϕ). (32)

In the next few subsections we derive the formulas for the computation of absorbed, extinction, scattered, and radiated power, along with a simple proof for their positivity.

### Iii-a Absorbed Power

The absorbed power is the power flowing into the body and is given by the integral of the inward-directed normal component of the Poynting vector. With the help of (30) and the divergence theorem, we derive the absorbed power in terms of volumetric quantities:

 (33)

where is the outward-directed surface normal. Now, inserting the current expansions (13) in the inner products of (33), we get

 (34)

where

 ^M=(CMχ)−1. (35)

In (34), we have used the definition of the equivalent polarization currents (8) in order to replace the total fields. The substitution is admissible for nonzero susceptibilities, otherwise the currents are identically zero and there is no contribution to the total work. Finally, the absorbed power takes the form

 Pabs=12Rex∗^M^Gx. (36)

### Iii-B Extinction Power

The extinction power is the total power removed from the incident field (the sum of the absorbed and the scattered powers) due to the presence of the scattering object , and is given by similar considerations as in the absorbed power computation, as follows:

 (37)

The computation of the work done by the incident fields on the polarization and magnetization currents can be simply expressed in terms of quantities from the linear system of JM-VIE solution. More specifically, the incident fields are related to the right hand side vector as shown in (19), and the associated inner products admit the following representation:

 ⟨j,beceχe⟩\tiny V+⟨m,bmcmχm⟩\tiny V=1ce∑αx∗e,α⟨pα,be⟩\tiny Vχe,αβ+1cm∑αx∗m,α⟨pα,bm⟩\tiny Vχm,αβ=1cex∗eM−1χebe+1cmx∗mM−1χmbm=x∗^Mb. (38)

The final formula for the extinction power reads

 Pext=12Rex∗^Mb. (39)

It is also useful to write as the real part of an analytic/causal function (see Appendix), both from a theoretical perspective (to get an analogue of the optical theorem as in [30]) and from a practical perspective (e.g. for transforming frequency averaging into a complex frequency, as in [31, 32]).

### Iii-C Scattered Power

The power scattered from an object is given by the real part of the integral of the outward-directed normal component of over 4:

 (40)

Obviously, the scattered power can be expressed in terms of quantities arising in JM-VIE linear system, as the difference between the extinction power and the absorbed power:

 Psca =Pext−Pabs (41)

In cases where the scattering mechanism is weak compared to absorption, formula (41) may be prone to numerical instabilities, i.e., computing a small number as the difference of two almost equal and possibly large approximate values. Therefore, it would be useful to derive some additional formulas for this case, that are immune to numerical instabilities. In doing so, we resort again to the conservation laws (30):

 (42)

where

 Wsca=⟨j,esca⟩\tiny V+⟨m,hsca⟩\tiny V. (43)

The first term of the right hand side is given by

 ⟨j,esca⟩\tiny V% =⟨j,1ce(Nj−j)⟩% \tiny V−⟨j,Km⟩\tiny V=1ce∑αβx∗e,α(⟨pα,Npβ⟩\tiny V−⟨pα,pβ⟩\tiny V)xe,β−∑αβx∗e,α⟨pα,Kpβ⟩\tiny Vxm,β=1cex∗e(N−G)xe−x∗eKxm (44)

and with similar considerations, the second term reads

 ⟨m,hsca⟩\tiny V% =1cmx∗m(N−G)xm+x∗mKxe. (45)

Hence, combining (43), (44) and (45) with (26), we get

 Wsca=x∗Ascax. (46)

Finally, we obtain an alternative formula for the computation of the scattered power:

 Psca=−12Rex∗Ascax. (47)

As mentioned above, there is a trade-off in numerical complexity for getting more stable formula: the matrix in (47) is dense, so the cost of the numerical evaluation scales like instead of the scaling of (41). Note, though, that JM-VIE formulations typically result in very large linear systems and fast solvers are employed for their numerical solution. In this case, the complexity of evaluating (47) scales like ([33, 34, 8, 22] among others), and the associated operators in (26) have been pre-computed in the actual numerical solution of the JM-VIE linear system.

Here we consider radiation from sources immersed in , in particular from elementary sources, i.e., electric and magnetic point Hertzian dipoles. The power radiated from is given by

 (48)

Obviously, the total fields generated in this case are singular at the location of the sources, hence we resort to the natural generalization of the divergence theorem, where the derivatives are taken in the weak/distribution sense,

 (49)

where

 Psup=−∫Ω⟨Wftot⟩=−12Re(⟨jf,etot⟩\tiny V+⟨mf,htot⟩% \tiny V) (50)

is the power supplied by the source and the absorbed power in , already defined. Note, that in the case of dissipative media, both supplied and absorbed powers are infinite [35] (a problem related also to the ill-defined local density of states [36, 37]). Nevertheless, the radiated power (i.e., power flowing from the surface of ) is still a finite quantity and represents the outward power flow from a dipole source with constant amplitude. Otherwise, the notion of the “insulated” dipole could be used (as in [35]), especially when the actual supplied power or the efficiency of the radiator is under scrutiny.

The supplied power can be easily derived with the help of the DD source,

 (51)

As with (41), the above formula is prone to catastrophic cancellations, especially considering that the values there could be very large. An alternative (and more intuitive) formula for the radiated power can be obtained by expanding the fields also for the case of the absorbed power,

where is defined in (23). Hence, by combining (51) and (53) the radiated power admits an elegant quadratic form:

Interestingly, the derived quadratic power formula has the same computation complexity as (52).

### Iii-E Positivity

In addition to their efficiency and simplicity, the above formulas manifest the positivity of power in passive media in a numerically stable fashion: positivity is preserved by discretization. As observed in [24], this is not always the case for SIE formulas. Similar behavior is also expected from the difference formulas (41) and (52) when the parts have almost equal values. While those formulas are analytically exact, their potential reliance on a large cancellation to leave a positive remainder makes them susceptible to numerical inaccuracy and a loss of positivity when they are approximated via a discretized basis. In contrast, we show here that our scattered (47) and radiated (52) power formulas are expressed in terms of manifestly positive-definite quadratic forms, and hence this positivity is preserved by any Galerkin discretization.

We begin with the notion of the Hermitian decomposition: Every matrix can be decomposed into the Hermitian () and the skew-Hermitian () components, as follows:

 B=BH+BSH (55)

where

 BH =(BH)∗=B+B∗2 (56) BSH =−(BSH)∗=B−B∗2.

Obviously, the quadratic forms of these components are purely real and imaginary, respectively. Since, in all formulas we are computing the real part of the associated quadratic forms, the positivity is guaranteed if the Hermitian component is positive- or negative-semidefinite, depending on the sign of the final formula. Hence, the sufficient conditions read

 Rex∗Bx=x∗BHx≥0iffBH⪰0. (57)

In the case of the absorbed power formula (36),

 BH=^G^MH (58)

where

 ^MH=^M+^M∗2=C−1M−1χ−(M−1χ)∗2⪰0 (59)

for any passive material ().

The positivity of the scattered (47) and radiated (54) power is guaranteed if

 BH=AHsca⪯0. (60)

The negative definiteness of carries over to , as explained in [26].

Finally, in the case of the extinction power, positivity is straightforward by definition: the extinction power is the sum of the absorbed and scattered power, which we showed above that are positive.

## Iv Computational Validation

In what follows, we validate the new formulas by using them to compute some representative test cases. The JM-VIE formulation (referred herein as VIE, for simplicity) is numerically solved by means of an in-house FFT-based fast solver [20]. More specifically, the unknown equivalent polarization and magnetization currents are approximated by a series of piecewise constant basis functions for each Cartesian component, with the support of each member of the discrete set being a voxel. A uniform grid of voxels is used for the discretization of the box that encloses the objects under study5. The arising 6-D singular Galerkin inner products in (16) are first reduced to 4-D singular (of lower order) integrals over the surfaces of the voxels [20], and then computed by means of DIRECTFN open-source package [39, 40]. The benchmark results are obtained with analytical formulas (Mie theory), and with a surface integral equation (SIE) method, and more specifically with the open-source package scuff-em [28, 41]. The default choice for the results regarding the VIE method are the boxed formulas presented above.

### Iv-a Plane Wave Excitation

We begin by computing efficiencies for scattering and absorption of a spherical particle of radius irradiated by a linearly-polarized -traveling plane wave with electric field:

 einc=e0e−ik0z^x.

Fig. 1 plots efficiencies as functions of the dimensionless “size parameter” . Efficiencies are obtained from cross sections , where , by dividing by the geometrical cross section () of the sphere, . The results are in good agreement with the efficiencies obtained by integrating the associated Poynting vectors by means of a SIE method [25]. For the computation of the scattering efficiency, we present additional results based on the difference formula (41), where one can identify the expected instabilities (blue +) discussed in Section IIIC. Note that the two missing data points assume negative values; there is no guarantee for the positivity of the difference formula, as is the case in the PM.

### Iv-B Dipole Excitation

We now consider the same sphere but irradiated by a Hertzian dipole (with moment equal to 1) directed along -axis, and located at the center of the sphere. Fig. 2 plots the radiated power flowing through the surface of the spherical particle, normalized with respect to the power radiated by the dipole in free-space (). We validate our VIE results by comparing against a reference analytical solution obtained via a Mie series[29],

 PMie=Z0c2a122k20√6π|eian3|2|D|2 (61)

where

 D=[n2(a3+a−i)−a+i]sin(na)+na[−in2(−1+a(a−i))+a−i]cos(na)

with being the speed of light, , and . Note that the power radiated by the dipole in free space .

We consider three different scenarios in Fig. 2, corresponding to different values of (with ). Specifically, we consider spheres with , , and , corresponding to free-space, lossless, and dissipative media, respectively. As evidenced by the results, the quadratic formula (54) is stable both for low and high frequencies, where convergence to the exact solution is attained as the mesh discretization gets finer. Note that due to the uniform mesh used in this work, low resolution meshes suffer from staircase approximation errors.

Next, we consider the same dipole radiation but using the difference formula of (52), and show that it leads to the aforementioned instabilities (Section IIID). Fig. 3 plots the radiated power as obtained from (52) (solid lines), for the case of the lossy dielectric sphere of Fig. 2. Also shown are the corresponding results from the Mie (black line) and VIE (54) (open circles) solutions. As discussed in Section IIID and [35], both the supplied and the absorbed power are infinite in this case. Nevertheless, the radiated power flowing through the surface of the sphere is the finite quantity presented in Fig. 2. As depicted in the inset of Fig. 3, the supplied and absorbed (not shown) powers diverge with the resolution of the mesh as . Consequently, the difference suffers from catastrophic cancellations that render the difference formula (52) practically useless. This result highlights the importance of the quadratic formula (54), which is remarkably stable and identically positive.

Finally, we consider the case of an inhomogeneous dielectric cube irradiated by a -directed dipole placed at the center. The continuous profile of the permittivity is given as follows:

 ϵr(z)=ϵr,l+z+R2R(ϵr,h−ϵr,l),z∈[−R,R] (62)

where and . The radiated power presented in Fig. 4 is computed using the quadratic formula (54) and converges as we refine the discretization. This is a particularly interesting example, since the continuous profile of the inhomogeneity rules out methods based on SIE formulations.

## V Conclusion

A collection of simple and stable formulas is presented for the computation of absorbed, scattered, extinction, and radiated power in VIE formulations. The proposed formulas (boxed equations in the manuscript) are accurate in a wide range of frequencies, and are based solely on volumetric quantities found in the associated linear system of equations. In addition, they preserve the positivity of the computed power, thus accurately capturing the physics of the problem. Thus, there is no need for significant post-processing, such as the evaluation of the fields and the integration of the Poynting vector along enclosing surfaces. By construction, the presented scheme is immune to the well-known instability issues that occur in Poynting’s method. The efficient and compact absorption/radiation formulas presented herein are expected to be especially useful—besides applications in classical scattering/radiation problems—in computations of EM fluctuation phenomena, including radiative heat transfer and Casimir forces between complex bodies. Our analysis is based on a current-based VIE formulation, but similar formulas may be easily derived, with only minor modifications, for the case of VIE formulations based on fields or fluxes.

## Vi Acknowledgments

This work was supported in part by grants from the Singapore-MIT programs in Computational Engineering and in Computational and Systems Biology, from the Skolkovo-MIT initiative in Computational Mathematics, and from the Army Research Office through the Institute for Soldier Nanotechnologies under Contract No. W911NF-07-D0004.

[On the analyticity of ]

## Appendix A On the analyticity of Pext

As mentioned in Section IIIB, it is important to write the formula for the extinction power (39) in a form where the analyticity (in the lower half of the complex- plane) is shown explicitly, so as to be able to exploit it both from a theoretical and a practical perspective. Indeed, (39) can be written with the help of (19) and (35) as follows:

 Pext=12Re(einchinc)∗x=12Re(einchinc)∗W(einchinc)

where is the matrix arising from the discretization of the operator relating the incident fields to the induced currents, i.e.,

 (jm)=W(einchinc).

Causality implies that is an analytic function in the lower half of the complex- plane [42]. Alternatively, it is sufficient to assume passivity rather than causality, since the former implies the latter in a time-invariant linear system [43].

Finally, we can eliminate the complex conjugation by exploiting the conjugate symmetry of the Fourier transform of any real incident field:

As desired, this is the real part of an analytic function in the lower-half complex- plane, as long as the incident fields are entire (everywhere-analytic) functions of (which is true for all typical incident fields, such as planewaves, gaussian pulses, or any pulse that is compactly supported in the time domain [44]).

### Footnotes

1. As explained in Section IIIC, an additional matrix vector product may be needed.
2. More on the use of operator can be found in [15, 19, 22].
3. Note that for a lossless ambient medium.
4. Note here the plus sign!
5. Of course, one could choose different schemes for the numerical solution of the VIE method, e.g. based on a tetrahedral mesh coupled with a FMM solver [38].

### References

1. D. H. Schaubert, D. R. Wilton, and A. W. Glisson, “A tetrahedral modeling method for electromagnetic scattering by arbitrarily shaped inhomogeneous dielectric bodies,” IEEE Trans. Antennas Propag., vol. 32, no. 1, pp. 77–85, Jan. 1984.
2. D. T. Borup and O. P. Gandhi, “Fast-Fourier transform method for calculation of SAR distributions in finely discretized inhomogeneous models of biological bodies,” IEEE Trans. Microw. Theory Tech., vol. 32, no. 3, pp. 355–360, Apr. 1984.
3. C. Y. Shen, K. J. Glover, M. I. Sancer, and A. D. Varvatsis, “The discrete Fourier transform method of solving differential-integral equations in scattering theory,” IEEE Trans. Antennas Propag., vol. 37, no. 8, pp. 1032–1041, Aug. 1989.
4. M. F. Cátedra, E. Gago, and L. Nuño, “A numerical scheme to obtain the RCS of three dimensional bodies of resonant size using the conjugate gradient method and the fast Fourier transform,” IEEE Trans. Antennas Propag., vol. 37, no. 5, pp. 528–537, May 1989.
5. P. Zwamborn and P. M. van den Berg, “The three-dimensional weak form of the conjugate gradient FFT method for solving scattering problems,” IEEE Trans. Microw. Theory Tech., vol. 40, no. 9, pp. 1757–1766, Sep. 1992.
6. H. Gan and W. C. Chew, “A discrete BCG-FFT algorithm for solving 3D inhomogeneous scatterer problems,” J. Electromag. Waves Applicat., vol. 9, no. 10, pp. 1339–1357, 1995.
7. S. A. de Carvalho and L. de Souza Mendes, “Scattering of EM waves by inhomogeneous dielectrics with the use of the method of moments and 3-D solenoidal basis functions,” Microwave Opt. Tech. Lett., vol. 23, pp. 42–46, Oct. 1999.
8. C.-. Lu, “A fast algorithm based on volume integral equation for analysis of arbitrarily shaped dielectric radomes,” IEEE Trans. Antennas Propag., vol. 51, no. 3, pp. 606–612, Mar. 2003.
9. M.-K. Li and W. C. Chew, “Applying divergence-free condition in solving the volume integral equation,” Progress in Electromagnetic Research, vol. 57, pp. 311–333, 2006.
10. G. Rubinacci and A. Tamburrino, “A broadband volume integral formulation based on edge-elements for full-wave analysis of lossy interconnects,” IEEE Trans. Antennas Propag., vol. 54, no. 10, pp. 2977–2989, Oct. 2006.
11. M. M. Botha, “Solving the volume integral equations of electromagnetic scattering,” Journal of Computational Physics, vol. 218, pp. 141 – 158, 2006.
12. M. I. Sancer, K. Sertel, J. L. Volakis, and P. V. Alstine, “On volume integral equations,” IEEE Trans. Antennas Propag., vol. 54, no. 5, pp. 1488–1495, May 2006.
13. N. A. Ozdemir and J.-F. Lee, “A nonconformal volume integral equation for electromagnetic scattering from anisotropic materials,” Proc. Antennas Propag, Soc. Int. Symp., pp. 2889–2892, 2006.
14. ——, “A nonconformal volume integral equation for electromagnetic scattering from penetrable objects,” IEEE Trans. Microw. Theory Tech., vol. 43, no. 4, pp. 1369–1372, Apr. 2007.
15. L. E. Sun and W. C. Chew, “A novel formulation of the volume integral equation for electromagnetic scattering,” Waves in Random and Complex Media, vol. 19, no. 1, pp. 162–180, 2009.
16. J. Markkanen, C.-. Lu, X. Cao, and P. Ylä-Oijala, “Analysis of volume integral equation formulations for scattering by high-contrast penetrable objects,” IEEE Trans. Antennas Propag., vol. 60, no. 5, pp. 2367–2374, May 2012.
17. M. C. van Beurden and S. J. L. van Eijndhoven, “Gaps in present discretization schemes for domain integral equations,” 2007 International Conference on Electromagnetics in Advanced Applications, ICEAA’07, pp. 673–675, 2007.
18. ——, “Well-posedness of domain integral equations for a dielectric object in homogeneous background,” Journal of Engineering Mathematics, vol. 62, no. 3, pp. 289–302, 2008.
19. J. Markkanen, P. Ylä-Oijala, and A. Sihvola, “Discretization of volume integral equation formulations for extremely anisotropic materials,” IEEE Trans. Antennas Propag., vol. 60, no. 11, pp. 5195–5202, Nov. 2012.
20. A. G. Polimeridis, J. F. Villena, L. Daniel, and J. K. White, “Robust J-EFVIE solvers based on purely surface integrals,” 2013 International Conference on Electromagnetics in Advanced Applications, ICEAA’13, pp. 379–381, 2013.
21. A. G. Polimeridis and J. K. White, “FFT-JVIE algorithm for computation of electromagnetic fields in inhomogeneous dielectric objects,” Proc. Antennas Propag, Soc. Int. Symp., pp. 254–255, 2013.
22. A. G. Polimeridis, J. F. Villena, L. Daniel, and J. K. White, “Stable FFT-JVIE solvers for fast analysis of highly inhomogeneous dielectric objects,” Journal of Computational Physics, vol. 227, no. 14, pp. 7052 – 7068, 2014.
23. A. F. Peterson, “Efficient solenoidal discretization of the volume EFIE for electromagnetic scattering from dielectric objects,” IEEE Trans. Antennas Propag., vol. 62, no. 3, pp. 1475–1478, 2014.
24. A. M. Kern and O. J. F. Martin, “Pitfalls in the determination of optical cros sections from surface integral equation simulations,” IEEE Trans. Antennas Propag., vol. 58, no. 6, pp. 2158–2161, Jun. 2010.
25. M. T. H. Reid and S. G. Johnson, “Efficient computation of power, force, and torque in BEM scattering calculations,” arXiv e-prints, 2014.
26. A. W. Rodriguez, M. T. H. Reid, and S. G. Johnson, “Fluctuating surface-current formulation of radiative heat transfer for arbitrary geometries: Theory and applications,” Physical Review B, vol. 88, no. 5, p. 054305, 2013.
27. M. T. H. Reid, A. W. Rodriguez, and S. G. Johnson, “Fluctuation-induced phenomena in nanoscale systems: harnessing the power of noise,” Proc. of the IEEE, vol. 101, no. 2, pp. 531 – 545, 2013.
28. M. T. H. Reid, A. W. Rodriguez, J. White, and S. G. Johnson, “Efficient computation of Casimir interactions between arbitrary 3d objects,” Phys. Rev. Lett., vol. 103, p. 040401, Jul 2009. [Online]. Available: http://link.aps.org/doi/10.1103/PhysRevLett.103.040401
29. J. D. Jackson, Classical Electtrodynamics.   3rd ed. John Wiley & Sons, 1999.
30. H. Hashemi, C. W. Qiu, A. P. McCauley, J. D. Joannopoulos, and S. G. Johnson, “Diameter-bandwidth product limitation of isolated-object cloaking,” Phys. Rev. A, vol. 86, p. 013804, 2012.
31. X. Liang and S. G. Johnson, “Formulation for scalable optimization of microcavities via the frequency-averaged local density of states,” Optics Express, vol. 21, no. 25, pp. 30 812–30 841, 2013.
32. O. D. Miller, C. W. Hsu, M. T. H. Reid, W. Qiu, B. G. Delacy, J. D. Joannopoulos, M. Soljacic, and S. G. Johnson, “Fundamental limits to extinction by metallic nanoparticles,” Physical Review Letters, vol. 112, no. 12, p. 123903, 2014.
33. V. Rokhlin, “Rapid solution of integral equations of scattering theory in two dimensions,” Journal of Computational Physics, vol. 86, pp. 414 – 439, 1990.
34. J. R. Philips and J. K. White, “A precorrected-FFT method for electrostatic analysis of complicated 3-D structures,” IEEE Trans. Comput.-Aided Design Integr. Circuits Syst., vol. 16, no. 10, pp. 1059–1072, Oct. 1997.
35. C. T. Tai and R. E. Collin, “Radiation of a Hertzian dipole immersed in a dissipative medium,” IEEE Trans. Antennas Propag., vol. 48, no. 10, pp. 1501–1506, Oct. 2000.
36. S. Scheel, L. Knoll, and D. G. Welsch, “Spontaneous decay of an excited atom in an absorbing dielectric,” Phys. Rev. A, vol. 60, no. 5, pp. 4094–4104, 2008.
37. C. V. Vlack and S. Hughes, “Finite-difference time-domain technique as an efficient tool for calculating the regularized green function: applications to the local-field problem in quantum optics for inhomogeneous lossy materials,” Optics Letters, vol. 37, no. 14, pp. 2880–2882, 2012.
38. S. Järvenpää, J. Markkanen, and P. Ylä-Oijala, “Broadband multilevel fast multipole algorithm for electric-magnetic current volume integral equation,” IEEE Trans. Antennas Propag., vol. 61, no. 8, pp. 4393–4397, Aug. 2013.
39. A. G. Polimeridis, F. Vipiana, J. R. Mosig, and D. R. Wilton, “DIRECTFN: Fully numerical algorithms for high precision computation of singular integrals in galerkin SIE methods,” IEEE Trans. Antennas Propag., vol. 61, no. 6, pp. 3112–3122, Jun. 2013.
40. “DIRECTFN package,” 2012. [Online]. Available: {http://web.mit.edu/thanos_p/www/Software}
41. http://homerreid.com/scuff-EM.
42. L. Landau, E. M. Lifshitz, and L. P. Pitaevskii, Electrodynamics of Continuous Media.   2nd ed., Butterworth-Heinemann, Oxford, 1984.
43. A. Welters, Y. Avniel, and S. G. Johnson, “Speed-of-light limitations in passive linear media,” arXiv.org e-prints, arXiv:1405.0238, 2014.
44. R. S. Strichartz, A Guide to Distribution Theory and Fourier Transforms.   World Scientific Pub Co Inc, 2003.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters