Optimized Analytical Computation of Partial Elements Using a Retarded Taylor Series Expansion

The aim of this article is to efficiently and accurately calculate the integrals of the full-wave (FW) partial element equivalent circuit (PEEC) method. The accuracy of the analytical formulas calculated by the standard precision can be compromised when using nonuniform mesh to properly model the high-frequency effects. The numerical errors can be avoided by using a high-precision arithmetic, i.e., higher number of digits, however, at the expense of significantly higher computation time. This article presents an analytical approach for calculating the FW-PEEC interaction integrals of two elementary volumes/surfaces based on the Taylor expansion, which allows a high computational speed preserving the accuracy with a relative error of less than 0.1%. The proposed solution is verified compared to the high-precision arithmetic and the standard Gaussian integration for two examples of strip lines. Moreover, it is shown that the accuracy of FW-PEEC integrals can affect the convergence of an iterative PEEC matrix solver.


I. INTRODUCTION
E LECTROMAGNETICS modeling has received increasing interest in recent decades due to its ability to predict the electromagnetic performance of an electronic system in an early design phase much before the realization of a physical prototype. For this purpose, consequently, the numerical solution of Maxwell's equations has acquired increasing importance and has become a powerful tool for designers. Many methods have developed over the years, such as the finite-difference time-domain technique [1], the finite-element method [2], the method of moments [3], and the partial element equivalent circuit (PEEC) method [4], [5], [6], [7].
Nowadays, the most challenging objective is the development of a method that is accurate from the very low frequencies, including dc, to the highest frequencies that can easily reach Manuscript  the tens of gigahertz. Among the methods listed above, the PEEC method is different as it can transform an electromagnetic problem into a circuit model, in which the circuit parameters model the dissipative phenomena and electromagnetic coupling. The latter are known as partial inductances and potential coefficients [7]. Since the PEEC method is based on the principle of volume equivalence, the currents and charges are assumed to radiate in the background medium, and therefore, free-space Green's function has to be considered. If the quasi-static assumption can be done, for electrically small problems (e.g., involving low frequencies and/or geometrically small structures), Green's function simplifies and the partial elements become frequency independent. If an orthogonal mesh is used, analytical formulas are available for both partial inductances L p and coefficients of potential P [7]. However, the existing analytical formulas are affected by significant numerical errors for certain PEEC structural mesh necessary to model the skin and proximity effects with higher accuracy. In [8], a systematic strategy to select a proper analytical formula depending on the dimensions and positions of two elementary volumes is proposed for the accurate computation of partial inductances under the quasi-static hypothesis.
When the quasi-static hypothesis is not valid longer, full-wave Green's function has to be considered. In this case, the calculation of the partial elements modeling the magnetic and electric field couplings, L p and P , respectively, is normally carried out by using numerical integration, which, however, is very computationally expensive. This limitation can be overcome by using Taylor expansion, as presented in [9]. Nevertheless, similar numerical issues can also be found in the full-wave calculation of the partial elements based on Taylor expansion. The aim of this article is to define a multifunction strategy that allows preserving the accuracy in the calculation of the coefficients of Taylor expansion of the partial elements as the frequency, size, and distance between the domains vary. An orthogonal mesh is assumed in this article.
The rest of this article is organized as follows. Section II summarizes the problems in the computation of integrals involved in the mutual partial inductance between two elementary parallelepipeds and the conditions to be used to set a correct strategy preserving the accuracy of the quasi-static partial inductance. The details of the computation of the Taylor series expansion coefficients, for different domains, e.g., volumes, surfaces, and lines, are given in Section III along with the strategy proposed to switch from one to another analytical formula in the different scenarios. The range of applicability is described in Section IV along with an example. Finally, Section V concludes this article.

