Efficient Characterization of Interconnects With Arbitrary Polygonal Cross Sections Using Fokas-Derived Dirichlet-to-Neumann Operators

A novel technique to accurately characterize interconnects with general, piecewise homogeneous material parameters and arbitrary polygonal cross sections is presented. To compute the per-unit-of-length (p.u.l.) complex inductance and capacitance matrices of the considered structures, we apply a boundary integral equation (BIE) framework, invoking a Dirichlet-to-Neumann (DtN) formalism to recast the problem at hand. The pertinent operators are constructed by means of the numerically fast Fokas method, leveraging fully analytical expressions for the pertinent matrix elements. Numerical examples of various multiconductor transmission lines demonstrate that our proposed scheme is flexible and precise. Since our method is not limited to rectangular cross sections, manufacturing effects such as etching can also be taken into account. Moreover, the examples are not restricted to resistance, inductance, conductance and capacitance (RLGC) data as we also consider signal attenuation, slow-wave factors, and crosstalk.


I. INTRODUCTION
M ODERN interconnect structures are becoming ever more miniaturized and complex, operating up to very high frequencies for which the corresponding wavelengths are comparable to their longitudinal dimension.However, as the dimensions of their cross section most often remain sufficiently small with respect to the smallest wavelength of the transported signals, it is still possible to describe these interconnects using frequency-dependent, per unit-of-length (p.u.l.), RLGC transmission line parameters [1] derived from a suitable quasi-transverse magnetic (quasi-TM) electromagnetic analysis, as thoroughly motivated and validated in [2] and [3].The quasi-TM regime implies that the longitudinal component of the magnetic field remains negligible.Moreover, the broadband nature of the signals, with frequency components ranging from dc to several tens of GHz, necessitates an analysis that accounts for current crowding, conductor losses, and the skin effect in the cross section as these phenomena change dramatically over the considered frequency range.
Furthermore, magnetic materials are rapidly gaining ground in emerging interconnect topologies, such as metaconductor configurations with alternating metal layers [4], while modern manufacturing techniques allow for intricately shaped structures [5].In view of these evolutions, it is imperative that state-of-the-art modeling schemes support these material properties and geometrical particularities to keep pace.
In the presence of good conductors at high operating frequencies, the skin effect imposes a very fine discretization and, consequently, an intractable number of unknowns for volumetric methods, while special care must be taken when evaluating the conductive medium Green's function integrals in BIE techniques.The differential surface admittance (DSA) formalism, first proposed in [15], mitigates these challenges by reformulating the problem via an exact, global operator to a situation where each conductor is replaced by an equivalent surface current density, residing in the background medium.While this technique was, at first, limited to nonmagnetic materials and circular or rectangular cross sections to fully exploit analytical expressions for the pertinent matrix elements, several extensions have been presented: semiconducting materials [3], triangular cross sections [16], magnetic contrast [17], and 3-D lossy interconnects [18], [19].The alternative construction of the DSA operator in [20], [21], and [22] expanded the applicability to arbitrary cross sections and layered media.However, this alternative again relies on Green's function integrals inside the (conductive) medium, which might occasionally lead to time-consuming numerical integration routines for intricate geometries as shown in this contribution.Interconnect modeling remains an actively researched topic.Very recently, for example, a potential-based modeling approach was presented in [23] and the capacitance problem was thoroughly discussed in [24].
In an earlier conference paper [25], we proposed a novel method to extract the p.u.l.resistance and inductance of a multiconductor transmission line with polygonal cross sections, applying a DSA operator [15] derived with the Fokas method [26], [27], which solves the relevant integral equation without relying on quadrature rules.In this article, we substantially extend our technique [25] to enable the characterization of these interconnects in terms of their p.u.l.capacitance and conductance [28] as well, representing the dielectric media by means of a dedicated Dirichlet-to-Neumann (DtN) operator that is once more constructed using the Fokas method.As such, a comprehensive framework is obtained, enabling the full signal integrity (SI) analysis of interconnects in terms of, for example, signal attenuation, slow-wave factors, and crosstalk, applicable to multiconductor configurations with polygonal cross sections and combined dielectric and magnetic material contrast.
The remainder of this article is organized as follows.In Section II, we formulate the foundations of our method, concisely reiterating the inductance problem and elaborating the solution of the capacitance problem.The construction of the novel, pertinent DtN operator for the capacitance problem, leveraging the Fokas method, is detailed in Section III, concluding the theoretical part.Various relevant examples of interconnect structures, including intricately shaped coplanar waveguides (CPWs), are treated in Section IV, where we compare the numerical results of our proposed method to references found in the literature.Finally, conclusions are provided in Section V.

