An Adaptive Rational Fitting Technique of Sommerfeld Integrals for the Efﬁcient MoM Analysis of Planar Structures

: The integral equation method is one of the most successful computational models for microwave devices or integrated circuits in planar layered media. However, the efﬁcient and accurate evaluation of the associated Green’s function consisting of Sommerfeld integrals (SIs) is still a remaining challenge. To mitigate this difﬁculty, this work proposes a spatial domain rational function ﬁtting technique (RFFT) for SIs so that the approximation accuracy is controllable. In conjunction with an adaptive sampling strategy, the proposed RFFT minimizes the orders of rational functions, and the resultant SI evaluation efﬁciency is optimized. In addition, we investigate the semi-analytical singularity treatment for the rational expression of SIs in method of moment (MoM) implementation. Extensive simulation of representative planar devices validates the correctness of the proposed method and demonstrates its superior performance over conventional SI approximation methods.


Introduction
Since their invention in the last century, microstrip circuits and antennas have become some of the most popular types of devices delivering electromagnetic energy and signals [1][2][3]. Operating from 100 MHz to 100 GHz, the original microstrip structure has now evolved into multiple layered planar configurations where thin metallization lies between the interface between two dielectric substrates [4]. In the era of 5G communication, such planar devices become more ubiquitous due to the inherent nature of integration and minimization for stratified structures [5].
Compared with volumetric discretization schemes such as the finite-difference timedomain (FDTD) [6] and finite element method (FEM) [7], the integral equation method discretizing only the metallization surfaces has had great success in modeling these planar microwave device because of its accuracy and efficiency [8][9][10][11]. With Green's function consisting of Sommerfeld integrals (SIs) in planar layered media available, the integral equations are established for the true or equivalent currents on metallization patterns only, and the approximated truncation boundary condition is avoided. In the past decades, considerable efforts have been made to improve the accuracy and efficiency of the method of moment (MoM) that solves these integral equations numerically. Generally, the efforts to fasciate the popularity of integral equations methods can be categorized into three kinds: (a) the integral equation formulations with higher accuracy or lower complexity in numerical implementation, (b) dedicated strategies to accelerate numerical discretization of integral equation, i.e., filling the impedance matrix, and (c) fast solvers for the matrix equation so electrically large or complicated structures can be efficiently simulated.
Unlike the free-space case, the optimal formulation for the integral equation in the planar layered medium is not unique because the scalar potential caused by horizontal and vertical dipoles are different [12]. To avoid the hypersingular behavior in the integral equation that cause difficulties in numerical implementation, three mixed-potential integral equation (MPIE) formulations (noted as A, B, and C) with singular weakly kernels were proposed [13], and formulation C was widely used in later studies because of its convenience in modeling targets penetrating dielectric interface. Later on, modified or extended MPIE formulations with improved numerical performance were developed for various applications at hand [9,14]. For instance, to overcome the low-frequency breakdown problem and support wideband simulation, the argument electric field integral equation (A-EFIE) [14] and current and charge integral equation (CCIE) [15] were developed. Most recently, a non-Galkerin quasi-MPIE formulation was proposed for efficient scattering analysis of homogeneous dielectric targets in a layered medium [16].
In a layered medium, Green's function does not have a closed-form expression and is written as a tensor consisting of infinite integrals, also known as Sommerfeld integrals (SIs). Due to the slow decaying and oscillating kernel of SIs, the partition and then summation numerical integration for SIs is time-consuming. This causes the discretization of integral equations (i.e., the efficient filling of impedance matrix) has been a long-term numerical challenge that limits the performance of MoM when simulating planar circuits or antennas. To accelerate the convergence of numerical integration, extrapolation techniques such as weighted averages methods and double exponent quadrature were developed [17,18]. On the other hand, the closed-form approximation to SIs provides a faster solution in terms of calculation efficiency, and representative methods of this kind include the discrete complex image methods [19,20] and rational function fitting technique [21]. With improved robustness and accuracy in recent years, the unpredictable accuracy and numerical instability issue of closed-form approximation methods still remain to be addressed [5]. Compromise strategies between evaluation efficiency and the accuracy of SIs include the use of interpolation techniques [22] and the tailored DCIM method [23], where the accuracy of SI evaluation is firstly refined before the MoM procedure. However, due to the singularities in the layered Green's function changing dramatically in the near-field region, a robust interpolation scheme with moderate sampling density is still challenging.
Besides the efforts with SI evaluation, another kind of approach to accelerate the filling of the impedance matrix is to adopt proper basis functions or integration methods to reduce the computational overload when evaluating the relevant reaction integrals. For example, based on DCIM, the analytical expression for evaluating double integrals for planar microstrip geometries is derived in [24] to avoid the double 2D integration. Furthermore, studies show that the double 2D integration can be transformed into a quasi-1D integration in a more general multilayer medium, which further decreases simulation time [10]. With the development of modern computation architecture, parallel computation techniques such as graphic processing units (GPU) are also considered to accelerate the MPIE simulations in recent years [11].
As any progress in SI evaluation techniques will directly result in faster and more accurate simulation tools for circuits and antennas in layered medium, in this work, we consider a closed-form rational function fitting technique (RFFT) for SIs. To minimize the order's RFFT expression with a controllable accuracy, an adaptive sampling procedure is proposed for a given spatial region of interest. In addition, a semi-analytical scheme to evaluate impedance matrix elements with a nearly singular kernel is presented to alleviate calculation difficulties in numerical integration.
The rest of the paper is organized as follows. In Section 2, after a brief introduction to MPIE formulation, the adaptive RFFT for SIs in Green's function is presented and verified. The Duffy transform to handle near-singular integration based on the rational expression is derived in the second subsection. Numerical results and comparisons for planar circuits and antennas in multilayer medium are presented in Section 3 to demonstrate the effectiveness of the proposed methods. Finally, conclusions and some remarks are given in Section 4.