II. MUTUAL PARTIAL INDUCTANCE COMPUTATION
In the quasi-static PEEC method, the magnetic field coupling between two elementary volumes i and j carrying uniform currents is described by the partial inductance, which requires the computation of the following double-folded volume integral: where R = ||R i − R j || is the distance between two points inside the two cells, and S i and S j are the cross sections normal to the current directions of two volumes i and j, respectively. When the quasi-static hypothesis is not satisfied, the full-wave computation is required. In this case, the coefficients (1) must consider the retardation through the exponential term, as shown in the following equation: where β = 2πf √ μ 0 0 . In [9], it is shown how it is possible to compute the coefficients (2) by resorting to the Taylor series expansion of the exponential term. In particular, the most efficient expansions reported in [9] work as described in the following. The double-folded integral (2) is rewritten into an equivalent form as follows: where R cc is the center-to-center distance between the two volumes. Then, the term e −jβR cc is approximated with its Taylor expansion, and it is moved out the integral, leading to By performing some trivial algebraic manipulations, it follows that where Integrals in (6) are frequency independent, and hence, they are computed only once and then used to compute a full-wave coefficient at any frequency. The analytical formulas for integrals (6), in the case of orthogonal volumes, are provided in [9] and [10]. The case of nonorthogonal volumes can be handled in two ways, in addition to the brute-force numerical integration.
1) The coefficients of Taylor expansion are computed once at the beginning and reused at different frequencies. 2) As a second option, a fine orthogonal mesh is used to approximate the nonorthogonal objects. Thus, the number of elementary unknowns can be very high. Such an approach is usually used in conjunction with an iterative method, in which the matrix-vector products are accelerated, e.g., via fast Fourier transform (FFT) [11], [12], [13], [14].

A. Problems in the Analytical Evaluation of Integrals (6)
Let us consider the following example, also presented in [8] for the quasi-static case. In particular, it considers two volumetric cells arranged as follows.

3) The distance between the cells in x, y, and z directions
is Δx = s x1 · h, where h ranges from 10 −2 to 10 5 , Δy = 0, and Δz = 0, i.e., the minimum distance between cells R min varies between 0.1 μm and 1 m. For all these geometrical configurations, the computation of integrals in (6) is performed in the standard double-precision arithmetic by using the analytical solution provided in [9] and [10], while the reference values are computed using the Symbolic MATLAB Toolbox [15] (with at least 32 digits). The result of such test is reported in Fig. 1 .
In particular, for orthogonal volumes, when integrals (6) are computed in double precision, some numerical errors occur [16], [17]. More recently, Kovacevic-Badstuebner et al. [8] have proposed a technique that solves definitively the numerical problems related to the analytical computation of I The behavior of this approach has been summarized in Appendix A for completeness.

III. IMPROVED INTEGRAL (6) EVALUATION
The novelty of this article is to extend the technique presented in [8] for the computation of integrals I (6). To this aim, the closed-form solutions for these integrals when suppressing one or more dimensions at the time are needed. Similarly to the definitions made in (13), for I  p−p , respectively. The closed-form solutions for all these integrals, that represent the major novelty introduced in this article, are provided in Appendix B for the sake of readability. v−v , that can be found in [8] and [17], defined in (13), has a similar counterpart in the solution of the integrals for I (o2) v−v , given in [9] and Figs. 12-15, in the sense that there are similar summations of terms having the same exponential values. This is also confirmed by the similar error map that can be observed in Fig. 1(a) and (b). Moreover, condition 1 (14a) in this case is the following: 3) Finally, for integral I Fig. 1(c), it can be observed that more relaxed thresholds of suppression can be used since its computation fails for higher values of sizes and distances. Experimentally, by running a very large set of test cases, the best thresholds found for the suppression strategy are the following: 1 = 10 −4 , 2 = 3 × 10 −4 , 3 = 10 −3 , and 4 = 3 × 10 −1 . Furthermore, condition 1 (14a) in this case reads Hence, by applying the proposed technique to the same example of Section II-A, with the reported choice of thresholds, we obtain the results summarized in Fig. 2, where, also in this case, the reference values are computed using the Symbolic MATLAB Toolbox [15] (with at least 32 digits).
As seen, applying the proposed strategy for the computation of integrals (6) leads to very low errors. This allows us to conclude that the evaluation of (5) is affected now only by the Taylor expansion truncation of the exponential term and not by the digits of precision used in the evaluation of integrals (6).

