Numerical method for solving electromagnetic wave scattering by one and many small perfectly conducting bodies

In this paper, we investigate the problem of electromagnetic (EM) wave scattering by one and many small perfectly conducting bodies and present a numerical method for solving it. For the case of one body, the problem is solved for a body of arbitrary shape, using the corresponding boundary integral equation. For the case of many bodies, the problem is solved asymptotically under the physical assumptions $a\ll d \ll \lambda$, where $a$ is the characteristic size of the bodies, $d$ is the minimal distance between neighboring bodies, $\lambda=2\pi/k$ is the wave length and $k$ is the wave number. Numerical results for the cases of one and many small bodies are presented. Error analysis for the numerical method are also provided.


Introduction
Many real-world electromagnetic (EM) problems like EM wave scattering, EM radiation, etc [17], cannot be solved analytically and exactly to get a solution in a closed form. Thus, numerical methods have been developed to tackle these problems approximately. Computational Electromagnetics (CEM) has evolved enormously in the past decades to a point that its methods can solve EM problems with extreme accuracy. These methods can be classified into two categories: Integral Equation (IE) method and Differential Equation (DE) method. Typical IE methods include: Method of Moment (MoM) developed by Roger F. Harrington (1968) [12], Fast Multipole Method (FMM) first introduced by Greengard and Rokhlin (1987) [11] and then applied to EM by Engheta et al (1992) [8], Partial Element Equivalent Circuit (PEEC) method [4], and Discrete Dipole Approximation [2]. Typical DE methods are: Finite Difference Time Domain (FDTD) developed by Kane Yee (1966) [3], Finite Element Method (FEM) [5], Finite Integration Technique (FIT) proposed by Thomas Weiland (1977) [16], Pseudospectral Time Domain (PSTD) [15], Pseudospectral Spatial Domain (PSSD) [18], and Transmission Line Matrix (TLM) [9]. Among these methods, FDTD has emerged as one of the most popular techniques for solving EM problems due to its simplicity and ability to provide animated display of the EM field. However, FDTD requires the entire computational domain be gridded [19], that results in very long solution times. Furthermore, as a DE method, it does not take into account the radiation condition in exact sense [13,14], which leads to certain error in the solution. On the other hand, spurious solutions might exist in DE methods [20][21][22]. Most importantly, most of DE methods are not suitable if the number of bodies is very large.
In [23-25, 27, 28, 30], A. G. Ramm has developed a theory of EM wave scattering by many small perfectly conducting and impedance bodies. In this theory, the EM wave scattering problem is solved asymptotically under the physical assumptions: a d λ, where a is the characteristic size of the bodies, d is the minimal distance between neighboring bodies, λ = 2π/k is the wave length and k is the wave number. In [29], a numerical method is developed for solving EM wave scattering by many small impedance bodies. In this paper, the problem of EM wave scattering by one and many small perfectly conducting bodies is considered. A numerical method for solving this problem asymptotically based on the above theory is presented. For the case of one body, the problem is solved for a body of arbitrary shape, using the corresponding boundary integral equation. For the case of many small bodies, the problem is solved under the basic assumptions a d λ and the assumption about the distribution of the small bodies N (∆) = 1 In Sections 2 and 3, the theory of EM wave scattering by one and many small perfectly conducting bodies is presented. The numerical methods for solving these problems are also described in details. Furthermore, error analysis for the numerical methods of solving the EM scattering problem are also provided. In Section 4 these methods are tested and numerical results are discussed.