Theory and Formulations
Before beginning the discussion of the proposed methods in this work, we first briefly introduce the problem addressed in this work, the associated MPIE formulation, and its MoM discretization. In this work, the efficient electromagnetic analysis of RF circuits or devices in a multilayer planar medium is considered, as illustrated in Figure 1, where the material discontinuity lies only in z axis of the Cartesian (x, y, z)-coordinate system. The dielectric layers from up to down are denoted as layer i = 1, 2, ..., N with a permittivity and permeability of ( i , µ i ). Here, we assume the e jωt time dependence, and that the metallic circuit or devices are thin perfectly conducting sheets lying between the interface between different media. Such structures can be found in various modern radio frequency devices such as antennas, filters, and couplers because of the inherent nature of the integration of planar structures. Simplifying the metallization patterns in the planar layered medium as perfectly conducting sheets and following the most popular MPIE formulation [13], the unknown electric current J on the metallic surfaces satisfy the following equation: where S m denotes the surface of S in the m th layer, andn m is the unit vector normal to S m . The magnetic vector potential A mi (r) and scalar potential φ mi (r) in m th layer caused by electrical currents or charges in i th layer are expressed as where q is the surface charge density and ∇ · J = −jωq. The vector potential Green's functionK mi A takes the formK The entries ofK mi A and scalar potential Green's function K mi φ are all Sommerfeld integrals and have a general form of where J 0 (·) is the zero-order Bessel function, and ρ is the horizontal distance between source point r and field point r. In (5),K mi is denoted as the spectral domain Green's function, and is calculated by using the transmission line theory [13].
Applying the Galerkin testing procedure, the MPIE can be discretized into a matrix equation ZI = V, V is the right-hand-side vector based on the excitation on circuit, and I is the unknown vector relating to the unknown currents by basis functions. The entries of impedance matrix Z are given by where g m and g n are the basis and testing functions defined in mesh unit s m and s n , respectively.