II. FORMULATION OF THE METHOD
Assume a quasi-TM, time-harmonic regime with an e ȷ ωt convention.We consider an interconnect structure consisting of piecewise homogeneous, polygonal components S n , with boundaries C n , surrounded by a homogeneous background medium (ϵ e , µ e ), as illustrated in Fig. 1.In particular, the polygonal components S n comprise a number of signal conductors and reference conductors, with material parameters (ϵ i , µ i , σ ) and (possibly lossy) dielectrics, characterized by their complex permittivity and permeability ϵ i and µ i .Their longitudinal dimension is aligned along the z-axis of the Cartesian coordinate system (x, y, z).
Recall that the quasi-TM assumption is valid provided the interconnect cross section remains small with respect to the wavelength.We briefly summarize the main results of [25] in Section II-A and demonstrate the novel, Fokas-based solution of the capacitance problem in Section II-B.

A. Inductance Problem
For the inductance subproblem, we need to model the good conductors in the considered interconnect.By means of a single-source equivalence theorem, we substitute the background medium for the component's material, keeping  the same exterior fields (e e , h e ) by imposing an equivalent surface current density j s = j s,z ẑ on the contour C. The original interior fields (e i , h i ), however, do change to the quantities (e ′ i , h ′ i ), bearing no physical meaning.For each of the conductors, the equivalent surface current density is given by We illustrate our approach for one conductor with a polygonal cross section in Fig. 2, denoting the M corner points (x m , y m ) in the transversal plane using complex numbers ζ m = x m + ȷ y m .The electric field e = e z ẑ at the boundary of the object is identical for both situations (a) and (b) depicted in Fig. 2. It is mapped to its normal derivative, and hence to the corresponding magnetic fields h i and h ′ i , by invoking DtN operators X and X ′ , respectively.Recalling (1), we obtain where Y is the desired DSA operator.In [27], we demonstrated in full detail how this operator can be constructed for polygonal structures by means of an extended Fokas method for the Helmholtz equation with complex wavenumber.For conciseness, we only restate the final result here, viz., the discretized version of (2), in terms of local, pulse-shaped basis functions b i where the vectors v and j collect the expansion coefficients of e and j s , while cc is the Gram matrix of the local basis functions.
Invoking the approach presented in [15] for a configuration with N conductors, we determine the N × N p.u.l.resistance Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
and inductance matrices R and L as follows: evaluating the matrix elements of A via where the pertinent 2-D static Green's function is defined as The second expression is used when an infinite, perfectly electrically conducting (PEC) ground plane is present.In this case, r ′′ is the mirror image of r ′ across this plane.Finally, the incidence matrix T is given by with ℓ i being the length of segment i in the mesh.

B. Capacitance Problem
To fully characterize the interconnect structure and to perform a SI analysis, the capacitance subproblem needs to be solved as well.Therefore, in this article, the approach developed for rectangular cross sections in [3] is extended to arbitrary polygonal cross sections.We first focus on the dielectric objects, replacing them with their contrast surface charge density with e n being the outward pointing normal component of the electric field.The charge densities ρ eq we introduced now reside in the homogeneous background medium with permittivity ϵ e .Consequently, the induced potential φ, which satisfies the Laplace equation ∇ 2 φ = 0 in source-free regions, is determined by where G is again the static Green's function (6).Next, we approximate e n by the normal derivative of the potential: −∂ n φ, which is valid in the quasi-TM regime [2], [3].Invoking the pertinent DtN operator D, this quantity is linked to the potential φ itself, on the boundary of the dielectric.We now propose to utilize the computationally efficient Fokas method to determine D, instead of relying on the eigenmode expansion in [3], seamlessly extending the applicability of this technique to polygons.The corresponding derivation is provided in Section III.Substitution of this result in (8) yields We employ the same discretization procedure with pulse-shaped basis functions as before, collecting all ρ eq and φ values for the dielectrics in vectors ρ d and φ d , respectively, ultimately leading to the following discretized version of (10): In this expression, dd is the Gram matrix of the pulse basis functions and D is the DtN matrix corresponding to D. Now, we also include the N conductors, with ρ eq and φ values collected in ρ c and φ c , and discretize (9) as follows: where the elements of A ′ are determined as in (5), but with µ e replaced by 1/ϵ e .The subblocks pertaining to the conductors (subindex c) and dielectrics (subindex d) are identified as Combining ( 11) and ( 12), we obtain Lastly, for each of the good conductors, we may assume a fixed potential that is constant over the entire boundary.We consider N situations, where only one of the conductors has a unit potential, while the others are equal to zero.Solving (15) for ρ c and integrating this quantity over the conductor boundaries, for each of the N excitations φ c , finally yields the (columns of the) sought-after N × N complex capacitance matrix C + G/ȷ ω, connecting the constant potential values on each conductor to their total charges.

