The Partial Elements Equivalent Circuit Method: The State of the Art

This year marks about half a century since the birth of the technique known as the partial element equivalent circuit modeling approach. This method was initially conceived to model the behavior of interconnect-type problems for computer-integrated circuits. An important industrial requirement was the computation of general inductances in integrated circuits and packages. Since then, the advances in methods and applications made it suitable for modeling a large class of electromagnetic problems, especially in the electromagnetic compatibility (EMC)/signal and power integrity (SI/PI) areas. The purpose of this article is to present an overview of all aspects of the method, from its beginning to the present day, with special attention to the developments that have made it suitable for EMC/SI/PI problems.

The Partial Elements Equivalent Circuit Method: The State of the Art Giulio Antonini , Senior Member, IEEE, Albert E. Ruehli, Life Fellow, IEEE, Daniele Romano , and Fabrizio Loreto Abstract-This year marks about half a century since the birth of the technique known as the partial element equivalent circuit modeling approach.This method was initially conceived to model the behavior of interconnect-type problems for computer-integrated circuits.An important industrial requirement was the computation of general inductances in integrated circuits and packages.Since then, the advances in methods and applications made it suitable for modeling a large class of electromagnetic problems, especially in the electromagnetic compatibility (EMC)/signal and power integrity (SI/PI) areas.The purpose of this article is to present an overview of all aspects of the method, from its beginning to the present day, with special attention to the developments that have made it suitable for EMC/SI/PI problems.

