035067-1-T HYBRID FINITE ELEMENT AND MOMENT METHOD SOFTWARE FOR THE SERAT ARRAY 1st Quarterly Report Sanders, a Lockheed Martin Co. P.O. Box 868 Nashua, NY 03061-0868 Naval Air Warfare Center Mission & Sensor Systems Division RF Sensors Branch Patuxent River, MI 20670 January 1997 35067-1-T = RL-2482

PROJECT INFORMATION PROJECT TITLE: REPORT TITLE: U-M REPORT No: CONTRACT START DATE: END DATE: DATE: Hybrid Finite Element Design Codes for the SERAT Array 1ST Quarterly Report Hybrid FE and MM Software for the SERAT Array #035067-1-T October 1996 SEPTEMBER 1998 January 9, 1997 (1ST Quarterly Report) Roland Gilbert SANDERS, INC, A Lockheed Martin Co. MER 24-1583 P.O. Box 868 Nashua, NH 03061-0868 Phone: (603) 885-5861 Email: RGILBERT@mailgw.sanders.lockheed.com SPONSOR: SPONSOR CONTRACT No. UM-PRINCIPAL INVESTIGATOR: John L. Volakis EECS Dept. University of Michigan 1301 Beal Ave. Ann Arbor, MI 48109-2122 Phone: (313) 764-0500 FAX: (313) 747-2106 volakis@umich.edu http://www-personal.engin.umich.edu/ volakis/ J. Gong (UM), D. Jackson (UH), J. Volakis (UM) and D. Wilton (UH) CONTRIBUTORS TO THE REPORT:

Contents 1 INTRODUCTION 2 2 MEETINGS 2 3 HIGHLIGHTS OF THE 1ST QUARTERLY ACTIVITY 2 3.1 Code Development and Integration Plan......................... 2 3.2 FSS and Simple Moment Method Code Development..................... 4 3.3 Small Array hybrid FEM code......................... 8 A ANTICIPATED TASKS AND SCHEDULE 11 B ATTRIBUTES OF PRISMATIC ELEMENTS 14 C POST PROCESSING PACKAGES 15 D Field Representations for Multi-Layered Periodic Media 18 E Extraction of Asymptotic Forms of Series 25 F Code Validation 32 G Acceleration of the Free-space Periodic Green's Function by the Method of Singh 33 H Acceleration of the Free-Space Periodic Green's Function by the Ewald Method 35 I Results Comparing Acceleration Methods 41 J Two-Dimensional Fourier Transform Properties for Periodic Structure Applications 46

PROJECT CHRONOLOGY April 1996 July 1996 August 1996 Sept. 20 1996 October 1996 October 1996 Nov. 15 1996 Proposal Submission Answers to Proposal Questions Contract Negotiations Kickoff meeting at Ann Arbor (attended by Sanders, UM and UH) Contract Signed between U-M and Sanders in Mid October Subcontract to the Univ of Houston (formalized in early November) SERAT Review meeting (at Nashua) 1

1 INTRODUCTION The goal of the SERAT project at the University of Michigan (with subcontract to Univ. Of Houston) is to develop a suite of software for the analysis of strip and slot dipoles on multilayered substrates backed by a frequency selective surface. The dipoles are equipped with photonic switches permitting variable electrical dipole lengths for broadband performance and the FSS is suitably designed to simulate a variable substrate thickness for optimal operation. A general view of the geometry is given in Figure 1. The UM/UH team proposed to construct a code which combines various computational modules interfaced with appropriate pre-processors and post-processors. The computational modules include: Stand-alone moment method simulation of the FSS with up to 10 layers with commensurate and non-commensurate periodicities. Simple moment method simulation of the antenna elements on the FSS panels Hybrid FEM simulation modules for small arrays, planar periodic arrays and curved arrays on FSS panels. Various options for modeling the FSS and for mesh truncation were proposed to provide a compromise between speed and accuracy. These are outlined in the proposal and summarized in the attached milestone chart (repeated from the proposal) 2 MEETINGS As noted above, the UM/UH SERAT activity began earnestly in October 1996 when the UM and Sanders contract was formally signed. Therefore, this is our first quarterly report covering the period from 1 October to 31 December 1996. The UH subcontract was initiated immediately afterwards and was formally in place by early November. The kickoff meeting took place on September 20, 1996 and was attended by Roland Gilbert (Sanders), Henry Karwacki (Sanders), J. Volakis (UM), J. Gong (UM), S. Bindiganavale (UM), D. Wilton (UH) and D. Jackson (UH). At the meeting, Dr. Gilbert presented an overview of the SERAT program and UM/UH presented the code development plan and related R&D effort at some detail. Follow-up discussions between the Sanders and UM/UH centered on the various parts of the codes, their feasibility and the capabilities of the Green's function for non-commensurate periodicities. UM gave a follow-up review of the SERAT analysis software at the Nov 15, 1996 review meeting of the SERAT program held at Nashua. 3 HIGHLIGHTS OF THE 1ST QUARTERLY ACTIVITY 3.1 Code Development and Integration Plan Before beginning with the implementation of the tasks outlined in the proposal (attached in section 4 of this report) and summarized in the milestone chart, it was necessary to designate a careful code development plan. This was done by Dr. Gong who is in charge of the overall code development plan. An overview of the proposed code integration plan is depicted in 2

SERAT panel:: singti s.' " tri' crossed o. it.i ifre nt slots FI am eletes FSS elements antenna panels Curved SERAT Figure 1: An Example of SERAT Configuration 3

Figure 2. Basically, the analysis engine (FEM, hybrid and moment method modules) will require a single geometrical interface, and the pre-processors will be used to translate and/or reformat data as needed for the various computational modules. The post-processors will include the necessary data interfaces for graphical display and parameter extraction. We investigated a number of Unix platform plotting and visualization packages in October and we selected two freeware packages. One of the them is the xmgr plotting package with ease to use graphical interface allowing for plotting and data display with X windows. This package reads in spreadsheet output datafile which can be also ported on personal computer (PC and Macs) for plotting and display. A second package was also investigated for less interactive capability but easy to integrate with our developed code and to be called directly from inside the code. The output will be high quality postscript format which can be displayed and printed on any platform. In addition to the above, a rather major effort during this quarter was devoted to the investigation of fast solver for sparse FEM system. Specifically, we investigated the convergence of new BCG, QMR and GMRES algorithms and concluded that BCG is still the most attractive of all these algorithms. However, a new fast solver was recently introduced for inverting sparse systems. This solver, referred to as CVSS, was developed by NASA funding and is available to U.S. companies, agencies and institutions. CVSS is an LU type solver specialized to sparse systems but because of its efficient memory use, speed and limited need for interface with user, it was shown to be competitive with the iterative methods. At this stage, we will continue the investigation of these solvers and we will report on them after their use in a specific SERAT analysis module. In accordance with the schedule shown in the milestone chart, two major activities were carried out in parallel: 3.2 FSS and Simple Moment Method Code Development This subproject is carried out at the Univ of Houston(UH) by Profs. Wilton and Jackson. As stated in the proposal, these moment method codes will serve for the evaluation of the multilayered FSS and simple antenna(strip and slots) element analysis in the presence of the FSS. A first version of the code(s) will be completed by the end of the 2nd quarter for commensurate periods with the more general code versions for non-commensurate periods to be available at the end of the first year. UH is on schedule for the delivery of these codes. At the start of the contract, UH had two codes available for target scattering and antenna radiation in the presence of multilayered media (See Figure 3). One of them was capable of modeling periodic FSS-type structures but was not equipped with the ease to use I/O interfaces. Moreover, the employed periodic Green's function was not optimized for fast convergence. The second code (referred to as EIGER) is equipped with a much improved I/O interface and its implementation is based on the more attractive and rigorous mixed potential integral equation(MPIE) formulation. However, the code is not specialized to periodic structures. The latter implies a necessary replacement of the Green's function and since neither code employed a fast convergent Green's function, UH decided to employ the EIGER code as the basis for the development of the new moment method SERAT code. 4