III. CONSTRUCTION OF THE DTN OPERATOR FOR THE CAPACITANCE PROBLEM
To determine the DtN operator D that was introduced in Section II-B, we need to solve the boundary value problem (BVP) governed by the Laplace equation, subject to the Dirichlet boundary conditions given by the known potential φ, on the bounded domain S of a conductor.
To construct the operator in an efficient way and for general polygonal cross sections, we leverage the Fokas method, or unified transform method [26].Therefore, we start by considering the known potential φ(x, y) and another function v(x, y), both satisfying the Laplace equation in the domain S ∇ 2 x y φ = 0 (16) Then, from Green's second identity, we get where C is the boundary of S. Recalling our definition of the complex variable ζ = x + ȷ y (see also Fig. 2), we find a particular solution of the Laplace equation ( 17) where λ is a parameter, which allows us to reformulate (18) as It should be mentioned that v = e ȷ λ ζ , where • denotes the complex conjugate, also satisfies the Laplace equation (17) and must be considered as well to obtain the complete final solution.Nonetheless, for conciseness, we will proceed with (19) only, since the procedure is entirely analogous for both v and v.
As mentioned earlier, in this work transforming (20) into where α m = −ȷ H m .To this end, we invoked the Fourier transform of the Legendre polynomials with α = −ȷ τ and I p+1/2 the modified Bessel function of the first kind and order p + 1/2.Subsequently, ( 24) is evaluated for well-chosen spectral collocation points λ ∈ C per side, with > P leading to an overdetermined system to be solved for D m p M × P times, with each time all C m p equal to zero, except for one pair ( p, m).As such, the DtN matrix D in the Legendre domain is determined.Note that the λ-values (26) are optimized for the Laplace equation [29] and differ from the set used for the inductance problem [25].More specifically, Fig. 3.
for each side m, the associated collocation points λ are chosen such that the terms pertaining to this particular side are dominant in (24).After a final transformation from the Legendre polynomial expansion to more practical local pulse basis functions, we obtain the DtN matrix D, enabling us to determine the unknown function ∂ n φ m at the boundary C, provided φ m is known.

IV. NUMERICAL RESULTS
Our formalism is now applied to a number of application examples.Note that, each time, one or more reference conductors are selected.Consequently, the p.u.l.parameter matrices, determined as elaborated in the previous sections, should still be reduced to reflect Sections II and III.For Z = R + ȷ ωL, we express all voltages with respect to u ref on the reference conductor(s), effectively subtracting the corresponding row(s) of Z from those pertaining to the signal conductors.Furthermore, we impose that the current flowing through the reference conductor(s) is minus the sum of the currents in the remaining signal conductors.The same procedure can be applied for , where the total charges are now taken into account, instead of the currents.
The number of collocation points per side and the number of Legendre polynomials P on each side are set to ( , P) = (40, 20) for the inductance problem and to ( , P) = (150, 75) for the capacitance problem.

A. Inverted Embedded Microstrip Lines
As a first example, we consider the pair of coupled inverted embedded microstrip (IEM) lines, shown schematically with annotated dimensions in Fig. 3.The signal conductors (indicated with 1 and 2) and the top reference plane (indicated with ref.) consist of aluminum (σ = 3.77 × 10 7 S/m) and the lines are embedded in silicon dioxide (SiO 2 ) with relative permittivity ϵ r = 3.9(1 − 0.001ȷ ).For the silicon substrate, we employ the following properties: ϵ r = 11.7 and σ = 10 S/m.
We determine the elements of the 2 × 2 p.u.l.parameter matrices R, L, G, and C, comparing the results of our proposed approach to [30], which invokes a DtN formalism for rectangles, in Figs. 4 and 5.An excellent agreement is observed.As is done in [30], we consider SiO 2 as the background medium.The remaining components are discretized Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.with segments of length ≈0.5 µm, amounting to a total of 964 unknowns.
The mesh segments are again ≈0.5 µm long, except for the semiconducting substrate, where the mesh density is ten times lower, leading to a total number of 1176 edges.For this structure with only one signal conductor, we determine the p.u.l.R, L, G, and C. The propagation constant γ is then given by γ ≜ β − ȷ α = (−(R + ȷ ωL)(G + ȷ ωC)) 1/2 .Furthermore, we define the slow-wave factor and the attenuation The slow-wave factor and attenuation obtained by means of the method in this work, with = 1 mm, are compared to those provided in [2], in Fig. 7. Once more, we observe that both approaches yield very similar results.