I. INTRODUCTION
T HE numerical solutions of Maxwell's equations were sig- nificantly enhanced between the end of the '60s to the '70s.Until then, Maxwell's equations were used to obtain analytical solutions.The numerical solutions were rapidly advanced with the first more powerful computers.Then, within a few years, we witnessed the birth of several fundamental numerical approaches which are in use today.Among them, we list the finite difference time domain (FDTD) method [1], the finite element method (FEM) [2], the method of moments (MoM) [3], [4], and the finite integration technique (FIT) [5].
The partial element equivalent circuit (PEEC) method also was conceived in the early '70s.The IBM mainframe computers which were the most advanced at that time, required modeling with more complicated inductive interconnect topologies.Specifically, the electronic system problem required relatively complex numerical modeling of inductive voltage drop in the Giulio Antonini, Daniele Romano, and Fabrizio Loreto are with the UAq EMC Laboratory, Department of Industrial and Information Engineering and Economics, University of L'Aquila, 67100 L'Aquila, Italy (e-mail: giulio.antonini@univaq.it;daniele.romano@univaq.it;fabrizio.loreto@graduate.univaq.it).
Albert E. Ruehli is with the EMC Laboratory, Missouri University of Science and Technology, Rolla, MO 65409 USA (e-mail: albert.ruehli@gmail.com).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TEMC.2023.3302700.
Digital Object Identifier 10.1109/TEMC.2023.3302700power and ground voltages connected to ICs.The concept of partial inductance-based complex circuit analysis was introduced in [6], where it was shown how to solve circuits based on multiple parts of conductors of finite length.This made complex inductive structure models solvable.Integrated circuit design and packaging also required the calculation with capacitive elements.The capacitance computation for 3-D conductors was presented in [7] using the concept of coefficient of potential.The two concepts of partial inductance and coefficient of potential were integrated into the integral equationbased formulation published in [8], which can be considered the foundation of the PEEC method.More than 20 years after the introduction of PEEC, the applications were mostly in the modeling of interconnections in printed circuit boards [9], [10] and packaged electronics [11], including dielectric as well [12], [13].
New applications of PEEC, other than that of interconnections date back to 1998.Specifically, new problems like the modeling of lightning protection systems and its coupling to coaxial cables were solved [14].
The number of papers on the PEEC method has been constantly increasing over the years, as shown in Fig. 1.
An important added value of the PEEC method is represented by the fact that, by providing a circuit interpretation of Maxwell's equations, it is totally compatible with circuit solvers, such as SPICE [61], [62], [63].A SPICE solver modified to take the propagation into account has been presented in [64].It is also worth mentioning that the PEEC method is listed among the numerical methods to solve Maxwell's equations listed in the IEEE Standard for Validation of Computational Electromagnetics Computer Modeling and Simulations (IEEE Std 1597.[65]. In this work, we retrace the 50 years' history of the method and the evolution that led it to be considered today among the most suitable ones to be used in the field of electromagnetic compatibility (EMC), signal, and power integrity (SI/PI).

II. BASIC PEEC FORMULATION
The PEEC method is based on the electric field integral equation (EFIE) and the continuity equation (CE).Differently from the MoM (Z-MoM) [3], the continuity equation is kept separate from the EFIE such that the solution includes both currents as well as potentials or charges.The basic derivation of PEEC equations is carried out in this section for a conductor in the Laplace domain.The modeling of different materials is described in Section V.
The electric field at a point r in a conductor is where σ is the conductor conductivity, E inc (r, s) is the external electric field impressed at point r at a Laplace frequency s.A(r, s) and Φ(r, s) are the magnetic vector and electric scalar potentials [66], respectively.They read and where τ = |r − r |/c 0 is the speed of light in the background medium.If we replace the vector potential A(r, s) in (1), we get which has J(r, s) and Φ(r, s) to be determined.In addition to the EFIE, the continuity equation is used that makes the problem well-posed.If the charges are assumed to be located on the surface of volumes, which is a reasonable hypothesis for good conductors, then the continuity equation has to be written as In this case, the electric scalar potential has to be redefined as Then, current and charge densities J(r, s) and q(r, s) are expanded in terms of pertinent basis functions b ∈ R 3 and p ∈ R where I n (s) and Q j (s) are the expansion weights that must be determined at each angular frequency s, N v and N s represent the number of the elementary volume and surface regions, respectively.Then, expansions (8) are substituted into (2)-( 4), and ( 6) and the so-called Galerkin's testing or weighting process [67] is used to generate a system of equations for the unknowns weights I n (s), n = 1, . . ., N v and Q m (s), m = 1, . . ., N s by enforcing the residuals of (3)-( 5) to be orthogonal to a set of weighting functions that are chosen coincidently with the basis functions.Then, two inner products are defined They are used to average ( 4) and ( 6) on elementary volumes and (7) on elementary surfaces.Importantly, the averaging process leads to double integration in the partial elements, e.g., (21) and (28), making symmetric the corresponding matrices.
The Galerkin's method applied to the EFIE (4) yields, in a matrix form where vectors Φ(s) and I collect the potentials to infinity and the currents flowing through the volumes, respectively, the matrix A is the usual circuit incidence matrix, R is a matrix accounting for conductor losses, L p is the partial inductance matrix that models the magnetic field/inductive coupling due to the currents, and vector V s (s) represents the voltages induced on elementary volumes by the incident electric field E inc (r, s) [68].
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Equation ( 10) can be regarded as the Kirchhoff voltage law (KVL) enforced to loops constituted by each ohmic-inductive branch and its closures at infinity.Applying the Galerkin's method to ( 6) and ( 7) yields and where matrix Y(s) models lumped elements and vector I s (s) accounts for eventual lumped current sources.Equation ( 12) is well suited to identify an equivalent circuit for the displacement currents in the background medium.Indeed, they read Some simple algebraic manipulations allows to rewrite (12) as where Ĩc = TI c and matrices D and T are N s denotes the number of elementary surfaces where the charge is localized.Equation ( 14) admits the circuit synthesis sketched in Fig. 2 for the case of N s = 3 patches.
Using (12) to eliminate the charges from (11) yields Equation ( 16) can be regarded as the Kirchhoff current law (KCL) enforced to each node corresponding to elementary volumes .
(17) Equation ( 17) is the modified nodal formulation (MNA) of the PEEC circuit [69] equations.It is also worth noticing that additional lumped elements can be easily added to the model just by stamping them into the MNA matrix [57].We note that the inversion of matrix P(s) in ( 17) can be time-consuming when the problem size is large.However, this matrix inversion can be avoided.As an example, the admittance part can be premultiplied by P(s).

III. MESH GENERATION AND BASIS FUNCTIONS
The choice of the basis functions is directly related to the mesh.Over the years, several different meshing techniques have been considered.

A. Orthogonal Meshing
Initially, the meshing for most PEEC solvers was based on orthogonal Manhattan meshing only, and piecewise constant basis functions were used exclusively [8].
Hence, the basis function for the current density is chosen as where tn is the tangential unit vector for the current direction in volume cell V n .With such a choice of the basis function, the corresponding weight represents the current flowing in the volume V n with orientation tn .Assuming that the free and bound charge densities are located on the surface of conductors, the basis functions used to expand the charge densities are chosen as follows: This choice implies that the corresponding weight Q m represents the charge on cell m and is assumed to be uniformly distributed on the surface.Fig. 3 shows a simple PEEC ohmic-inductive cell structure for a metal sheet.It represents the layout of the nodes with resistive and inductive cells.Please note that in this example, half-width cells are used at the edge of the conductors so that the cell nodes can be connected at the edge, by assembling the substructures which may consist of multiple sheets.We notice that the inductive and resistive cells are rectangular of zero or a finite thickness.An example of uniform mesh for a finite thickness cell is shown in Fig. 4.
Today, for most problems, the high-frequency skin effect can be included in the model by subdividing the cross-section such that the redistribution of the high-frequency current to the surface can be modeled.Thus, when skin-effect has to be caught at high frequencies, the standard PEEC formulation requires a fine mesh of the volume.Such a model can be very inefficient if uniform subsections for the current flow are used, as shown in Fig. 4.However, several different solutions can be used like highly nonuniform meshing [70].Fig. 5 shows an example of a nonuniform meshing of the cross-section of a conductor.
Theoretically, a surface integral formulation of the PEEC method (e.g., the PMCHWT method) can give a rigorous evaluation of the skin-effect losses for piecewise homogeneous conductors [71], [72].However, problems may occur in the numerical implementation due to the high contrast between the conductor and the surrounding dielectric, which makes the resultant linear system highly unbalanced.In [73], the interaction integrals involving the curl of the magnetic and electric vector potentials are computed through the Taylor series expansion of the full-wave Green's function, leading to analytical forms that are rigorously derived, thus reducing the numerical issues.When the skin-effect is well-developed, a valuable alternative option is represented by the use of the surface impedance accounting for the nonuniform distribution of the current density [74], [75].More recently, it has been coupled to voxellized fast Fourier transform (FFT)-based PEEC models [19].

B. Nonorthogonal Meshing
The need to model more complex geometries has led to the development of the PEEC method which is also able to handle problems with non-Manhattan geometry.
An unstructured PEEC formulation by dual discretization for 2-D geometries has been proposed in [85] and [86] and then used in [87] to analyze two-port TEM cells for VHF applications.Then, it has been extended to volume, dielectric, and magnetic media [88], [89], [90], [91].A 3-D PEEC formulation for structured and unstructured PEEC models that is based on the cell method [92] has been presented in [93].It guarantees a circuit interpretation of Maxwell's equations also in magnetic materials and a reduction in the number of degrees of freedom (DoFs) (i.e., required memory) compared to the full-wave unstructured PEEC formulation previously proposed in [89] for electric and magnetic media.
Triangular cells have been applied for a long time in the MoM approach using the Rao-Wilton-Glisson (RWG) basis functions [94], [95].Triangular meshing has also been considered for the PEEC modeling of conductors in [74], [75], [96], [97], and [98].Triangular cells have been applied for a long time in the MoM approach using RWG basis functions [94], [95].Triangular Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
More recently, an isogeometric analysis (IgA) is proposed for the PEEC method for electrostatic problems in [99].It is shown that using the spline-based geometry concepts from IgA allows for extracting circuit elements without an explicit meshing step, with a lower number of DoFs and faster convergence compared to the standard PEEC approach.

IV. COMPUTATION OF THE PARTIAL ELEMENTS
Once this decomposition takes place analyzing a given EM problem, it is straightforward to place these effects in a circuit context where electric and magnetic phenomena are always concentrated and well separated, and where each single EM interaction is described by a "partial element."Therefore, the so-called partial elements are subdivided into two categories.
1) Partial Inductances: Describing the mutual or self effects due to electric currents flowing through volumes.2) Coefficients of Potential: Describing the mutual or self effects due to electric charges localized on surfaces.3) Resistances describing losses.In this section, the partial elements' rigorous mathematical definitions as well as their approximations typically employed in a PEEC model are described and discussed.