EM wave scattering by one perfectly conducting body
Let D be a bounded perfectly conducting body, a = 1 2 diamD, S be its C 2 -smooth boundary, and D := R 3 \ D. Let and µ be the dielectric permittivity and magnetic permeability constants of the medium in D . Let E and H denote the electric and magnetic fields, respectively, E 0 be the incident field and v E be the scattered field. The problem of electromagnetic wave scattering by one perfectly conducting body can be stated as follows 3) where ω > 0 is the frequency, k = 2π/λ = ω µ is the wave number, ka 1, λ is the wave length, E is a constant vector, and α is a unit vector that indicates the direction of the incident wave E 0 . This incident wave satisfies the relation ∇ · E 0 = 0. The scattered field v E satisfies the radiation condition (2.6). Here, N is the unit normal vector to the surface S, pointing out of D. By [·, ·] the vector product is denoted and α · x is the scalar product of two vectors.
The solution to problem (2.1)-(2.6) can be found in the form see [23]. Here, E is a vector in R 3 and ∇ × E is a pseudo-vector, that is a vectorlike object which changes sign under reflection of its coordinate axes. E 0 is the incident plane wave defined in (2.5) and J is an unknown pseudo-vector that is to be found. J is assumed to be tangential to S and continuous. J can be found by applying the boundary condition (2.3), or equivalently [N , E ] = 0, to (2.7) and solving the resulting boundary integral equation If D is sufficiently small, then equation (2.9) is uniquely solvable in C (S) and its solution J is tangential to S, see [23]. The asymptotic formula for E when the radius a of the body D tends to zero can be derived as follows, see [23]. Rewrite equation (2.7) as where x 1 ∈ D, an arbitrary point inside the small body D, and the second term in (2.11) is much greater than the last term Then, the asymptotic formula for E when a tends to zero is Thus, when D is sufficiently small, instead of finding J , we can just find one pseudovector Q. The analytical formula for Q is derived as follows, see [30]. By integrating both sides of (2.8) over S, one gets This is equivalent to When a → 0, this equation becomes (2.20) where in the second term, summations over the repeated indices are understood, e p , 1 ≤ p ≤ 3, are the orthogonal unit vectors in R 3 , |D| is the volume of D, |D| = c D a 3 , and in the third term we use this estimate then equation (2.20) can be rewritten as follows where Γ is a 3 × 3 constant matrix and it is defined by (2.25) in which summations are understood over the repeated indices. Thus, Q can be written as where I := I 3 , the 3 × 3 identity matrix. This formula is asymptotically exact as a → 0.

Numerical method for solving EM wave scattering by one perfectly conducting spherical body
In this section, we consider the EM wave scattering problem by a small perfectly conducting spherical body. Instead of solving the problem (2.1)-(2.6) directly, we will solve its corresponding boundary integral equation Then the solution E to the EM wave scattering problem by one perfectly conducting body can be computed by either the exact formula (2.7) or the asymptotic formula (2.17). Scattering by a sphere has been discussed in many papers, for example [1] in which Mie solves the EM wave scattering problem by separation of variables. The EM field, scattered by a small body, is proportional to O(a 3 ).
Suppose S is a smooth surface of a spherical body. Let S be partitioned into P non-intersecting subdomains S i j , 1 ≤ i ≤ m θ , 1 ≤ j ≤ m φ , using spherical coordinates, where m θ is the number of intervals of θ between 0 and 2π and m φ defines the number of intervals of φ between 0 and π. Then P = m θ m φ + 2, which includes the two poles of the sphere. m θ is defined in this way: m θ = m φ + |φ − π 2 |6m φ . This means the closer it is to the poles of the sphere, the more intervals for θ are used. Then the point (θ i , φ j ) in S i j is chosen as follows Note that there are many different ways to distribute collocation points. However, the one that we describe here will guarantee convergence to the solution to (2.27) with fewer collocation points used from our experiment. Furthermore, one should be careful when choosing the distribution of collocation points on a sphere. If one chooses φ j = j π m φ , 1 ≤ j ≤ m φ , then when j = m φ , φ j = π and thus there is only one point for this φ regardless of the value of θ as shown in (2.30). The position of a point in each S i j can be computed by (2.30) and the outward-pointing unit normal vector N to S at this point is For a star-shaped body with a different shape, only the normal vector N needs to be recomputed. Rewrite the integral equation (2.27) as This integral equation can be discretized as follows , and ∆ j is the surface area of the subdomain j . This is a linear system with unknowns This linear system can be rewritten as follows where by the subscripts x, y, z the corresponding coordinates are denoted, e.g.

Error analysis
Recall the boundary integral equation (2.8) Integrate both sides of this equation over S and get see Section 2. Once J is found from solving (2.46), Q can be computed by Q = S J (t )d t . Then one can validate the values of J and Q by checking the following things • Is J tangential to S as shown in Section 2? One needs to check J (s) · N s .
The relative error of Q can be computed as follows where R H S := −|D|∇×E 0 . This will give the error of the numerical method for the case of one body.
Furthermore, one can also compare the value of the asymptotic Q a in formula (2.26) with the exact Q e defined in (2.12) by

