The Spectral Integral Method (SIM) for the Scattering from an Arbitrary Number of Circular PEC Cylinders

We present an accurate spectral integral method (SIM) for the analyses of scattering from multiple circular perfect electric conductor (PEC) cylinders. It solves the coupled surface integral equations by using the Fourier series and addition theorem to decouple the system. The SIM has exponential convergence so that the error decreases exponentially with the sample density on the surfaces, and requires only about 2-3 points per wavelength (PPW) to reach engineering accuracy with less than 1% error. Numerical results demonstrate that the SIM is much more accurate and efficient than the method of moments (MoM), and thus can be potentially used as the exact radiation boundary condition in the finite element and spectral element methods.


I. INTRODUCTION
T HE scattering of multiple circular perfect electric con- ductor (PEC) cylinders is of significant scientific and engineering interest.In this work, we develop a fast spectral integral method (SIM) for simulating the scattering from an arbitrary number of non-overlapping circular PEC cylinders as shown in Figure 1.The radii and center positions of all cylinders can be arbitrary.
Such a scattering problem is classical and can be treated by analytical [1]- [5] and numerical methods [6]- [8].The analytical solutions use cylindrical harmonics to expand the fields everywhere in terms of cylindrical harmonics and solve for the expansion coefficients.Traditional numerical methods with the low-order method of moments (MoM) [9], finite element method and finite difference method can also be applied to solve this problem, but are time consuming, although the fast multipole method [10] and multilevel fast multipole algorithm [11] have been applied successfully to accelerate the MoM.
Here we present a new perspective of this classical problem through fast surface integral equations for the unknown surface current density (rather than the fields everywhere).We develop a spectral integral method (SIM) for this scattering problem through the fast Fourier transform (FFT) algorithm and singularity subtraction to speed up the convergence of numerical solutions for the surface integral equations.Although the SIM was first named in [12], its elemental ideas go back to the Fig. 1.Configuration of P non-overlapping circular PEC cylinders of radius a p centered at O p = (ρ p , φ p ) (p = 1, 2, • • • , P ).Any point on the p-th cylinder can be conveniently written in its local polar coordinates as ρ = (a p , φ p ).The virtual outer circle with radius a 0 centered at the origin encloses all PEC cylinders.The z direction is perpendicular to the page.original work by Bojarski using FFT for a single cylinder [13], with the important singularity subtraction of the Green's function [14].The SIM was further developed for objects in a layered medium [15].Both [12] and [15] use an iterative method to solve the system equation through FFT, thus the computational complexity is O(KN log N ) for N surface unknowns, where K is the number of iterations.The FFT has also been used for the scattering of multiple cylindrical layers [16], although without treating the singularity of the Green's function and hence with slow numerical convergence.Recently, Zhu et.al. developed a direct SIM solution for the scattering of a circular cylinder with the singularity subtraction of the Green's function where no iteration is required [17]; it was further extended by Guan et.al. to multiple cylindrical layers [18].To our knowledge, however, the SIM has not been available for multiple objects.
In this work, we extend the SIM to analyze the scattering of multiple circular PEC cylinders.(Note that a similar problem in elastic waves has been treated by using the surface integral equation and addition theorem in [21]; it may also use the SIM developed in this work.)We first derive the coupled surface integral equations for P non-overlapping circular PEC cylinders.Through the addition theorem and Fourier series, we can rewrite the coupled integral equations into a set of algebraic equations in the spectral domain, which is then solved easily by the stabilized biconjugate gradient (BiCGSTAB) method [19] with an effective preconditioner.We show that this SIM needs a sampling density (SD) of only 2-3 points per wavelength (PPW) to reach the 99% engineering accuracy for the electric current density on the surfaces of the PEC cylinders.The new contribution of this work is the extension of the spectral integral method to an arbitrary number of circular PEC cylinders by using the fast Fourier transform and singularity subtraction, thus allowing a low sampling density close to the Nyquist rate of two points per wavelength.This new contribution paves the way to the future application of SIM to multiple cylinders with inhomogeneous media through its combination with the finite element method (FEM) or spectral element method (SEM) for the interior regions, as shown for the multiple layers of circular cylinders with inhomogeneities in [20].
The organization of the paper is as follows: In Section II, we present the detail formulations of the SIM.Numerical results are shown in Section III.Finally, concluding remarks are given in Section IV.