A. Partial Inductance for Both Domains
Depending on the highest frequency of interest, the elements are differently formulated.The inductances for quasi-static or low-frequency solutions are where tm and tn represent the unit vectors in the direction of the current flow and S m and S n are the current conductor crosssection areas for the current.Further, V m and V n represent the threefold integrals for the conductor volumes.The inductances are called to be partial since a single partial inductance does not result in a meaningful inductance result.The analytical evaluation of the partial inductances ( 20) is possible for rectangular Manhattan shapes [6].This approach has several advantages since it results in nonsingular partial selfinductance and it can handle large aspect ratio shapes which would be very costly using numerical techniques.
The most complete (full wave, full spectrum) solution requires that the exponential term is included under the integral as where β = ω/c 0 The exponential e −jβR m,n is the factor accounting for the propagation delay effects taking place between the two volumes.Rigorously, the delay exponential function is fully involved in the integration process, since the propagation effects are connected to the whole points in the volumes V m , V n .We observe that using (21) requires the recomputation for the partial inductances at each frequency sample, which is expensive.The numerical evaluation of ( 21) can be performed using numerical integration using methods like the Gauss quadrature.Of course, the computational effort increases linearly with the number of frequency samples required and, as usual, the singularities for the self-term have to be taken care of.Unfortunately, a full wave solution can be obtained only for limited cases of (21).
A direct improvement of ( 20) can be obtained by considering a unique time delay between the volumes, to extract the exponential delay from the integral.Usually, the delay 21), the center-tocenter (CC) approximation of the partial inductance L F W p m,n is obtained In the recent work [100], it has been shown that the Taylor series approximation can be used to compute the full-wave partial elements.Using the Taylor series of the exponential around zero (24) A better approximation for far-apart cells is obtained by taking the center distance R m,n /c 0 delay between the two cell centers outside of the integral and then expanding the exponential e −jβ(R−R m,n ) .This leads to the following approximation of the partial inductance L F W p m,n : (25) For TD analyses, (20) can be used as is since it is frequency independent.The TD counterpart of ( 21) can be recovered via the inverse Fourier transform (IFT) but is computationally expensive since it requires the computation of the partial element at many frequency samples over a wide frequency range.The CC approximation (22) can be easily translated to the time domain as where δ(t) denotes the delta Dirac function.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.More recently, rigorous TD expressions for partial inductances and coefficients of potentials have been proposed in [101], [102], and [103] assuming rectangular Manhattan shapes.An effective way to recover TD partial inductances for nonorthogonal shapes of the volumes is to use the modified numerical inversion of the Laplace transform [104] that uses the Cauchy's theorem.

B. Coefficients of Potential
There is a great similarity between the coefficient of potential partial elements and the partial inductance in Section IV-A.
In the case of good conductors, it is assumed that the charge is on their surface because it moves to the surfaces at a fast rate [57].Also in the case of homogeneous dielectrics, the polarization charge can be assumed to be only on their surface.Under the quasi-static approximation, the basic form of the coefficient of potential is where, in this case, S m,n represents the area surface of the cells.The full-wave coefficient of potential is In lossy dielectrics, such as silicon, due to the reduced value of their conductivity, the charge is not restricted to its surface, e.g., the relaxation time is not as small as in standard conductors.In this case, the coefficients of potential have to be computed as a double-folded volume integral In [105], analytical formulas are presented for (28) under the hypothesis of orthogonal geometries.

V. MODELING OF CONDUCTORS, DIELECTRICS, AND MAGNETIC MATERIALS
Over the years, the PEEC method has had many improvements concerning the possibility of modeling materials of different types.This section briefly revised the models which have been introduced since the beginning.

A. Conductors
In its original 1974 version [8], the PEEC method only considered conductors.Fig. 7 shows an example of the equivalent circuit of an elementary volume of a conductor.Conductive dispersive materials are in use today like graphene.PEEC modeling of graphene interconnects has been presented in [106] and [107].