General method for solving EM wave scattering by one perfectly conducting body
In this section, we present a general method for solving the EM wave scattering problem by one perfectly conducting body, whose surface is parametrized by

EM wave scattering by many small perfectly conducting bodies
Consider a bounded domain Ω containing M small bodies D m , 1 ≤ m ≤ M , and S m are their corresponding smooth boundaries. Let D := M m=1 D m ⊂ Ω and D be the complement of D in R 3 . We assume that S = M m=1 S m is C 2 -smooth. is the dielectric permittivity constant and µ is the magnetic permeability constant of the medium. Let E and H denote the electric and magnetic fields, respectively. E 0 is the incident field and v is the scattered field. The problem of electromagnetic wave scattering by many small perfectly conducting bodies involves solving the following system 3) where v satisfies the radiation condition (2.6), ω > 0 is the frequency, k = 2π/λ is the wave number, ka 1, a := 1 2 max m diamD m , and α is a unit vector that indicates the direction of the incident wave E 0 . Furthermore, Assume that the distribution of small bodies D m , 1 ≤ m ≤ M , in Ω satisfies the following formula where N (∆) is the number of small bodies in ∆, ∆ is an arbitrary open subset of Ω, and N (x) is the distribution function Note that E solves this equation if µ =const. Once we have E , then H can be found from this relation From (3.10) and (3.9), one can get (3.2). Thus, we need to find only E which satisfies the boundary condition (3.3). It was proved in [23] that under the radiation condition and the assumptions a d λ, the problem (3.1)-(3.5) has a unique solution and its solution is of the form where J m are unknown continuous functions that can be found from the boundary condition. Let When a → 0, the asymptotic solution for the electric field is given by Therefore, instead of finding J m (t ), ∀t ∈ S, 1 ≤ m ≤ M , to get the solution E , one can just find Q m . This allows one to solve the EM scattering problem with a very large number of small bodies which is impossible to do before. The analytic formula for Q m can be derived by using formula (2.26) and replacing E 0 in this formula by the effective field E e (x m ) acting on the m-th body where the effective field acting on the m-th body is defined as When a → 0, the effective field E e (x) is asymptotically equal to the field E (x) in (3.13) as proved in [23] and [30]. Let E em := E e (x m ), where x m is a point in D m . From (3.14), and (3.15), one gets (3.16)

Numerical method for solving EM wave scattering by many small perfectly conducting bodies
For finding the solution to EM wave scattering in the case of many small perfectly conducting bodies, we need to find E em in (3.16). Apply the operator (I +Γ) −1 ∇× to both sides of (3.16) and let A m := (I + Γ) −1 ∇ × E em . Then Solving this system yields A m , for 1 ≤ m ≤ M . Then E can be computed by where Equation (3.17) can be rewritten as follows  (3.23) in which by the subscripts x, y, z the corresponding coordinates are denoted, e.g.
for i = j ; when i = j : a i i = 0, b i i = 1, and c i i = 0;

Error analysis
The error of the solution to the EM wave scattering problem by many small perfectly conducting bodies can be estimated as follows. From the solution E of the electromagnetic scattering problem by many small bodies given in (3.11) we can rewrite it as (3.34) Comparing this with the asymptotic formula for E when a → 0 given in (3.13) we have the error of this asymptotic formula is

EM wave scattering by one perfectly conducting spherical body
To illustrate the idea of the numerical method, we use the following physical parameters to solve the EM wave scattering problem by one small perfectly conducting sphere, i.e solving the linear system (2.34)-(2.36) • Speed of wave, c = (3.0E + 10) cm/sec.
• The body is a sphere of radius a, centered at the origin.
Recall that We choose a coordinate system centered at the center of the sphere such that t = (0, 0, a) and s = aN . Then Therefore, matrix Γ is For example, Table 1 shows the exact and asymptotic vector Q when the radius of the body is a = (1.0E − 09) cm and the number of collocation points used to solve the integral equation (2.27) is P = 766. Note that a = (1.0E − 09) cm satisfies ka 1. The point x 1 in (2.17) is taken at the center of the body, the origin. Table 2 and 3 show the exact and asymptotic vector E = (E x , E y , E z ), the electric field, at the point x outside of the body, respectively. The distance |x − x 1 | is measured in cm in these tables. 1.0E-21 * Q e 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.3925i Q a 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.3760i Table 2: Vector E e for one perfectly conducting body with a = (1.0E −09) cm and P = 766 collocation points.
1.73E-08 1.0000 + 0.0010i 0.0000 + 0.0000i 0.0000 + 0.0000i 1.73E-07 0.9999 + 0.0105i 0.0000 + 0.0000i 0.0000 + 0.0000i 1.73E-06 0.9945 + 0.1045i 0.0000 + 0.0000i 0.0000 + 0.0000i  Table 4.  Table 5 compares the asymptotic Q a versus exact Q e and asymptotic E a versus exact E e , when P = 1386 collocation points, |x − x 1 | = 1.73E − 05 cm, and with various a. The errors shown in this table are relative errors, see the error analysis in Section 2.2. As one can see from this table, the smaller the radius a is, compared to the distance from the point of interest to the center of the body, the more precise the asymptotic formulas of E and Q are.
Furthermore, the numerical results also depend on the number of collocation points used. The more collocation points used, the more accurate the results is.