II. THEORY
Consider P circular PEC cylinders of radii a p and surface S p centered at O p = (ρ p , φ p ) (p = 1, 2, • • • , P ) in a global polar coordinate system with its origin at O 0 .The infinite background domain outside the PEC surfaces is filled with a homogeneous medium with permittivity ǫ = ǫ r ǫ 0 and permeability µ = µ r µ 0 .The objective of this work is to develop a fast method to solve for the scattered electromagnetic fields everywhere in space.For a TM z incident electric field E inc = ẑE inc z , the scattered electric field from these PEC cylinders can be written as where ρp × H| S p is the induced electric current density on the PEC cylinder surfaces, ρp is the unit outward radial vector in the p-th cylinder's local coordinates, H is the magnetic field, and is the two-dimensional Green's function of the background medium with the wavenumber k = ω √ µǫ.As the tangential electric field E z = E inc z + E sct z on the surfaces of all PEC cylinders must vanish, we have the electric field integral equation where R = |R|, R = ρ − ρ ′ .According to (1), equation (4) consists of the following P coupled electric field integral equations: To solve these P coupled integral equations, we focus on the scattered electric field produced on the q-th cylinder by the current on the p-th cylinder There are two situations depending on the values of p and q: the self interaction p = q and mutual interactions p = q.
Below we discuss their treatments within the SIM.

A. Self Interaction in the p-th Local Coordinate System
For p = q, we have the self interaction in (6) that describes the scattered field on the p-th cylinder caused by its own induced current.Naturally, it is more convenient to use the pth local polar coordinate system to express the above equation ( 6).This local coordinate system is centered at O p of the p-th cylinder, so we define the local polar position vector for the source on S p as ρ ′p = (ρ ′p = a p , φ ′p ), and that for the receiver on S p as ρ p = (a p , φ p ). Thus the self-induced scattered field on the p-th cylinder is 0 (kR pp )J p z (a p , φ ′p )dφ ′p (7) where for both the source and receiver on the p-th cylinder and the integral kernel H (2) 0 (kR pp ) is weakly singular because R pp can be zero.The treatment of this weakly singular problem using FFT has been solved by the previous SIM for a single PEC cylinder in [17], so we will skip its detail and only summarize the final result below.
Define the Fourier series and its Fourier coefficient for the induced electric current, and z (a q , φ q )e −jnφ q dφ q (10) for the scattered field, respectively.Then, from the well-known convolution theorem, the scattered field on S p due to J p z can be written explicitly in the spectral domain (i.e., Fourier domain) from ( 7) as ẽsct,pp where hpp n is the Fourier coefficient of H (2) 0 (kR pp ), which is given in detail in [17] with the singularity subtraction technique for treating the weak singularity of the Hankel function; such a singularity subtraction [14] is essential for the fast convergence of the SIM.Note that the weak singularity exists only for the self interaction, and not for the mutual interactions discussed below because R qp = 0 when q = p.

B. Mutual Interactions and Coordinate Translation
Now we focus on the mutual interaction case where p = q.In this case we define ρ > = max(ρ, ρ ′ ) and ρ < = min(ρ, ρ ′ ), so from the addition theorem (28) in Appendix we can rewrite (6) as in the global coordinate system.Alternatively, it is more convenient to write equation ( 12) in the p-th local coordinate system centered at O p as because for this case we know ρ p > ρ ′p as long as the cylinders do not overlap.Note that this is particularly convenient for the integral over S p as ρ ′p = (a p , φ ′p ) has a constant radius a p .On the other hand, E sct,qp z (ρ p ) written in the p-th local coordinates in (13) is not very convenient for an observation point on the q-th cylinder surface as |ρ p | is not a constant when q = p.However, we can use the q-th local coordinate system ρ q to express the observation point by noting that where ρ pq = O q − O p = (ρ pq , φ pq ) is the constant translation (displacement) vector from O p to O q (see Fig. 1).Furthermore, as ρ q < ρ pq because the P cylinders do not overlap, from the addition theorem in (33) in Appendix we can rewrite H (2) m (kρ p ) in (13) as m−n (kρ pq )e j(m−n)φpq (15) Therefore, in the q-th coordinate system equation ( 13) becomes m−n (kρ pq )e jnφ q e −jmφ ′p J p z (φ ′p )dφ ′p (16) With the Fourier series, equation ( 16) can be written in the spectral domain as This completes the spectral domain formulation for the mutual interactions between the p-th and q-th cylinders (p = q).