B. Dielectric Materials
The equivalent circuit model of lossless ideal dielectrics was introduced in 1992 making it possible to model realistic interconnects [12].The polarization of the dielectric is modeled by the so-called excess capacitance which, assuming a Manhattantype mesh, for a cell of length and cross-section S, can be written as where ε r is the relative permittivity of the dielectric.It is to be noticed that polarization currents are bound to the dielectric and generate magnetic fields and interact exactly as electrical currents.The equivalent circuit of an elementary volume of a dielectric is sketched in Fig. 8.With the increase of the frequencies of interest, it has become more and more important to model the dispersive and lossy behavior of dielectrics.PEEC models of dispersive dielectrics of Debye and Lorentz type have been presented in [108].Fig. 9 shows an example of the equivalent circuit of a dielectric elementary volume of Lorentz type.
PEEC models of lossy dielectrics, like silicon, and anisotropic dielectrics have been presented in [105] and [109], respectively.Finite-sized piecewise homogeneous dielectrics easily require massive equivalent circuits.Novel compact models of finite-size dielectrics that are based on the surface equivalence principle and the quasistatic assumption are presented in [110] and [111].

C. Magnetic Materials
The modeling of magnetic materials has been addressed first by considering the magnetization and the constitutive which, being the flux density divergence-less, can be rewritten as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Enforcing (32) in each direction of the Cartesian reference yields an additional set of equations Furthermore, the magnetic vector potential has an additional contribution due to the magnetization resulting in a modified KVL −A T Φ(s) + RI(s) + sL p (s)I(s) + sL m (s)M(s) = V s (s).
(34) Hence, the set of equations to be solved reads This formulation was proposed for the first time in [112].Under the quasi-static hypothesis, rigorous analytical formulas are proposed in [113] for integrals accounting for flux density due to current and magnetization densities for the Manhattan-type discretizations.The formulation presented in [112] has been extended in [114] for modeling 3-D magnetic plates for lowfrequency magnetic shielding problems.An augmented PEEC formulation has been proposed in [106] to model dispersive and lossy linear magnetic materials again under the quasi-static hypothesis.A novel PEEC formulation that is able to handle 3-D nonlinear problems has been proposed in [115].A behavioral magnetostatic hysteresis model has been implemented in a PEEC environment and presented in [116].Fig. 10 shows an example of a closed magnetic circuit that is excited by a loop driven by an enforced current.The time evolution of the B and H z-components, probed at the center of the excited limb, are presented in Fig. 11 along with the corresponding B-H curves.

VI. PEEC SOLVERS
The equivalent circuits generated by the PEEC method are well suited to be analyzed in both the frequency and the time domains.

A. Frequency Domain Solvers
Frequency domain solvers solve the MNA set of equations ( 17) over a set of frequency samples s k = jω k , k = 1, . . ., N f , in the frequency bandwidth of interest.
For small problems, with a number of DoFs not exceeding 30 000, the solution can be addressed using direct solvers, e.g., the LU decomposition.The direct solution can be accelerated by resorting to the multiscale block decomposition and the adaptive cross approximation (ACA) technique [117], as described in [118], [119], [120], and [121].
Typically, the filling of partial inductances and coefficients of potential matrices, L p (jω) and P(jω), respectively, is computationally heavy since they must be computed at each frequency.Their computation can be accelerated by resorting to the fast multipole method (FMM) [122].Indeed, these matrix entries fluctuate slowly with frequency and have a polynomial behavior that can be efficiently interpolated by the Lagrange interpolation scheme, as was shown in [123].In [124], an effective methodology for the interpolation in space of the partial inductances required by the PEEC method is developed.It makes use of the cubic spline interpolation method, that, under the hypothesis of cubic meshed regions, guarantees always a significant speed-up without loss of accuracy.
For large problems, iterative solvers become mandatory because matrices become so large that they cannot be stored.In this case, matrix-vector products must be computed.An approach that is based on the compression of the partial inductance matrix utilizing the QR decomposition of the far coefficients submatrices is proposed in [125].To speed up the matrix-vector products, the translational invariance of elementary magnetic and electric field interactions can be exploited.In particular, in [126] a method is proposed that allows one to exploit providing significant memory saving and an excellent improvement of the computation performances.This method was then applied to inductance extraction in [127].A Tuckerenhanced and FFT-accelerated version of the method has been proposed for capacitance and inductance extraction in [128] and [129], respectively.Such an approach has been applied to PEEC models of power electronics applications in [18] and [19].
The ACA coupled with hierarchical matrix (H-matrix) arithmetics [130] can also provide an effective method to increase the size of the largest solvable problems.This approach has been investigated in the framework of the PEEC method in [131] and [132].
1) dc and Low-Frequency Solution: A rigorous dc solution of PEEC models of conductors has been proposed in [133] using a two-step process.First, a magnetostatic problem is analyzed by solving a purely resistive network providing the currents in the circuit.Then, an electrostatic problem is solved by enforcing the neutrality of each conductor, which is completely disconnected from the others, and using the voltage drops known from the magnetostatic problem as a boundary condition.This methodology has then been extended in [134] to PEEC models consisting of conductors but also of ideal dielectrics and magnetic materials.Fig. 12 shows the geometry of a loop inductor and the corresponding distribution of node dc potentials.
2) PEEC Models With Dyadic Green's Function: The toutcourt application of the PEEC method to problems involving large layered media, as they occur in PCB modeling or in problems involving transmission lines over or buried in the ground, would require the dielectric to be discretized using a 3-D grid.Thus, it would end up in a large system of algebraic equations with a computational cost that easily becomes prohibitive.A valuable alternative for such problems is represented by the use of the dyadic Green's functions for layered media (DGFLM).This formulation requires the discretization only of conductors, while the Green's functions include the features of the multilayer substrate.A novel numerical solution for the mixed potential integral equation for layered media using the PEEC method has been presented in [135].The same concept has been used for modeling patch antennas [136], PI problems, multilayered PCBs [137], [138], [139], [140], and 3-D IC/packaging analysis [141], [142].
When the PEEC method is adopted to analyze lightning transients in wire/plate structures in the air and lossy-ground space, it is necessary to evaluate the Green's functions for layered media for both vector and scalar potentials for source and field points in any of these two layers.This approach has been applied to address both transient current and voltage in a variety of structures with arbitrarily oriented lines and plates [143], [144], [145], [146], [147], [148], [149], [150], [151], [152], [153].
Time domain solvers: Full-wave PEEC models can be analyzed in the time domain using three different approaches which are as follows.

