A simple numerical experiment of GREEN’s function expansion in the Fast Multipole Method

In this paper the theoretical foundation of the fast multipole method (FMM) applied to electromagnetic scattering problems is briefly presented, the truncation of the GREEN’s function expansion is revisited, and the well established truncation criteria, in terms of the relative accuracy of the solutions of the electric field integral equation, is revised from a numerical experiment. From this numerical procedure an interesting result for the number L of poles is reported. In FMM L is the number of terms in the GREEN’s function expansion and it determines the precision of such an expansion. In our experiment a lesser value of L is obtained compared to previous studies.


Introduction
It is well known that entries of the system matrix [Z], generated in electromagnetic scattering and radiation problems solved via the Galerkin's form of Method of Moments (MoM), are reactions [1]. In the ordinary implementation of MoM reactions (interactions) are calculated individually, consequently the bigger the number of unknowns the greater the time expended in filling the matrix [Z]. If the Fast Multipole Method (FMM) is used inside the MoM, interactions are divided in nearby and distant interactions by grouping the source sub-domains in, say, M cubic clusters. A cubic cluster is usually defined as the set of sources subdomains within a cube of diagonal d = √ 3w, where w is a fraction of wavelength λ, typically w = λ/2 or w = λ/4 [2]. A nearby interaction is defined as the interaction between sources in the same cluster, or between sources belonging to adjacent clusters. A distant interaction is defined as the interaction among sources belonging to groups separated by one or more groups. In applying the FMM, at first, the matrix [Z] is partially filled with nearby reactions, which are computed in the classical pairwise fashion. By the way, entries of [Z] corresponding to distant interactions are filled with zeroes; therefore, an intermediate sparse matrix [Z near ] is obtained. On the other hand, distant interactions are calculated by factoring the scalar GREEN's function by means of a multipole expansion. In this manner, the scalar GREEN's function is expanded in a product of three factors which are pre-calculated, stored, and repeatedly used for all distant interactions. A considerable amount of time is saved in this way as compared with the time needed in calculating these interactions in the ordinary way. Then, the contribution of distant interactions are incorporated automatically inside a iterative solver avoiding their explicit computation in [Z].
In the literature, there are surprisingly few papers on numerical experiments regarding the truncation parameter of the multipole expansion. There are some, as [3,4] where the topic is treated in general form. In [5,6,7] the FMM truncation criteria is discussed for the case of Multi Level FMM.
In its original conception, and by definition [3], the FMM is based on the expansion of the kernel inside the integral equation, i.e. the GREEN's function. Indeed, this original approach is used in all FMM implementations [2] getting very accurate results.
In [3,4] heuristic formulas for determining the value L of terms, or poles, in the expansion of the GREEN's function G are given. These formulas are based on a given error of such an expansion in comparison to G itself. In this article the value of L is determined based on the error of estimated weights I n calculated using such an expansion of G in comparison to values of I n calculated using G with no expansion.
The main objective of this paper is to investigate the truncation criteria for the GREEN's function expansion by carrying out a simple numerical experiment based on the fundamentals of the single level FMM. This paper is organized as follows: In Section 2 the basics of the multipole expansion in electromagnetic scattering problems are briefly described. In Section 3, we perform a numerical experiment to determine a proper value for the number L of poles in the multipole expansion and a comparison with the values of L derived from common semi-empirical formulas proposed in the literature is made. The conclusions are presented in Section 4.

The multipole expansion in the FMM
The key idea in the FMM resides in the factorization of the GREEN's function G(r, r ′ ) for distant interactions in the form [3,2] α G(r, where is the field at the observation point r within the cluster with center at a due to a source located at r ′ α within the cluster with center at b as depicted in Figure 1. The expansion (1) includes the so called aggregation, transmission, reception, radiation, and disaggregation functions.
Let us review these functions in more detail. The sources in the cluster with center b produce a field at r given by The term r − r ′ α in (2), where α = 1, ..., N , can be expressed as (see Fig. 1) can be further expanded as follows [1]: where d = |d| and D = |D|, with d < D, and j ℓ (κd) and h (2) ℓ (κD) are the spherical BESSEL and HANKEL functions of the first and second kind, respectively, whereas P ℓ (d·D) is the LEGENDRE polynomial of degree ℓ, andd andD are normalized vectors along the corresponding directions. The identity in (3) expresses the field radiated from the source point b in terms of waves originating from the point a.
Since the function j ℓ (κd)P ℓ (d ·D) can be expanded as an integral of properly weighted plane waves in the form [8]: (4) whereκ = sin θ cos φx + sin θ sin φŷ + cos θẑ, and the integration domain is a unit sphere. Thus, the expansion (3) can be rewritten as where ‚ 4π ( )dΩ ≡´2 π 0´π 0 ( ) sin θdθdφ. For numerical purposes the summation in (5) needs to be truncated at some point, say ℓ = L, and the integral must be evaluated numerically; values for L and cubature rules for the integral are suggested in [3]. Such a truncation enables us to interchange the summation and integration. Hence, we obtain where the integral has been approximated by the cubature rule ‚ 4π f (κ) dΩ ≈ N q=1 f (κ q ) Ω q , whereκ q denote the sampling/quadrature nodes. As suggested in [3], 2L azimuthal samples and L polar samples, for a total of N = 2L 2 , have been taken over the unit sphere. The function in brackets is commonly known as transfer function, viz.
and it plays the role of f 2 (a, b) in (1), mediating the interaction from b to a. The function e κκq·d in (6) can be factorized in the form where e κκq·rra is a sort of reception function, which yields the disaggregation process by the cubature rule, whereas e −κκq·r r ′ α b is a sort of radiation function; hence, the coherent addition of all the radiation functions belonging to sources within the cluster centered at b yields the process of aggregation.
Finally, by grouping all previous results, we arrive at [3]: where we see that the radiation function e −κκq·r r ′ α b and the reception function e κκq·rra play the roles of f 1 (r, a) and f 3 (b, r ′ α ) in (1), respectively.