C. System of Equations and Its Solution
Using (11) and (17) in the integral equations (4) in the spectral domain, we obtain the following algebraic equations for the unknowns { jp m } in the spectral domain where and ẽinc,q n is the Fourier series coefficient of E inc z on the q-th cylinder surface.After we choose a sampling density, we can calculate the total number of sampling points on every circle.However, this number isn't always an integer, and therefore, we could find the closest odd to the number as the real total number of sampling points of that circle M 0 .And M = (M 0 − 1)/2.
Equation ( 18) has an infinite number of coupled equations, so its exact solution is difficult.However, numerically it is unnecessary to carry out this infinite sum; usually the summation over m is truncated to −M ≤ m ≤ M (i.e., M multipoles), and only N = 2M + 1 equations are kept for index n.Then the truncated equation ( 18) consists of P N unknowns { jq m } in the system.The SIM solves the truncated P N equations in (18) by using the stabilized biconjugate-gradient (BiCGSTAB) in the spectral domain, where { jp m } are solved.Here we can use the induced current densities in the isolated single PEC cylinders (i.e., no mutual interactions among the cylinders) as the initial solution, where q p m = −2/[πωµa p hpp m ].Furthermore, to speed up the convergence of the iterative solution, the diagonal matrix equivalent to the solution of isolated single PEC cylinders is used as the preconditioner.Thus, the truncated equation ( 18) is now solved by where the preconditioner is Q pq nm = q p m δ mn δ pq , j denotes the vector of jp m , and ẽinc denotes the vector of ẽinc,p m .For all the numerical examples in the next section, the number of iterations is only less than 8 for the residual error to decrease to less than 10 −6 .Then the FFT algorithm is used to convert the spectral domain current density to the spatial domain.If M multipoles are used in (19), the computational complexity is O(KP 2 N 2 ), where K (usually less than 8) is the number of iterations.

D. Far-Field Solution
It is often required to evaluate the far-zone field distributions.For this purpose, we denote O 0 = (0, 0) corresponding to q = 0 as the origin of the global coordinate system; thus ρ = ρ 0 .By virtue of equation (34) in Appendix, on a circle with a radius a 0 large enough (a 0 ≥ ρ p0 for all p) to enclose all the P circular cylinders, the addition theorem reads m−n (kρ 0 )e j(m−n)φ 0 (22) Thus, from (13) we have the scattered field given by the sum of all contributions from these cylinders: When a 0 = |ρ 0 | → ∞, we can use the asymptotic form of Hankel functions in (23) to arrive at where The radar cross section is thus given by where the incident electric field is that of a plane wave.Note that the near field in ( 23) is written in a double summation over m and n, while the far-field in (24) and the radar cross section are simplified to only two single summations A p and B p through (25).

III. RESULTS AND DISCUSSIONS
In this section, we show three numerical experiments to verify the accuracy and efficiency of the SIM.The first two examples compare the SIM and MoM results for the induced current density on the PEC surfaces; they verify that the developed SIM is more efficient and accurate for the EM scattering of multiple circular PEC cylinders.The third example shows the accuracy of the RCS solution from the SIM using the simplified far-field expression.The background medium of these three models is set as vacuum.The general configuration of our problem is shown in Fig. 1.Structures corresponding to every figure can be obtained by substituting the center positions and radii described in the corresponding caption into Fig. 1.A 100 MHz TM z plane wave of unit magnitude is incident along the +x direction with its phase being zero at the origin, i.e., E inc z = exp(−j 2π 3 x).The computation is performed by MATLAB R2020a on a computer with the system Microsoft Windows 10 Professional X64, 2019.The processor is Intel(R) Core(TM) i5-1035G4 CPU @1.10GHz, 1498 MHz and the memory is 16 GB.

A. Example 1: Three Identical PEC Cylinders
The first model consists of three circular PEC cylinders with the same radius a 1 = a 2 = a 3 = 5 m and centered at O 1 = (0, 0), O 2 = (0, 20) and O 3 = (35, 21) m, respectively, in Cartesian coordinates.SIM and MoM are used to calculate the induced current densities on these cylinders.The MoM used the 3-point Gaussian quadrature to calculate the integration of every segment.In order to obtain a solution with roughly the same accuracy, the SD is set as 30 PPW in MoM and 3 PPW in SIM (thus M = 15 for all cylinders), as shown in Fig. 2(a).Because the current density in both methods is calculated at the sampling points, a different SD leads to different observation points.Therefore, in order to compare the two results, we adopt the sampling points of SIM as our common points for comparison and then use cubic spline interpolation to obtain the corresponding values of these points in the results calculated by MoM.The relative difference of the current density errJ z,SM is defined as where J z,S and J z,M denote the current density computed by SIM and MoM, respectively.In this simulation, The 2-D MoM solver in FEKO requires 139.524 seconds of CPU time and 410.175MB memory, while the cost of SIM is 0.196 s and 24.03 MB, respectively.The relative difference between these two methods is 0.3%, indicating that SIM is more effective and less memory consuming than MoM, and their results are in good agreement with each other.Furthermore, the traditional low-order MoM requires a much higher SD to reach good accuracy.Based on the data above, we choose the SIM result at 20 PPW as the baseline reference to study the error convergence in SIM.Fig. 2(b) shows how the relative error of the SIM current density changes with the sampling density.It can be seen that it is an exponential convergence, and the SIM requires less than 3 PPW to reach 99% accuracy.
To further demonstrate the effects of the preconditioning with the initial single-cylinder solutions, we show the error convergence curves for the cases with and without the preconditioning in Fig. 2(c).It is shown that it takes only 6 iterations for the preconditioned solution to converge to a residual error of 10 −6 , compared to 32 iterations without preconditioning.