3) Time Stepping Methods:
The enforcement of the Kirchhoff voltage and current laws to the PEEC model yields the following set of neutral delayed differential equations (NDDE) [57] The vector of the unknowns x(t) ∈ n u ×1 is given as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where i(t) are the branch currents, φ sr (t) are the scalar potentials for surface nodes, φ i (t) are the scalar potentials for internal nodes, v d (t) are the excess capacitance voltages for dielectric branches, and q s (t) represent the surface charges.Furthermore, the state space matrices C(t), G(t), and B are where * represents the convolution operator, C d is the excess capacitance matrix [13], R is the branches resistance matrix, A s is the incidence matrix for the surface nodes, A i is the incidence matrix for the internal nodes, Γ is the dielectric region selection matrix, M is the surface-to-node reduction matrix, and G le is the load conductance matrix (assuming by now for simplicity that only resistive lumped elements are connected to the PEEC model).Also, the time-dependent partial inductance matrix L p (t) and the coefficient of potential matrix P(t) are considered as impulsive, that is where the superscripts DL and D stand for delay-less and delayed, respectively; N L p is the number of significant delays between elementary volumes while N P is the number of significant delays between elementary surfaces; τ cc,i = R cc,i /c 0 , i = 1, . . ., N L p and τ cc,q = R cc,q /c 0 , q = 1, . . ., N P denote the delays between the centers, identified by R cc,i and R cc,q , respectively, of the spatial supports of the basis functions of currents and charges; c 0 is the speed of the light in the background medium.Finally, the source vector u(t) is given as where v s (t) and i s (t) are the voltage and current sources, which are applied to branches and nodes, respectively.The solution of (36) has been addressed in [154] experimenting with different temporal basis functions.Fig. 14.Computed and measured ESD coupling waveform (see [159] for the details).Time domain PEEC modeling has been combined with a thermal solver to perform an electrothermal analysis in [105].Fig. 13 shows the surface temperature of a multilevel chip interconnect [155].
In [156], the analysis of delayed PEEC models has been accelerated by using the waveform relaxation technique.An innovative time domain digital wave PEEC solver has been presented in [157] under the quasi-static hypothesis.
PEEC time domain solvers have been used to analyze ESD problems in [158] and [159].Fig. 14 shows the calculated and measured voltage waveforms induced by a 2-kV ESD event.
a) Time domain stability issues: Stability has for a long time been a challenge for time domain methods [160].It is an issue also for PEEC models like (36).
The stability analysis of solutions consists of two basic aspects.
First, the stability of the model, and thus, of its differential equation system, is to be considered.If the PEEC model is unstable, its numerical solution is generally unstable.However, sometimes, one can obtain a stable numerical solution by using a numerical scheme whose region of absolute stability allows one to obtain a stable solution even for unstable models, but it is not granted that this solution represents well the physical response of the system.
In [161], it is shown that the discretization of PEEC models can lead to differential equations with unstable solutions.A stabilization scheme is introduced using so-called split cells.Damping resistors introduced in parallel to self-inductances in the standard PEEC have been suggested as an effective tool for the stabilization of the solution [162].This modification, although it does not increase the number of unknowns, since it increases the damping, may affect the accuracy of the results.
The accuracy in the computation of partial elements may have an impact on the stability.We will distinguish two cases: 1) quasi-static (QS) PEEC models (propagation delays are neglected); PEEC models are described by RLC equivalent circuits and, Kirchhoff laws lead to ODEs; 2) full-wave PEEC models (propagation delays are considered; they are described by RLC circuits with delayed sources and Kirchhoff laws lead to NDDEs.In both cases, PEEC model stability is significantly affected by the mesh.Indeed, the aspect ratios of the elementary volumetric and surface regions impact the accuracy of the computation of partial elements.It happens when they are computed using analytical formulas, which is possible under the quasi-static hypothesis and using orthogonal meshes.Double precision floating point (IEEE Standard 754), introduces a truncation error in their computation.It is known that for cells with extreme form factors (the ratios between the sizes of the cells) or for extreme ratio distance/size, the truncation error might be so large that the computation of the partial elements is completely dominated by the numerical error.In [19], guidelines for using the analytical formulas are provided depending on the sizes of the elementary domains and their distances.
Then, as pointed out in [163], [164], and [165], the CC approximation is recognized as one further cause driving the model to instability.In [163], a low pass filter (LPF) is used for each partial element, and then the state-space description of the resulting delay-free system is constructed and the stability is assessed by checking the eigenvalues of the system matrix.However, to match the behavior at sufficiently high frequencies, the order of each LPF has to be as high as 38, resulting in quite a significant increase in the number of DoFs.This limitation is mitigated in [166], where time domain models of integrals of Green's functions are proposed leading to full-spectrum convolution macromodeling.
All these types of approximations may lead to right-half-plane unstable poles causing the so-called late time instability.It is important to observe that the unstable poles are very weak in amplitude and typically have a very large imaginary part, so they do not impact significantly the low frequencies.Fig. 15 shows an example of poles of a zero thickness conductor PEEC model, 10-cm long and 2-cm wide.In particular, the poles have been computed first neglecting all the delays (Quasi-static) and then using a first-order Taylor expansion to approximate the exponentials in a polynomial form.
To establish a priori the stability of a PEEC model, it would be necessary to compute explicitly its poles.Due to the theoretical complexity of a time-delayed system like PEEC, however, the existing stability tests are all either sufficient or necessary, but not both.The analytic computation of poles is possible only in a few simple cases or when delays are neglected.On the other hand, sufficient conditions hold under strong assumptions, considering only a reduced number of delays, and with a significant computational cost, making their utility questionable [167], [168], Fig. 15.Stable and unstable poles of delayed PEEC models.[169], [170].An effective pole-based stability test is proposed in [171].
Second, an unstable solution may be caused by the numerical method itself.Different numerical methods have different regions of absolute stability.They can be described as A-stable, L-stable [160].
The influence of the numerical method features on the stability of full-wave PEEC models has been investigated in [172], where the P-stability has been defined as a criterion for choosing an appropriate numerical method.Bellen et al. [167] studied a PEEC model with one single delay, and propose a sufficientcondition test based on checking some matrix-norm inequalities.Extension of this method to models with more delays, however, remains unclear.In general, the backward Euler (BE) and Lobatto III-C have also been proposed as suitable numerical methods for the PEEC method.
a) The Inverse Fourier Transform: The easiest way, although not the most efficient one, to obtain TD responses of PEEC models is represented by the use of the FT applied to the FD results made available by an FD solver of the equations in (17) over a pertinent set of frequency samples s k = jω k , k = 0, . . ., N f .This approach has some disadvantages as follows.
1) The frequency response must be computed over the entire bandwidth of interest in order to evaluate the time domain counterpart; the maximum frequency for which the model is accurate based on some criteria is decided by the mesh; frequency responses at higher frequencies are less accurate.This will only impact the accuracy of the IFT computed waveforms, not the stability because the IFT entails a weighted sum of multiple contributions and it is not affected by the stability problems of time-stepping methods.2) If the sources to the system are waveforms with sharp edges (such as a pulse with sharp rise/fall times), a very large number of frequency points would be required in order to avoid aliasing problems [173].3) Nonlinear materials or lumped elements cannot be easily handled.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
In addition, in order to achieve accurate results, it is necessary to let the energy vanish at the end of the simulation.This means that the signal exciting the system as well as the associated response have to be band-limited in frequency and time.This poses a serious limitation in computing the step response since the FT of a step function is not band-limited.As an approximation, FT of a square (more likely trapezoidal) pulse going back to zero once the steady-state output is reached can be used.Additional simulation time is required to ensure that the system response goes to zero as well.This results in a larger number of frequency points where the solution of the system equations is needed.Moreover, the IFT needs the dc solution of the system which is a well-known issue for EM models (see Section VI-A1).
a) The Numerical Inversion of the Laplace transform: The numerical inversion of the Laplace transform (NILT) technique represents a valid alternative to time-stepping solvers.This technique, introduced in the early '70s of the last century [174], has gained more and more interest in the TD characterization of linear-time-invariant (LTI) systems.Essentially, given the state vector X(s) of an LTI system expressed in the complex frequency domain in the form as it appears in ( 17), the well-known Laplace inverse transform can express its TD counterpart as [175] x where α > e(p k ) ∀p k , where p k are the poles of X(s).The Laplace inverse integral can be approximated by replacing st with z and e z with its [N/M ] Padé approximant [176] ξ N,M (z), finally employing the Cauchy theorem of residues [177].Assuming M an even integer, the NILT approximation for the TD state vector results in being z i and K i the Padé poles and residues, respectively.It is evident from (45) that the knowledge of the system's poles is not required to build the transient waveforms associated with the corresponding Laplace domain quantities.Moreover, the method's accuracy is not bounded by the time-step choice and, consequently, does not depend on the number of time samples employed.Furthermore, also the solution's stability of the solution is not affected by the time-step.Hence, the responses of LTI PEEC models remain always stable when reproduced by NILT.For illustrative purposes, Fig. 16 shows two different sets of poles z i /t corresponding, respectively, to the choices M = 6 and M = 12, with N = M − 2. The two sets are computed over the same time window [1 − 5] ns.
It is evident that, as time t increases, each set of poles approaches the origin of the complex plane, causing a loss of accuracy with the evaluation time.To mitigate this limitation, a modified NILT, known as NILTn, has been recently proposed in the framework of TD simulation of multiconductor transmission lines (MTL) [178].The MNA representation of full-wave PEEC Fig. 16.Two Padé poles sets moving in the complex plane as the evaluation time increases.models in the Laplace domain (43), being the NILT scheme easily applicable.The NILT technique has been applied for the first time to retarded PEEC models in [179], where its stability properties are emphasized when compared to time-stepping methods applied to the resolution of retarded PEEC models.For illustrative purposes, in Fig. 17 the port voltage response of a power divider, a typical microstrip device, is depicted, assuming a trapezoidal waveform at its feeding port.The EM model of the structure has been built through the PEEC model considering the propagation delays, giving rise to an LTI retarded system.The NILT results agree with those of the IFFT, while those obtained through the BD2 time-stepping scheme exhibit clear late-time stability issues (t > 1 ns).
Subsequently, in [180], the modified NILTn has been successfully employed for retarded PEEC models showing accuracy improvements and reduced computational costs.
It is worth noticing that the NILT approach does not compute the response based on the singularities of the PEEC model, but rather on the relevant Padé matrix.Hence, the eventual righthand side (RHS) unstable poles of the PEEC model are automatically filtered and do not impact the accuracy of the solution.