A. Evaluation of the Range of Applicability
The range of applicability of the proposed formulas when using the double-or quadruple-precision arithmetic is analyzed by calculating the integral (5) for different volumes sizes and distances. In particular, the distance between two cubic volumes is varied in the range from 0 to 1 m leading to 30 geometrical configurations. Such distance is calculated between the two nearest edges of the two volumes. In addition, the tests are performed for the edge size of the two volumes ranging from 10 μm to 1 mm. In the presented analysis, the integrals are computed based on five different approaches.
1) Numerical: (2) is calculated through the adaptive numerical integration of MATLAB with an error threshold set to 10 −9 . In particular, the integration is performed with the Gauss-Kronrod quadrature formula [18], where the adaptive strategy is implemented with seven evaluation Gauss points and 15 evaluation Kronrod points.
2) Double Taylor: The integral I (FW) v−v in computed as reported in (5) when the following condition is satisfied: In (10), c 0 denotes the speed of light in free space, f denotes the frequency, and s max denotes the largest value between the lengths, the widths, and the thickness of the two bars, for which I (FW) v−v is computed. Clearly, the integral I (9) is computed by adopting the suppression strategy given in [8]. Finally, the validity of the approximation (9), along with the condition (10), has been proven by running a large number of numerical tests. The CPU times required for these methods are summarized in Table I along with the details of the frequency range and frequency samples for each case. Their computation has been carried out in MATLAB on a PC equipped with two physical processors operating at 3.46 GHz. From Table I, it is evident that the use of Taylor expansion with the double-precision arithmetic is extremely advantageous in terms of performances, while an accurate numerical integration is prohibitive.The accuracy of the Double Taylor, Quad Taylor, Suppression Taylor, and Proposed Taylor approaches is defined by a relative error computed as where the same notation of (10) has been used. As can be clearly seen from Figs. 3(d), 4(d), and 5(d), the proposed strategy always guarantees an error below 0.1%. Furthermore, it is evident that the Quad Taylor approach can be used as a reference method if the λ/30 criterion (12) is satisfied.
The following examples demonstrate how the calculation of L p and P matrices influences the calculated vectors of currents and potentials. These coefficients are computed in C++, and more specifically, they are computed in parallel on seven threads in order to accelerate all the simulations.

B. Example 1: Coplanar Stripline Example
As the first example, a coplanar stripline structure described in Fig. 6 is modeled. The conductive ground, the lines, and the plate above the structure are made of copper. In addition, the    lines are embedded into a dielectric with relative permittivity r = 4. All the geometrical details are given in Fig. 6.
The PEEC analysis has been performed from 100 MHz to 5 GHz by using a nonuniform mesh (with minimum mesh size of 0.4 μm and maximum mesh size of 1.5 mm in each direction) resulting in 25 680 PEEC inductive branches and 8688 capacitive surfaces.
The CPU time and relative error for the double-folded integrals used to fill L p and P matrices (2)-(6) are summarized in Table II showing that the use of quadruple-precision computation leads to a high computation time, while the proposed solution shows a drastic reduction in the computation times. In addition, the Proposed Taylor approach is also faster than the Double Taylor, since, when several integration dimensions are suppressed, the derived closed formula evaluation requires less operations than the standard volume-volume integral solution.
The choice of the integrals used to fill L p and P matrices by the Proposed Taylor method, between the developed closed formulas given in Appendix B, is summarized in Table III. In  TABLE II  CPU TIME AND RELATIVE ERROR FOR THE DOUBLE-FOLDED INTEGRALS USED TO FILL L p AND P MATRICES (2)-(6) FOR THE COPLANAR STRIPLINE EXAMPLE this example, the approximation (9) is practically unused since less than 1% of coefficients are corrected with it, and hence, the accuracy is reached by only applying the suppression strategy introduced in Section III.
In the next step, the accuracy in the computation of currents and voltages between a node of the structures and all the other nodes is evaluated for the Double Taylor and Proposed Taylor approaches. The errors are computed by using (11) and presented in Fig. 7. In this case, Fig. 7 demonstrates how the errors in L p matrix influence the calculation of the current distribution, while the error on voltage distribution is less influent for the same reasoning given in [8].
In any case, the errors of the current distributions have a significant impact on the postprocessing step to calculate electromagnetic fields.