Proposed Adaptive RFFT
In the previous studies, the RFF methods are usually first performed in the spectral domain to approximate the SI integrands, i.e.,K mi in Equation (5), then the closed-form expressions using integral identities are obtained. Unfortunately, spatial accuracy of such treatment is somehow difficult to evaluate or control in practical applications [25].
To achieve an accuracy-controllable closed-form approximation, we considered the rational fitting function (RFF) directly in the spatial domain. After extracting poles and quasi-static terms in original SIs, for a given z and z, the remaining SI K mi p is approximated by a piecewise rational function, and given as and ρ 0 is a threshold horizontal distance separating the singular and regular behavior of SIs, and it takes the value of one dielectric wavelength from our empirical studies. The reason for choosing RFF here is based on the following considerations. First, the RF approximation is well known for producing an excellent local approximation to arbitrary functions. Secondly, evaluating rational functions is generally more efficient than the existing closed-form approximation scheme to a Sommerfeld integral that contains the Hankel function or other special functions. In addition, we find that RFF expression facilitates the efficient evaluation of MoM impedance entries, as explained later in this paper.
Although an unified RFF model can be established to approximate K mi p with satisfying accuracy, we found that the piecewise approximation in (7) is much faster in terms of numerical calculation. As the metallization pattern lies only on the interface between different dielectric layers, only a finite amount of discrete z and z values need to be evaluated in MoM simulations.
Both K mi,near pr and K mi,far pr in (7) are written by the rational function as follows: where P n is a polynomial function with order n. Here, superscript in the above notation is omitted for simplicity. To find the lowest order of polynomials and their associated unknown coefficients in (8) with a preset accuracy, an iterative procedure is illustrated in Algorithm 1, shown below. The fundamental idea here is to use two RFF models to approximate the SI in the spatial domain, and the additional sampling point ρ next for improving the approximation accuracy is determined by the discrepancy between these two models This iterative procedure continues until the relative error between these two models is lower than the preset threshold ε 0 . The proposed algorithm is summarized in Algorithm 1, and the rational function orders updating paths for the two models are presented in Figure 2. Clearly, the order updating paths shown in this diagram cover larger subsets of m, n combinations and hence yields a better fitting performance. In this algorithm, the Bulirsch-Stoer (BS) algorithm [26] is recommended for the robust and efficient evaluating or updating of unknown coefficients in an RFF model (8). Combing the Richardson extrapolation and modified midpoint method, the BS algorithm finds the numerical solution of RFF or ordinary differential equations with a high accuracy and relatively low computational complexity, and its numerical implementation can be found in the well-known literature (chapter 3.4) [27].

Algorithm 1: Adaptive RFF for SIs
input : preset accuracy ε 0 , maximum iteration number I itr , spatial region z, z , and ρ ∈ [ρ 1 , ρ 2 ] output : RFF model K pr calculate K mi p with randomly selected samples evaluate model 1 K 1 pr with order m 1 ,n 1 evaluate model 2 K 2 pr with order m 2 ,n 2 I cnt = 0 while ||K 1 pr − K 2 pr ||/||K 1 pr || ≥ ε 0 and I cnt < I itr do find next sample point ρ next based on (7) evaluate K mi p (ρ next ) using numerical integration I cnt = I cnt + 1 increase m 1 , n 1 and m 2 , n 2 according to the diagram shown in Figure  To verify the performance of the proposed adaptive RFFT, the CPU time to obtain approximation parameters, SIs evaluation time, and the approximation accuracy are presented in Tables 1 and 2 for two configurations of layered media. Although there are more advanced DCIM schemes available [25], here we use the classic two-level DCIM scheme [20] as the benchmark for comparison as the resultant number of complex images is moderate and approximation accuracy is decent for the examined spatial region after removing quasi-static terms. The approximation accuracy here is calculated by comparing the results between closed-form schemes and the numerical integration method [17]. In Table 1, the SIs in K 11 xx and K 11 φ are evaluated and compared for a microstrip structure with a dielectric constant of 9.6 and a thickness of 0.00254 λ 0 , while in Table 2, a two-layer PEC-backed dielectric substrate with a thickness of 0.05 λ 0 , 0.1 λ 0 is considered, and their dielectric constants are 10 − 0.1 j and 2.2 − 0.2 j, respectively. From these two tables, we can see that the calculation time to establish closed-form expressions is a few times larger than the counterpart of DCIM. Nevertheless, it is still a neglectable portion compared to the numerical evaluation of SIs. The typical convergence history of the proposed adaptive RFFT is shown in Figure 3, where the preset accuracy ε 0 is set to be 0.0001. Although the orders of polynomials in RFFs are comparable to the number of complex images in two-level DCIM, the proposed method provides a speed increase of ten to twenty times, and more importantly it has a higher approximation accuracy.   Relative error 1 0 < ρ < λ 2.6 × 10 −3 3.8 × 10 −3 2.7 × 10 −3 3.8 × 10 −3 λ < ρ < 2λ 4.9 × 10 −4 6.9 × 10 −6 7.2 × 10 −3 8.9 × 10 −5 1 Compared to numerical integration with Euler transformation [17].