VII. COMPARISON WITH THE METHOD OF MOMENTS
Both the MoM and PEEC integral equation-based approaches are using what is known as the weighted residual method (MWR or WRM) [181], [182] to discretize the formulation into numerical solutions.Originally, the mathematical MoM notation was part of fundamental WRM techniques.As pointed out in [181], WRM applies to the solution of both differential as well as integral equation-based problems.Hence, WRM can be used for both approaches we consider here.
To give a unique name to the MoM considered here, we call it the Z-MoM method since the formulation leads to an impedance formulation.Of course, the Z-MoM name can also be used for the time domain versions since the basic equations are the same.Further, the Z-MoM method was devised in a time before PEEC [4].The key difference between the two methods is in the fundamental derivation of the solution which results in different unknowns.The derivation of the PEEC method is given in Section II, where it is clear that the unknowns are both the potential Φ and the current I, which leads the equations in the MNA form (17) which can be implemented with other SPICE type circuit domain models.
The Z-MoM method is based on the same fundamental integral equation ( 4) as is PEEC.However, the charge or potential variable is eliminated using the continuity equation (5).This results in the impedance formulation of the Z-MoM method.
Another difference exists between Z-MoM and PEEC meshing approaches.The main issue is based on the fact that Z-MoM has mesh cells only for the current variables, while PEEC requires appropriate, separate cells for both the potential as well as the currents, as is obvious from the formulation in Section II.It is clear that the inductive and resistive partial elements connect between nodes while the capacitive cells are connected to the nodes.This applies to all types of cells including nonorthogonal meshing.In the original PEEC method, the meshing was mostly based on orthogonal bricks or hexahedron volumes and rectangular cells or quadrilaterals to mesh the surfaces.The Z-MoM mostly uses triangular meshes for surfaces [94] and tetrahedral meshes for volumes [95].This has led to the development of dedicated basic functions, known as the RWG basis functions for triangles [94] and the Schaubert-Wilton-Glisson (SWG) basis functions for tetrahedral volumes [95].
The use of triangular/tetrahedral cells allows for an accurate model of complex geometries.However, easily results in a large number of cells and unknowns when compared to quadrilateral meshes.
An important issue for the Z-MoM formulation is the so-called low-frequency breakdown.This is because the low-frequency response including the dc [183] solution is poor or missing.This high-pass behavior is an important issue for the solution of some EMC, SI, and PI problems while it is not an issue for the solution of high-frequency antenna problems.The Z-MoM formulation has been improved by resorting to augmented formulations which also assume the charges as unknowns, and also by adopting suitable scalings to improve the condition of the system [184], [185].In [186], it has been shown that the MNA form of PEEC models attenuates the low-frequency breakdown even for problems that do not have closed loops which also results in dc solutions.

