A Diagonalized Multilevel Fast Multipole Method With Spherical Harmonics Expansion of the k -Space IntegralsThomas F. Eibert, Member, IEEEAbstract—Diagonalization of the fast multipole method (FMM) for the Helmholtz equation is usually achieved by expanding the multipole representation in propagating plane waves. The resulting -space integral over the Ewald sphere is numerically evaluated. Storing the -space quadrature samples of the method of moments (MoM) basis functions constitutes a large portion of the overall memory requirements of the resulting algorithm for solving the integral equations of scattering and radiation problems. In this paper, it is proposed to expand the -space representation of the basis functions by spherical harmonics in order to reduce the sampling redundancy introduced by numerical quadrature rules. Aggregations, plane wave translations, and disaggregations in the realized multilevel fast multipole method (MLFMM) are carried out using the -space samples of a numerical quadrature rule. However, the incoming plane waves on the ?nest MLFMM level are expanded in spherical harmonics again. Thus, due to the orthonormality of spherical harmonics, the testing integrals for the individual testing functions are simpli?ed into series over products of spherical harmonics expansion coef?cients. Overall, the resulting MLFMM can save a considerable amount of memory without compromising accuracy and numerical speed. Index Terms—Fast integral methods, fast multipole methods, integral equations, iterative methods.

I. INTRODUCTION

M

ETHOD OF moments (MoM) solutions of integral equations are among the most successful numerical methods for the solution of electromagnetics radiation and scattering problems. Most important is the robustness and insensitivity against dispersion errors. The drawback of large computation complexity due to the fully populated system matrices has been overcome since the introduction of fast integral methods such as the fast multipole method (FMM) and its multilevel versions [1], [2]. In this paper, we will concentrate on diagonalized high-frequency versions of FMM working with propagating plane waves and associated multilevel fast multipole methods (MLFMMs) [1], [3], [4]. Computational implementation of diagonalized FMM algorithms using propagating plane waves is based on numerical integration of the corresponding -space integrals over the Ewald sphere. The quadrature sampling rate is dictated by the spectral content of the diagonalized translation operator and this results in oversampling of the -space representations of the basis/testing functions of the MoM procedure employedManuscript received October 16, 2003; revised August 5, 2004. The author is with Institute for High-Frequency Physics and Radar Techniques, FGAN e.V, 53343 Wachtberg, Germany. Digital Object Identi?er 10.1109/TAP.2004.841310