Singularity Treatment
Besides the efficient SI evaluation, the fast calculation of the impedance matrix elements based on the proposed closed-form expression is also considered. For the sake of efficient calculation and convenience of singularity analysis, by finding the roots of the denominator, the rational function expression in (7) is rewritten as the sum of the simple rational functions that follow: Substituting (10) into (6), there are two elementary integrals to be evaluated where r 0 is the free vertex of the RWG basis function. Clearly, the magnitude b i is nearly zero when the associated SI is nearly singular, hence the numerical integration of (6) can be time-consuming when the horizontal distance ρ between r and r is close to zero. To overcome this difficulty in numerical integration, we consider the Duffy transformation technique [28]. Let s n be a triangle associated with a basis function g n with tree vertices r 1 , r 2 , and r 3 , as shown in Figure 4a, and r p is the vector of field point r project to plane s n . The integration domain s n can be divided into three sub-triangles s 1 , s 2 and s 3 using the addition vertex r p . In the Duffy transformation, take s 1 as an example. It can be mapped into a transformed domain, as shown in Figure 4b,c, and given by where scalars α and β are In this way, the integrals in (11) can be rewritten as the summation of three singularityfree integrals on s 1 , s 2 and s 3 , i.e., I 1 = I s 1 1 + I s 2 1 + I s 3 1 and I 2 = I s 1 2 + I s 2 2 + I s 3 2 . For the integral on s 1 we have  Observing the integrals on sub-triangles shown in (14), we find that only three types of integration are required to be evaluated here.
Fortunately, analytical expressions are available for double integral (15b) and (15c) and given in Appendix A. The integral in (15a) can be converted to a one-dimensional integral by eliminating u The first two terms in (16) also have analytical expressions and are presented in Appendix A. The third term in (16) generally does not have an analytical expression, depends on the position of r p , and the evaluation strategy may differ.

•
If r p lies on the direction vector l 1 , i.e., l 1 × l 2 = 0, we can define scalar p l 2 /l 1 , and the integral I (1) can be analytically evaluated. The detailed expressions given in Appendix B depend on the value of p.
• If r p is far from vector l 1 , the above integral has to be evaluated numerically using adaptive Gaussian integration. • If r p is close to vector l 1 but not lies on it, to avoid the numerical integration for almost singular integrand as |l 2 + vl 1 | approx to zero, here we use the approximated condition then derive the approximated analytical expression for I (1) as shown in Appendix B.
In addition, for the special case where r p coincide with r 1 or r 2 , the analytical expressions for above three integrals in (15) are derived in Appendix C.
To verify the validity of the semi-analytical strategy to evaluate integral (15a), here we examined the performance of the proposed method in two extreme cases. As shown in Figure 5, the required number of Gaussian quadrature points to achieve the relative error is less than 10 −5 . In this figure, lengths are normalized by the average length of the source triangle edge length, and the triangle edges are illustrated by black lines in xy−plane. The vector r p relating to the field point varies in the range x ∈ [−1, 1] and y ∈ [−1, 1], and the number of required points are represented by the color of each pixel. Obviously, around two or three Gaussian points are sufficient to achieve good accuracy when r p is far from the source triangle, even though b i is quite small. More Gaussian points are required when r p is located near the triangle edges. Nevertheless, the use of approximation (17) makes the maximum number of Gaussian points still below 20.

Numerical Validation
In this section, based on the MoM analysis of planar structures in microwave engineering, numerical examples to demonstrate the calculation efficiency of the proposed method are presented through comparisons.