B. Example 2: Five Different PEC Cylinders
The second model consists of five PEC cylinders with different radii of 30, 18, 24, 12 and 36 m, and centered at positions Similar to the first model, the SIM and MoM are both used to calculate the induced current densities.The sampling density is set as 3 PPW for SIM and 30 PPW for MoM, respectively to maintain consistency in accuracy (thus M = 94, 56, 75, 37, 113 for the five cylinders, respectively).For clarity, only the current density on three of the five cylinders is shown in Fig. 3(a).The 2-D MoM solver in FEKO requires a much larger memory than that on our desktop computer.Therefore, we just used our own MoM code, which requires 40187.89s of CPU time and 8.38 GB memory, while the cost of SIM is 0.847 s and 0.20 GB, respectively.The relative difference between the two methods is only 0.3%.This result indicates that SIM is much more efficient than MoM and requires a much lower memory with approximately the same accuracy.
Based on the data above, we still choose the SIM result at 20 PPW as the baseline reference for studying the convergence of SIM.The error convergence curve of SIM is also exponential, as shown in Fig. 3(b).It can be seen that good accuracy of better than 99% can be achieved with a SD between 2 and 3 PPW.This indicates that the convergence rate of SIM is very fast.Through this study, we can see that SIM saves 90% sampling points compared with MoM when computing current densities.Moreover, benefiting from its much fewer sampling points, the semi-analytical method and characteristics of SIM itself, the proposed SIM is much more efficient than MoM.

C. Example 3: RCS of Five Identical PEC Cylinders
The above near field obtained by the SIM can be used to calculate the field everywhere, including the far field and radar cross section (RCS).Here we use the simplified formulas in (26) to obtain the RCS.The model used in this simulation consists of five PEC cylinders with the same radius a   3D) is used as a reference.Note that FEKO apparently is unable to solve such a special 3-D problem with the z-direction periodic boundary condition with the multilevel fast multipole method.The SD is set as 3 PPW for SIM (i.e., M = 19 for all cylinders), and the mesh is set as "fine" in FEKO.It is shown in Fig. 4(a) that the current densities computed by SIM and FEKO are in good agreement: the relative difference between them is only 0.75%.The bistatic RCS curves are obtained over the 360 degrees equally divided by 720 observation angles, as shown in Fig. 4(b).It can be observed that good agreement has been achieved, and the relative difference between the two solvers is only 0.924% for RCS.In this case, the SIM takes 0.555 s, while FEKO takes 389.789 s, so the SIM is 701 times faster than the 2-D MoM in FEKO.

IV. CONCLUSIONS
We have developed a spectral integral method for the scattering of an arbitrary number of non-overlapping circular PEC cylinders through the fast Fourier transform and addition theorem.It is an extension of the previous SIM for a single object.It is shown that this SIM can achieve very accurate results with a sampling density of only slightly above 2 points per wavelength.Numerically, to achieve an accuracy of 99%, it can be 100 times more efficient and 6 times less memoryconsuming than the MoM in the simulated cases.Thus, it provides an efficient and highly accurate technique to solve the scattering problem of multiple circular PEC cylinders.Furthermore, this SIM can be further extended to the cases of multiple dielectric cylinders.More importantly, using the SIM as an exact radiation boundary condition, we can combine it with the finite element method and spectral element method to obtain highly efficient solutions for multiple objects with arbitrary inhomogeneity and anisotropy.These will be part of our future work.

Fig. 2 .
Fig. 2. Example 1: Scattering of three circular PEC cylinders of the same radius of 5 m and centered at (0, 0), (0, 20) and (35, 21) m, respectively, in Cartesian coordinates.(a) |J z | on the PEC surface obtained by the SIM at 3 PPW and MoM at 30 PPW.(b) Relative error of J z versus the sampling density in SIM.(c) Comparison of the relative residual error versus iteration between the methods with and without preconditioner.