Revision of truncation criteria for the expansion of the GREEN's function
The GREEN's function factorization described in the previous section is a key aspect of FMM. This factorization is achieved by expanding G(r, r ′ ) = e −κ|r−r ′ | /|r − r ′ | in the form where we have coined the GF E abbreviation to such an expansion. GF E stands for GREEN's Function Expanded and it is a function of the number L of poles in the expansion. The error ǫ G (L) of the estimation is defined as follows: where G is the exact value of the GREEN's function for a given set of source and test points. The order of ǫ G is the accuracy of the estimation GF E(L), and it depends on the value of L. It is well known that the term T (κ,κ, D) in the expansion of Eq. (9) exhibits a kind of oscillation for L greater than κD, causing inaccuracies in the value of GF E(r, r ′ ) regarding the exact value G(r, r ′ ). On the other hand, L must be greater than κd for convergence of Eq. (9). In [3] two semi-empirical formulas for setting a proper value of L for a given precision are suggested: where α = 5 is used for single precision (∼ 10 −6 ), and α = 10 is used for double precision. Further, CHEW and SONG in [4] stated that Eq. (11) yields a good approximation of L for κd up to 40, but the formula where β is the number of digits of accuracy, is always a good approximation. It must be noted that the terms accuracy and precision are used here equivalently. Equation (12) is due to ROKHLIN as well. In [2] it is stated that the minimum among L (from any of the two semi-empirical formulas given previously) and κD must be chosen as the actual number of multipoles to be used in the expansion given in Eq. (9): In order to validate the value of L suggested in Eqs. (11) (probably α = 5) and (12) (probably β = 6) for an accuracy of 6 digits (or 10 −6 ), we have carried out a numerical procedure for measuring the relative accuracy of the solutions of the Electrical Field Integral Equation (EFIE) in a problem of radar cross section (RCS) estimation. In this procedure an array of two perfect electric conductor (PEC) squares, with dimensions λ/2×λ/2, which have been separated by distances D = λ, 2λ, 3λ, and 5λ successively, has been used as the benchmark scatterer (see Fig. 2) leading to four cases of study. The Method of Moments was applied to solve the EFIE, and the Rao-Wilton-Glisson (RWG) basis functions were used as basis and weighting functions. An incident frontal plane wave was used as excitation.
For the purpose of this paper the error of the estimation has been defined as follows: where I G n is the magnitude of the n th electrical current weight of the n th RWG basis function computed from the application of the ordinary Method of Moments, and I GF E n (L) is the magnitude of the same coefficient computed by the same procedure but expanding through Eq. (9) the GREEN's function for far interactions. Then, we have defined the order of error ǫ I (L) as the relative accuracy of the estimated weights I GF E n (L). The computation of weights I n was accomplished from the direct solution of the equations where [V ] contains samples of the incident plane wave impinging frontally to the scatterer. The entries Z mn = Lf n , f m of matrices [Z] in Eq. (15) and [Z near ] in Eq. (16), where f n.m is a RWG function, were computed numerically as described in [9]: The coefficients of matrix [Z far ] were computed from (17) by substitution of Eq. (9) in it [2]: matrix I is the identity dyadic, and For all the simulated scenarios L α , L β , κd and κD were also computed, and the results are shown in Table 1.  Table 1 and for the case D = λ the actual value of L to be used is 7, which is less than L β = 13 (for β = 6), so there would not be guarantee of achieving the desired accuracy of 6 digits if such small separation of groups is used. For the rest of cases, that is D = 2λ, 3λ and 5λ, the actual value of L that would be used in accordance with the formula (12) is 13. Despite those values of L, we have found that from our numerical procedure (see Fig. 3) a value of L = 6 is enough for this specific problem to achieve a relative accuracy of 6 digits in the EFIE's solutions.

Conclusions
In this article we have revisited the theoretical foundation of the Fast Multipole Method and established criteria for truncating GREEN's function. For the truncation of the GREEN's function it is well known that the number L of multipoles must satisfy the condition κd < L < κD, and to find an optimum value of L some semi-empirical formulas have been proposed in the literature. In the usual approach the accuracy of the solution is given in terms of the accuracy of the GREEN's function expansion. In the proposed experiment, we have set the value of L based on in the accuracy of EFIE's solutions. As a result, for a given accuracy in the EFIE's solutions a lower value of L is necessary in comparison with that required for the same accuracy in the GREEN's function. This result could be explained by considering that, as a summation is involved in the calculation of weights, far interactions are so weak in comparison with the near ones, that their contribution is greatly minimized.