A Microstrip Low-Pass Filter
We first considered an MoM simulation of a microstrip-design low-pass filter to verify the correctness and efficiency of the proposed SI evaluation technique and the singularity treatment strategy. The microstrip structure is identical to the aforementioned one whose SIs evaluation performance is compared in Table 1. The metallization pattern is printed on top of the substrate surface and shown in Figure 6, with two ports at the front and end of this structure, the simulated S11 and S12 parameters are presented in Figure 7a,b, respectively. This low-pass filter has a length of around 46 mm, a width of around 16.93 mm, and the dielectric substrate length is 0.254 mm. Due to the reciprocity of this device, S22 and S21 parameters are omitted here for simplicity. In addition, here we also presented the S-parameter results evaluated by MoM using DCIM and by a commercial MoM tool FEKO in the same plot. In Figure 7, a good agreement between these three methods can be observed, except for S12 at the stopband where the value is lower than −50 dB, which is extremely sensitive to round-off errors in numerical calculation. To verify this point, at 2.9 GHz where the S12 has around 10 dB discrepancy between different methods, the surface current intensity of the LPF circuit by excitation port 2 with a voltage source is presented in Figure 8. The current intensity near port 1 almost vanished to zero in all three simulations, as the operation frequency is located in the stopband of this LPF. The electric current intensity plots are also almost identical and indicate the correctness and accuracy of the proposed method, despite the fact that the S21 error is large here.   In addition to the S-parameter curves, the MoM calculation CPU time using different SI evaluation schemes is presented in Table 3 for comparison. In particular, we also provide the CPU time for filling the impedance matrix in the MoM procedure in the second and fourth columns of this table. Clearly, without losing simulation accuracy, compared to twolevel DCIM, the proposed method has an increase in speed of approximately in terms of total calculation time and a larger increase in speed in terms of filling the impedance matrix.

A Stacked Antenna
The following example considers a stacked planar antenna simulation using MoM and layered Green's function. The layered medium configuration is shown in Figure 9a, where an additional infinite thin metallization (shown in blue) sheet is inserted into the dielectric substrate and connected to the feeding pin; we treated this layered medium as a four-layered structure. The layout of this antenna design is shown in Figure 9b, where the four rectangular sheets shown in red in this figure are printed on top of the dielectric substrate and act as radiation elements in this antenna design so that an omnidirectional radiation pattern can be achieved. Using the microstrip edge port as excitation, the input impedance and radiation pattern of this antenna with infinitely large ground and substrate can be efficiently simulated. Similar to the previous numerical example, the input impedance of this antenna is calculated by MoM with different SI evaluation schemes, i.e., the proposed methods and the conventional two-level DCIM. The input impedance values in the frequency band from 1.6 GHz to 3.2 GHz with 50 MHz interval are presented in Figure 10, where the resistance and reactance are plotted in Figure 10a,b. Again, a very good agreement between the different methods in calculating the input impedance of this antenna can be observed in the examined frequency band. The CPU time comparisons for total computational time and time taken for filling the impedance matrix are also provided in Table 3. Again, an approximately seven times increase in speed of can be seen in this example.

Conclusions
Aiming for a more efficient MoM simulation for planar circuits and antennas, the fast closed-form evaluation of SIs in the layered medium is considered in this work. Using the rational function approximation to SIs directly in the spatial domain, the proposed adaptive RFFT achieves fast evaluation of SIs with a preset approximation accuracy. In addition, the semi-analytical expressions for the near-singular SIs in MoM implementation are derived and verified based on the proposed RFFT. Representative MoM simulations for planar devices in microwave engineering are studied to demonstrate the correctness and efficiency of the proposed methods and an increase in speed of several times can be achieved when compared to classic DCIM. As an adaptive fitting technique for SIs, the proposed method is robust and flexible in the approximation of Green's function in an arbitrary layered medium without considering the vertical location of source and field points. Hence, it can be integrated into existing MoM tools easily to accelerate the electromagnetic analysis of planar circuits and antennas.

Appendix B
Substituting l 2 = pl 1 into the integrand of (15a), we have then analytical expressions for I (1) can be derived base on the value of p. If p > 0 or p < −1, then

Appendix C
If r p coincides with r 1 or r 2 , the analytical expressions in Appendices A and B are rewritten as follows and If r p = r 1 , then if r p = r 2 ,