EM wave scattering by one perfectly conducting ellipsoid body
In this section, we consider the EM wave scattering problem by a small perfectly conducting ellipsoid body. The method for solving the problem in this setting is the same as that of Section 2.1 except that one needs to recompute the unit normal vector N .
To get the solution of this problem, one can follow the steps in Section 2.1. In particular, one needs to solve the linear system (2.34)-(2.36).
To illustrate the idea, we use the same physical parameters as described in Section 4.1, except that the body now is an ellipsoid. Let S be its smooth surface. The way we partition S into many subdomains S i j is the same as the way we partition a spherical body as described in section 2.1. Then, the position of the collocation point in each subdomain S i j is defined by (x, y, z) i j = (a cos θ i sin φ j , b sin θ i sin φ j , c cos φ j ), (4.10) where a, b, and c are the lengths of the semi-principal axes of the ellipsoid. The outward-pointing normal vector n to S at this point is 11) and the corresponding unit normal vector N is For example, Tables 6 and 7 show the exact and asymptotic vector E = (E x , E y , E z ), the electric field, got from solving this EM wave scattering problem with one perfectly conducting ellipsoid body, when the semi-principle axes of the body are a = (1.0E − 08) cm, b = (1.0E − 09) cm, c = (1.0E − 09) cm, and Table 6: Vector E e for one perfectly conducting ellipsoid body with a = (1.0E − 08) cm, b = (1.0E − 09) cm, c = (1.0E − 09) cm, and P = 1052 collocation points.
As for the case of one body, we need to verify the following things: a) Is J tangential to S? In fact, this vector J is tangential to the surface S of the body, J · N s = O(10 −13 ). b) Are Q and J correct? We check the relative error described in Section 2.2, Error = |Q+ΓQ−R H S| |R H S| = 14%. The more collocation points used, the smaller this error is, for example, with P = 1762 collocation points, this error is only 3.6%. c) How accurate is the asymptotic formula (2.17) for E ? The accuracy of the asymptotic formula for E in (2.17) can be checked by comparing it with the exact formula (2.7) at several points x outside of the body, max(a, b, c) where x 1 is the center of the body. The relative errors are given in Table 8. Table 9 shows the relative errors between the asymptotic E versus exact E , when P = 1052 collocation points, |x − x 1 | = 1.73E − 07 cm, and with various semi-principle axes a, b, and c. As one can see from this table, the smaller the semi-principle axes are, compared to the distance from the point of interest to Table 8: Relative errors between the asymptotic and exact formulas for E when P = 1052 collocation points, a = (1.0E − 08) cm, b = (1.0E − 09) cm, and c = (1.0E − 09) cm.
|x − x 1 | E e vs E a 1.01E-07 1.73E-04 1.01E-06 1.73E-07 1.01E-05 1.73E-10 Table 9: Relative errors of the asymptotic E when P = 1052 collocation points. the center of the body, the more accurate the asymptotic formulas of E is.