Code Development Plan Planar or curved Planar periodic FSS only single element and od code small finite array a tcoe code option code option p Geometry t0 Hardwired 7 - User C meshes/elements selections Cd Free choices C)1 ff EM solution 5 Output a t ----- I ----- I f

Milestone Chart for EM Model Development (Tasks 1 and 2) Quarterly Progress Task FSS GreenCcsfnctin io: de:: (U of Houston) Mesh GeneratorforAntenna Elements Mesh generator forFSS elements v Single Element and Small Array v vv v! -!!P!!. ~ -. PlaIar and -Curvi IBC... ''. I .. 1. I.. I -.111. I- I. - -- I.-., Planar-FEM/Moment Method * s, Curved-FEM for antenna and FSS - II Planar Periodic A rray - FEM with IBCs 2nd 3rd 4th 5th 6th 7th Q. Q. Q. Q. Q. Q. ------- -~ =---— lo. ---- - --- -- --- ~ — - c Simple Moment Method code.IL FEM and Moment Method for FSS FEM for antenna and FSS Curved Array Cylindrical. Approximate Doubly Curved A. 1 f Doubly Curved with fast integral algorithms for mesh truncations I ll I Software Integration and I/O Displays Tr 1,. Validation Software Support I A A 6

M+1 'g M+1 Z = bM+l I EM' IM = bM Z-{'M - r -- -- ~ — r ~- -- 8M4' 11M-1 z b-l cl I 0 0. * - b3 Z- 2 e1,, z =b -'h * ---- -- — lr -- -~ — 0o 4o Figure 3: Obstacle or antenna radiator in the presence of multilayered medium 7

In addition to upgrading EIGER, during the past three months, UH emphasized the development and implementation of acceleration techniques for a fast implementation of the periodic multilayered Green's function. This is the crux of the moment method code implementation and its efficiency is therefore essential in delivering a useful stand-alone and hybrid FEM code. Instead of using the Singh et.al. approach, UH employed an acceleration technique introduce by Ewald. This technique replaces the periodic (Floquet expansion) Green's function (a surn of spectral integrals) by sums of sampled values of the spectra on a reciprocal lattice. Basically, the summations of the Floquet expansion are replaced by a fast converging semi-definite integrals which can be cast in terms of complimentary error functions. Using these replacements, the resulting Floquet sums are rapidly converging. Initial results based on the free standing periodic structure show that speed-ups on the order of 50 to 100 can be achieved using the Ewald approach. Actual series convergence curves are given in Figures 4 and 5 for the two sums (spectral and spatial) of the Floquet representation. It is seen that the new series converges using 2 to 3 terms, whereas the older approach (based on Singh's method) requires 50 to 200 terms to achieve convergence using the same error criterion. A short description of UH's acceleration technique and MPIE formulation is given in the attached section to this quarterly report. More details on the technical aspects of the implementation will be given after code completion. In the next quarter, UH plans to implement the new acceleration technique for the multilayered periodic media and to complete the moment method code for commensurate periods. The code will then be delivered to UM for incorporation into the FEM hybrid code and for stand-alone analysis. 3.3 Small Array hybrid FEM code UM began the development of this module in early November 1996. According to our schedule, the plan is to complete a first version of the small array and planar periodic codes by the end of 4th quarter. Over the past two months, we developed the formulation of the FEM code using second order prismatic elements. Instead of proceeding to modify one of our existing FEM codes, we opted to develop a totally new code to accommodate the planned mesh truncation options, IBC simulations of the FSS and specialization to the FSS configurations. Developing a new code from scratch is certainly a more time consuming task but will result in a better overall code. Moreover, this approach will avoid potential difficulties typically associated with older codes written for different applications. A key feature of the code under development is the use of higher order prismatic elements which will permit a straightforward application of the code to doubly curved arrays and FSS panels. We anticipate that the testing of the code will begin toward the end of January. At that time we will also test the use of higher order IBCs for modeling FSS panels. Attached to this quarterly report in the Appendices are the Code Development Flow Chart, the feature descriptions of the chosen FEM elements, and the post processing package. 8

0.........................................................................................Z................... MEN* <LU o~go cI~cI0 40U elect I I O( U) N 0 0 OC 0 z 0 0 So 0 'm - T,r- T M z L0 D I,.. cO 00 0 0 0 0 0 o ~ o o O O o Vo o o o W o 0 LU W- ITl L (LN33g3d) 4 0'Sd3 Figure 4: Singh's method 9

E=Eop=3.54 CTI 0t p 0 w 0 Lu 100 10 I 0.81 001o 0.0001 0.0001 I1E-005 I Ew006 I E-007 I E-008 I E-00O9 I E-OlO 0 1 2 3 4 5 N

A ANTICIPATED TASKS AND SCHEDULE The following anticipated 'ready-for-delivery' time is referred with respect to the project starting time from October 1996. Single Element or Small Array Code: * Scheduled Ready-for-Delivery Time: sixth months * Purpose: to model elementary parts of the SERAT panel and acts as an FEM kernel * Features: - Modular (FE-BI, FE-PML, FE-IBC, etc.); - Capable for fast solutions; - Capable of modeling antenna/FSS elementary coupling; - Simplified feed models; - Rendering reconfigurable design of SERAT platform (thickness) * Challenging Efforts: - Development of a high-order p-type element FEM code - Fast solution scheme (robust iterative/direct solvers) - Mesh automation/semi-automation (triangularization) - Speed-up PML convergence; improve BI efficiency (whenever used) - User-friendly environment design Periodic Array Code: * Scheduled Ready-for-Delivery Time: one year * Purpose: to model the periodic infinite SERAT panel * Features: - Able to model infinitely large SERAT panel (in contrast to MoM) - Capable of incorporating inhomogeneity, feedlines, etc. - Fast and efficient - Appropriate for reconfigurable design * Challenging Efforts: - Period FEM development (periodic BC's) - Periodic BI development - Periodic IBC development - Combination with BI/IBC 11

- Feed model incorporation - Fast solver testing - User interface design - Reconfigurable design ability Full Hybrid FEM Code * Scheduled Ready-for-Delivery Time: 16-18th months * Purpose: to accurately model the entire non-commensurate SERAT * Features: - Considers the entire couplings of FSS-FSS, antenna-antenna, FSS-antenna, feedline, etc. - Considers the finite sized platform (edge effects, etc.) - Considers the curved platform * Challenging Efforts: - Mesh generation - User-friendly interface design - Curved platform incorporation - Speed-up solution time as much as possible!!! Package Integration and Validation * Universal file format design/modification * Modular environment design and integration * Software validation * Software documentation 12

Code Development Flow Chart I r ' ^ V ----- -— I- -- ------- 9 1 ^ '*^ S Si-.; ---------------..........:b.Pi I 1^^ I i^.Rt s ^ ^ ^ ^ ^.................. ^;M M M ~m Wam. U$ xt I"$M X ~M W W %" ~ M Am W*0 " *W 'Ma so*M M Wm*W M - 4* OM k Ak a 1^-^-r. — " —j bf A Code tngine 1 I 1 3 FE Methods - - Moment Methods 3; - Hybrid Methods I.0 i n i y t9 T ara r I ^ ^ ^ ^ ^ ^ ^, ^^^^ 8.. _,.........................'i.... "lk w w" *w = Xw s0. 1%..w 00 w Wm. Wow

B ATTRIBUTES OF PRISMATIC ELEMENTS * Suited for multi-layered structures * With minimum number of unknowns for large thicknessslot-ratio * Good system condition * More accurate with higher order version of basis functions * Meshing needed only on a surface * Extendable to doubly-curved structures 14

C POST PROCESSING PACKAGES A FEW ATTRIBUTES OF XMGR PLOTTING PACKAGE * User defined scaling, tick marks, labels, symbols, line styles, colors. * Batch mode for unattended plotting. * Read and write parameters used during a session. * Polynomial regression, splines, running averages, DFT/FFT, cross/auto-correlation. * Hardcopy support for PostScript, HP-GL, and FrameMaker.mif format. Figure 6 is the sample layout by the XMGR package. ABOUT PGPLOT PLOTTING PACKAGE PGPLOT is a Fortran subroutine package for drawing scientific graphs on various graphics display devices. It is intended for the Fortran programmer who wishes to write a program generating graphical output. For most applications, the program can be device-independent, and the output can be directed to the appropriate device at run time. Figure 7 shows the graphic sample plots using the PGPLOT package. 15

:i I} C) b0 * -+-D *PI... C14 bJO 0 6 - I I 16

PGPLDr Exomple IS; PMOVE and PGOFlW 0 60 120 1 B0 I I I I I i2 I.1I i I i 2 i -KI 0 - 1 -2 -( *.*/ \ * \ / I I I I 0 -1 -2 90 160 270 360 450 54-0 630 Fig ure 7: Sample Graphic Plots Using PGPLOT Package

D Field Representations for Multi-Layered Periodic Media The electric field due to equivalent electric and magnetic currents is written in the mixed potential integral equation (MPIE) representation as Es(r) = -jwA-V - 1V x F, (1) with the corresponding magnetic field expressed as 1 Hs(r) = -jwF- VI + -V x A. (2) Al The magnetic vector potential is defined as A(r) = A(r, r'). J(r') dS' (3) while the electric vector potential is F(r)= sF(r, r') M(r') dS'. (4) The electric scalar potential is >(r) = V -A(r) = jf V. gA(r, r') J(r) dS' = Vs J(r)KD(r rT) dS'+ z J(r')P,(r,r') dS', (5) wher the last equality follows from using the relation j: 2 GA(r,r') = -V'K (r,r')+ zPz(r,r'), (6) 18

followed by an integration by parts. Michalski shows that such a splitting of the scalar potential is always possible and permits the order of the singularity in the scalar potential kernel to be reduced by transferring a derivative from the kernel onto the current. The corresponding magnetic scalar potential is @(r) = -WV. F(r) = fv vgF (r, r).M(r') dS' = f V'. (r') I (r, r') dS' + j M(r')QQzr, r') dS', (7) which can be written as k2V. F(, r')= -V' (r, r') + ZQ'(r, r'). (8) The vector potential dyadic Green's function gA is not unique, but the most convenient form for numerical work has the form GA(, r) o0 g(TAr ) = 0 GAy(r,r') O GA(r,r') G(r, r') Gz(r, r') CO 00 = E E (z)z-jkpq. (p-p') (9) q(A z') (9) p=-oo q=-oo where Aq constitutes a set of sampled values of the spectral representation of the magnetic vector potential for a point source in a layered medium. The sampled values are on the reciprocal lattice of the wavenumber domain and correspond to the grating lobe wavenumbers of the periodic structure. The series may be said to be a Floquet modal or grating lobe series representation for the magnetic vector potential. The elements of the dyad, in turn, may be 19

written as pq (Z/ Zt) = 1 Vih (z z')) 3w 0 0 pk[pqIh (z z') - I (z, Z') J ppq 0 0 I e(z, Z'), ijw., C (10) where, for example, V/(z, z') is the voltage at z due to a unit current source at z' on the cascaded transmission lines representing the layered media for TEz (h-) polarized fields with transverse wavenumbers ktpq. Similarly, I{(z, z') is the current at z due to a unit voltage source at z' for the corresponding TMz (e-) polarized fields, and so on. The corresponding scalar potential grating lobe series representations are 1 00 00 - jK': (r, r') = -' E E - P-P')i p=-oo q=-oo (11) where ) Vh(z, z') - Ve(z, ') ppq (12) and P(r, r') = 1 00 00 kp(PP') 1 E E Pzpq(ZZe- (P-P') p=-oo Q=-oo (13) where Pzpq(Zz') k2 [vV (IZ) z -Ve (z,') ppq (14) 20

The electric vector potential is completely dual to the magnetic vector potential above and is given by gF (r, r') - 0 GF (r, 'r') 0 1 GF (r, r') GF r, ')G F(r, r') II I 0ii ~ ~O g(z ZI) 6,Jkpq'(P-PI) Ap=-oo q=-oo (15) where gpq(Z ' Ive(z~z) 0 c kxpq [Ve(ZZI)_-Vh(ZZ,) 0 1 IV6(z~z) 3w 1 jk2 v ZZq Vt'(, zJ 0 0 c IVh,(z Zl) 3W/I (16) Similarly, the corresponding scalar potentials dual to those above are given as follows: K"'(r, r') -: ~ k",(z, ZI)eikC iP' p=-oo q=-oo where (17) (18) "~v t)= vI (z, z') - I,'(z, z') pg kZZ I kpq and Q, Qr, r') I- 0 0 Qzqzz) e-3kpq (P-P') Ap=-oo q=-coo 21 (19)

where Qzpq(Z I z) k2p Ii z ') - Ii,( Z, Z' (20) Various derivatives of the potentials are needed to complete the field descriptions. In particular, the curl of the magnetic vector potential is V xA (r) = Is v x gA(r, r') -J(r') dS' (21) ay aG A~ aG A a9y ax V x>< A(r, r) a(z aGcA aGA a~x atyy ax (22) aGA a9y 0 - 1X j p >< A I' -ikpq' (fPfJ') where the spectral form of the nabla operator is defined as - + A a -AA+ A ( j ktpq+ j9z kxpqjkypq +Z, Thus, we have V O> A ~(Z',Z) vpq xp iikx kyg[ 1h _Ie ] P k-pq [ LhIe IIIh ppq kypq V~h w) i - kxP kYI[h - fte],Ih ppq u kxpq k ijh r -I - - V. La I WE1 V WEI V 0 1 (23) 22

where the transmission line equations are used to replace any z-derivatives with voltage or current quantities again. Similarly, the gradients of the electric scalar potential functions are defined as V K(r, r') 1 00 00 1 i E VS pqkpq,(z, z') -3kpq( Ap=- 0 q=- P p=-oo g=-oo (24) where V pqk7(z z') -jk ( h- ppq jkzpq (kzpq e + Z W I ppq kzp Ih ) zpq / (25) and VP,(r, r') 1 00 00 ( - E E V zpq(q) e —kpq p=-oo q=-oo (26) where VpqPzpq(z,) ') ktp pq L ppq w u' kzpq - v) + z ppq ( Wi/h kze kzpq WC6 J (27) By duality, we obtain the curl of the electric vector potential, V x F(r) = V x F F(r, r') M(r') dS, JS (28) where aGF 9y V X F(r, r') OGF (9z aG9x ) ax ( acy acz Oy az aGF ax a9G ax aGF zz ay aGF ax 0 (29) oaG ay 23

-1 ~X F z' -kpq- (P-P') and where Jvpq xg (pz' 6kxpq kypg e _ Vh] kppq k2BK _Vh] _ ve ppq ~ jf w; Iv — Ve- _ h] + cVe k2 [ - ppq - kxp ve w v - WAkIpg yh1 0 ji (30) Similarly, the gradients, of the scalar potentials are given by V K',(r, r') Ap=-oo q=-oo (31) where - T I lvpq Kpq (z, z Ap=-oo q=-oo h - 'ktpq IV6 - ivtl I k2 I- ppq -.jkzpq ikzqV + z ZP2 h kppq W/I -WC ve1 kzpq V, (32) and V Q, r, r') 1 00 00kp-( P' Ei Z Z VpqQzpq(Zi Z') Cikp3PP' Ap=-oo q=-coo (33) where Vpq Qzpq (zI ') - WCtp f A c zpq ( C ye p pq h w~' ZPVh) ktpq, ) kZ~pq (34) 24

E Extraction of Asymptotic Forms of Series When z and z' are well-separated, the convergence of the grating lobe series is, in principal, exponential. When z or z' are in the same or in adjacent layers, however, we attempt to enhance the convergence of the grating lobe series representations. The procedure for doing this may be summarized as follows: * The asymptotic forms of the transmission line quantities are determined for large spectral values. These forms contain terms of the form e-jkz and e- kzl( —(-z 2jKz 2jkz which may be identified as so-called direct and quasi-static image terms, respectively. They appear to arise as a direct or homogeneous medium contribution from the actual source located at z = z' or from its image at z = 2d - z'. The direct contribution is exact; the image contribution is only asymptotically correct, but becomes exact as u -+ 0. Hence it is often said to be a the quasi-static image. * Terms from the asymptotic series series are removed from the grating lobe series (Kummer's transformation) to enhance its convergence. * The contributions from terms removed from the series must be replaced. It may be shown, however, that they constitute contributions from an array of phased point sources (spatial representation), and may be obtained by using Poisson's transformation on the removed series. We thus obtain a hybrid spectral representation in which most of the layered media effects are in the spectral (grating lobe) series, while effects of the direct and first-bounce images of the source in adjacent layers are primarily contained in the form of spatial series. * The spatial series summations may be accelerated using Singh's method, Ewald's method, or other acceleration methods (discussed later). By standard transmission line theory, the asymptotic forms of the transmission line voltages and currents are obtained as follows: 25

V (z, z') TS jkzIz+z'-2bmi,+lI +-'- /i" ekzlz+z'-2bmIt = k,,u~i-'s) e+kzzz'I m =m -II. ~ v)ekzIz-z'1,r -~[(m'+ Mmm+( +=)+Smm- M+F8)&Uzz' + +mm (i$ e-3kz Iz-z'2bI +i +Tne k-zZbmI) (35) Iis(z, z') I( nz-z)6k, z-z'l 4 g + + '- bi )ekzlz+Z'-2bm1+iI + Fg( z) sgn(z + z' - 2bmi,kt+z2m1),) m~ y1- T'k) sgnzzz- z2emzI ~ PmS' F sgn(z + z' - 2bm,+)e~kc- '2m11 + ~sgn(z~2 - 2bm)eik-zz'bm1) = M f - (36) 26

[IS eikzzlz'2bmI2bI FI 1 (.klz-z'I eIskkzlzz'zb'fbm IIs) iz- Ke - 00c - 00I I5( z) 1 ZIS (1 - Is 3k, Iz-z'l 6-, 00 11 m =M'+ 1, m - 1 /5 - 'kz Jz-zl I rl - 8mm (EI smm/j I- Folo) + sm,ml-i (I - 0) 6-kI-' 6jkzIz+z-2bm1,+iI +?s e~zzz2mI) (37) 1 2 (sgn(z - z )6kzzzI_ F,~ sgn(z + z' - 2bmi,+i) ejkz3 '2m+u - 7P sgn(z + z' - 2bmi)e73kzzz2m11 (1~ F,) sgn(z - e zIzl VVK5( z, z') -1( 2k\ 1+"sgn(z - zI)e-k~zI m=m I 1 2 [(m,mi + sm,m'+i (I + / 8.g-~ + F~)sgn(z - ZI)63ze z -mm'( sgn(z + z' - 2bmt,+)6-3kzz ~ m+ + 'sgn(z + z' - 2bmt)e-k3z '2m/) (38) where 8mm' is defined as 8mm' 0 m $~m 2 7

and z is in layer m, z' is in layer in'. The spectrally asymptotic values of the layer interface reflection coefficients are given by C-M + c'ErbI~ and,h / ilm'~1 -. I 00 Am'+1 + /1im Using these asymptotic expressions, one the potential quantities: determines the following asymptotic forms for OAq(Zz Z') r~j0AO (Z - 0AO2p 01 0 AO 0 yypq [ OA 0AO q AOAJ zxpq ZYPq ZzpqJ (39) where gAO Oxxpq - AO - A yypq 2j'kzpq [ Jm,m' + Sm,m'+i (t + ij~) + SMIMI(1 (I+ F0 e-3,kpq mzz l)],,i th 'kzpq iz-zl I ) e-3 ~ 8 hm,kz.pqlz+z'-2b 1+iI + (40) gAO Ozxpq _ iikxpq jk2~ /C smm/+j (FIZ - JI) + mm'M-i(WOO 100)~) sgn(z - z) Cizp-3 ' 6 m~ml((T'~ - F~) sgn(z + z -F,~)sgn(z +z' -2bm - 2b)Ml+i) 6j3kzpq lz+z'-2bm/ +II,) e..3kpqIz+zI-2bmu I) 28

(41) zyPq - j2 K m,m+l (rF - r) +.mm'-i( _( F 6,rnh - OT)) sgn(z -z) e-kzp -,rn,m' (PIe - rF ) sgn(z + z' - 2bmi+i) e ikzpql+'mi+ - 4 — +(IP - F'h) sgn(z + z' - 2bm) e-kz z+z'-2bmr) (42) gAO zzpq j= [(Sm,m' + Sm,m'+l(1 - Foo) + Sm,m'-l(1 - Lie)) e-jkzp -2 kzpq - - tb, k k-zp lz+z'-2bmI/ kqm' k Iz+z'- m+1 (~m f e me MMF(ZeI (43) Aq(zz') p Kq (z,z') 2= 27k [\(8mm' + -m,m'+i(l+ 1) + memi(l+ )) e kPqlz-'l mml ( kzpq kzz+ 4e ~-jkzlz^zt-2bm,+11 _4._ pc e-Jkpqlz+z'-2bm/l "(mm' lZ (44) 29

gpq zz) O pq ( z ~z) F 1?xp 0 0 1 pq 0FO q 0I yypq gFO ~FO ~FoJ Lzxpq ozypq Ozzpq (45) where gFO!yyp~q 2j'kzpq mm(IK [ ( ml,m+ 6m m/+i (I 6jkzpqlz+z'-2bmt~l I + - Jjo"-) + sm,mt-i1 F0,Jkzpq lz+z'-2bm/ I) 6hJkzpq iz-zII (46) gFO Ozxpq __ kxpq Jk2, [ ( f~ 6 F( -44-4 - T;~) + Sm,m'-i (FOO -F,)) sgn(z z') 6jzp 3.' - m,rm I (F(F - F,~) sgn(z + z' - 2bmi+i) cizpI+Z2mIi + (FZ - F~~) sgn(z + Z' - 2bm) eizPIzz2bII) (47) gFO Ozypq __ kypq jk2p /+( 6(m,m - F4h~ + 8mm(F l -,)sgn(z - z' e-kzqlz' - 8mr,mi (F(F - F,) sgn(z -i z' - 2bmi+i) e-ikzpqlz+z'2bmiI + ( F(' - F,~) sgn(z + z' - 2bm) Cez- I~~bI) 30

(48) zzpq 2= 2 [ m + m,m+ (l 0 + rm+ -i(l + rc)) e 2] j zpq + (' e-jkzpqlz+z'-2bm, + r e-kpql+'-2bm + (49) (49) 1= [(65mm + Sm,m'+(1- Fe) + m,m-l(1 — )) + -jkpqiz' 2w 'kAzpq, L -- m ' (F p-jkkzpqlz+z'-2bm,+1\I + 00 -jkpqlz+Z'-26 ^mm l +c c' -2 m Io (50) Similar asymptotic forms exist and may be obtained for the potentials Pz and Qz as well as for all the derivatives of the potentials, but this task has not yet been completed. Instead, focus has been on the acceleration of the spatial domain representations which arise from applying the Poisson transformation to these terms. A summary of two-dimensional Fourier transform definitions and properties, as well as a brief derivation of the Poisson transformation for periodic structures on skewed lattices is contained in the appendix J to this report. 31

F Code Validation The preceding formulas have been validated by the following tests: * A Fortan 90 code has been developed which handles the multilayered medium problem, but which does not yet handle periodicity. This code is being modified to handle the periodic problem using the acceleration procedures outlined. Its principal differences are 1) the grating lobe sums are actually integrals over the spectral variables, and 2) the direct and image contributions are not an array of contributions, but single point source contributions. Using this code, the following items have been validated: - Electric fields and associated potential representations, all associated transmission line calculations, and incident field expressions for multi-layered media by solving for currents on conductors in multi-layered media using the electric field integral equation (EFIE). - Correctness of direct and quasi-static image extraction terms from boundaries both above and below the source point. - Correctness of magnetic field and associated potential representations by obtaining electric currents on conductors by duality. * A separate FORTRAN 77 code exists which implements the formulation presented here for conductors in multi-layered material only. This code does not have the more advanced acceleration procedures begin developed here, nor the ability to handle magnetic currents or fields. Furthermore the associated code is less flexible and maintainable than the Fortran 90 code. Nevertheless, it serves as a benchmark for validating many of the new code's capabilities. 32

G Acceleration of the Free-space Periodic Green's Function by the Method of Singh The free-space periodic Green's function has the spatial representation 00 -jkRmn G(r,r')= - e-jmk-a -jnkyOb (51) m=-oo n=-oo 47 where Rmn = ^/( -- '- ma)2 + (y- y'- nb)2 + (z - ')2. The wavenumbers ko and kyo are given by ko = k cos bo sin 0o and kyo = k sin o0 sin 0o, where (0o, o) either specify the angles to which the array is scanned or the direction of propagation of a plane wave exciting the array. An alternative, spectral-domain representation of the periodic free-space Green's function is 1 00 00 1 G(r,r') - 1 -jk.n(p-p'l)e-Jkzmn Iz-z' ab?n=n=oo n=-oo Jzmn where k = z(k + y ), p' = x + y y, (tmn — kx 2m7 n -i k/ 2ny0 -7la b and the wavenumbers kzmn are given by kzmn - /k2 - k2 ZmnXm Yn I where Re(kzmn) > 0, Im(kzmn) < 0. 33 (52)

In the method of Singh, acceleration of the spectral form is achieved by adding and subtracting the periodic Green's function for a medium with an imaginary wavenumber k =-ju (which will be called the "modified Green's function"). The spectral form of the modified Green's function is subtracted from the spectral form of the free-space Green's function, and then the spatial form of the modified Green's function is added back on. This results in a hybrid spectral/spatial form of the free-space Green's function, given as G(r r') - - E E ekj tmn'(PP') [2km I- - ] m=-o n=:-00 I Zmn?mn +- E" E e-jmk, a -jnkyo b ---— 3) 00 -- Ul-Rmn m=-oo n=-oo 4Rmn where kZmn in the medium with the imaginary wavenumber has been denoted as rimn, given by Nimn -= U2k- - ^ (54) XCM Yn 7 where Re(Kmn) > O, Im(Kmn) < O. The free-space periodic Green's function is now given as the sum of two series, one denoted as a "spectral" series (the first one) and one denoted as a "spatial" series (the second one). (Of course, these series are different than the spectral and spatial series defined previously, which are alternative representations of the total Green's function.) The variable u is called a "smoothing" or "acceleration" parameter. As u is increased the spatial series converges faster, while the spectral series converges slower. Ideally, u should be chosen to obtain an optimum balance for the convergence rates of the two series, to obtain an accurate calculation with the fewest total terms. The optimum value of u will depend on the error criterion specified for convergence - the smaller the error criterion, the smaller the value of u becomes. This is because the spatial series converges exponentially, while the spectral series converges algebraically. An alternative acceleration method, which is also a hybrid method, is the Ewald method. This method results in two series that each have exponential (gaussian) type of convergence. This method is discussed next. 34

H Acceleration of the Free-Space Periodic Green's Function by the Ewald Method We begin with the spatial domain form of the Green's function, G(r,r') = Z - ~e-jktoo Pmn m-=-oo n=-oo 47Rmn (55) and make use of the identity - 9 2 k e-jkRmn 2 j o-R2 s2 ds Rmn a O (56) where s is a complex variable and the path is chosen such that the integrand remains bounded as s -+ 0, arg(s) E [arg(k) + q, arg(k) + (, (57) (58) and decays as s -4 oo, arg(s) E [-,. (59) (60) Observing that -- < arg(k) < 0, it is convenient to restrict s to the intersection of these two regions, i.e. the angular sector arg(s) E [arg(k) +,. 4[4 (61) (62) Next (56) is substituted into (55), and the parameter E is introduced to split the integral into two terms, as G(r, r') = G1(r, r') + G2(r, r') (63) 35

(64) where G, (r, r') G2(r, r") C-j too'Pn 2~ c — RnS42+ ds, rn=-oo n=-oo C e too'Pmn j;EO e-RmnnsT-2+ ds, m=-oo n=-oo /7 (65) (66) Using the identity 2 Oo 2 2+ k 2 N/7F J eRn 472 dsI 1 2Rmn [- km erfc(RmniE 1k) + eacmlerfc (RmnE + 2 (67) G2(, ')can be written as G2(r, r') ejktoo mn [ ekRmn m=-o~o n=-00 7rm erfc (RmnE ~ e~k erfc(RmnE + )] Poisson transformation of (55) (c.f. Appendix) now yields (68) m=o =oo 4lkZmnekt (P ) x [eukzn3 zerfc (ik Iz - - z'lE) + e ikzmn,~I z-Z' Ierfc (Ik +F lz -zIJE) (69) 36

Each of the new series Gl(r, r') and G2(r, r') generally converge at a much faster rate than the original series G(r, r'). The series G2(r, r') is a modified spatial series, and is therefore denoted as the "spatial" part of the Ewald representation, while G1 (r, r') is denoted as the "spectral" part. 37

H.1 Choice of the optimum E parameter The two series G1 and G2 both converge exponentially. The parameter E controls the convergence rate. As E' becomes larger the the "spatial" series G2 will converge faster, while the "spectral" series G, will converge slower. The optimum parameter is that which makes the two series converge at the same rate, so that equal number of terms are required in the calculation of both series (this assumes that the calculation time for each term in the two series is the same). To prove this assertion, assume that the Ewald parameter has first been chosen so that both series converge at the same rate (E = Eopt), so that each of the series are summed from -N to +N. The total computation time due to the summation of both series is Topt - CN2 + CN2 = 2CN2 (70) where C is some constant. Now assume that a different Ewald parameter has been chosen, and denote r = E/Eopt. In order to obtain the same accuracy in each summation, the summation limits must be changed. The limit for the spatial sum is inversely proportional to E, while the limit for the spectral sum is proportional to E. Hence the new computation time is T= CN r2 + -). (71) The ratio of computation time to that with the optimum Ewald parameter, in order to obtain the same accuracy in the summations, is is TI/Tpt = r2 + (72) The minimum value of the above function occurs at r = 1. Hence, the selection of the optimum Ewald parameter will yield the least number of total terms in the two series, for a given accuracy criterion. To derive a simple formula for the optimum Ewald parameter Eopt, we enforce the condition that the asymptotic rate of convergence is the same for the series GC and G2. For simplicity, the optimum Ewald parameter is derived here only for the case z = z'. However, the value of z - z' does not affect the asymptotic rate of convergence of the series. 38

The two series are written in the form 00 G1= E A(, (73) m=-(nn oo G2= 5 A4^. (74) m t = n m= —,x We then enforce the condition that Am ~ A2) (75) as m, n -+ oc. Note that the m and n indices are switched in the subscripts of the two series in the above equations. This is is because the indices m and n are multiplied by the lengths a and b to construct Rmn in the spatial sum G2, whereas they are divided by these lengths to construct the wavenumber term kzmn in the spectral sum G1. The asymptotic form of the summands is calculated by using the asymptotic form of the complimentary error function, -z2 erfc(z) -~ Equation (75) then yields e- 2E e-RnmE = -- 4abvr (r2E) nEmn 87 Vr(RnmE) Rnm where / 27r M2 277n 2 mn aV -- + (b ) The exponentials in this equation will converge at the same rate provided amn R-nmE. 2E 39 (76) (77) (78) (79)

Hence, E= \RV 2Rnm This yields E ( a b ) (n7r)2 + (nb)2 Multiplying and dividing by /ab yields the result E = b. ab (80) (81) (82) The above choice of E makes the exponential terms in the two series converge at the same rate. It is not completely clear yet if two series will converge at exactly the same rate however, due to the presence of the other terms in Eq. (77). Using Eq. (80), the ratio R of the LHS to the RHS terms in this equation is nE2R2 R = n a2 ab mn (83) In view of Eqs. (80) and (82), it follows that R = 1. Hence, the optimum Ewald parameter is Eopt - (84) 40

I Results Comparing Acceleration Methods During the first quarter of this contract, the Ewald method has been extensively tested and compared with a previous acceleration scheme, the method of Singh et al. The Ewald method has been tested for a wide variety of lattice parameters, including lattice spacing, interelement phase shift, and location of the observation point. Only a few sample results will be presented here. The conclusion is that the Ewald method is a very robust method, requiring only a few terms in each series to obtain very accurate results. The convergence rate is almost always significantly better than that of the Singh method, at least for the case z z', which is of most interest since this is when the original series converge the slowest. The results presented here are for a typical case with a = b = 0.5Ao, x = y = 0.25Ao, and kxo = kyo = 0. Figure 8 shows the convergence rate of the Singh method, using the Singh acceleration parameter u = 1.0 (which is a typical value). The spatial and spectral curves denote the alternative spatial and spectral forms of the periodic Green's function, with no acceleration. The hybrid curve is the combined spatial/spectral representation used in the Singh method. The integer N is the summation limit in the double sums. As can be seen, the Singh method provides significantly better convergence than using no acceleration at all. However, Fig. 9, which shows the same result on a different scale, shows that N must typically be 10 or more before a small error is obtained. Figure 10 shows the percent error in the spatial and spectral parts of the Singh (hybrid) method as a function of N (In the Singh method, the final result is the sum of these two parts - the designation "spatial" and "spectral" therefore has a different meaning in this figure than in the previous figures). In order to obtain very small errors, on the order of 10-5, N needs to be on the order of 30. Figure 11 shows the same type of error plot as Fig. 10, but with the Ewald method, using the optimum acceleration parameter Eopt. The curves labeled "spatial" and "spectral" now denote the series G2 and G1 (which are in essence spatial and spectral series modified by the complimentary error function terms). It is quite remarkable that fairly accurate results (less than 0.1% error) are obtained with only N = 1. And extremely accurate results (percent error less than 10-7) are obtained using N = 2. Many results similar to Fig. 11 have been generated to examine the convergence rate of the Ewald method as the various parameters of the lattice are changed. In all of the cases where the lattice is square (a = b), the conclusion has been the same as for the sample case shown in Fig. 11. With an oblong lattice (a 57 b), the convergence may be slower. This aspect is still under investigation. Also, when z ~= z', faster convergence may be obtained by using an Ewald parameter that is slightly greater than that given by the Eopt formula. This is also under investigation at the present time. 41

.....................................................................................................................................i.. 0 *: * ' _ ___ 0.S \:: of c r r: 0 *. I::L CD CC Figure 8: Comparison of convergence rates for the spatial form of the periodic Green's I:.: ~: u.: _ I. i I:: y ' me T n shon fx y = = 0..?rrrrrr,6f....:....... '~ '~....................... k, ~o - ko-O. 42

I\ ~' - w o * -. \ \*. Fo - /.- X: ~....... - I ~: Figure 9: The same comparison of convergence rates shown in Fig. 8, plotted on a different scale. scale/ It) CD 0 ~~.. ' ( C4 ~ scale 43

0 II ~ f a: I, i I i r r I * <i i Ii U * I I 0 u, o0 0 C) IO tO z 0 0 10 O) T% 0C rr CM 0 Co 0 0 O TO so I w VW (O 0 To I UJm O CD w Wm (N33M3d) Mo0M3 Figure 10: Percent error in the spatial and spectral series that form the hybrid method of Singh, plotted versus the summation limit in the series. The geometry is the same as in Fig. 8. The acceleration parameter u has been chosen as 1.0, which is a typical value. 44

iD t..t t -- --.o i,..... - I:, |, = I $: n I I l —.- / ',. -. * — t. '.: /i:' *. I Iii Li i L Li L L Figure 1: Percent error in the spatial and spectral series that form the Ewald method, plotted versus the summation limit in the series. The geometry is the same as in Fig. 8. The Ewald parameter E has been chosen to be the optimum value E,. ~?: - 2, i i $:... i: /! i I I / / | L.L.....L.i.....-..... o LiI CI LI W I LI Ewald parameter E has been chosen to be the optimum value Eopt. 45

J Two-Dimensional Fourier Transform Properties for Periodic Structure Applications This appendix collects a number of two-dimensional Fourier transform relations useful in periodic structure analysis. Though most of these relations are well known, they are somewhat scattered throughout the literature and often notationally inconsistent. Furthermore, they are often specialized to rectangular lattices, or have other limitations. For reference, therefore, we collect here the most basic definitions and provide simple proofs for the most useful of these formulas. Use of vector notation throughout not only results in more compact expressions, but also emphasizes the coordinate-free nature of the results. J.1 The Two-Dimensional Fourier Transform The two-dimensional Fourier transform of a function f(p) = f(x, y) is defined as F(k) = f (p) e Pdx dy (85) -oo -oo and the inverse transform as f (P) = (27r)2 J F(k) e-kp dk dky. (86) We use the notation f(p) ~ F(k) (87) to denote the correspondence between f(p) and its transform F(k). The exponential factor k.p is often written as k p = k~x + kyy (88) and this form is more convenient when one-dimensional transforms, whose familiarity is assumed, are applied in succession to obtain properties of two-dimensional transforms. However, not only is this form less compact, but it is usually less convenient for analyzing skewed 46

lattices, where lattice points in the spectral and spatial domains (k and p domains, respectively) do not necessarily fall on lines of constant coordinates. Using the Fourier transform definition (85), it is easily proved that f(-p) e kP dx dy = (-k) (89) - ) (-oo -0O OO and (p) e k'P dx dy = F*(-k), (90) -00oo -oo00 where the asterisk denotes complex conjugate. These results are summarized as f(-p) = F(-k) (91) f*(p) F*(-k). (92) J.2 Transform of the Two-Dimensional Delta-Function The two-dimensional delta function may be defined as S(p-p') = (xr- x')8(y- y'), (93) and hence its Fourier transform is given by roo roo ~ i S (p- p')ejkPdxd dy eP' (94) Fourier inverse transforming yields the useful identity (1 2r eik(P'-Pdkdk = (p -p'). (95) Thus, (p - p'),= ejkP (96) 47

J.3 The Two-Dimensional Shift and Multiplication Theorems Consider the Fourier transform of the "shifted" function f(p- p'): 00 -00 J 00J ejk.p F(k), (97) where the substitution p" = p - p' is used. Thus, the shift theorem states that shifting in the spatial domain is equivalent to multiplying by a linearly varying exponential in the spectral domain. When the shift is in the spectral domain, the exponential multiplier appears in the spatial domain, a result which is often called the multiplication theorem: 1 00 f00 (27)2 / F(k - k') e-jP dk, dky = e-k'pf(p (27)2 J2 O J oo These results are summarized as f(P- ') < ejkP' F(k) e-jk"p f(p): F(k-k'). (98) (99) (100) J.4 The Two-Dimensional Derivative Theorems Define V = a- + y Ax 9y (101) and V = X a + Ya Vkx dky (102) 48

Then Vf(p) = (22 (-jk) (k) e-kp dkxdky, (103) so that Vf(p) (-jk)F(k). (104) The rate of change of a function f (p) along the direction of the unit vector i, i I is given by df() = V.f(p) (105) and hence t(p). V f(p) (-jk.t) F(k). (106) Similarly, pf(p) < -jVF(k) (107) and (p- ) f(p) -3 - d ) = -t VF(k). (108) J.5 Two-Dimensional Convolution Theorems The two-dimensional convolution between two functions f(p) and g(p) is defined as f f(p')g(p - p') dx' dy'. (109) 49

With the substitution p" = p - p' and subsequent replacement of the dummy coordinate p" by p', one finds that 00 roo roo roo f(p') g(p -p') dx' dy' g(p') f(p - p') dx' dy', (110) -CO -00 -00 -00 or using the symbol '*' to denote convolution, f(p) *g(p) = g(p) *f(p). (111) The Fourier transform of the convolution integral is 0 f(p) *g(p)e3kPdx dy (112) -00 J-o00 J J /J f(p') [(p - P') eikP dx dy dxdy' - 00 0 L0 - J f(p') G(k)ejk'p' dx'dy' = F(k)G(k), (113) evaluated with the aid of the shift theorem. Thus we have f(p) * g(p) - F(k) G(k). (114) In a similar manner, it is shown that (2 2 J L F(k) * G(k)e-jk' d dkx = (27r)2f(p)g(p) (115) or, equivalently, f(p)g(p) 2 2 F(k)*G(k). (116) A number of alternative forms of the convolution theorem may be obtained by combining (115) or (116) with (89) and (90). 50

J.6 The Power Theorems The power density at a point p is often proportional to the quantity f(p)g*(p), and hence the total power crossing the transform plane is given by _ of(p)g*(p) dxdy -00 00 = J 0 f(p)g*(p)eI"*P dxdy k'= -oo -oo k=0 (27')2 k]=O 1 {CO o - - ( J J F(k)G (k) F kdkx dkyI where (115) has been used. This result, JCO o -00 J-00 f(p)g*(p)dxdy (2) oo J oo (27r)2 -0O 0 F (k) G* (k) dkx dky, (117) is known as the power theorem. If both f(p) and g(p) are real, the conjugation symbol on the left hand side may be removed, resulting in a form often called Parseval's theorem: -o oo f(p)g(p)dxdy = (2-2 i F(k)G*(k) dk dky, (27 )2 _-0 _-0 (f,g real). (118) Combining this with (89) and (90), we have alternatively f J f(p)g(-p) dxdy (2 - 2 F(k)G(k) dk dky, (f,g real). -oo J-oo (r-)2 0oo0 -oo (119) When f(p) = g(p), the power theorem becomes Rayleigh's theorem: / f(p)j2dxdy = IF (k)l2 dk dky. - J -o (22 51 51 (120)

J.7 The Two-Dimensional Sampling Function The one-dimensional sampling or replicating function comb (x) is a periodic train of delta functions defined as comb(x) = E 6(x-n) n=-oo (121) and which has the Fourier transform 0o oo N comb (x) e kdx - ek = lim ekn -o no N —)oo (-1n J- 00 72=-00 7=-N - lim (2N + l1)sn(2N )2 N-oo (2VN + 1) sin 2 ={ 0, k ~ 2m7r co, k = 2mr. (122) These properties suggest that the transform is also a periodic train of delta functions located at k = 2m7r for m an integer. To determine the magnitude of each delta function, we calculate the area under each as follows: I(2m+1)7r I N i lim e jkn dk = 27r d(2m-1 )r \ n=-N (123) where the orthogonal properties of the exponential functions are invoked. Thus, we assert f00 00 ] comb(x) ejkdx = 27r 5 (k - 2m7r) = -00 m=-o = comb (k ) 00 - 6( -m) m=-oo (124) so that the comb function is found to be its own Fourier transform. Comparing (122) and (124), we also establish the important identity comb ( ) = 2wr n=-oo Th 00 (125) 52

Now consider the two-dimensional version of the sampling or function E E ((P-Pmn), (126) m=-o n=-oo where Pmn = ms1 + ns2 and s, s2 are the lattice vectors of a skewed cell of a periodic structure. Eq. (126) represents, for example, an array of point sources corresponding to the scalar Green's function of a periodic structure. The Fourier transform of (126) is 00 00 00 00 ikP dx dy / E E (p -Pmn)) ejPdxdy mx m=-oo n=-=-oo 00 co = h ek ^Pmn m=- o) n=-oo 00 CO m= —Do n=-oo - > 3 e Jmk slejnk.s2 m=- o n=-oo00 comb (k s) comb (k ), (127) where (125) is used. The product of the two comb functions on the left hand side of (127) forms another sampling or replicating function on the so-called reciprocal lattice in the transform domain. The principal planes of the reciprocal lattice are perpendicular to si and s2 and spaced a distance of 2i and 27r respectively. Hence if the spatial domain lattice is skewed, so is the reciprocal lattice. The reciprocal lattice vectors are easily found to be 27r(s2 X z) A kl = — A — 27 r( x Si) k2 = A — A (128) where A = Isl x s21 is the area of a unit cell in the spatial lattice, and the lattice vectors satisfy the biorthogonality conditions ki sj = 27rSi, (129) 53

where &ij is the Kronecker delta. Eq. (127) implies that the transform of the two-dimensional sampling function is also a sampling function with two-dimensional delta functions located at the lattice plane intersections where neither comb function in (127) vanishes. However, the amplitude of each of these delta functions depends on the lattice skewness, and it may be determined by integrating over a unit cell in the reciprocal lattice plane. For this purpose, we introduce the normalized unit cell coordinates k s1 - S2 62 (130) 27r Since this transformation of coordinates from (k,, ky) to (&1, 2) is linear, the Jacobian of the transformation is constant and is easily found to be ( merely by equating areas of the reciprocal lattice unit cell computed both in lattice and normalized coordinates: 27r(s2 x i) 27r(, x 1i) [unit dkx dky = (k1 xk2). = 27r( ) 2 x s X J cell A A (22 (x ). = (27)2: (81 xs2). A A2 A (2) d (131) Thus since unit comb ( ) comb (k) dkx dky (132) cell (27r )2 I.1 b (27rw)2 = ( comb (1) comb (2)dld = A- - (133) each two-dimensional delta function of (127) has amplitude (2-) and therefore the sampling function and its transform are related as (2oo7 c2 oo 00 E E 6(p-Pmn) A 2 AE ] 6(k kpq), (134) m=-oo n=-oo p=-oo q=-oo where kpq = pkl + qk2. 54

A similar relationship useful in periodic structure analysis is obtained using (134) with the shift theorem: E S(P - Pmn) e-jkooPmn A E (k-kpq) (135) _ rA m=-oo n=- o p=-oo q=-oo where now kpq = koo + pkl + qk2. J.8 The Poisson Transformation Choosing g(p) to be the two-dimensional sampling function with linear progressive phase shift, (135), in the power theorem yields the so-called Poisson transformation: 00 C 00 CO \ * ) (P E E S(p-Pm) e-koo dxdy (136) - 00 - 00 nq=- -oo f (Pmn) e i P m=-oo n=-oo cxoo / (97r}2 00 00 (27r)2 F(k) (A Z (k- kpq) dkx dk 1 00 co = E E F(kpq) (137) p=-00 q=-00 The resulting transformation, E E f(pmn)ejkoo.pn = F(kpq) (138) m=-oo n=-oo p=-oo q=-oo converts a spatial summation into a spectral sum, often with improved convergence properties. 55