—An efﬁcient Ping-Pong algorithm is presented to deal with electromagnetic scattering from planar reﬂector antennas of arbitrary geometric shapes. The proposed Computational Electromagnetics (CEM) algorithm does not require formation, storage and inverse of the impedance matrix used by the Method of Moments (MoM). Instead, a rigorous formula is obtained for computation of the surface current from the scattering electric ﬁeld. Besides, the algorithm only contains convolution operations with the scalar Green’s function and second deriva- tives operations to iteratively update the surface currents and scattering electric ﬁeld. Also, hypersingular integration arising from the singularities of the surface current and scattering electric ﬁeld is properly regularized. Furthermore, the algorithm converges fast and satisfactory results are obtained after a few iterations. Besides, the ﬁrst iteration gives the solution for reﬂector antennas on perfect absorption substrates. Numerical solutions are obtained for some typical antennas and comparisons are made to those in the literature and from the MoM. Index Terms —Electromagnetic scattering, reﬂector anten- nas, convolution, singularity, hypersingular integration, scalar Green’s function, Method of Moments (MoM), Computational ElectroMagnetics (CEM). I. I NTRODUCTION PLANAR reﬂector antennas and arrays have many impor- tant applications such as smart phased array electromag- netic wave steering [1] and programmable metasurfaces for radiation pattern control [2]. Careful designs of the reﬂector antennas and arrays are required to achieve the speciﬁc perfor- mance such as low Radar Cross Section (RCS [3]) and better antenna radiation patterns [4]. However, the accurate simulation of the electromagnetic scattering of the planar reﬂector antennas and arrays is chal- lenging, partly due to the boundary edge effect of the reﬂector antennas, which usually generates singularities in surface cur- rent and scattering electromagnetic ﬁeld [5], [6]. So efﬁcient Computational ElectroMagnetics (CEM) methods are much needed. For electromagnetic scattering from metallic objects, the Method of Moments (MoM [7]) has been extensively used for computation of electromagnetic scattering of metallic objects. The MoM solves the Maxwell’s equations based on the the Electric Field Integral Equation (EFIE [8] ), one of the Surface Integral Equations (SIEs [9], [10]). For objects other than metals such as dielectric objects, the Magnetic Field Integral Equation (MFIE) and the combination of the EFIE and the MFIE, i.e., the Combined Field Integral Equation (CFIE) Shaolin Liao (Corresponding Author: liaoshlin@mail.sysu.edu.cn) is with School of Electronics and Information Technology, Sun Yat-sen University, Guangzhou, Guangdong Province, China. Lu Ou is with College of Computer Science and Electronic Engineering, Hunan University, Changsha, Hunan, China 410082. She was partly supported by China Postdoctoral Science Foundation No. 2020M672488.[11]. The MoM surface current solution of the EFIE can be obtained through inverse of the MoM impedance matrix that is constructed as follows: 1) surface meshing — dividing the metallic object surface into many electrically small unit cells; 2) surface current basis representation — then the surface current within each unit cell is expressed in terms of the sum of some basis functions like Rao-Wilton-Glisson (RWG) basis [12], [13]; 3) electromagnetic ﬁeld computation — electric ﬁeld on all surface unit cells due to the surface current is computed; 4) test of the MoM boundary condition — test basis functions are used to obtain the MoM matrix equation that consists of the left-side product of the surface current vector under its basis function and the impedance matrix constructed through the inner product the scattering electric ﬁeld and the test basis on the left and the surface current basis function on the right; as well as the right-side incident electromagnetic electric ﬁeld vector under the test basis; and 5) solving of the unknown surface current vector — perform inverse of the MoM impedance matrix to obtain the unknown surface current vector under its basis function through some numerical matrix solver such as the conjugate gradient method. The MoM CEM method has its own advantages compared with other methods such as the Finite Elements Method (FEM) and the Finite Difference Time Domain (FDTD) method. Firstly, the MoM only requires calculation of the electromag- netic ﬁeld and storage of the surface current on the boundary surface, greatly reducing the computational effort and stor- age resource. Secondly, fast electromagnetic ﬁeld propagation methods such as the Fast Multipole Method (FMM [14]) and the Fast Fourier Transform (FFT [10], [15]) method are readily available to accelerate the MoM impedance matrix computation. At last but not the least, efﬁcient matrix inverse algorithms [16], [17], [18] such as the iterative Generalized Minimal REsidual Method (GMRES [18]) can be used to solve the MoM equation. Even though the above advantages of the MoM, it has its own issues: ﬁrst, the surface has to be discretized into complex mesh unit cells such as RWG triangles. then the careful choice of the basis function is required for accurate representation of the problem. also, the large impedance matrix has to be stored during the computation, whose size grows asN2, withNbeing the number of unit cells. What’s more important, inverse of the impedance matrix is required — even with efﬁcient iterative algorithm such as GMRES, careful pre-conditioning and many iterations are required to obtain accurate solution. Finally, MoM doesn’t work so well when the surface current contains singularities, which is the case of the electromagnetic scattering from the reﬂector antennas. This is because, around the singularities of the surface current, much SUBMITTED TO IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATIONS, 2020 2 denser mesh is required. This means that the number of unit cellNbecomes too large, and the size of the impedance matrix N2grows so fast that the inverse of the impedance matrix becomes increasingly intractable. From these aspects, it would be really helpful if no impedance matrix forming and inverse are involved to obtain the solution of the surface current, which is the focus of this paper: to develop some electromagnetic algorithm without impedance matrix. The main contribution of this paper is summarized as follows, 1) A Ping-Pong algorithm without impedance matrix form- ing, storage and inverse: both the surface current and the scattering electromagnetic ﬁeld are updated alternatively and iteratively. 2) Rigorous formulas of the forward propagation operator and its inverse operator: the forward propagation opera- tor is to obtain the scattering electromagnetic ﬁeld from the surface current and the inverse operator is to obtain the surface current from the scattering electromagnetic ﬁeld. 3) Simple computation: only convolution with the scalar Green’s function and the following second derivatives are required to iteratively update the surface current and scattering electromagnetic ﬁeld. 4) Hypersingular integration regularization: the hypersin- gular integration arising from the singularities of the surface current and the scattering electric ﬁeld around the antenna boundary edges is properly regularized. 5) Fast convergence: the ﬁrst iteration of the Ping-Pong algorithm gives the solutions of reﬂector antennas on perfect absorption substrates and it usually takes a few iterations to converge to satisfactory results. The organization of this paper is as follows. Section II formulate the problem studied. Followed by Section III, where the electromagnetic ﬁeld is represented in the 2D Fourier spectral domain. Then, Section IV presents the Ping-Pong algorithm. After that, numerical validation is shown in Section V, followed by discussion in Section VI. Finally, Section VII concludes the paper. II. P ROBLEM FORMULATION The electromagnetic scattering problem is shown in Fig. 1. The electromagnetic source ﬁeld is incident upon some planar metallic reﬂector antenna (the square, disc, star and triangle here) and the goal is to solve for the unknown surface current. Imposing zero total electric ﬁeld boundary condition on the metallic surface, the EFIE is obtained, as given Eq (1), h Es ==(x;y) +Ei ==(x;y)i (x;y) = 0; (1) where is the binary mask that deﬁnes the geometry area of the reﬂector antenna; the superscripts sandidenote the scattering and incident electric ﬁeld, respectively. Also, Es ==(x;y)andEi ==(x;y)are the scattering electric ﬁeld and incident electric ﬁeld on xyplane without the ^zcomponent, respectively, Es;i ==(x;y) = ^xEs;i x(x;y) + ^yEs;i y(x;y): Figure 1. Electromagnetic scattering of planar reﬂector antennas of arbitrary shapes (clockwise): a square antenna of length Lx=Ly= 1, a circular disc antenna of radius R= 5=(2)orkR= 5, a star antenna of radius R=, and an equilateral triangle antenna of radius R= 1. Dashed lines denote where the surface currents are evaluated in the numerical simulation. The scattering electric ﬁeld Escan be expressed in terms of convolution of the surface current Jand the electric dyadic Green’s function Ge, as shown in Eq (2), Es(x;y) j!=Ge(x;y)~J(x;y) (2) =1Z 11Z 1Ge(xx0;yy0)J(x0;y0)dx0dy0; where ~denotes the 2D convolution operation; !is the angular frequency; is the permeability; and the electric dyadic Green’s function is given as Eq (3), Ge(r) =g(r)I+1 k2rrg(r); (3) whereIis the identity matrix; kis the magnitude of the wave vectork= [kx;ky;kz]; andg(r)is the scalar Green’s function at the observation point vector r= [x;y;z ], g(r) =ejkr 4r; r=jrj;k=jkj=!p: III. E LECTROMAGNETIC FIELD IN THE 2D F OURIER SPECTRAL DOMAIN Due to the planar geometry of the reﬂector antennas, it is convenient to express the electromagnetic ﬁeld in terms of plane wave superposition in the Fourier spectral domain [19], [20]. A. Spectrum of the Scalar Green’s Function According to Weyl’s identity [21], [22], [23], the scalar Green’s function can be expressed in terms of plane wave SUBMITTED TO IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATIONS, 2020 3 superposition as follows, g(rr0) =ejkjrr0j 4jrr0j=1Z 11Z 1ejkr0 j4kzejkrdkxdky;(4) which is valid when the observation point ris above the source point r0,i.e.,z > z0; alsokzis the wavevector along zdirection given by, kz=8 < :q k2k2xk2y forq k2x+k2yk, jq k2x+k2yk2forq k2x+k2y>k. From Eq. (4), the 2D Fourier transform pair is then obtained as follows, g(x;y)()G (kx;ky) =j 4kzej(kxx0+kyy0+kzz0);(5) where G(kx;ky) =1 2Z1 1Z1 1g(x;y)ej(kxx+kyy)dxdy; g(x;y) =1 2Z1 1Z1 1G(kx;ky)ej(kxx+kyy)dkxdky: Similarly, when the observation point ris beneath the source point r0,i.e.,z < z0, the 2D Fourier spectrum of the scalar Green’s function is given by g(x;y)()G (kx;ky) =j 4kzej(kxx0+kyy0kzz0);(6) In particular, when the source point is set to the origin r0= 0, the following is obtained, g(x;y)()G (kx;ky) =j 4kz: (7) B. Spectra of Derivatives of the Scalar Green’s Function With the 2D Fourier spectrum of the scalar Green’s function given in Eq. (5), the spectra of the derivatives of the scalar Green’s function are readily obtained. Lemma 1. The 2D Fourier transform pairs of the derivatives of the scalar Green’s function are @n @xngx(x;y)()G(n) x(kx;ky) =(j)n+1 4kn x kz; @n @yngy(x;y)()G(n) y(kx;ky) =(j)n+1 4kn y kz; @n@n0 @xn@yn0g(x;y)()G(n;n0) x;y (kx;ky) =(j)n+n0+1 4kn xkn0 y kz; where the superscript ndenotes the order of derivative and the subscript x;ydenotes the variable the derivative takes on. Proof: Taking the derivatives on the observation point r of Eq. (4), Lemma 1 can be proved.C. Spectrum of the Electric Dyadic Green’s Function The 2D Fourier spectrum of the electric dyadic Green’s function can be obtained from its deﬁnition in Eq. (3) and Lemma 1 . Lemma 2. The 2D Fourier transform pair of the electric dyadic Green’s function is Ge(x;y)()Ge(kx;ky) =2(z) k2^z^zj 4kz I^k^k ; where ^kis the unit wave vector ^k=k=k; and ^k^kcan be considered as a matrix working on an input current spectrum vector on its right side (the dot product of the right k) and producing an output electric ﬁeld spectrum vector (the left k). Proof: The ﬁrst term of the electric dyadic Green’s functionGein Eq. (3) is the product of the scalar Green’s function with a constant identity matrix, whose 2D Fourier spectrum is justG(kx;ky)I. The second term of the electric dyadic Green’s function Gecan be expanded into the sum of all second derivatives whose 2D Fourier spectra can be obtained from Lemma 1 . At last, the ^z^zterm is due to the change of electromagnetic wave propagation direction from +^zto the^z, as shown in Eq. (5) and Eq. (6), respectively. D. Spectrum of the Scattering Electric Field The 2D Fourier spectrum of the scattering electric ﬁeld Es in Eq. (2) is the product of the spectrum of the surface current Jand the spectrum of the electric dyadic Green’s function. Lemma 3. The 2D Fourier spectrum of the scattering ﬁeld is given by Es(kx;ky) =! 2kz I^k^k J(kx;ky): (8) Proof: According to the convolution property of the Fourier transform, the spectrum of the scattering ﬁeld is the product of the spectrum of the surface current and the spectrum of the dyadic Green’s function given in Lemma 2 . IV. T HEPING-PONG ALGORITHM OF THE ELECTROMAGNETIC SCATTERING PROBLEM In the proposed Ping-Pong algorithm, the solution of the scattering electromagnetic ﬁeld can be obtained ﬁrst in the 2D Fourier spectral domain and then transformed back to the spatial domain. A. Scattering Electric Field Radiated by Surface Currents The scattering electric ﬁeld can be expressed as the convo- lution of the surface current with the dyadic Green’s function in Eq. (2). It contains higher-order singularities due to the derivatives on the scalar Green’s function [26], which can be reduced to only convolutions with the scalar Green’s function. SUBMITTED TO IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATIONS, 2020 4 Lemma 4. The scattering electric ﬁeld can be expressed in terms of only convolution with the scalar Green’s function kernel and second derivatives after that, Es(x;y) =L J(x;y) (9) =! j8 >>< >>:I+1 k22 664@2 @x2@2 @x@y@2 @x@z @2 @x@y@2 @y2@2 @y@z @2 @x@z@2 @y@z@2 @z23 7759 >>= >>; g(x;y)~J(x;y) ; whereLis the forward propagation operator. Proof: The differential operators are on the source point coordinater0of the scalar Green’s function in Eq. (2), which can be ﬁrst transferred to the observation point coordinate r, Es(x;y) j!=Ig(x;y)~J(x;y)+ (10) 1 k21Z 11Z 1rrg(xx0;yy0)J(x0;y0)dx0dy0: where the following relations have been used, r0g(xx0;yy0) =rg(xx0;yy0); r0r0g(xx0;yy0) =rrg(xx0;yy0): withr0denoting the derivative on the source point r0. Then the derivatives on the observation point coordinate r in Eq. (10) can be further transferred outside the integral and Lemma 4 is proved. In particular, for the problem of 2D electromagnetic scat- tering from planar reﬂector antennas, Eq. (9) in Lemma 4 reduces to the following, Es(x;y) =L J(x;y) (11) =! j8 >>< >>:I+1 k22 664@2 @x2@2 @x@y @2 @x@y@2 @y2 @2 @z@x@2 @z@y3 7759 >>= >>; g(x;y)~J(x;y) : B. Surface Current from the Scattering Electric Field The the surface current can be obtained from the scattering electric ﬁeld Es, which is the inverse process of Lemma 4 . Lemma 5. The surface current can be expressed in terms of the scattering electric ﬁeld as J(x;y) =L1n Es(x;y)o (12) =4! j( I+1 k2"@2 @y2@2 @x@y @2 @xy@2 @x2#)h g(x;y)~Es ==(x;y)i ; whereL1is the inverse operator of the forward propagation operator Lgiven in Eq. (9) of Lemma 4 ;is the permittivity. Proof: With the help of Lemma 3 , the 2D Fourier spectrum of the surface current J(kx;ky)can be solved inthe spectral domain of the EFIE equation in Eq. (1), Jx(kx;ky) =2 !" k2k2 y kzEs x(kx;ky) +kxky kzEs y(kx;ky)# ; Jy(kx;ky) =2 !k2k2 x kzEs y(kx;ky) +kxky kzEs x(kx;ky) ; from which it can be seen that the 2D Fourier spectrum of the surface current can be re-written in terms of the product of the 2D Fourier spectrum of the scalar Green’s function and the 2D Fourier spectrum of the incident electric ﬁeld, as well as thekm xkn y;m= 1;2terms, Jx(kx;ky) =j8 ! k2k2 y G(kx;ky)Es x(kx;ky)+ kxkyG(kx;ky)Es y(kx;ky) ; Jy(kx;ky) =j8 ! k2k2 x G(kx;ky)Es y(kx;ky)+ kxkyG(kx;ky)Es x(kx;ky)]: Then the 2D inverse Fourier transform in the spatial domain gives the convolution of the scalar Green’s function with the incident electric ﬁeld, 2G(kx;ky)Es u(kx;ky)()g(x;y)~Es u(x;y); u=x;y: Finally thekm xkn y;m= 1;2terms correspond to the deriva- tives of the convolution in the spatial domain and Lemma 5 is proved. Eq. (12) in Lemma 5 is believed to appear for the ﬁrst time and is the main result of this paper that leads to the following Ping-Pong algorithm. C. The Ping-Pong Algorithm With the help from Lemma 4 andLemma 5 , the solutions of the surface current and scattering electric ﬁeld can be obtained according to the EFIE, as shown in Eq. (1). Theorem 6. The surface current and scattering electric ﬁeld can be updated iteratively as follows, Jn+1(x;y) =L1n Es ==;n(x;y) Ei ==(x;y)o ;(13) Es n+1(x;y) =L Jn+1(x;y) ; (14) where = 1 is the complementary area of the reﬂector antenna. Proof: There are two constraints according to the EFIE equation in Eq. (1): 1) the scattering electric ﬁeld cancel the incident electric ﬁeld on the surface of the reﬂector antenna, i.e.,Es ==(x;y) =Ei ==(x;y); and 2) the surface current only exists on the surface of the reﬂector antenna, i.e., J(x;y) = 0. Eq. (13) is obtained from constraint 1) according to Eq. (12) of Lemma 5 and Eq. (14) is obtained from constraint 2) according to Eq. (9) of Lemma 4 . In particular, the ﬁrst iteration of Theorem 6 gives solution for the perfect absorption substrate: the scattering electric ﬁeld SUBMITTED TO IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATIONS, 2020 5 only exists on the reﬂector antenna surface and disappear outside the reﬂector antenna surface, Es ==(x;y) = Ei ==(x;y); Es ==(x;y) = 0; with the solution of the surface current given by, J1(x;y) =L1n Ei ==(x;y)o : (15) D. Singularities Regularization Besides the regular singularity of the scalar Green’s func- tion, the surface current and scattering electric ﬁeld could contains singularities around the antenna boundary edges due to the the derivatives operations after the convolution with the scalar Green’s function, as shown in the forward propagation operator in Lemma 4 and the inverse operator in Lemma 5 . These singularities have to be regularized before integration [24], [25], [26], [27]. Two types of singular integration are met during the iteration of the Ping-Pong algorithm: 1) Regular singular integration: for observation point rthat is not at the singularities of the surface current or the scattering electric ﬁeld around the antenna boundary edges@ , the integration can be readily regularized [27]. Because there are two types of regular singularities, one from the scalar Green’s function and the other from the surface current or scattering electric ﬁeld, the integration domain can be divided into two sub-domains around these two singularities, besides the rest of the sub- domain. For example, consider the point-type singularities of the surface current around the antenna corners in the convolution integration of the surface current with the scalar Green’s function of Lemma 4 : First, in the sub- domain around the singularity of the scalar Green’s function, denoted as g, I g ~=g(x;y)~J(x;y) g =ZZ gejkr0 4J(r0; 0)dr0d 0; where the observation point has been set to the origin (0;0). Then, in the sub-domain around the singularity of the surface current, denoted as J, I J ~=ZZ Jejkr0 4jr0roj r0J(r0; 0) dr0d 0; whereI J ~is deﬁned similarly as I g ~; also the origin of the polar coordinate has been set to the singularity point of the surface current and rois the observation point in this local coordinate. Now the regularized surface currentr0J(r0; 0)should contain no singularities and the integration can be readily evaluated. Similarly, the line-type singularities of the surface cur- rent around the edges of the antenna boundary edges can be done similarly [27] in the local Cartesian coordinate.Finally, the integration in the sub-domain other than g and Jcan be done through quadrature methods. 2) Hypersingular integration: for observation point rthat is at the singularities of the surface current or the scattering electric ﬁeld around the antenna boundary edges@ , both the regular singularity of the scalar Green’s function and the singularities of the surface current or the scattering electric ﬁeld exist. This is the hypersingular integration [28] that can be done in the following steps: 1) similar to the regular singular inte- gration above, ﬁrst remove the regular singularity of the scalar Green’s function in the polar coordinate; 2) then, divide the integration domain into two sub-domains: a small disc sub-domain around the singularity with radius r, dentoed as , and the rest sub-domain; 3) after that, the radius of the small disc sub-domain is chosen such that the rest of the sub-domain can be accurately evaluated with the conventional quadrature integration method; 4) ﬁnally, the hypersingular integration [28] on the small disc sub-domain is performed. It is noted that only the singularity with 1=(r0)p;p1generates hypersingular equation; also, it is required that p<2or otherwise the solution will diverge everywhere. So only the hypersingular integration with 1p<2is relevant, I  ~=ZZ ejkr0 4J(r0; 0)dr0d 0; =ZZ ejkr0 4J0(r0; 0) (r0)pdr0d 0; whereI  ~is deﬁned similarly as I g ~and1p < 2is the order of singularity such that J0(r0; 0) = J(r0; 0)(r0)pis not singular. Then, the integration can be evaluated through Taylor expansion [28] of J0(r0; 0)ejkr0, which contains the following singular part, rZ r0=!02Z 0=0J0(0;0) 41 r0dr0d 0=Cln (r)Cln (); forp= 1 withC=J0(0;0)=2; and rZ r0=!02Z 0=0J0(0;0) 41 (r0)pdr0d 0=C0 (r)p1C0 p1; for1< p < 2withC0=J0(0;0)=[2(1p)]. The second term is the scattering electric ﬁeld singularity when!0:ln()forp= 1 and1=p1forp>1. Similarly, the hypersingular integration due to the line- type singularities of the surface current around the antenna edges can be evaluated in the local Cartesian coordinate. At last, similar treatment can be done for the hypersingular convolution integration of the singular scattering electric ﬁeld with the scalar Green’s function in Eq. (12) of Lemma 5 . SUBMITTED TO IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATIONS, 2020 6 Figure 2. A square reﬂector antenna of length Lx=Ly= 1: surface currentsJxandJy. Figure 3. A square reﬂector antenna of length Lx=Ly= 1: line plots of surface current JyJxalong ^x=^ydirections and comparison with results from the literature [5]. V. N UMERICAL SIMULATIONS Numerical solutions of four typical reﬂector antennas are obtained with the Ping-Pong algorithm: a square reﬂector antenna of length Lx=Ly= 1, a circular reﬂector disc antenna of radius R= 5=(2)orkR= 5, an equilateral triangle reﬂector antenna of radius R= 1and a star reﬂector antenna of radius R= 1, as shown in Fig. 1. Comparisons of the solutions with those in the literature and from the MoM are also made. Also, without loss of generality, normal plane wave incidence with ^ypolarization is used. Firstly, Fig. 2 shows the surface currents JxandJyfor the square reﬂector antenna of length Lx=Ly= 1, obtained with ﬁve iterations. It can see that the the surface current Jy(normalized to the incident magnetic ﬁeld Hi) disappears on the antenna boundary edges that are perpendicular to ^y direction, which is due to the ^ypolarization of the incidence electromagnetic wave. Also, the surface current Jyhas line- type singularities around the antenna boundary edges that are perpendicular to ^xdirection and the the surface current Jx has point-like singularities around the antenna corners. Then, Fig. 3 shows the comparison of results obtained from the Ping-Pong algorithm and those from the literature [5] using the MoM with a mesh of 90 rectangular patches: dominant surface current Jyis shown along both the ^xand^ydirections, as shown in Fig. 1, the red dotted line and blue squares are for Figure 4. A disc reﬂector antenna of radius R= 5=(2)orkR= 5: surface currents JxandJy. Figure 5. A disc reﬂector antenna of radius R= 5=(2)orkR= 5: line plots of surface current Jyalong ^x=^ydirections and comparison with results from the literature [6]. the^xdirection and the purple dotted line and green squares are for the ^ydirection. It can be seen that good agreement has been achieved. Then, in Fig. 4, the surface currents JxandJyare shown for the circular disc antenna of radius R= 5=(2)orkR= 5, obtained with three iterations. It can be seen that singular surface current Jyhappens on the antenna boundary edges that are perpendicular to the ^xdirection, due to the ^y-polarized incidence ﬁeld. The surface current Jyshows the resonant characteristics distribution: two surface current peaks are on each side of the disc reﬂector antenna. Also, Fig. 5 shows the comparison of results obtained from the Ping-Pong algorithm and those from the literature [6]: dominant surface current Jy is shown along both the ^xand^ydirections. As shown in Fig. 1, the red dotted line and blue squares are for the ^xdirection and the purple dotted line and green squares are for the ^y direction. It can be seen that all features of the surface current are obtained by the Ping-Pong algorithm. Also, Fig. 6 shows surface current JxandJyfor the triangle reﬂector antenna of radius R= 1, obtained with four iterations. Similarly, it can be seen that the surface current disappear on the antenna boundary edge that is perpendicular to the electric ﬁeld polarization direction ^y. Also, line-type singularities are observed along the edges that are not perpen- dicular to the electric ﬁeld polarization direction ^yand point- type singularities are observed around the antenna vertices. SUBMITTED TO IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATIONS, 2020 7 Figure 6. A triangle reﬂector antenna of radius R= 1: surface currents JxandJy. Figure 7. A triangle reﬂector antenna of radius R= 1: line plots of surface currentJyalong ^x=^ydirections and comparison with results of the MoM. In particular, the surface current is the strongest at the top vortex of the triangle reﬂector antenna, which is the joining point of the surface current ﬂows along both edges of the triangle reﬂector antenna. To show the validity of the Ping- Pong algorithm, the comparison with the results obtained from the MoM has been made in Fig. 7 for the dominant surface currentJyalong both ^xand^ydirections (see Fig. 1): the red dotted line and blue squares are for ^xdirection and the purple dotted line and green squares are for the ^ydirection. It can be seen that the solution of the Ping-Pong algorithm can gives better result around the singularity point. This is mainly due to the ﬁnite size of the mesh used in the MoM, limiting its performance around the singularities. Finally, Fig. 8 shows the surface currents JxandJyfor the star reﬂector antenna of radius R= 1, obtained with 3 iterations. Again, it can be seen that line-type singularities happen around the antenna boundary edges and point-type singularities at the vertices of the star reﬂector antenna. Also, Fig. 9 shows the comparison of the results between the Ping-Pong algorithm and the MoM, which again shows good agreements in dominant surface current Jyin both ^xand^y directions (see Fig. 1): the red dotted line and blue squares are for^xdirection and the purple dotted line and green squares are for the ^ydirection. At last, to show the convergence and the efﬁciency of the Ping-Pong algorithm, Fig. 10 compare the convergence Figure 8. A star reﬂector antenna of radius R= 1: surface currents Jx andJy. Figure 9. A star reﬂector antenna of radius R= 1: line plots of surface currentJyalong ^x=^ydirections and comparison with results of the MoM. performance of both algorithms, which clearly shows that the Ping-Pong algorithm outperforms the MoM: it converges in around 5 iterations, while the MoM converges in around 50 iterations, meaning 10 times as fast. VI. D ISCUSSION The proposed Ping-Pong algorithm is very efﬁcient for the electromagnetic scattering of 2D metallic structures: it works in a straight-forward iteration way to sequentially update the surface current and the scattering electric ﬁeld, with each other in the last update as input. Also, usually only a few iterations give satisfactory results. Besides, the ﬁrst iteration of Theorem 6is the solution of the perfect absorption substrate, as shown in Eq. (15), which provides good initial solution of the Ping- Pong algorithm. In addition, the numerical simulation conﬁrms the singu- larities of the surface current around the antenna boundary edges that are not perpendicular to the polarization direction of the incidence electric ﬁeld, as well as around the vertices of the antenna. These singularities generate hypersingular integration during the iteration of the Ping-Pong algorithm, whose treatment has been discussed in Section VI. More advanced hypersingular integration methods [29] would deﬁ- nitely improve the algorithm. What’s more, fast algorithms such as FFT [10], [15] and FMM [14] could be used to accelerate the computation of the SUBMITTED TO IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATIONS, 2020 8 Figure 10. Convergence of the Ping-Pong algorithm and the MoM shows the efﬁciency of the Ping-Pong algorithm. convolution with the scalar Green’s function, further improv- ing the efﬁciency of the Ping-Pong algorithm. At last, although numerical simulations are only performed for individual reﬂector antennas in this paper, the Ping-Pong algorithm is expected to work equivalently well for 2D metal- lic antenna arrays and metasurfaces. VII. C ONCLUSION The efﬁcient Ping-Pong algorithm has been studied for the electromagnetic scattering problem of planar reﬂector antennas of arbitrary shapes. The algorithm does not compute, store and perform inverse of the impedance matrix required in the MoM. Instead, rigorous formulas have been obtained for the forward propagation operator for computation of the scattering electric ﬁeld from the surface current and its inverse operator for computation of the surface current from the scattering electric ﬁeld. Then, the algorithm interactively updates the surface current and scattering electric ﬁeld with the forward propagation operator and its inverse, which consists of only convolutions with the scalar Green’s function and the fol- lowing simple second-order derivatives. Also, only regular singularities of the scalar Green’s function is involved in the integral, greatly simplifying the integration. In addition, hypersingular integration that arises from the singularities of the surface current and the scattering ﬁeld has been properly regularized. What’s more, the Ping-Pong algorithm converges very fast and usually takes only a few iterations to obtain satisfactory results. In addition, the ﬁrst iteration gives the solution of reﬂector antennas on perfect absorption substrates. Numerical simulation has been performed for four typical antennas, i.e., a square, a disc, a triangle and a star. The comparisons to those of the literature and from the MoM show good agreements. At last, the algorithm is universal and very promising in simulation of other 2D metallic structures, including antenna arrays and metasurfaces.