EM wave scattering by one perfectly conducting cubic body
In this section, we consider the EM wave scattering problem by a small perfectly conducting cubic body. Again, the method for solving the problem in this setting is the same as that of Section 2.1 except that one needs to recompute the unit normal vector N .
One can follow the steps outlined in Section 2.1 to solve this problem. That means, one needs to solve the linear system (2.34)-(2.36).
For illustration purpose, we use the same physical parameters as described in Section 4.1, except that the body now is a cube. Suppose the cube is placed in the first octant where the origin is one of its vertices, one can use the standard unit vectors in R 3 as the unit normal vectors to the surfaces of the cube.
For example, Tables 10 and 11 show the exact and asymptotic vector E = (E x , E y , E z ), the electric field, got from solving this EM wave scattering problem with one perfectly conducting cubic body, when the half side of the body is a = (1.0E − 07) cm and the number of collocation points is M = 600. Note that a = 1.0E − 07 cm satisfies ka 1. The point x 1 in (2.17) is taken at the center of the cubic body. Each row in Tables 10 and 11 shows the exact and asymptotic Table 10: Vector E e for one perfectly conducting body with a = (1.0E − 07) cm and M = 600 collocation points.
As before, for the case of one body, we need to verify the following things: a) Is J tangential to S? In fact, this vector J is tangential to the surface S of the body, J · N s = O(10 −13 ). b) How accurate is the asymptotic formula (2.26) for Q? We check the relative error described in Section 2.2, Error = |Q+ΓQ−R H S| |R H S| = 1.13%. c) How accurate is the asymptotic formula (2.17) for E ? The accuracy of the asymptotic formula for E in (2.17) can be checked by comparing it with the exact formula (2.7) at several points x outside of the body, |x − x 1 | a where x 1 is the center of the body. The relative errors are given in Table 12. Table 13 shows the relative errors between the asymptotic E versus exact E , when M = 600 collocation points, |x − x 1 | = (1.73E − 06) cm, and with various a. From this table, we can see that the smaller the side of the cube is, compared to the distance from the point of interest to the center of the body, the more accurate the asymptotic formula of E is.

EM wave scattering by many small perfectly conducting bodies
To illustrate the idea, consider a domain Ω as a unit cube placed in the first octant such that the origin is one of its vertices. This domain Ω contains M small bodies. Suppose these small bodies are particles. We use GMRES iterative method, see [10], to solve the linear system (3.21)-(3.23). The following physical parameters are used to solve the EM wave scattering problem • Speed of wave, c = (3.0E + 10) cm/sec.
• Volume of the domain Ω that contains all the particles, |Ω| = 1 cm 3 .
Note that the distance d satisfies the assumption d λ. The radius a of the particles is chosen variously so that it satisfies the assumption ka 1. For illustration purpose, the problem of EM wave scattering by many small perfectly conducting bodies is solved with M = 27 and 1000 particles.  Table 15 shows the relative errors of E when there are M = 27 particles in the cube, the distance between neighboring particles is d = (1.0E − 07) cm, and with various radius a. Figure 1 shows the relative error of the asymptotic E . From this figure, one can see that when the ratio a/d decreases from 1.0E − 01 to 1.0E − 04, the error of the asymptotic solution decreases linearly and rapidly from 8.16E − 06 to about 8.16E − 18. The smaller the ratio a/d is, the better the asymptotic formula (3.13) approximates E .   Figure 2 show the results of solving the problem with M = 1000 particles, when the distance between neighboring particles is d = (1.0E − 07)  cm, and with different radius a. From these table and figure, one can see that the relative error of the asymptotic solution in this case is also very small, less than 3.02E − 04, when the ratio a/d < 1.0E − 01. In this case, the error of the asymptotic E is greater than that of the previous case when M = 27. However, this time, the error is also decreasing quickly and linearly when the ratio a/d decreases from 1.0E − 01 to 1.0E − 04. Therefore, the asymptotic formula (3.13) for the solution E is applicable when a d .

Conclusions
In this paper, we present a numerical method for solving the EM wave scattering by one and many small perfectly conducting bodies. One of the advantages of this method is that it is relatively easy to implement. Furthermore, one can get an asymptotically exact solution to the problem when the characteristic size of the bodies tends to zero. To illustrate the applicability and efficiency of the method, we use it to solve the EM wave scattering problem by one and many small perfectly conducting bodies. Numerical results of these experiments are presented and error analysis of the asymptotic solutions for the case of one and many bodies are also discussed. For the case of one small body, one can always find the exact solution using the described method. For the case of many small bodies, the accuracy of our method is high if a d λ. The problem of EM wave scattering is much harder to treat, compared to scalar wave scattering [6,7,26,31]. For scalar wave scattering problem, a fast algorithm is developed in [31] to deal with billions of particles. This is still open to EM wave scattering. One might consider this as future research.