VIII. COMPUTATIONAL COMPLEXITY OF PEEC MODELS
Being an integral equation-based method, similar to MoM, the solution of PEEC models in both the frequency and time domains for quasi-static models requires the solution of a blockdense linear system.In the frequency domain, small problems, with a number of unknowns N < 30 000, can be solved with direct solvers.They involve a factorization step followed by a solve step.The factorization step comprises an LU factorization or QR factorization, which is generally computationally very expensive.Further, direct solvers are advantageous when one is interested in multiple RHSs.A naive direct solver costs O(N 3 ), which is prohibitively resource-demanding for large system sizes.To reduce the computational complexity, fast methods are used.Many dense matrices arising out of N -body problems possess a hierarchical low-rank structure.This low-rank structure is exploited to construct hierarchical matrices and hierarchical matrices-based fast direct solvers [130], [187], [188].In the framework of the PEEC method, multiscale compressed techniques exploiting the low-rank nature of magnetic and electric field interactions have been presented in [118] and [189].In [190], it has been shown that the storage requirements and computation time of these types of techniques scale as O(N 3/2 ) and O(N 2 ), respectively, for asymptotically high frequency.
Larger-dimensional problems require iterative solvers that require repeated vector-matrix products with O(N 2 ) complexity.It is possible to accelerate such vector-matrix products by resorting to methods, such as the FFM [193] and Barnes-Hut [194].Using a multilevel approach, a O(N logN ) computational complexity can be achieved asymptotically.Further, for fast convergence in problems with high-condition numbers, iterative solvers are coupled with a preconditioner.More recently, the FFT-based technique presented in [126] has been applied to the PEEC method in [18] and [19] allowing obtaining excellent memory and CPU time savings.
In the time domain, the system matrix is typically very sparse.This enables the efficient use of multifrontal techniques to perform the LU decomposition even for very large problems [195], [196], [197].Nevertheless, when using time-stepping techniques, the RHS vector has to be updated at each time step requiring many matrix-vector products.In this case, FFT-based techniques can also be used to accelerate these matrix-vector products [198].