C. Example 2: Microstrip Line Example
As the second example, a microstrip line structure described in Fig. 8 is modeled. The conductors are made of copper (conductivity σ = 5.8 × 10 7 ), while the dielectric substrate is made of the standard FR4 material ( r = 4.2 and tanδ = 0.02).
The PEEC analysis has been performed from 500 MHz to 5 GHz by using an orthogonal mesh approach, in which all the structures are discretized into voxels (that are practically rectangular bars) having all the same sizes. The sizes of the voxels along the x, y, and z axes are s x = 84μm, s y = 384μm, Fig. 7.
-errors of the current and voltage vectors for the coplanar stripline example. and s z = 10μm, respectively. In particular, the mesh is achieved through a rectangular grid having N x = 20, N y = 140, and N z = 26, where N x , N y , and N z are the number of voxels along the x, y, and z axes, respectively, resulting in 186 288 inductive branches, 18 864 capacitive surfaces, and 61 600 nodes. Since the mesh is made by a regular 3-D grid in which all the voxels have the same size, we can use an iterative solver in which the matrix-vector products can be accelerated through the FFT [11], [12], [13], [14]. In particular, we use the generalized minimal residual method (GMRES) iterative solver, as described in detail in [11], with the threshold of convergence set to 10 −4 and with the maximum number of iterations set to 3000. Since the matrix-vector products are performed via FFT, it follows that just one matrix row needs to be computed and stored to build the circulant tensors required for the matrix-vector products involving the matrices L p and P [12]. The CPU time and relative error for the double-folded integrals used to fill the rows of L p and P matrices (2)-(6) are summarized in Table IV showing that the use of the quadrupleprecision computation leads to a high computation time, while the double-precision computations are extremely fast (since, in this example, just few matrix rows are computed to build the circulant tensors).
The choice of the integrals used to fill L p and P matrices by the Proposed Taylor method is summarized in Table V, where the notation is the same as used for Table III. Also, in this example, the approximation (9) is practically unused since less than 1% of coefficients are fixed by (9).
In addition, the scattering parameters are shown in Fig. 9. As can be clearly seen, the proposed solution shows a very good agreement with the Quad Taylor method, while the Double Taylor method leads to less accurate results at high frequency.
Finally, the number of GMRES iterations is shown in Fig. 10, in which it can be noted that the Quad Taylor method and the proposed method perform the same number of iterations, while the Double Taylor method exhibits a higher number of iterations reaching also the limit of 3000 iterations at high frequency without satisfying the 10 −4 threshold used for the convergence criterion.  Clearly, from this results, it is evident that an inaccurate subset of L p and P coefficients has a negative impact on the convergence of an iterative method, which worsens with increasing frequency. This is a strong limitation because the use of an iterative method is strictly needed to solve problems with a large number of unknowns that appear for practical examples of 3-D circuit layouts.

V. CONCLUSION
In this article, the accuracy of PEEC electromagnetic models at high frequency has been significantly improved by exploiting and extending strategy for selecting the right analytical formula depending on the dimensions and positions of the two PEEC elementary volumes/surfaces. This allows the computation of the interaction integrals with a relative error of less than 1%, avoiding the numerical integration and the quadruple-precision arithmetic. A provided set of analytical formulas, required by the technique, can also be useful to any other numerical method requiring double-folded volume or surface integrations of freespace Green's function over rectangular domains. Finally, the numerical results show that significant speedup is achieved while preserving the accuracy of the solution. In addition, the last example also highlights that inaccurate mutual coefficient computation can compromise the convergence of an iterative solver required to solve large computational problems.
In [8], starting from the six-variable-folded integral I The algorithm needed to switch between the various formulas in (13) is given in [8] and reported in the following for completeness. By defining the following -condition: condition A = condition 1 and (condition 2 or condition 3 ) the algorithm reported in Fig. 11 can be applied to understand what dimension must be suppressed for a volume.  12-15. It is important to underline that only the closed-form solutions are given, while all the algebraic steps performed to achieve them are skipped due to the lack of space.