to solve the integral equation [4]. To obtain a computationally ef?cient algorithm, the -space samples of the basis/testing functions should be precomputed and stored in core memory of the computer prior to starting the iteration procedure. Due to the aforementioned oversampling, a large amount of core memory is wasted and considerable memory savings can be achieved by using a more adequate -space representation of the basis/testing functions. For both electric ?eld integral equation (EFIE) and magnetic ?eld integral equation (MFIE), FMM can be formulated using basis/testing function components tangential to the Ewald sphere only. Therefore, it is suf?cient to store -space representations for these two vector components. However, the spectral content of Cartesian vector components is smaller than that of vector components tangential to the Ewald sphere. Therefore, storage of -space representations of three Cartesian vector components can be more memory ef?cient than storage of two vector components tangential to the Ewald sphere. Moreover, Cartesian vector components are continuous functions on the Ewald sphere allowing for an ef?cient -space representation using spherical harmonics without any problems arising from Gibb’s phenomenon. In this paper, the concept of spherical harmonics expansions of the -space integrals motivated by an ef?cient -space representation of the basis/testing functions is worked out. It is shown how a computationally ef?cient MLFMM can be realized leading to considerable memory reductions as compared to previous implementations. Spherical harmonics expansions are well known in formal FMM derivations [4], but in contrast to this paper practical implementations work with discretized -space representations only. Discussions on memory requirements and error analysis are given and various results are shown. II. FORMULATION ) surWe consider a time harmonic (time dependence face integral equation formulation using the EFIE for metallic objects. However, the presented ideas can be applied to MFIE and combined ?eld integral equation (CFIE) formulations for metallic and dielectric objects as well. According to [1], [5], a Galerkin-type MoM equation system can be derived as (1) where cients, are the unknown surface current expansion coef?are the excitation vector elements due to an incident

0018-926X/$20.00 ? 2005 IEEE

EIBERT: DIAGONALIZED MLFMM WITH SPHERICAL HARMONICS

815

wave or delta-gap voltage sources, and is the number of unknowns. The MoM matrix elements are given by

A. Multilevel Fast Multipole Method In the following discussion of an MLFMM implementation, we only describe the steps that are different from standard implementations such as described in [1] and [2]. For an ef?cient algorithm, the number of spherical harmonics must be considerably smaller than the number of quadrature points used to numerically compute the integral in (3). This assumption is typically valid if the spherical harmonics expansion is employed for Cartesian vector components of the basis functions. In the initialization step, the spherical harmonics expansion coef?cients are computed and stored for all basis functions using (5) and (7). This can be done by applying Gauss Legendre or other quadrature rules [4]. This step adds some computation time to the initialization time of MLFMM but is typically negligible as compared to other initialization steps. For the evaluation of matrix-vector products within the iteration loop of the iterative solver, the summation of the outgoing waves of all basis functions in each MLFMM group on the ?nest level is carried out according to

(2) is the unit dyad and is the wavenumber in the considered medium. Also, and are Rao–Wilton–Glisson expansion and testing functions, respectively [5], and is the surface of the scattering or radiation object. A diagonalized FMM representation of (2) as basis for an ef?cient evaluation of matrix-vector is found to be [1] products involving (3) with the translation operator (4) and the -space representation of the basis functions (5)

(10) where is the index of the considered group, is the number of basis functions in this group, and is the expansion coef?cient for the whole group. These computations are done using Cartesian vector components and typically need less computer time than the corresponding summations of quadrature samples in standard MLFMM. However, some additional computer time is required to subsequently obtain the outgoing waves at the quadrature points for all groups on the ?nest MLFMM level using (11) with the quadrature locations , . The Cartesian vector com-

and indicate the group centers of the FMM groups conand , respectively, and is the unit vector in the taining direction of . Also, is the second kind spherical Hankel is the Legendre polynomial of degree , function of order , and the denotes complex conjugation. In order to achieve a memory ef?cient representation of , we employ a spherical harmonics expansion according to (6) with orthonormalized spherical harmonics sion coef?cients are obtained from [6]. The expan-

(7) If we further introduce (8) and utilize the orthonormality of the spherical harmonics, (3) simpli?es to the series (9)

are transformed into spherical vector ponents of components and only the and components are retained. The translations of outgoing waves into incoming waves as well as the aggregations and disaggregations between different MLFMM levels including interpolation and anterpolation are performed as in standard MLFMM using the numerical quadrature samples [1], [2]. However, when all incoming waves are collected in a certain group on the ?nest level, a spherical harmonics representation is computed by virtue of [see (8)]

(12) employing the numerical quadrature rule, where is the index of the receiver group on the ?nest level. Prior to the evalution of the integral, it is important to transform the spherical into Cartesian vector components in vector components of order to avoid Gibb’s phenomenon due to the discontinuities in

the spherical vector components and to guarantee that only few spherical harmonics are required. This step adds again computation time as compared to standard MLFMM, but this is compensated by the reduction of computation time in the ?nal testing step according to [see (9)] (13)

B. Error Analysis and Memory Requirements Detailed error analyses of MLFMM have been given in various papers such as [1], [3], and [4] and are beyond the scope of this paper. Also, it is clear that the error of the introduced spherical harmonics expansion of degree is controllable once we assume that the error of numerical integration in conventional MLFMM is controllable. The chosen value of is basically dependent on the spectral content of the radiation patterns of the basis functions as given by (5). It is found that can be one degree smaller if the spherical harmonics expansion is performed for Cartesian vector components instead of spherical vector components. This does not seem to be substantial but is signi?cant if is small. Also, it is noted that spherical harmonics expansions for spherical vector components exhibit bad convergence behavior due to Gibb’s phenomenon caused by the and . discontinuities at Due to the translation operator, the spherical harmonics exthan that in (6). However, pansion in (8) requires a larger due to the orthogonality of the spherical harmonics, the higher order contributions of (8) will be strongly attenuated by the corresponding expansion coef?cients of (6) in the evaluation of (9). Therefore, we may argue that the employed for both expansions is basically determined by the value required for the expansion of (5). In [4], it is stated that the spherical harmonics (10) representing the outgoing waves of a group are rapidly decaying for [ : order of multipole expansion; see (4)] in an appropriately designed MLFMM. Since this statement refers to spherin our imical vector components, we may choose plementation with Cartesian vector components. For the storage according to (7), we need of all Byte, where is the number of basis functions, in case of EFIE and in case of CFIE [2], and for single precision and 16 for double precision. The leading “3” in this expression indicates that we need to store three Cartesian vector components. In contrast to this, we need Byte if we store the -space samples 1 Gauss Leat the quadrature points. According to [4], 2 gendre quadrature points were assumed and only two spherical vector components tangential to the unit sphere must be stored. In both cases, the symmetry of the basis function radiation patterns was utilized [2]. Assuming a 1 Mio unknown problem by solution with CFIE, the memory requirement with for storing the outgoing waves of all basis functions is reduced MByte to MByte (equivafrom lent to a factor of 5.2) using the spherical harmonics expansion for Cartesian vector components of the basis functions. Further considerable memory reduction is achieved by representing the

Fig. 1. Bistatic RCS of metallic stealth ?ight object, f = 1000 MHz, VV polarization, emitter at 160 . More information about object and measurements is found in [7]. In contrast to [7], we worked with the full-size model. TABLE I COMPARISON OF CPU TIMES AND MEMORY REQUIREMENTS (TOTAL MEMORY AND MEMORY OF THE RADIATION PATTERNS OF THE BASIS FUNCTIONS) OF CONVENTIONAL MLFMM (CON) AND SPHERICAL HARMONICS EXPANSION MLFMM (SEP ) FOR THE PROBLEM IN FIG. 1 [L: SEE (4), P : SEE (6) AND (8)]

outgoing and/or incoming waves of the individual groups on the ?nest MLFMM level in spherical harmonics. III. RESULTS We consider the stealth ?ight object illustrated in Fig. 1 discretized using 79 173 unknowns and discussed in detail in [7]. CFIE computations for varying values of the multipole order [see (4)] are compared for conventional MLFMM and MLFMM with spherical harmonics expansion (SE-MLFMM) as introduced in this paper. The numbers of numerical quadrature points 1 [4] in both cases and the were chosen according to 2 group size on the ?nest level was (six levels in total). A biconjugate gradient (BiCG) solver requiring two matrix-vector products per iteration was used to solve the linear equation systems (residual error: 1 10 ). The radar cross-section (RCS) results in Fig. 1 compare conventional and SE-MLFMM reand as well as . For , sults for the SE-MLFMM results are almost identical with the convencauses deviations on tional results, whereas choosing the order of 2 dB, which is typically not acceptable. Convenstarts to show similar deviations tional MLFMM using results. Table I compares memory but is closer to the and central processing unit (CPU) time requirements for various parameter combinations. For small , the memory requirements (single precision for memory intensive data) are dominated by the near-coupling matrix (185 MByte) and the savings by using SE-MLFMM are little. However, for larger

EIBERT: DIAGONALIZED MLFMM WITH SPHERICAL HARMONICS

817

the memory savings increase since the choice of is not recommended in practical comnot affected by ( putations). Using SE-MLFMM, high solution accuracy can be achieved with very little additional memory since the memory demand for the coarser MLFMM groups and the translation, aggregation, and disaggregation operators is almost negligible. From the CPU times per iteration for an AMD Athlon 1900+ given in Table I, it is clearly seen that the SE-MLFMM is even more ef?cient than the conventional MLFMM for our speci?c code implementations. This is because the CPU time for the additional steps in SE-MLFMM [see (11) and (12)] is compensated by reduced CPU time requirements for (10) and (13). The largest CFIE problem we did on a PC-level computer was a metallic sphere with a diameter of 40 represented by a triangular surface model with 874 179 electric current unknowns. Seven MLFMM levels were used, where the group size on the and the degree of the multipole ex?nest level was together with . The pansion on the ?nest level was 10 ). BiCG solver required 55 iterations (residual error: The total CPU time was 25 920 s, where 18 832 s were consumed for the iterations. The MLFMM setup took about 591 s, and about 6497 s were needed for explicitly computing the near-coupling matrix elements. The totally consumed memory was about 1562 MByte, including 760 MByte for the near-coupling matrix, 252 MByte for the radiation patterns of the basis functions, 290 MByte for the incoming waves of the groups in the various levels, and 57 MByte for the translation operators. 1 quadrature points Using conventional MLFMM with 2 would change the memory for the radiation patterns of the basis functions from 252 to 960 MByte and the memory for the incoming waves of the various groups from 290 to 430 MByte. IV. CONCLUSION This paper presented a methodology for optimizing the memory requirements of the diagonalized multilevel fast multipole method without compromising the speed and accuracy of the algorithm, where an ef?cient storage of the -space representations of the individual basis functions is achieved by an expansion in spherical harmonics. For a variety of numerical problems, it was found that a spherical harmonics expansion with degree of two (equivalent to six independent functions)

Thomas F. Eibert (S’93–M’97) received the Dipl.-Ing. (FH) degree from Fachhochschule Nürnberg, Nürnberg, Germany, in 1989, the Dipl.-Ing. degree from the University of Bochum, Bochum, Germany, in 1992, and the Dr.-Ing. degree from the University of Wuppertal, Wuppertal, Germany, in 1997, all in electrical engineering. From 1997 to 1998, he was with the Radiation Laboratory, Electrical Engineering and Computer Science Department, University of Michigan, Ann Arbor. From 1998 to 2002, he was with Deutsche Telekom, Darmstadt, Germany. In 2002, he joined the Institute for High-Frequency Physics and Radar Techniques of FGAN e.V., Wachtberg, Germany, where he is Head of the Department of Antennas and Scattering. His major areas of interest are numerical and analytical techniques for electromagnetic and terrestrial wave propagation problems from low frequencies up to millimeter waves and all kinds of antenna technologies for radar applications. Dr. Eibert is member of URSI Commission B and the German VDE/ITG.