IX. MODEL ORDER REDUCTION OF PEEC CIRCUITS
The equivalent circuits generated by the PEEC modeling can be very large and easily require prohibitive storage and significant computing resources.Model order reduction (MOR) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
techniques [199], [200] can be adopted to compress the C, G, B matrices in (36) preserving the accuracy of the reduced model at the ports.Equation ( 36) can be rewritten in the Laplace domain as Several researchers have applied MOR techniques to examples of quasi-static PEEC models without retardation.See, for example [201], [202], and [203].Procedures for reducing PEEC models with retardation were first proposed in [204] and [205].Both of these references construct reduced-order models for the frequency domain transfer function.For PEEC models with retardation the transformed system matrix Q(s) contains many elements with factors of the form e −sτ for some τ corresponding to a time delay (retardation) in the circuit.All of those papers expand each delay term in an infinite Taylor series.The authors in [204] used a single expansion point and applies asymptotic waveform (AWE) model reduction [206] to the resulting infinite order linear system.In that article, an example is included for which good approximations were obtained for low frequencies up to 1 GHz.Chiprout et al. [205] used complex frequency hopping (CFH), considered several expansion points, applied AWE to the infinite linear system obtained at each such point, and then combined the pole and residue information to obtain an approximation to the transfer function of the PEEC system.It includes examples showing that this method is accurate to approximate the transfer function of a system up to 4.8 GHz.A similar approach has been applied to a PEEC model in conjunction with an FFT grid representation in [207].Furthermore, in [207], an equivalent first-order system is computed by means of a single-point Taylor expansion, and a corresponding orthogonal projection matrix is computed by means of a block Arnoldi algorithm [208], [209].However, the NDDE formulation is not preserved.In [210], a readily parallelizable procedure is proposed for generating reducedorder frequency-domain models from general full-wave PEEC systems.Multiple expansion points and piecemeal construction of pole-residue approximations are adopted to approximate the transfer functions of the PEEC systems through an Arnoldi recursion.The reduction of equivalent first-order systems become computationally expensive and sometimes not feasible when large delays are involved since exponential terms with large delays need many terms in the Taylor expansion to be accurately approximated.The multipoint expansion [211] addresses this issue and is able to accurately reduce NDDE systems with large delays since a small expansion Taylor order can be used for each expansion point and the accuracy of the reduced model is increased by adding new expansion points.In [212], an adaptive algorithm is proposed to choose the expansion points, assuming that the order of the Taylor expansion is fixed for each expansion point.
When applying MOR techniques, in the simplest case, the reduced-order model can be obtained by performing a singlepoint expansion reduction.Denoting the orthogonal basis as K ∈ n u ×n r , the reduced order system is given by sC where the following congruence transformations are used: G r,i = K T G i K, i = 0, . . ., n τ (48b) and X r (s) is a vector containing the state variables in the reduced domain.To obtain a compact and accurate reduced order model over a wide frequency range it is important to construct the orthogonal basis K. To this purpose, the Arnoldi algorithm [208] can be adopted because of its numerical reliability and robustness.Indeed, the Arnoldi algorithm is able to provide the orthogonal basis for the transfer function moments without computing these moments explicitly.Unfortunately, adapting the Arnoldi algorithm for NDDE systems is not straightforward, and, thus, the original NDDE system must be transformed into a suitable form for standard Arnoldi-based reduction.To address this issue, by expanding the exponential factors e −sτ k in a Taylor series form and using a companion form [213], an equivalent first-order system is computed.Then, the Arnoldi algorithm is applied to the first-order equivalent system, allowing the computation of the corresponding orthogonal projection matrix K for the original NDDE system, and a reduced NDDE system is finally obtained [213].This single point expansion-based MOR algorithm for NDDE systems was proposed assuming s = 0 as an expansion point.If any other expansion point s = s 1 , s 1 = 0 is selected, setting s = s 1 + σ, where s 1 is a frequency shift and σ is the new Laplace variable, the NDDE system (46a)-(46b) reads where and the algorithm described in [213] is applied.While the single-point MOR approach [213] is able to preserve the NDDE formulation, it may be not able to reduce NDDE systems with large delays since the reduction of equivalent first-order systems becomes computationally expensive and sometimes unfeasible because many terms in the Taylor expansion of the exponential terms with large delays are required to preserve the accuracy.A more valuable alternative is represented by the multipoint expansion [211] that is able to accurately reduce NDDE systems with large delays since a reduced order of the Taylor expansion can be used for each expansion point while the accuracy of the reduced model is achieved by adding new expansion points.If the order of the Taylor expansion is fixed for each expansion point, an adaptive algorithm can be used to identify the expansion points.As described in [213], the NDDE formulation is preserved in the reduced model.At each expansion point, the MOR algorithm described in [213] is applied and the corresponding projection matrix K i , i = 1, . . ., n points is computed, where n points denotes the number of expansion points.The final projection matrix K is based on the orthogonalization of the stack column collection of all single expansion point projection matrices.
Parameterized reduced order models for sensitivity analysis of PEEC circuits have been presented in [218] and [219].
A different approach is presented in [220] to perform the reduction of quasi-static PEEC models.By introducing a general circuit transformation that can be directly applied to the circuit configuration of the PEEC model, the proposed MOR method progressively reduces the order of the original problem by at least one order of magnitude without involving any matrix operations.Furthermore, another attractive feature of this method is that its computational overhead is dominated by the operation of outer products in combining processes, which takes more than 95% of overall computing time.This feature allows the MOR process to be significantly accelerated by multicore parallel computation using a massive GPU acceleration technology.

X. CONCLUSION
The purpose of this work is to summarize the progress made in PEEC methods over the years.The key uniqueness of the approach is the connection between electromagnetic and circuit theory.This covers the important issue areas of combined systems with integrated circuits and EM problems.Over the course of fifty years, the PEEC method has had many developments and improvements regarding the types of materials that can be treated, the meshers, and the solution methodologies, both in the time and frequency domains.The number of scientific articles dealing with it has steadily increased over time as has the number of applications using it.In particular, its full-wave nature combined with circuit interpretation and rigorous dc solution has made it popular with engineers dealing with EMC/SI/PI problems.A lot has been done, the best is yet to come.

Fig. 1 .
Fig. 1.Number of articles on the PEEC method from 1972 to present.

Fig. 2 .
Fig. 2. Equivalent circuit modeling the displacement currents in the background medium.

( 23 )
leads to the expanded version of the partial inductance L F W p m,n

Fig. 11 .
Fig. 11.Closed magnetic circuit.(a) Time evolution of B z and H z at the center of the excited limb.(b) B x versus H x at the same point.

Fig. 13 .
Fig. 13.Transient temperature map on the surface of a multilevel on-chip interconnect.