C. Overlay and Elevated Coplanar Waveguides
As a third example, we investigate the overlay coplanar waveguide (OCPW) and elevated coplanar waveguide (ECPW) structures of Fig. 8.In both cases, the central signal conductor is elevated by h = 15 µm, either overlaying the ground planes by a distance o (OCPW), or separated from the equally elevated ground planes by a distance g (ECPW).All conductors consist of 3-µm-thick gold (σ = 4.1 × 10 7 S/m) and are positioned on top of a dielectric with permittivity ϵ 1 = 4.6(1 − 0.02ȷ )ϵ 0 and a height of 520 µm.The lower horizontal and sloped sections of the signal conductor are w = 28 µm and t = 40 µm wide, respectively.Finally, the ground-to-ground spacing w + 2s is fixed to 170 µm, hence s = 71 µm.
We discretize the conductors under study with segments of length 1 µm, and, once again, mesh the substrate ten times coarser.For the OCPW, the overlay o is varied between −30 and 50 µm, where a negative value indicates a separation between the signal and reference conductors.Likewise, for Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the ECPW, the gap g is changed from 3 to 27 µm.As such, the largest number of segments for the OCPW amounts to 1172, while the ECPW required at most 980 edges.Moreover, we consider a situation with and without dielectric, that is, in the latter case ϵ 1 = ϵ 0 , and thus the structures are positioned in free space.For all of the aforementioned instances, we show the attenuation as a function of Re[Z c ] at 50 GHz in Fig. 9, where Z c is the characteristic impedance of the line Comparing the results obtained by means of the proposed technique to the reference taken from [11], we remark a very close match.

D. Ridge-Type CPW
The structure presented in Fig. 10   We apply the same discretization strategy as, for example, in Section IV-C, leading to a problem size of 1410 unknowns.
To analyze the characteristic behavior of the interconnect, we extract its p.u.l.parameters and determine the slowwave factor (27), the attenuation (28), and the characteristic impedance (29).The corresponding plots are shown in Figs.11 and 12, respectively, compared to the reference provided by [31].All results nearly coincide, except for the slow-wave factor.This discrepancy can be attributed to the fact that LiNbO 3 is actually an anisotropic medium, with ϵ r = 28 in the direction perpendicular to the LiNbO 3 -SiO 2 interface, while we assumed ϵ r = 43 isotropically in the entire cross section.the effective relative permittivity is overestimated, explaining a slightly slower propagation predicted our method.

E. Crosstalk in a Multiconductor Transmission Line
As a final example, we revisit the multiconductor transmission line [32], we considered earlier in [25], but now add a lossy dielectric slab (ϵ r = 4, tan δ = 0.01).The trapezoidal signal lines and the rectangular reference conductor consist of metal with conductivity σ = 3.57 × 10 7 S/m and relative permeability µ r ∈ {1, 10}.A schematic representation, annotated with dimensions in micrometer, is provided in Fig. 13.We will study three different signal conductor cross sections: the default trapezoid, an inverted trapezoid (with the lengths of the short and long bases interchanged), and a 1.5-µm-wide rectangle, leaving all other dimensions unaltered.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Multiconductor transmission line (σ = 3.57 × 10 7 S/m, µ r ∈ {1, 10}) with three trapezoidal signal lines, a finite rectangular reference conductor, and a lossy dielectric substrate (ϵ r = 4, tan δ = 0.01).All dimensions are in µm.Fig. 14.Far-end crosstalk between lines 2 and 1 for a 1-mm-long line with a cross section as shown in Fig. 13, for various signal conductor shapes and µ r ∈ {1, 10}.
The p.u.l.R, L, G, and C matrices are determined for all of the aforementioned configurations and utilized in a time-domain simulation of a 1-mm-long transmission line.We excite line 2 with a step signal, rising from 0 to 1 V in a rise time of 25 ps, by means of a voltage with internal impedance 1 .The other lines are quiet and terminated at the near end with a 1-resistor.At the far end, all lines are loaded with a 1-k resistor, in parallel with a 0.1-pF capacitance.The far-end crosstalk on victim line 1 is shown in Fig. 14.We clearly observe the influence of the conductor shape and the relative permeability on the maximum crosstalk level.The delay is also clearly affected by the relative permeability.This was confirmed through the application of method described in [20] and [21], which showed excellent agreement for all six curves in Fig. 14.

