Efficient Spectral Domain MoM for the Design of Circularly Polarized Reflectarray Antennas Made of Split Rings

1 Efficient spectral domain MoM for the design of circularly polarized reflectarray antennas made of split rings Rafael Florencio, Rafael R. Boix, and José Encinar, Abstract—The Method of Moments (MoM) in the spectral domain is used for the analysis of the scattering of a plane wave by a multilayered periodic structure containing conducting concentric split rings in the unit cell. Basis functions accounting for edge singularities are used in the approximation of the current density on the split rings, which makes it possible a fast convergence of MoM with respect to the number of basis functions. Since the 2-D Fourier transforms of the basis functions cannot be obtained in closed-form, judicious tricks (controlled truncation of infinite summations, interpolations, etc.) are used for the efficient numerical determination of these Fourier transforms. The implemented spectral domain MoM software has been used in the design of a circularly polarized reflectarray antenna based on split rings under the local periodicity condition. The antenna has been analyzed with our spectral domain MoM software, with CST and with HFSS, and good agreement has been found among all sets of results. Our software has proven to be around 27 times faster than CST and HFSS. Index Terms—Multilayered media, moment methods, periodic structures, reflectarrays, circular polarization. I. Introduction REFLECTARRAY antennas are an interesting alternative to reflector antennas and microstrip arrays owing to their versatile radiation performance, low profile, light weight, ease of fabrication, simplified feeding system, etc. [1]. In the particular case of circular polarization (CP) reflectarrays, the variable rotation technique (VRT) introduced by Huang and Pogorzelski is a design approach that has proven to be very successful [2]. In the frame of this technique, each element of the reflectarray is rotated a different angle to achieve the adequate phase shift in the impinging circularly polarized wave that makes it possible to generate the prescribed radiation pattern after reflection. The VRT requires that the reflection coefficients of the two orthogonal linear components of the impiging wave electric field are of equal magnitude and 180◦ out of phase in order to keep the sense of circular polarization (LHCP to LHCP or RHCP to RHCP) after reflection [2], [3]. One reflectarray element that has turned out to be especially suitable for CP applications is that based on split rings. Han et al. used the VRT and single split rings to design a Kaband RHCP reflectarray at 31.75 GHz with measured crosspolarization level below 40 dB at broadside [4]. Stacked split rings of different size were used in [5] to design a dual frequency RHCP reflectarray operating in the C and Ka bands. Smith et al. [6] designed a single layer dual frequency CP reflectarray by using an element with two concentric split rings. The inner rings and the VRT were used to adjust the elements phase at 29.75 GHz for RHCP radiation, while the outer rings and the VRT were used at 19.95 GHz for LHCP radiation. An extension of this work was presented in [7] where the dual frequency CP reflectarray was combined with an FSS made of non-split rings to combine the operation of the reflectarray with that of a second CP microstrip array antenna working at L band. Zhao et al. used the single layer element based on two concentric rings and the VRT to design a broadband CP reflectarray at a center frequency of 10 GHz with a bandwidth larger than 30% for 1-dB gain variation, and a bandwidth larger than 40% for an axial ratio smaller than 3-dB [8]. In all previous papers dealing with CP applications, the elements of the reflectarrays consist of symmetric split rings with two gaps and one axis of mirror symmetry. Zhang et al. proposed in [9] the use of non-symmetric single gap split rings which can be used as phasing elements for both linear polarization (LP) and CP reflectarray operation. By using this element, a LP reflectarray was designed at 20 GHz with a bandwidth of 23% for a 1-dB gain variation and an aperture efficiency of 60%. In the design of a CP reflectarray antenna made of split rings, the use of the VRT requires to adjust the size of the gaps and their orientation for every single element in order to achieve the required phase shift in the impinging circularly polarized wave, and in order to ensure a phase difference of 180◦in the reflection coefficients of its two orthogonal linear components. When estimating the gap size and their orientation in each element, it is customary to assume that the elements are located in a periodic environment, which is called the local periodicity assumption [10], [11]. The validity of this assumption is justified by the fact that leads to reflectarray designs in which the simulated performance agrees very well with the measured performance [5], [7], [8]. The design of a midsize CP reflectarray antenna made of split rings under the local periodicity assumption requires the analysis of a large number of different multilayered periodic structures owing to the wide variety of gap sizes and orientations that have to be adjusted in the split rings of the different reflectarray elements for phase adjustment [7]. Owing to this, a powerful numerical tool is needed for the analysis of these multilayered periodic structures. In this paper we apply the Method of Moments (MoM) in 2 the spectral domain to the analysis of the scattering of a plane wave by a multilayered periodic structure including concentric split rings in the unit cell. The implemented MoM software is subsequently used in the design of a circularly polarized reflectarray antenna made of split rings. When applying the spectral domain MoM, the 2-D Fourier transforms of the basis functions are expressed as infinite summations of Hankel transforms. This Hankel transform approach was used in the past in conjunction with the spectral domain MoM for the determination of the resonant frequencies of microstrip ring resonators [12], the input impedance and radiation patterns of single and stacked annular-ring microstrip antennas [13], [14], the resonant frequencies of shielded single and coupled microstrip ring resonators [15], [16], and the radar crosssection of concentric annular microstrip rings [17]. However, to the authors’ knowledge, it has never been used for microstrip split rings and for multilayered structures in periodic environments. We have only found one paper dealing with the analysis of microstrip split ring resonators in which a magnetic wall approximate model is used for the determination of the resonant frequencies, and where effective dimensions are introduced to account for edge effects [18]. Whereas most papers dealing with the analysis of microstrip ring structures use either magnetic wall cavity mode functions [13], [19] or subsectional piecewise functions [15], [17] as basis functions for the electric current density on the rings, in this paper we have used basis functions that account for edge singularities [16], [19] since these basis functions provide a fast convergence of MoM with respect to the number of basis functions, which implies that only small MoM matrices have to be inverted. Since the Hankel transforms of the edge singularity basis functions cannot be obtained in closed form, we have used especially tailored quadrature rules for their numerical computation, which provide very accurate results with a small number of quadrature points. Also, we have introduced judicious stop criteria for the summation of the infinite series leading to the 2-D Fourier transforms of the basis functions, and we have finally carried out Chebyshev interpolations of these Fourier transforms as a function of the spectral radial variable. As a result of all these strategies, we have developed a very efficient spectral domain MoM approach for the analysis of multilayered periodic structures containing split rings. The spectral domain MoM software has been used in the design of a circularly polarized reflectarray antenna and the results obtained have been compared with results provided by commercial softwares CSTrand HFSSr. Good agreement has been found between the three sets of results, our software being between one and two orders of magnitude faster than CSTrand HFSSr. II. Numerical procedure Figs. 1(a) and (b) show a multilayered periodic structure backed by a ground plane. Concentric conducting split rings are located in each unit cell, and are printed on the upper layer. The center of the rings has been chosen to be coincident with the geometrical center of the unit cell in the z=0 plane. The conducting rings and the ground plane will be assumed to be perfect electric conductors (PEC). In a generic case, there would be Rconcentric rings in the unit cell (R=2 in Fig.1(b)). Also, each ring would be split into Sarcs (S=2 in Fig.1(b)) so that the total number of arcs in the unit cell would be L=RS (L=4 in Fig. 1(b)). Let ρ1land ρ2lbe the inner and outer radius of the l-th arc respectively (l=1,...,L), and let wl=ρ2l−ρ1lbe the width of the l-th arc. The two ends of the l-th arc are characterized by the angles ϕ1land ϕ2l, such that ϕ1l<ϕ2l. The angle ϕkl (k=1,2; l=1,...,L) is the angle subtended between the k−th end of the l-th arc and the semiinfinite line directed along the positive xaxis with origin at the center of the ring. The angle ϕkl is taken as positive if, starting from the x-directed semi-infinite line, it is drawn in the counterclockwise sense, and is taken as negative if it is drawn in the clockwise sense (e.g., in Fig. 1(b), ϕ21,ϕ12,ϕ22, ϕ13 and ϕ23 are positive, but ϕ11,ϕ14 and ϕ24 are negative). z xy qinc em 10 em 20 , ,h2 h1 b bbbbb (a) b a y zx a j 22 2r 23 2r j 13 j 23 j 12 a2 a1 21 j 21 j 11 j 24 j 14 (b) Fig. 1. Two-layered periodic structure. The unit cell contains two concentric split rings with different axes of mirror symmetry. A plane wave impinges on the multilayered periodic structure. (a) Side view. (b) Top view. In the particular cases treated in this paper, we will assume that each individual split ring of the periodic structure has an axis of mirror symmetry going through the center of the ring (in the case of Fig. 1(b), this means that ϕ21 −ϕ11=ϕ22 −ϕ12 and ϕ12 −ϕ21=360◦+ϕ11 −ϕ22, and that ϕ23 −ϕ13=ϕ24 −ϕ14 and 360◦+ϕ14 −ϕ23=ϕ13 −ϕ24). For the r-th ring (r=1,...,R), this mirror symmetry axis makes an angle αr(−180◦≤αr<180◦) with the semi-infinite line directed along the positive yaxis with origin at the center of the ring (see Fig. 1(b)). The angle αris taken as positive if, starting from the y-directed semi-infinite line, it is drawn in the counterclockwise sense, and is taken as negative if it is drawn in the clockwise sense (e.g., in Fig. 1(b), α1and α2are both positive). Please note that the restriction introduced by these mirror symmetry axes is not necessary in the mathematical derivations presented in the rest of this section. The concentric split rings of the periodic structure are printed on a multilayered substrate containing Nllayers (Nl= 3 2 in Fig.1(a)). The i-th layer has a thickness hiand a complex permittivity εi=ε0εr,i(1 −j tan δi) (i=1,...,Nl). In the following, a time dependence of the type ejωtwill be assumed and suppressed throughout. As commented in the introduction, the unit cell of the multilayered periodic structure of Figs. 1(a) and (b) has been used as constituent element of several CP reflectarray antennas [6]–[8]. In the design of a CP reflectarray antenna under the local periodicity assumption, each cell is characterized by means of a CP 2 ×2 complex reflection matrix, RCP, which relates the RHCP and LHCP components of the reflected and incident electric fields when a plane wave impinges on the cell surrounded by a periodic environment. In this paper, this impinging plane wave is assumed to propagate in an arbitrary incidence direction given by the angular spherical coordinates θinc and ϕinc (see Fig.1(a)). If Eref RHCP and Eref LHCP are the RHCP and LHCP complex components of the reflected electric field, and Einc RHCP and Einc LHCP are the RHCP and LHCP complex components of the incident electric field, then Eref RHCP Eref LHCP !=RCP · Einc RHCP Einc LHCP !(1) where RCP = RRHCP,RHCP RRHCP,LHCP RLHCP,RHCP RLHCP,LHCP !.(2) When writing (1), we assume that the plane wave impinging on the periodic structure does not generate grating lobes after reflection, which requires that the condition max(a,b)< λ0/(1 +sin θinc) (λ0=2π/ω √µ0ε0) is fulfilled, aand bbeing the periods of the unit cell in the xand ydirections respectively (see Fig. 1(b)). A. Spectral domain MoM formulation In order to obtain the matrix RCP of (2) for the multilayered periodic structure of Figs. 1(a) and (b), we need to obtain the scattered electric field for impinging plane waves with both RHCP and LHCP polarizations. These scattered electric fields can be obtained in terms of current density excited at the metallized interface z=0, J(x,y), by the impinging waves. Let Ems(x,y,z) be the electric field generated in all space by a plane wave impinging on the multilayered substrate of Figs. 1(a) and (b) in the absence of the conducting split rings. The current density J(x,y) induced on the split rings will be the solution of the following electric field integral equation (EFIE) [20] ˆ z×"Ems(x,y,z=0) + +∞ X m=−∞ +∞ X n=−∞ZSmn GE(x−x0,y−y0,z=0,z0=0) ·J(x0,y0)dx0dy0#=0(3) (x,y)∈S00 where Smn (m,n=. . . , −1,0,1, . . .) is the metallized portion of the z=0 plane within the mn-th periodic unit cell, and GEis the non–periodic dyadic Green’s function of the multilayered substrate [21]. Since J(x,y) is a Floquet-periodic function of xand y, in order to solve the EFIE of (3), we only need to determine J(x,y) within one unit cell, e. g., the cell C00 covering the rectangular domain {0≤x≤a; 0 ≤y≤b}. For that purpose, we expand J(x,y) in C00 in terms of known basis functions Jlq(x,y) (l=1,...,L;q=1,...,Nb) as shown below J(x,y)= L X l=1 Nb X q=1 clqJlq(x,y) (4) where Jlq(x,y) (q=1,...,Nb) is the set of Nbbasis functions used to approximate the current density on the l-th arc of the C00 unit cell. When (4) is introduced in (3) and Galerkin’s version of MoM is applied, the following system of linear equations is obtained for the unknown coefficients clq [20] L X l=1 Nb X q=1 Γkl,pqclq =ekp (k=1,...,L;p=1,...,Nb).(5) If we invoke Parseval’s identity for 2-D Fourier transforms, the MoM matrix entries Γkl,pq of (5) can be expressed in the spectral domain as double infinite summations given by [22] Γkl,pq =ab +∞ X m=−∞ +∞ X n=−∞he Jd kp(kxm,kyn)∗it ·e G E,c (kx=kxm,ky=kyn,z=0,z0=0) ·e Jd lq(kxm,kyn) (6) where e G E,c (kx,ky,z=0,z0=0) is the continuous 2-D Fourier transform of GE(x,y,z=0,z0=0) [21], kxm = k0sin θinc cos ϕinc +2πm/a,kyn =k0sin θinc sin ϕinc +2πn/b, k0=2π/λ0and e Jd lq(kxm,kyn) is the discrete 2-D Fourier transform of Jlq(x,y), which is given by e Jd lq(kxm,kyn)=1 ab ZS00 Jlq(x,y) e−j(kxm x+kyny)dxdy (7) Finally, the coefficients ekp of the system of equations (5) can be obtained in the spectral domain as ekp =−ab he Jd kp(kx0,ky0)∗it·Ems(x,y,z=0) ×e−jk0(sin θinc cosϕinc x+sin θinc sinϕincy)(8) where the factor e−jk0(sin θinc cosϕinc x+sin θinc sinϕincy)has been explicitly included in (8) to absorb the dependence of Ems(x,y,z=0) on xand y. Equations (5) to (8) provide the spectral domain MoM formulation for the determination of the vector function J(x,y) of (3) and (4). The problem with this formulation is that the discrete 2-D Fourier transforms e Jd lq(kxm,kyn) of (6) (defined in (7)) cannot be obtained in closed form, even for the simplest choice of basis functions in (4) (i.e., a constant value for the two components of Jlq(x,y) on the surface of the split rings of Figs. 1(a) and (b)). In the following subsection, we will see how to tackle this problem. B. Basis functions and 2-D Fourier transforms Fig. 2 shows the l-th arc (l=1,...,L=4) of the unit cell C00 of Fig. 1(b). In Fig. 2 we have defined a shifted system of coordinates {x0,y0,z0}with origin at the center of the unit 4 b a y z x r 1 l z’ y’ x’ r 2l a/2 b/2 j 1l j 2l Fig. 2. Top view of the l-th arc in the left unit cell of Fig. 1(b). Radial and angular dimensions of the l-th arc are shown. cell in the plane z=0. We are going to introduce shifted polar coordinates ρ0and ϕ0, which are related with the original coordinates xand yas x=a 2+ρ0cos ϕ0(9) y=b 2+ρ0sin ϕ0(10) The unit vectors linked to the shifted polar coordinates are given by ˆρ0=cos ϕ0ˆ x+sin ϕ0ˆ y(11) ˆϕ0=−sin ϕ0ˆ x+cos ϕ0ˆ y(12) In the split rings used in CP reflectarray antennas, the condition wl=ρ2l−ρ1l(ρ1l+ρ2l)/2 (see Fig. 2) is usually fulfilled (i.e., the width of the rings is usually much smaller than their mean radius) [6]–[8]. Assuming this condition holds, we are going to neglect the radial component of the basis functions Jlq(x,y) of (4), Jlq(x,y)·ˆρ0, by comparison with their azimuthal component, Jlq(x,y)·ˆϕ0. Also, we will assume that the azimuthal component of Jlq(x,y) can be factored in terms of independent functions of ρ0and ϕ0, i.e., we will assume that the functions Jlq(x,y) can all be mathematically written as Jlq(x,y)=fl(ρ0)glq(ϕ0)ˆϕ0(13) (ρ1l<ρ0<ρ2l;ϕ1l<ϕ0<ϕ2l) As we will see in Section III, the assumption shown in (13) is justified by the fact that leads to values of the matrix RCP in (2) that match the values provided by the commercial software CSTr. Although the radial component, Jlq(x,y)·ˆρ0, of Jlq(x,y) has been neglected in (13), we would like to point out that this component can be easily accommodated in the mathematical formulation presented in the rest of this subsection. Since the functions glq(ϕ0)ˆϕ0of (13) are periodic functions of ϕ0with period 2π, these functions can be expanded as Fourier series of ϕ0as shown below glq(ϕ0)ˆϕ0= +∞ X i=−∞e gi lq ejiϕ0(ϕ1l<ϕ0<ϕ2l) (14) where e gi lq =1 2πZϕ2l ϕ1l glq(ϕ0)ˆϕ0e−jiϕ0dϕ0(15) Now, let us introduce the polar discrete spectral coordinates kρ,mn and kϕ,mn given by kρ,mn =q(kxm)2+(kyn)2(16) kϕ,mn =arctan kyn kxm !(17) If (9) to (14), (16) and (17) are introduced in (7), after some manipulations, it is possible to express e Jd lq(kxm,kyn) as e Jd lq(kxm,kyn)=2πe−jkxma+kynb 2 ab ×(e g0 lq ˜ f0 l(kρ=kρ,mn)+ +∞ X i=1 e−jiπ/2˜ fi l(kρ=kρ,mn) ejikϕ,mne gi lq +e−jikϕ,mn e gi lq∗)(18) where ˜ fi l(kρ) (i=0,1, . . .) are Hankel transforms of order iof the functions fl(ρ0) of (13), which can be expressed as ˜ fi l(kρ)=Zρ2l ρ1l Ji(kρρ0)ρ0fl(ρ0)dρ0,(19) and where the function Ji(·) of (19) is a Bessel functions of first kind and order i. Although the expression (18) for the computation of e Jd lq(kxm,kyn) looks simpler than (7), it presents two drawbacks. First, the functions ˜ fi l(kρ) cannot be obtained in closed form, even for the simplest choice of fl(ρ0) (e.g., a constant value in the interval ρ1l<ρ0<ρ2l). Second, (18) involves the determination of an infinite series, which has a deleterious impact on the computation of e Jd lq(kxm,kyn) by means of (18). In this paper, the functions chosen for Jlq(x,y) are entire domain basis functions which account for the singularities of J(x,y) at the edges of the arcs of Fig. 1(b). It is well known that these functions ensure a fast convergence of the spectral domain MoM with respect to the number of basis functions (LNbin (4)) [16], [19], [20], as will be demonstrated in Section III. The particular functions chosen for fl(ρ0) and glq(ϕ0) can be written as fl(ρ0)=1 r1−2 ρ2l−ρ1lhρ0−ρ2l+ρ1l 2i2 (20) (l=1,...,L) glq(ϕ0)=s1− 2 ϕ2l−ϕ1lϕ0−ϕ2l+ϕ1l 2!2 ×Uq−1 2 ϕ2l−ϕ1lϕ0−ϕ2l+ϕ1l 2!(21) (l=1,...,L;q=1,...,Nb) where Uq−1(·) is a Chebyshev polynomial of second kind and degree q−1. The integrals of (15) can be obtained in closed form for the functions glq(ϕ0) of (21), and their expression is given 5 in the Appendix A. However, the integrals of (19) cannot be obtained in closed form for the functions fl(ρ0) of (20). Fortunately, we have checked these latter integrals can be numerically obtained within very reasonable CPU times by means of Gauss-Chebyshev quadrature rules. In particular, when the functions fl(ρ0) of (20) are introduced in (19), the integrals can be rewritten as ˜ fi l(kρ)=Z+1 −1 hi l(kρ,t) √1−t2dt (22) where hi l(kρ,t)=ρ2l−ρ1l 2Jikρρ2l−ρ1l 2t+ρ2l+ρ1l 2 ×ρ2l−ρ1l 2t+ρ2l+ρ1l 2 (23) The integrals of (22) are amenable to be computed by means of Gauss-Chebyshev quadrature rules [23, Eqn. 25.4.38]. In particular, the integrals should be approximately obtained by means of the closed-form formula ˜ fi l(kρ)≈ Nqp X j=1 ηjhi l(kρ,t=tj) (24) where ηj=π/Nqp and tj=cos (2 j−1)π/2Nqp(j= 1,...,Nqp) ( j=1,...,Nqp). Since the functions hi l(kρ,t) are non-singular smooth functions of tin the interval −1≤t≤+1, a low number of quadrature points, Nqp, should suffice to obtain ˜ fi l(kρ) with a large accuracy. This has been confirmed by numerical simulations, as will be shown in Section III. Concerning the infinite series of (18), its convergence is strongly dependent on the value of kρ,mn. The larger the value of kρ,mn, the larger the number of terms that has to be retained in the series of (18) for an accurate determination of e Jd lq(kxm,kyn). Let us see how the convergence of this infinite series depends on kρ,mn. If we assume that a≈b(which holds in practical cases for most reflectarray antennas), by virtue of the mean value theorem for definite integrals, we can write that (see [23, Eqn. 9.3.1]) |˜ fi l(kρ=kρ,mn)||i| ∝Jikρ,mn ρ2l+ρ1l 2|i| ≈Ji kρ,mna 2!|i|∝ ekρ,mna 4i!i =eiln(ekρ,mna/4i)(25) which indicates that the series of (18) has an exponential convergence provided ekρ,mna/4i≤1, i. e., provided i≥ekρ,mna/4. Numerical simulations have shown that the infinite series of (18) has to be added in the interval 1 ≤ i≤ekρ,mna/4+8 for an accurate estimation of e Jd lq(kxm,kyn). Therefore, the larger kρ,mn, the larger the number of terms to be retained in the series. In practice, the computation of Γkl,pq requires to truncate the infinite summations of (6) within the intervals −Nmax ⩽ m,n⩽+Nmax. If we assume that a≈b, the maximum value of kρ,mn we will have to use in (6) after the truncation will be (see (16)) kmax ρ,mn =q(kmax xm )2+(kmax yn )2≈√22πNmax a(26) which means the maximum number of terms that will be necessary to retain in the infinite series of (18) will be imax =ekmax ρ,mna 4+8≈πeNmax √2 +8 (27) For the cases treated in the results section, we have used Nmax =50, which leads to imax ≈310 when using (27). For our purposes, this is the worst case scenario in the evaluation of (18), i.e., any evaluation of e Jd lq(kxm,kyn) within −Nmax ⩽ m,n⩽+Nmax will require to retain at most 310 terms in the computation of the infinite series of (18) when Nmax =50. Apart from the analytical stop criterium established above for the summation of the series of (18) (1 ≤i≤ekρ,mna/4+8), one additional fact that helps to reduce the CPU time required in the computation of this series is a well known property of Bessel functions, according to which Ji(x)≈0 if x<xth =Cthi when i>Nth. Numerical simulations show this property holds for Nth around 10 and Cth around 1. According to (19), this property of Bessel functions implies that ˜ fi l(kρ=kρ,mn)≈0 in (18) when i>Nth and (kρ,mna)/2<Cthi, which makes it unnecessary the evaluation of a large number of summation terms in the numerical computation of the series leading to e Jd lq(kxm,kyn) by means of (18). In Section III we will present a numerical method that makes it possible to easily estimate two new quantities dth(l) and bth(l) (which are a function of ρ1land ρ2l) for which ˜ fi l(kρ)≈0 if kρa<dth(l)(i−bth(l)) and i≥bth(l). According to this, it will only be necessary to compute ˜ fi l(kρ=kρ,mn) in the interval kmin ρ,mn(i)=dth(l)(i−bth(l)) a≤kρ≤kmax ρ,mn (28) when i≥bth(l). For i<bth(l), we will have to compute ˜ fi l(kρ= kρ,mn) in the interval kmin ρ,mn(i)=0≤kρ≤kmax ρ,mn. The integrals ˜ fi l(kρ) of (19) are a function of kρ,ρ1land ρ2lfor each value of i. In most applications of reflectarray antennas the distances ρ1land ρ2l(l=1,...,L) remain fixed from element to element and the angles ϕkl (k=1,...,S;l= 1,...,L) and αr(r=1,...,R) are varied to obtain the required phase shifts. So, a large amount of CPU time can be saved in the design of a large reflectarray antenna if the functions ˜ fi l(kρ) (l=1,...,L;i=0,...,imax) are interpolated as a function of kρfor fixed values of ρ1land ρ2l, and the interpolations are subsequently used when adjusting the angles ϕkl and αr in each element. The interpolations have to be performed in the interval kmin ρ,mn(i)≤kρ≤kmax ρ,mn for each value of l and i. In practice, we have divided the required interpolation intervals into Nint subintervals, and we have used Chebyshev polynomials up to seventh-degree in the interpolation in each subinterval. In Section III we will see that a maximum number of Nint =80 subintervals suffices to obtain interpolated values of ˜ fi l(kρ) within four significant figures, and that the number of required subintervals Nint decreases as iincreases, which is a consequence of the increase of kmin ρ,mn(i) (and therefore, of the reduction of the width of the original interpolation interval) as iincreases. III. Numerical results and validations Figs. 3(a) and (b) show the relative errors made in the evaluation of ˜ f0 l(kρ) and ˜ f100 l(kρ) by means of (24) as a function 6 (a) (b) Fig. 3. Magnitude of the relative errors in the computation of (a) ˜ f0 l(kρ) and (b) ˜ f100 l(kρ) by means of (24). The relative errors are plotted as a function of the number of quadrature points Nqp employed in (24). Parameters: a=b=5 mm; ρ1l=1.85 mm; ρ2l=2.05 mm. of kρ, and as a function of the number of Gauss-Chebyshev quadrature points Nqp employed in (24). The reference values used for ˜ f0 l(kρ) and ˜ f100 l(kρ) have been computed by means of the double exponential quadrature rule with 203 integration points (level of the quadrature rule M=5), since this is a numerical procedure which is particularly accurate when used in the integration of functions with integrable singularities at the end points of the integration interval [24]. Please note that a value of Nqp equal to 10 suffices to provide 4 significant figures in the values of ˜ f0 l(kρ) and ˜ f100 l(kρ) in nearly all cases, which indicates that (24) provides values that are accurate enough for ˜ fi l(kρ) within a very low CPU time consumption. In the rest of the results of this Section, we will always use Nqp =10 in (24). In Figs. 4(a) and (b) we plot the threshold values kth ρ(i) as a function of i, where kth ρ(i) are the values of kρfor which |˜ fi l(kρ)/˜ f0 l(kρ)|<10−4if kρ<kth ρ(i) and |˜ fi l(kρ)/˜ f0 l(kρ)|>10−4 if kρ>kth ρ(i). Therefore, for every value of i, the threshold value kth ρ(i) indicates the value of kρbelow which ˜ fi l(kρ) can be considered to be vanishing. Please note the plot of kth ρ(i) as a function of ifits very well to a straight line. So we have carried out a linear least squares fitting of kth ρ(i) as a function of igiven by kth ρ(i)a≈dth(l)(i−bth(l)) (29) which should be valid for i≥bth(l). For the least squares fitting, we have used Nsam samples of iin ascending order, ik(k=1,...,Nsam) the first sample being i1=20. In i 0 20 40 60 80 100 120 140 160 180 200 kρa 0 100 200 300 400 500 kρ th(i)a (a) i 0 20 40 60 80 100 120 140 kρa 0 100 200 300 400 500 kρ th(i)a (b) Fig. 4. The solid line kth ρ(i)astands for the threshold values of kρafor which ˜ fi l(kρ) is negligible if kρa<kth ρ(i)a, and ˜ fi l(kρ) is non-negligible if kρa>kth ρ(i)a.kth ρ(i)ais plotted as a function of i. Parameters: a=b=5 mm; ρ1l=1.85 mm and ρ2l=2.05 mm in (a); ρ1l=1.20 mm and ρ2l=1.40 mm in (b). the case of Fig. 4(a), it turns out that dth(l)≈2.343 and bth(l)≈9.484 with a coefficient of determination r2=0.9999. In the case of Fig. 4(b), it turns out that dth(l)≈3.388 and bth(l)≈8.371 with a coefficient of determination r2=0.9998. These results indicate that the values of dth(l) and bth(l) are slightly dependent on the values of the inner and outer radii of the rings, ρ1land ρ2l, which are the endpoints of the integration interval of (19). Since the fitting of kth ρ(i) to a straight line for i≥bth(l) is very good, in practice a small value of the number of samples of kth ρ(i), Nsam, suffices for the determination of dth(l) and bth(l). Once these values of dth(l) and bth(l) are known, one can use (28) to estimate the interval of values of kρin which the functions ˜ fi l(kρ) present non-negligible values with a view to computing e Jd lq(kxm,kyn) via (18). Figs. 5(a) and (b) show the errors made in the interpolation of ˜ f0 l(kρ) and ˜ f100 l(kρ) in the interval kmin ρ,mn(i)≤kρ≤kmax ρ,mn when this interval is divided into Nint subintervals of equal width, and Chebyshev polynomials up to seventh-degree are subsequently used to carry out the interpolation in each subinterval. Please note that whereas 80 subintervals are required for the determination of ˜ f0 l(kρ) with an accuracy of 4 significant figures in the whole interval kmin ρ,mn(i)≤kρ≤kmax ρ,mn, only 40 subintervals are required for the determination of ˜ f100 l(kρ) with the same accuracy. The explanation for this is that whereas the whole interpolation interval for ˜ f0 l(kρ) is roughly 0≤kρa.450, the interpolation interval for ˜ f100 l(kρ) is much smaller, 200 .kρa.450, which is a consequence of the fact that ˜ f100 l(kρ) is negligible for kρa.200 as shown in Fig.4(a). Fig. 6(b) shows the phase curves required for the design 7 (a) (b) Fig. 5. Magnitude of the relative errors in the interpolation of (a) ˜ f0 l(kρ) and (b) ˜ f100 l(kρ) as a function of the number of subintervals Nint used in the interpolation. Parameters: a=b=5 mm; ρ1l=1.85 mm; ρ2l=2.05 mm. of a CP reflectarray antenna (see (2)) made of single split rings (see Fig.6(a)) at a frequency of 19.95 GHz. By virtue of the VRT of [2], a linear phase variation is obtained for both ∠RRHCP,RHCP and ∠RLHCP,LHCP as a function of α2, which is the rotation angle of the mirror symmetry axis of the split ring. In order to keep the sense of CP after reflection in a reflectarray antenna, we need a phase difference of 180◦ between the reflection coefficient of the component of the electric field along the mirror symmetry axis of the split ring and the reflection coefficient of its orthogonal component [2], [3]. To ensure this phase difference is maintained, the length of the arcs in the reflectarray element has to be slightly adjusted as α2is varied [6]. This adjustment of the length of the arcs as a function of α2is shown in Fig. 6(c). In Fig. 6(b) we show the convergence of ∠RRHCP,RHCP and ∠RLHCP,LHCP as a function of the number of basis functions Nbused in the two arcs of the split rings. It can be seen that just three basis functions per arc suffice to achieve convergence (in fact, the results obtained with the three basis functions corresponding to q=1,2,3 in (21) overlap those obtained when Nb=7 and 1≤q≤7), which means that only 6×6 MoM matrices have to be inverted when analyzing the multilayered periodic structure containing single split rings in the unit cell. In order to validate the results obtained with the spectral MoM software described in Section II, these results are compared with results provided by the commercial software CSTr. Excellent agreement is found between both sets of results. In Fig. 7(b) we plot the phase curves for a dual-frequency CP reflectarray antenna made of concentric split rings (see b a y z x y’ x’ r 23 j 23 j 13 a 2 r 13 (a) α2 (deg) 20 40 60 80 100 120 140 160 180 RRHCP,RHCP and RLHCP,LHCP (deg) -300 -200 -100 0 100 200 300 RLHCP,LHCP q=1,2 RLHCP,LHCP q=1,2,3 RLHCP,LHCP q=1,...,7 RLHCP,LHCP CST RRHCP,RHCP q=1,2 RRHCP,RHCP q=1,2,3 RRHCP,RHCP q=1,...,7 RRHCP,RHCP CST (b) α2 (deg) 20 40 60 80 100 120 140 160 180 ϕ23 −ϕ13(deg) 147 147.5 148 148.5 149 149.5 150 150.5 (c) Fig. 6. (a) Single split ring. (b) Phase of RRHCP,RHCP and RLHCP,LHCP for the reflectarray element of (a) as a function of the rotation angle of the ring α2. Convergence is studied with respect to the number of basis functions (q=1,...,7) used in the two arcs of the split rings. Our results obtained with the spectral domain MoM are compared with results obtained with CSTr. (c) Variation of the length of the arcs of the split ring with the rotation angle α2to ensure a phase difference of 180◦in the reflection coefficients of the two orthogonal linear components of the CP impinging waves. Parameters: a=b=5 mm; ρ13=1.85 mm; ρ23=2.05 mm; h1=0.787 mm, εr,1=2.2, tan δ1=0.0009; h2=0 mm; θinc =30◦,ϕinc=0◦;f=19.95 GHz. Fig. 7(a)). While the outer rings are intended to adjust the phase at 19.95 GHz, the inner rings are intended to adjust the phase at 29.75 GHz as in [6]. In Fig.7(b) we plot the curves for phase adjustment at 29.75 GHz. Since this phase adjustment is controlled by the inner rings, in the application of the VRT ∠RRHCP,RHCP and ∠RLHCP,LHCP are varied by rotating the mirror symmetry angle of the inner split ring, α1, while keeping fixed the mirror symmetry axis of the outer ring α2. Note that whereas ∠RRHCP,RHCP and ∠RLHCP,LHCP are linear functions of α2in Fig. 6(b), ∠RRHCP,RHCP and ∠RLHCP,LHCP are not any longer linear functions of α1in Fig. 7(b). This is due to the fact that whereas in Fig.6(b) the whole reflectarray element is 8 b a y z x r 13 y’ x’ j 13 j 21 j 11 r 21 r 11 a1 a 2 j 23 r 23 (a) α1 (deg) 20 40 60 80 100 120 140 160 180 RRHCP,RHCP and RLHCP,LHCP (deg) -300 -200 -100 0 100 200 RLHCP,LHCP q=1,2 RLHCP,LHCP q=1,2,3 RLHCP,LHCP q=1,...,7 RLHCP,LHCP CST RRHCP,RHCP q=1,2 RRHCP,RHCP q=1,2,3 RRHCP,RHCP q=1,...,7 RRHCP,RHCP CST (b) α1 (deg) 20 40 60 80 100 120 140 160 180 ϕ21 −ϕ11(deg) 145 150 155 160 165 (c) Fig. 7. (a) Two concentric split rings. (b) Phase of RRHCP,RHCP and RLHCP,LHCP for the reflectarray element of (a) as a function of the rotation angle of the inner ring α1. Convergence is studied with respect to the number of basis functions (q=1,...,7) used in the four arcs of the two split rings. Our results obtained with the spectral domain MoM are compared with results obtained with CSTr. (c) Variation of the length of the arcs of the inner ring with the rotation angle α1to ensure a phase difference of 180◦in the reflection coefficients of the two orthogonal linear components of the CP impinging waves. Parameters: a=b=5 mm; ρ13=1.85 mm; ρ23=2.05 mm; ϕ13=−75.2◦mm; ϕ23=75.2◦mm; α2=0◦;ρ11=1.20 mm; ρ21=1.40 mm; h1=0.787 mm, εr,1=2.2, tan δ1=0.0009; h2=0 mm; θinc =30◦,ϕinc=0◦; f=29.75 GHz. rotated for phase adjustment, in Fig. 7(b) only the inner part of the element is rotated for phase adjustment, and therefore the proportionality between phase and rotation angle proven in [2] is lost. Fortunately, the rotation of α1in Fig. 7(b) provides enough phase range for reflectarray design. In Fig.7(c) we plot the adjustments in the length of the inner arcs that are needed to keep a phase difference of 180◦between the reflection coefficients of the two orthogonal components of the electric field of the impinging CP waves. These adjustments are larger than in the case of Fig. 7(b), which is attributed to the fact TABLE I CPU times (seconds)required for the generation of 31 points in the curves of Figs.6(b)and 7(b). The CPU times are for our MoMsoftware and CSTr. Figure TMoM CPU TCST CPU Fig. 6(b) 89.8 2383 Fig. 7(b) 152.9 2486 TABLE II CPU times (seconds)contributions (computation of Hankel transforms, interpolation of Hankel transforms,calculation of 2-D Fourier transforms of the basis functions and evaluation of matrix entries)to the total CPU time required by our MoMsoftware for the generation of one point in Figs.6(b)and 7(b). Figure THT TINT T2DFT TME TMoM TOTAL Fig. 6(b) 0.398 0.0005 2.868 0.007 3.290 Fig. 7(b) 0.570 0.0008 4.871 0.020 5.477 that the outer arcs of the element are not modified during the process of phase adjustment. As in Fig. 6(b), in Fig. 7(b) convergence of the spectral domain MoM software is achieved with just three basis functions per arc (which requires the inversion of 12×12 MoM matrices), and excellent agreement is found between the spectral domain MoM software and CSTr. We have computed the CPU time required by CSTrto generate 31 points of the curves of Figs. 6(b) and 7(b), TCST CPU, and the CPU time required by the spectral domain MoM software when 3 basis functions per arc are employed, TMoM CPU . When running the spectral MoM software, we have used the interpolated versions of ˜ fi l(kρ), and we have included in TMoM CPU the CPU time required to compute kth ρ(i) and interpolate the functions ˜ fi l(kρ). The results obtained for TCST CPU and TMoM CPU are shown in Table I. These CPU times have been obtained in a laptop computer with processor Intel Core i7-6700HQ at 2.6 GHz with four cores and 32 GB of RAM memory. The MoM code has been written in FORTRAN language. According to the results of Table I, the spectral MoM software is typically 27 times faster than CSTrin the analysis of the periodic structures studied in Fig. 6(b), and around 16 times faster than CSTrin the analysis of the periodic structures of Fig. 7(b). When the spectral domain MoM described in Section II is applied to the analysis of one single periodic structure containing split rings, there are different numerical steps having different CPU time contributions to the total CPU time required by MoM, TMoM TOTAL. These CPU time contributions include the CPU time required to compute the Hankel transforms ˜ fi l(kρ) that are needed for the interpolations by means of (24), THT, the CPU time required for the interpolations of the Hankel transforms by means of Chebyshev polynomials, TINT, the CPU time required for the evaluation of the 2-D Fourier transforms of the basis functions by means of (18), T2DFT, and the CPU time required for the determination of the MoM 9 matrix entries by means of (6), TME. Table II shows these different CPU time contributions to TMoM TOTAL (the CPU times have been obtained with the same computer used for Table I). In particular, the CPU times presented in Table II are for the periodic structure analyzed in Fig. 6(b) when α2=92.9◦ (ϕ23 −ϕ13=150.4◦in Fig. 6(c)), and for the periodic structure analyzed in Fig. 7(b) when α1=92.9◦(ϕ21 −ϕ11=150.0◦ in Fig. 7(c)). Note that the most important contribution to the total MoM CPU time is T2DFT, which is around 88% of TMoM TOTAL. This justifies the efforts carried out in Section II to compute the Hankel transforms ˜ fi l(kρ) appearing in (18) in an efficient way. If we had not optimized the computation of ˜ fi l(kρ), the method presented in this paper would not have been competitive by comparison with commercial software. Once the functions e Jd lq(kxm,kyn) of (18) have all been computed, the CPU times required to compute the double summations of (6) and to solve the system of equations of (5) are negligible (note that the maximum size of the MoM matrices to be inverted is 12 ×12). The next important contribution to the total MoM CPU time is THT, which represents around 11% of TMoM TOTAL. Table II shows the contribution of the time required for the interpolation of the Hankel transforms to TMoM TOTAL is negligible. At this point, we should remember that the computation and interpolation of the Hankel tranforms has to be carried out only once in the design of a reflectarray antenna where the inner and outer radii of the rings ρ1land ρ2l(l=1,...,L) remain constant in all the elements of the antenna. Therefore, we can neglect the contribution of THT to the CPU time required in the design of a whole reflectarray antenna. In our study of the computational performance of the spectral domain MoM software as a function of the different input variables involved, we have finally analyzed the computational complexity of TMoM TOTAL as a function of the number of quadrature points used in (24), Nqp, as a function of the number of terms retained in infinite series of (18), imax, and as a function of the number of basis functions used in the application of MoM, Nb. We have found that the effect of Nqp on TMoM TOTAL is irrelevant for 3≤Nqp ≤15, that the spectral MoM is roughly of complexity O((imax)2) for 48 ≤imax ≤200, and that the spectral MoM is roughly of complexity O(Nb) for 1 ≤Nb≤7. Based on the curves of Figs. 6(b), 6(c), 7(b), and 7(c), a dual-frequency pencil beam CP reflectarray antenna has been designed. The antenna radiates LHCP at 19.95 GHz and RHCP at 29.75 GHz as in the case of the antenna designed in [6]. The designed antenna is circular and consists of 5024 elements arranged in a 80 ×80 grid with cell size 5 mm ×5 mm (the diameter of the antenna is 400 mm). The antenna is intended to produce a focused beam in the direction θb=30◦and ϕb=0◦ (see [25, Fig. 5.a]). It is illuminated by a corrugated circular feed-horn with its phase center located at the coordinates x=−150 mm, y=0 mm, z=259.8 mm with respect to a coordinate system with origin at the center of the reflectarray (see [25, Fig. 5.a]). The horn is assumed to radiate LHCP waves at 19.95 GHz and RHCP waves at 29.75 GHz. The radiation pattern of the horn is modelled as a function cos7(θ), which provides an illumination level at the reflectarray edges 12 dB below the maximum. The antenna elements are dual concentric split rings as those shown in Fig.7(a). In the design Azimuth angle (deg) -50 -40 -30 -20 -10 0 10 20 30 40 50 Gain (dB) -20 -10 0 10 20 30 40 Copolar LHCP q=1,2,3 Copolar LHCP CST Copolar LHCP HFSS Crosspolar LHCP q=1,2,3 Crosspolar LHCP CST Crosspolar LHCP HFSS (a) Elevation angle (deg) -50 -40 -30 -20 -10 0 10 20 30 40 50 Gain (dB) -20 -10 0 10 20 30 40 Copolar LHCP q=1,2,3 Copolar LHCP CST Copolar LHCP HFSS Crosspolar LHCP q=1,2,3 Crosspolar LHCP CST Crosspolar LHCP HFSS (b) Fig. 8. Radiation patterns in (a) the azimuth plane and (b) the elevation plane for a dual-frequency CP reflectarray based on an element with two concentric split rings. Plots are presented for the analysis of the reflectarray under the local periodicity assumption with both our spectral domain MoM software (solid and dash-dotted lines), CSTr(dashed and dotted lines) and HFSSr(× and +). The results plotted are for LHCP radiation at 19.95 GHz. of the antenna at 19.95 GHz, the dimensions and orientations of the outer rings have been adjusted in accordance with Figs. 6(b) and 6(c), while assuming the inner rings were not present. However, when the antenna has been designed at 29.75 GHz, both the outer and the inner rings have been taken into account. The dimensions and orientations of the outer rings have been fixed at the values obtained for the design at 19.95 GHz, and the dimensions and orientations of the inner rings have been adjusted in accordance with Figs. 7(b) and 7(c). The design of the antenna has been carried out by means of the spectral domain MoM software described in Section II under the local periodicity condition, while using the exact incidence angle in each reflectarray element. Once the reflectarray antenna has been designed, it has been analyzed with the spectral domain MoM software (using three basis functions per arc with q=1,2,3 in (21)), and with the two commercial softwares CSTrand HFSSr, under the local periodicity assumption. In the analysis, the exact incidence angle has been considered in each reflectarray element. Figs.8 and 9 show the radiation patterns obtained in the azimuth and elevation planes at 19.95 GHz and 29.75 GHz respectively. The resulting gain is 35.88 dBi at 19.95 GHz and 39.42 dBi at 29.75 GHz. And the antenna efficiency is of 56% at 19.95 GHz and 57% at 29.75 GHz. In Fig. 8 the agreement between our MoM results, CSTrand HFSSris excellent. In Fig. 9 the agreement is excellent for the copolar component of the