F. CPU Timing
To conclude this section, in Table I, we provide the total CPU time required to solve the inductance (ind) and capacitance (cap) problems in Sections IV-A-IV-E, indicating the number of frequency points N f and the number of edges N e in the discretization and including the time to fill the pertinent Green's matrix A and the DSA or DtN matrix.All calculations were performed on a system with a 1.9-GHz CPU and 16 GB of RAM.For the configurations in Sections IV-C-IV-E, the same metrics are reported for implementation of the contour integral method realization of the DSA operator [20], [21] in Table II.
For both instances of the DSA method, the main portion of the computation time is dedicated to the construction of the Green's matrix of the outside problem (5), taking into account the fact that the DSA/DtN matrix has to be computed for every frequency point while the Green's matrix only needs to be filled once.Although the contour integral method [20], [21] is realized with fast point matching, the construction of the Fokas matrix, leveraging accurate Galerkin weighting, is comparatively fast.This is due to the fully analytical expression of its elements, sidestepping the more time-expensive use of a quadrature routine, which can sharply increase the computation time for larger structures, which require a fine discretization, as is the case for the inductance problem of Fig. 10.

V. CONCLUSION
We presented a comprehensive formalism to determine the p.u.l.resistance, inductance, conductance, and capacitance (RLGC) matrices of multiconductor configurations with arbitrary polygonal cross sections and general, piecewise homogeneous material properties.As such, our method is applicable to emerging interconnect topologies with magnetic materials and intricately shaped cross sections in modern manufacturing techniques.We solved the corresponding inductance and capacitance problems leveraging a DtN BIE formalism.To this end, we constructed the pertinent operators by means of the Fokas method for the Helmholtz equation with complex wavenumber and the Laplace equation.Our novel method was validated with various numerical examples.Important design and SI analysis data were presented, such as the characteristic impedance, slow-wave factor, attenuation and crosstalk, demonstrating the method's versatility, accuracy, and efficiency.

Fig. 1 .
Fig. 1.General multiconductor transmission line cross section considered in this work, including polygonal signal and reference conductors, along with (lossy) dielectrics.

Fig. 2 .
Fig. 2. Geometry of a component, illustrating the equivalence theorem, with (a) the original and (b) the equivalent situation.
, we limit S to a polygonal domain with M corners ζ m , m ∈ {0, 1, . . ., M − 1}, numbered counterclockwise, where ζ M = ζ 0 .Now, we introduce the following parameter representation on the side connecting ζ m and ζ m+1 : λh m , andB m = |h m |.To proceed, we project the known φ m and their unknown normal derivative ∂ n φ m for each polygon side m onto a basis of P Legendre polynomials P pφ m ∂ n φ m = 2 (α m ) A m C m p − B m D m p = 0 (24)

Fig. 4 .
Fig. 4. P.u.l.resistance and inductance matrix elements as a function of for the IEM structure of Fig. 3.

Fig. 12 .
Fig. 12. Characteristic impedance Z c of the CPW in Fig. 10, as a function of frequency.

TABLE I CPU
TIMING OF THE ADVOCATED METHOD FOR THE EXAMPLES STUDIED IN THIS SECTION.REPORTED TIMES INCLUDE THE COMPUTATION OF THE GREEN'S MATRIX ELEMENTS, THE CONSTRUCTION OF THE DSA/DTN MATRIX, AND THE TOTAL SOLUTION TIME FOR ALL FREQUENCY POINTS COMBINED TABLE II CPU TIMING OF THE IMPLEMENTED METHOD REPORTED IN [21] FOR THE EXAMPLES STUDIED IN THIS SECTION.REPORTED TIMES INCLUDE THE COMPUTATION OF THE GREEN'S MATRIX ELEMENTS, THE CONSTRUCTION OF THE DSA/DTN MATRIX, AND THE TOTAL SOLUTION TIME FOR ALL FREQUENCY POINTS COMBINED