Methods for modeling interactions between massively coupled multiple vias in multilayered electronic packaging structures

[0001] This application is based on a prior copending provisional application Ser. No. 60/294,657, filed on May 30, 2001, the benefit of the filing date of which is hereby claimed under 35 U.S.C. §119(e).

[0002] This invention generally pertains to the numerical simulation of reflected, transmitted, and coupled signals in respect to vias in structures having multilayered geometry, and more specifically, to an adaptation of the Foldy-Lax multiple scattering formula for modeling signal distortion and coupling between massively coupled multiple vias in multilayered chip packages, chip substrates, and printed circuit boards.

[0003] Components and devices utilizing electrical circuits involving integrated circuit (IC) chips, multi-chip modules, printed circuit boards, connectors, back planes, and cards are generally referred to as electronic packaging. The terms “electronic packaging structure,” “electronic package,” and “package” are used interchangeably throughout this disclosure to generally refer to any and all of the aforementioned devices and structures. Also included in this definition is the specific case of multilayered packaging wherein multiple substantially parallel conducting ground and power planes are separated by dielectric substrates.

[0004] It is desirable and has long been the practice in this technology to employ simulation techniques to establish performance baselines of electronic packages prior to finalizing hardware design and prior to building prototype electronic packages in order to rapidly and inexpensively predict system behavior. While such simulations were simply useful when developing discrete and low density circuits, the preponderance of large scale integration (LSI), very large scale integration (VLSI), and ultra large scale integration (ULSI) implementations has made simulation a necessity. The construction of physical prototypes and long design cycles are undesirable, and hence, rigorous, accurate, and fast simulation of packages is required.

[0005] Via structures are predominantly used in multilayered packages and printed circuit boards to connect together traces residing in/on different layers. These packages typically include at least two separated conducting planes that are metallic, and multiple vias extending between the two or more planes for connecting signal traces disposed on the different planes. These packages are used extensively, primarily because of high speed, high densities, and routing complexity. From electromagnetic theory, voltages at the active via ports are represented as magnetic ring current sources, which in turn generate electromagnetic modes inside the plane structure. Substantial electromagnetic coupling occurs, even in well separated vias. The electromagnetic modes propagate along the plane structure and are scattered by other cylindrical vias, induce voltages at idle ports and currents on idle vias, and in turn, affect the active signal-carrying via. Such parallel plate waveguide effects are induced by the multilayered geometry. As a higher order effect, the affected vias can in turn influence the original signal. Because of the waveguide modes, such coupling is not necessarily localized in space. This lack of localization poses considerable design problems, particularly at high frequencies, in achieving reliability, high speed, and in implementing a simulation. Such coupling can even cause unreliable behavior or complete signal failure, along with signal integrity loss, longer delays, and inappropriate switching of signals.

[0006] When logic gates alter their states, transient noise is generated on vias. Due to coupling between vias, spurious signals then occur throughout the multilayered structure, primarily through cylindrical wave modes coupling to parallel plate modes. If such signals reach terminal ports, the noise can affect input and output devices through spurious switching signals, thereby degrading system performance.

[0007] Several different methods have been developed to attempt to model the crosstalk or coupling mechanism through equivalent circuits. For instance, lumped inductance or capacitance models have been developed. Unfortunately, these models only function well at very low frequencies. In particular, propagation and high-frequency scattering effects are not modeled accurately with this method. Therefore, this method is not suitable for modeling high speed electronic packages.

[0008] High accuracy numerical techniques have also been developed to address this problem. The Method of Moments (MoM) is an integral equation method formulated in the frequency domain, while the finite difference time domain (FDTD) is a differential equation method formulated in the time domain. Both methods can theoretically be used to accurately solve the coupling problem. Unfortunately, these techniques require massive computational time and memory overheads and in practice, have therefore been limited to modeling only small portions of an entire electronic package. Accordingly, such techniques are impractical for modeling a complete multi-conductor, multilayered package structure. Indeed, the efficiency of these methods is dramatically worse than the efficiencies of analytic or semi-analytic approaches.

[0009] One semi-analytic approach employed to account for coupling noise between coupled vias is disclosed in “Coupled Noise Analysis for Adjacent Vias in Multilayered Digital Circuits,” by Qizheng Gu, et al., (IEEE Transactions on Circuits and Systems, Vol. 41, No. 12, December 1994). In this reference, coupling between two adjacent vias is analyzed using equivalent magnetic frill array models. These models are combined with even-odd mode decomposition. This method assumes symmetry of the two vias and neglects the influence of other vias. In reality, it is necessary to model large-scale distributed coupling effects among many vias due to large-scale mode coupling and crosstalk. Hence, this technique is not suited to use with multiple coupled vias in a multilayered environment.

[0010] Another semi-analytic approach is employed in the Speed2000 modeling application developed by SIGRITY Inc., of Santa Clara, Calif. The Speed2000 modeling application employs an impedance transformation technique. The technique assumes that the package has an inherent input impedance associated with a via that is essentially dependent upon the radius of the via. A numerical model is created for a portion of the physical electronic packaging structure containing the via. The effective input impedance of the numerical model is then determined and transformed to essentially match the input impedance of the physical structure. A time domain technique is subsequently used for the propagation problem. This method, although fast, may provide results having excessive errors. Errors arise because the input impedance of a via is dependent not only on the radius of the via, but also on the radius of the anti-pad, and because the effective impedance can only model the via inside the numerical model and its first order interaction with other portion of the physical structure, but may not account for the multiple scattering effect of all the vias inside the physical structure.

[0011] It would therefore be desirable to develop an electronic package modeling technique that can accurately model multiple coupled vias in a multilayered environment. Such a technique should account for the multiple scattering effect of all the vias inside the physical structure, as well as accurately modeling both interior and exterior conditions.

[0012] In accordance with the present invention, a method is provided for analyzing signal reflection, transmission, and coupling between vias in multilayered electronic packaging structures. A physical electronic packaging structure having at least two separated metallic planes and a via is assumed. Each via passes through apertures in the ground plane and in the power plane. These two planes form a parallel-plate waveguide. This package has an inherent admittance and scattering matrix associated with the via and the layer and is independent of the exterior trace.

[0013] By using the equivalence principle, each aperture can be replaced by a perfect electric conductor, with equal and opposite magnetic currents on the two sides of the aperture. Thus, a model of the electronic structure can be separated into an exterior problem related to traces with magnetic surface currents, and an interior problem related to magnetic surface currents; each problem can be solved separately. After the solutions of the exterior and interior problems are computed, the two results are related by imposing the condition of equal and opposite magnetic currents.

[0014] A numerical model of a portion of the physical electronic packaging structure containing the via is constructed. The admittance matrix is determined by solving the interior problem, and a scattering matrix is obtained by solving the exterior problem. The admittance matrix permits the exterior structure to be coupled into the model with the interior structure, by relating the port voltage and current, while the scattering matrix enables the exterior structure to be coupled into the simulation by relating the port incident and outgoing voltage wave.

[0015] For a single via, the interior problem is solved by first applying the dyadic Green's function, expressed in terms of vector cylindrical waves and waveguide modes, to determine the exciting waves. Then, the Foldy-Lax equation is used to relate the exciting waves, the incident waves, and the scattering waves. The waveguide modes are decoupled in the Foldy-Lax equations. This approach enables the admittance matrix of the entire interior geometry to be obtained. The currents on the via surfaces, essential for circuit simulation, are obtained by integrating the magnetic field over the cylinder surface. The equivalent circuit can be integrated into prior art circuit solvers, such as a Simulation Program with Integrated Circuit Emphasis (SPICE) solver. The cylindrical wave expansion of Green's function is derived directly from Maxwell's Equations and includes all the modes of the given structure, thereby facilitating accurate simulations.

[0016] The exterior problem is solved by first coupling the exterior problem into the numerical model using the admittance matrix to relate the port voltage and current. Next, the MoM and the matrix pencil method are applied. The solution of the exterior problem is the scattering matrix.

[0017] Multilayered structures can be simulated by using the technique of transfer matrix. The numerical model is then further developed to account for more complicated structures that contain from two to hundreds (or thousands of vias). The present invention preferably uses the Foldy-Lax equations for solving the multiscattering problem. This set of equations relates the incident wave to the exciting wave from the source and to all the scattered waves from other cylinders, thus accurately accounting for interactions among all of the vias. The simulation paradigm developed in the present invention yields accurate results that are essentially independent of either the number of vias involved or the effective radii of via pins in the package. The invention can be used to rapidly and accurately simulate the coupling effects of hundreds of massively-coupled vias in a densely packaged circuit, and in printed circuit boards.

[0018] The admittance matrix solution from the interior problem and the scattering matrix from the external problem are used to represent the electromagnetic interactions in this multilayered, multi-via geometry. These matrices can easily be modified to integrate directly into commercial electronic package simulation tools. Also, the admittance and scattering matrices for an n-port circuit network can be easily integrated into a system simulation by using microwave network theory.

[0019] For electronic packages having vias with small radii, only a single set of harmonics of the vias are coupled to each other, significantly reducing the complexity of the Foldy-Lax equation portion of the model. Additional simplifications to the Foldy-Lax equation portion of the model can be made if the separation distance between layers is small, thereby reducing the number of possible modes propagating from one cylindrical via to another cylindrical via.

Overview of the Present Invention

[0042] FIG. 1 schematically illustrates a multilayered electronic package 10 with multiple parallel, or nearly parallel, conducting ground and power planes that are separated by dielectric substrates. A plurality of via structures 12 are clearly shown connecting together traces 14 residing in/on different layers or planes. Each via passes through a corresponding aperture 16. The layered geometry includes matched ground and power planes 18a and 18b, defined by metallic layers separated by dielectric layers.

[0043] Consider the problem of modeling the dynamics of traces 14 connected to vias 12. Each via 12 passes through a corresponding aperture 16 in ground/power planes 18a and 18b of FIG. 1. For each matched pair of ground/power planes, the paired planes will form a parallel-plate waveguide (i.e., ground/power plane 18a forms a first parallel-plate waveguide, and ground/power plane 18b forms a second parallel-plate waveguide). By using the equivalence principle, each aperture can be replaced by a perfect electric conductor (PEC), with equal and opposite magnetic currents on the two sides of the aperture. Note that because the apertures pass through each ground/power plane, and each ground/power plane includes a conducting plane, by extension each ground/power plane can be treated as a PEC, which results in an exterior problem related to traces with magnetic surface currents, as indicated in FIG. 2, and an interior problem with magnetic surface currents, as indicated FIG. 3. The surface currents correspond to the portions of each ground/power plane through which apertures are formed and enable vias to pass. In each of FIGS. 2 and 3, the apertures have been replaced by a PEC with an equivalent magnetic current.

[0044] FIG. 2 shows trace 14, a PEC 20, and ground/power plane 18a, while FIG. 3 shows ground/power planes 18a and 18b, and multiple PECs 20. FIG. 4 illustrates the voltages and currents of the interior and exterior problems, showing an equal and opposite magnetic current for each problem. Thus, with the apertures replaced by PECs, a model having parallel-plate waveguides with cylindrical vias is achieved. The exterior problem and the interior problem can be solved separately, and after the solutions of the exterior and interior problems are completed, the two results are related by imposing the condition of equal and opposite magnetic currents. FIG. 3 also illustrates parameters to be used when executing the method of the present invention, including a via inner diameter (indicated as 3a), a via outer diameter (indicated as 3b), and a separation distance between the parallel layers (indicated as d). It should be noted that the dielectric constant of the material of the electronic package (i.e., of the material from which each parallel plane is fabricated) is preferably included as one of the parameters used in executing the method of the present invention. The dielectric constant is factored into the following equations because the value for k (the wave number) utilized in the following equations is dependent on the dielectric constant.

[0045] FIG. 5 is a flow chart showing the sequence of logical steps employed by the present invention to analyze an electronic package like that just described. A detailed mathematical description of the analysis, as well as mathematical examples of each step, follow a broad description of the process set forth below. It should be noted that a key element of the present invention is the solution of the interior problem by manipulating a via parameter, via distribution parameters, and the layering parameter, and simulating the multiple coupling effects among a large number of vias.

[0046] Referring to FIG. 5, in a block 30, the parameters relating to the vias and the parallel planes are determined and input into the model used for evaluating via interactions in accord with the present invention. For the vias, the parameters that are input include the inner and outer diameters, while for the two parallel planes, the parameters include the separation distance of the planes, and material properties of the planes that affect wave propagations. As described above, the present invention enables signal reflection, signal transmission, and signal coupling between vias in multilayered electronic packaging structures to be determined and evaluated. By employing the equivalence principle, the coupling among traces with many vias is decomposed into interior and exterior problems. A block 31, which includes blocks 32 and 34, corresponds to the steps employed to solve the interior problem, while a block 33, which includes blocks 40 and 42, corresponds to the steps employed to solve the exterior problem.

[0047] Referring now to block 31, the first component of solving the interior problem is to applying the dyadic Green's function, expressed in terms of vector cylindrical waves and waveguide modes, to get the exciting waves, as indicated in a block 32. The cylindrical shape of vias, as illustrated in FIG. 6, enables the present invention to utilize the vector cylindrical wave expansion of dyadic Green's function, which makes the simulation fast and easy to use. This characteristic is one key benefit of the present invention.

[0048] Next, as indicated in a block 34, the Foldy-Lax equation is used to relate the exciting waves, the incident waves, and the scattering waves for the interior problem (i.e., the parallel plate waveguide). Thus, the waveguide modes are decoupled in the Foldy-Lax equation, which enables the admittance matrix of the entire interior geometry to be obtained. Indeed, the solution to the interior problem is the admittance matrix.

[0049] Note that the currents on the via surfaces, essential for circuit simulation, can be obtained by integrating the magnetic field over the cylinder surface. As indicated by an optional block 38, the equivalent circuit can be input into one of the SPICE solvers if desired. As those of ordinary skill in the art will recognize, SPICE is one of the most widely used computer programs for analog circuit analysis and is a general-purpose circuit simulation program for nonlinear direct current (DC), nonlinear transient, and linear alternating current (AC) analyses. Circuits that can be analyzed with SPICE may contain resistors, capacitors, inductors, mutual inductors, independent voltage and current sources, different types of dependent sources, lossless and lossy transmission lines, switches, and common semiconductor devices, including diodes, bipolar junction transistors (BJT), junction field effect transistors (JFET), metal semiconductor field effect transistors (MESFET), and metal oxide semiconductor field effect transistors (MOSFET). SPICE was originally developed by the University of California at Berkeley, and many commercial circuit analysis software packages derived from SPICE are available. One embodiment of the present invention uses the Foldy-Lax equations to solve the interior problem, which enables equivalent circuits to be developed for use with a SPICE solver.

[0050] Referring once again to FIG. 5, as noted above, block 33 corresponds to the steps required to solve the exterior problem, and includes blocks 40 and 42. The exterior problem is solved by constructing a numerical model of a portion of the physical electronic packaging structure containing the via. The admittance matrix, determined in block 34, permits the exterior structure to be coupled into the simulation with the interior structure, by equating the port voltage and current of the exterior and interior problems, as indicated in block 40. The exterior problem of trace and bent trace is then solved using the MoM and matrix pencil method, as indicated in block 42. The solution of the exterior problem is the scattering matrix. Note that the scattering matrix is not the same as the scattering waves solved for using the Foldy-Lax equation, as described above with respect to the interior problem. The scattering waves for the interior problem refer only to the waves within the parallel plate waveguide. With respect to the scattering matrix (i.e., the solution of the exterior problem), the waves correspond to those waves along the trace, through the via and out at the opposite end of the via.

[0051] The admittance matrix solution from the interior problem, and the scattering matrix from the external problem, are used to represent the electromagnetic interactions in this multilayered multi-via geometry. These matrices can easily be modified for integration directly into commercial electronic package simulation tools. The admittance and scattering matrices for an n-port circuit network can also easily be integrated into a system simulation by using microwave network theory, as indicated in an optional block 44.

[0052] Multilayered structures can also be simulated by using the technique of transfer matrix. The numerical model can be further developed to account for more complicated structures that contain from two vias to hundreds (even thousands) of vias. An important technique of the present invention is the use of the Foldy-Lax equations for solving the multiscattering problem. This set of equations relates the incident wave to the exciting wave from the source and to the scattered waves from other cylinders, thus accurately accounting for interactions among all of the vias. The simulation paradigm developed in the present invention yields accurate results that are essentially independent of either the number of vias involved or the effective radius of via pins in the package. The invention can be used to rapidly and accurately simulate the coupling effects of hundreds of massively-coupled vias in a densely packaged circuit.

Cylindrical-Wave Expansion of Green's Function in Parallel Wavegide Modes

[0053] FIG. 7 provides additional detail about how the model of the present invention employs the dyadic Green's function, as mentioned above. Two PECs, at z=d/2 and z=−d/2 are shown, corresponding to power/ground planes 18a and 18b. A source 50 and a field point 52 are defined as shown, with the source corresponding to the voltage at aperture 16, which can be represented by a magnetic current ring source. Accordingly, the source point will be either at z=d/2 or z=−d/2. As illustrated, source 50 is at z=d/2. Field point 52 can be located anywhere inside the waveguide. For the dyadic Green's function of a parallel plate waveguide, the boundary condition is satisfied on the metal plate by forcing the tangential electric field to equal zero, which results in waveguide modes. For the dyadic Green's function with cylinder post, the previously obtained parallel plate waveguide Green's function will serve as the primary wave, and a boundary condition is satisfied on the cylindrical surface by forcing the tangential electric field to equal zero.

[0054] The excitation source of via 12 is the voltage at a port 54 (i.e., ports correspond to apertures in power/ground planes through which vias are directed), as in FIG. 8. The voltage is expressed in magnetic ring source and generates the magnetic field inside the structure. By integrating the resulting magnetic field over the cylinder corresponding to via 12, the current induced by the voltage source is obtained. The admittance/impedance matrix of the single via geometry can then be calculated and used in the circuit simulation, as shown in FIG. 9. From the admittance/impedance matrix, the equivalent circuit of the structure can be extracted, as shown in FIG. 10. Referring to FIG. 10, Y1=Y11+Y12;Y2/2=−Y12; −Y12=−Y21; where Y11, Y12, Y21, Y22 are admittance matrix elements. These admittance matrix elements can be calculated by using the admittance matrix of the interior problem.

[0055] FIG. 11 shows the propagation characteristics of a through-hole via, single-layered via, and multilayered via. The reflection increases with the layering because electromagnetic energy is coupled to the waveguide modes.

[0056] A more exhaustive mathematical description of the use of the dyadic Green's function in the present invention follows, with reference to the multi-via model illustrated in FIG. 3. Consider two perfect electric conductors at z=d/2 and z=−d/2 (noting that the apertures normally present in power/ground planes 18a and 18b have been replaced with a PEC, such that each power/ground plane can be considered a PEC). For the multi-via problem, the equivalence principle is invoked, with equivalent magnetic sources at z′=±d/2. These sources correspond to the original locations of the apertures in each power/ground plane, and the sources are disposed at each PEC 20 identified in FIG. 3.

[0057] Consider the general case of any magnetic source at z′ between the two PECs corresponding to power/ground planes 18a and 18b (i.e., between any PEC 20 of FIG. 3). As noted above, the present invention utilizes the dyadic Green's function to represent the magnetic field. The primary Green's function {double overscore (G)}_P is the direct wave from the source with multiple reflections by the parallel plates. The solution is expressed in terms of waveguide modes and also in terms of vector cylindrical waves for adapting to via scattering later. The notations employed herein are based on those previously employed in publications relating to scattering of electromagnetic waves (The Theory of Microwave Remote Sensing, Wiley Interscience, 1985; L. Tsang, J. A. Kong and R. T. Shin; and Scattering of Electromagnetic Waves, Volume 1; Wiley Interscience, 2000, L. Tsang, J. A. Kong and K. H. Ding and Scattering of Electromagnetic Waves, Volume 2; Wiley Interscience, 2001, L. Tsang, J. A. Kong, K. H. Ding and C. Ao), with the exp(-iωt) dependence in the three references changed to exp(jωt) dependence herein. Let

[0058] ${\stackrel{\_}{\stackrel{\_}{I}}}_t=\hat{x}\hat{x}+\hat{y}\hat{y}$
be the transverse dyad. Then

[0059] $\begin{array}{cc}\stackrel{\_}{H}=-j\omega \varepsilon \int \int d{x}^{\prime }d{y}^{\prime }\stackrel{\_}{\stackrel{\_}{G}}\left(\stackrel{\_}{r},{\stackrel{\_}{r}}^{\prime }\right)·{\stackrel{\_}{M}}_s\left({\stackrel{\_}{r}}^{\prime }\right)& \left(1\right)\end{array}$
where {double overscore (G)} is the dyadic Green's function and {overscore (M)}_s is the magnetic surface current.

[0060] The waveguide mode solutions are

[0061] $\begin{array}{cc}{k}_{\mathrm{zl}}=\frac{l\pi }{d}& \left(2\right)\end{array}$
with l=0, 1, 2, . . . for TM modes and l=1, 2, . . . for TE modes.

[0062] Magnetic TE and TM modal solutions are defined as follows:

[0063] $\begin{array}{cc}\mathrm{Rg}{\stackrel{\_}{H}}_n^{\mathrm{TE}}\left(k_{\rho },k_z,\stackrel{\_}{\rho },z\right)=\frac{e^{-jn\varphi }}{\eta }\left\{\begin{array}{c}-\hat{\rho }\frac{j{k}_{\rho }{k}_2}{k}{J}_n^{\prime }\left(k_{\rho }\rho \right)\cos \left(k_{z}z\right)\\ -\hat{\varphi }\frac{{\mathrm{nk}}_z}{k\rho }{J}_n\left(k_{p}p\right)\cos \left(k_{z}z\right)-\hat{z}\frac{k_{\rho }^2}{k}{J}_n\left(k_{\rho }\rho \right)j\sin \left(k_{z}z\right)\end{array}\right\}& \left(3\right)\\ \mathrm{Rg}{\stackrel{\_}{H}}_n^{\mathrm{TM}}\left(k_{\rho },k_z,\stackrel{\_}{\rho },z\right)=\frac{j}{\eta }\left\{-\hat{\rho }\frac{jn}{\rho }{J}_n\left(k_{\rho }\rho \right)-\hat{\varphi }{k}_{\rho }{J}_n^{\prime }\left(k_{\rho }\rho \right)\right\}{e}^{-jn\varphi }\cos \left(k_{z}z\right)& \left(4\right)\end{array}$
where Rg stands for regular, with Bessel functions being used. Without the Rg, the corresponding functions are Hankel functions of the second kind. Also, k_p=√{square root over (k²−k_z²)}. The two vector cylindrical wave functions are

[0064] $\begin{array}{cc}{\stackrel{\_}{m}}_n\left(k_{\rho },k_z,\stackrel{\_}{\rho }\right)=-\hat{\rho }\frac{jn}{\rho }{H}_n^{\left(2\right)}\left(k_{\rho }\rho \right)-\hat{\varphi }{k}_{\rho }{H}_n^{{\left(2\right)}^{\prime }}\left(k_{\rho }\rho \right)\mathrm{and}{\stackrel{\_}{n}}_n\left(k_{\rho },k_z,\stackrel{\_}{\rho }\right)=-\hat{\rho }\frac{j{k}_{\rho }{k}_z}{k}{H}_n^{{\left(2\right)}^{\prime }}\left(k_{\rho }\rho \right)& \left(5\right)\\ -\hat{\varphi }\frac{{\mathrm{nk}}_z}{k\rho }{H}_n^{\left(2\right)}\left(k_{\rho }\rho \right)+\hat{z}\frac{k_{\rho }^2}{k}{H}_n^{\left(2\right)}\left(k_{\rho }\rho \right).& \left(6\right)\end{array}$

[0065] In terms of the vector cylindrical waves expressed in modal solutions, the primary dyadic Green's functions are expressed below. Note that the upper sign is employed for ż>z′, and the lower sign for z<z′.

[0066] For ρ<ρ′,

[0067] $\begin{array}{cc}\begin{array}{c}{\stackrel{\stackrel{\_}{\_}}{G}}_P\left(\stackrel{\_}{r},{\stackrel{\_}{r}}^{\prime }\right)·{\stackrel{\stackrel{\_}{\_}}{I}}_t=-\frac{\eta }{4d}\sum _{n,l}\frac{{\left(-1\right)}^{n+1}}{k_{\rho l}^2}{f}_l\left(1+e^{\mp 2j{k}_{\mathrm{zl}}\left(z^{\prime }\pm d/2\right)}\right)e^{\pm 2j{k}_{\mathrm{zl}}\left(z^{\prime }\pm d/2\right)}\times \\ \mathrm{Rg}{\stackrel{\_}{H}}_n^{\mathrm{TM}}\left(k_{\rho l},k_{\mathrm{zl}},\stackrel{\_}{\rho },z\mp d/2\right){\stackrel{\_}{m}}_{-n}\left(k_{\rho l},k_{\mathrm{zl}},{\stackrel{\_}{\rho }}^{\prime }\right)e^{jn{\varphi }^{\prime }}·{\stackrel{\stackrel{\_}{\_}}{I}}_t+\\ \frac{j\eta }{4d}\sum _{n,l}\frac{{\left(-1\right)}^{n+1}}{k_{\rho l}^2}{f}_l\left(1+e^{\mp 2j{k}_{\mathrm{zl}}\left(z^{\prime }\pm d/2\right)}\right)e^{\pm j{k}_{\mathrm{zl}}\left(z^{\prime }\pm d/2\right)}\times \\ \mathrm{Rg}{\stackrel{\_}{H}}_n^{\mathrm{TE}}\left(k_{\rho l},k_{\mathrm{zl}},\stackrel{\_}{\rho },z\mp d/2\right){\stackrel{\_}{n}}_{-n}\left(k_{\rho l},k_{\mathrm{zl}},{\stackrel{\_}{\rho }}^{\prime }\right)e^{jn{\varphi }^{\prime }}·{\stackrel{\stackrel{\_}{\_}}{I}}_t\end{array}\mathrm{and}\mathrm{for}\rho >{\rho }^{\prime },& \left(7\right)\\ \begin{array}{c}{\stackrel{\stackrel{\_}{\_}}{G}}_P\left(\stackrel{\_}{r},{\stackrel{\_}{r}}^{\prime }\right)·{\stackrel{\stackrel{\_}{\_}}{I}}_t=-\frac{\eta }{4d}\sum _{n,l}\frac{{\left(-1\right)}^{n+1}}{k_{\rho l}^2}{f}_l\left(1+e^{\mp 2j{k}_{\mathrm{zl}}\left(z^{\prime }\pm d/2\right)}\right)e^{\pm j{k}_{\mathrm{zl}}\left(z^{\prime }\pm d/2\right)}\times \\ {\stackrel{\_}{H}}_n^{\mathrm{TM}}\left(k_{\rho l},k_{\mathrm{zl}},\stackrel{\_}{\rho },z\mp d/2\right)\mathrm{Rg}{\stackrel{\_}{m}}_{-n}\left(k_{\rho l},k_{\mathrm{zl}},{\stackrel{\_}{\rho }}^{\prime }\right)e^{jn{\varphi }^{\prime }}·{\stackrel{\stackrel{\_}{\_}}{I}}_t+\\ \frac{j\eta }{4d}\sum _{n,l}\frac{{\left(-1\right)}^{n+1}}{k_{\rho l}^2}{f}_l\left(1+e^{\mp 2j{k}_{\mathrm{zl}}\left(z^{\prime }\pm d/2\right)}\right)e^{\pm j{k}_{\mathrm{zl}}\left(z^{\prime }\pm d/2\right)}\times \\ {\stackrel{\_}{H}}_n^{\mathrm{TE}}\left(k_{\rho l},k_{\mathrm{zl}},\stackrel{\_}{\rho },z\mp d/2\right)\mathrm{Rg}{\stackrel{\_}{n}}_{-n}\left(k_{\rho l},k_{\mathrm{zl}},{\stackrel{\_}{\rho }}^{\prime }\right)e^{jn{\varphi }^{\prime }}·{\stackrel{\stackrel{\_}{\_}}{I}}_t\end{array}& \left(8\right)\\ \begin{array}{c}{f}_l=\frac{1}{2},\mathrm{for}l=0\\ =1,\mathrm{for}l=1,2,\dots \end{array}& \left(9\right)\end{array}$
Foldy-Lax Equations for Scattering with Multiple Cylinders

[0068] Once the dyadic Green's functions have been employed to determine the exciting wave as described above, the next step corresponds to block 34 of FIG. 5, in which the Foldy-Lax equations are employed to determine the admittance matrix. Note that block 34 represents the second key step of the interior problem solution (i.e., block 31 of FIG. 5) and simulates the coupling effect among a large number of vias in layered geometry by using a Foldy-Lax multiple scattering formulation. Such multiple scatterings are schematically illustrated in FIG. 12. Each aperture 16 in paired parallel plates corresponding to power/ground planes 18a and 18b is associated with a port, such that ports 58a–58d are disposed adjacent to power/ground planes 18a, and ports 59a–59d are disposed adjacent to power/ground planes 18b. Four vias 12 also shown, along with schematic representations of multiple scatterings among such vias, but it will be understood that many more vias might be included.

[0069] The Foldy-Lax multiple scattering equations are derived using vector cylindrical waves and waveguide modes as basis functions. In active via ports, the voltage sources generate a vector cylindrical wave, which propagates along the parallel plate waveguide, and such waves are scattered by other vias. There is strong multiple scattering among the vias with slow decay over distance, because of the waveguide effect. The Foldy-Lax equations relate the excited, incident, and scattered waves of each cylinder. This set of equations describes the coupling effect among a large number of vias. It is important to recognize that the Foldy-Lax equations expressed in vector cylindrical waves are decoupled among the waveguide modes. There is no coupling between Transverse Electric (TE) and Transverse Magnetic (TM) waves, because the cylinders are assumed to be perfect electric conductors (under the model employed in the present invention). In terms of the matrix equation, this decoupling significantly reduces the dimension of the matrix equation, as one only needs to solve the Foldy-Lax equation for each waveguide mode separately.

[0070] This decoupling is useful when multiple waveguide modes are excited for waveguides of moderate thickness. For such cases, if discretization in the z direction is done using sub-sectional basis functions, a large number of basis functions are required, and they are coupled in multiple scattering.

[0071] FIG. 13 shows the equivalent model of multiple via geometry, it being understood that there can be multiple excitations at different vias. Note the similarities between FIGS. 12 and 13, particularly that FIG. 13 also illustrates paired parallel plates corresponding to power/ground planes 18a and 18b, ports 58a–58d and ports 59a–59d, and four vias 12.

[0072] Using the Foldy-Lax equations, the magnetic fields can be calculated. By integrating a field over the cylinder, the induced current on each cylinder is obtained. The entire structure can be modeled as a multi-port microwave network, as is schematically illustrated in FIG. 14. The currents and voltages of the interior problem are related as follows:

[0073] $\begin{array}{cc}{I}_{\mathrm{iu}}=\sum _{j=1}^N{Y}_{\mathrm{ij}}^{\mathrm{UU}}{V}_{\mathrm{ju}}+\sum _{j=1}^N{Y}_{\mathrm{ij}}^{\mathrm{UB}}{V}_{\mathrm{jb}}{I}_{\mathrm{ib}}=\sum _{j=1}^N{Y}_{\mathrm{ij}}^{\mathrm{BU}}{V}_{\mathrm{ju}}+\sum _{j=1}^N{Y}_{\mathrm{ij}}^{\mathrm{BB}}{V}_{\mathrm{jb}}& \left(10\right)\end{array}$

[0074] FIG. 15 graphically illustrates the coupling effect between two vias. In the prior art, a coupling model based on a capacitor plate antenna has been developed. FIG. 15 compares the results obtained using the techniques of the present invention (referred to in FIG. 15 as the Foldy-Lax Approach) with those obtained using the capacitor plate antenna model of the prior art (referred to in FIG. 15 as the Even-Odd Mode Approach). The results are in good agreement, indicating that the techniques of the present invention are at least as accurate as prior art techniques and can be more easily executed.

[0075] FIG. 16 is a plan view of power/ground plane 18a, showing the geometry of a large number of adjacent vias, some of which are active vias 62, and some of which are passive vias 60. The objective is to obtain the scattering matrix of the two active vias in the presence of many passive vias. FIG. 17 graphically illustrates the coupling effects of passive vias. The solid lines represent the coupling effects as determined using the capacitor plate antenna model of the prior art, which ignores the effects introduced by passive vias. The dash lines represent the coupling effects introduced by passive vias, as determined using the techniques of the present invention. The graph of FIG. 17 is based on the coupling effect of two vias among 256 adjacent vias. The large difference in the coupling effects is due to the ability of the present invention to take into account the scattered waves from the adjacent vias.

[0076] A more exhaustive mathematical description of the use of the Foldy-Lax equations in solving the interior problem in accord with the present invention will now be provided. Consider N cylinders between the two parallel plates centered at {overscore (ρ)}₁, {overscore (ρ)}₂ . . . {overscore (ρ)}_N and magnetic surface current densities {overscore (M)}_s={overscore (M)}_su at ({overscore (ρ)}′, z′=d/2) and {overscore (M)}_s={overscore (M)}_sb at ({overscore (ρ)}′, z′=−d/2). The multiple scattering can be formulated in terms of Foldy-Lax multiple scattering equations. The upper sign is used for results with a source {overscore (M)}_su at ({overscore (ρ)}′, z′=d/2), and the lower sign for results with a source {overscore (M)}_sb at ({overscore (92 )}′, z′=−d/2).

[0077] The Foldy-Lax multiple scattering equations are formulated as follows. After multiple scattering, the final exciting field of cylinder p is:

[0078] $\begin{array}{cc}\begin{array}{c}{\stackrel{\_}{H}}_{\mathrm{ex}}^{\left(p\right)}=\sum _{n,l}{w}_{\mathrm{lm}}^{\mathrm{TM}\left(p\right)}\mathrm{Rg}{\stackrel{\_}{H}}_m^{\mathrm{TM}}\left(k_{\rho l},k_{\mathrm{zl}},\stackrel{\_}{\rho }-{\stackrel{\_}{\rho }}_p,z\pm d/2\right)+\\ \sum _{n,l}{w}_{\mathrm{lm}}^{\mathrm{TE}\left(p\right)}\mathrm{Rg}{\stackrel{\_}{H}}_m^{\mathrm{TE}}\left(k_{\rho l},k_{\mathrm{zl}},\stackrel{\_}{\rho }-{\stackrel{\_}{\rho }}_p,z\pm d/2\right)\end{array}& \left(11\right)\end{array}$
where w_lm^TM(p) and w_lm^TE(p) are exciting field coefficients to be solved by using the Foldy-Lax equations.

[0079] The T-matrix coefficients for perfect conducting cylinders are

[0080] $\begin{array}{cc}{T}_n^{\left(M\right)}=T_{-n}^{\left(M\right)}=\frac{J_n^{\prime }\left(k_{\rho l}a\right)}{H_n^{{\left(2\right)}^{\prime }}\left(k_{\rho l}a\right)}& \left(12\right)\\ T_n^{\left(N\right)}=T_{-n}^{\left(N\right)}=\frac{J_n\left(k_{\rho l}a\right)}{H_n^{\left(2\right)}\left(k_{\rho l}a\right)}.& \left(13\right)\end{array}$
The Foldy-Lax equations state that the final exciting field of cylinder q is equal to the incident field plus scattered fields from all cylinders, except the scattered field from itself. The scattered field from cylinder p incident on cylinder q can be re-expressed by using the translation addition theorem for the modal solution {overscore (H)}_n^TE and {overscore (H)}_n^TM.

[0081] The Foldy-Lax multiple scattering equations are in the following form:

[0082] $\begin{array}{cc}\begin{array}{c}{w}_{\ln }^{\mathrm{TM}\left(q\right)}=a_{\ln }^{\mathrm{TM}\left(q\right)}+\sum _{\underset{p\ne q}{p=1}}^N\sum _{m=-\infty }^{\infty }{H}_{n-m}^{\left(2\right)}\left(k_{\rho l}{\stackrel{\_}{\rho }}_p-{\stackrel{\_}{\rho }}_q\right)\times \\ e^{j\left(n-m\right)\varphi \stackrel{\_}{{\rho }_p}{\rho }_q}{T}_m^{\left(N\right)}{w}_{\mathrm{lm}}^{\mathrm{TM}\left(p\right)}\end{array}\begin{array}{c}{w}_{\ln }^{\mathrm{TE}\left(q\right)}=a_{\ln }^{\mathrm{TE}\left(q\right)}+\sum _{\underset{p\ne q}{p=1}}^N\sum _{m=-\infty }^{\infty }{H}_{n-m}^{\left(2\right)}\left(k_{\rho l}{\stackrel{\_}{\rho }}_p-{\stackrel{\_}{\rho }}_q\right)\times \\ e^{j\left(n-m\right)\varphi \stackrel{\_}{{\rho }_p}{\rho }_q}{T}_m^{\left(N\right)}{w}_{\mathrm{lm}}^{\mathrm{TE}\left(p\right)}\end{array}& \left(14\right)\end{array}$
where a_ln^TM(q), a_ln^TE(q) are the fields of the current source incident on cylinder q and can be calculated by using the primary Green's function of the previous section. In the Foldy-Lax equations, there is no coupling between different l because each l corresponds to a specific k_z. Neither is there coupling between the TE and TM waves, because the cylinders are perfectly conducting. Note that the primary Green's function is translationally invariant in the horizontal plane.

[0083] Using the translational addition theorem and translational invariance of {double overscore (G)}_P({overscore (r)}, {overscore (r)}′) in the horizontal direction, the field from the magnetic surface current source {overscore (M)}_s incident on cylinder q is:

[0084] $\begin{array}{cc}\begin{array}{c}-j\omega \varepsilon \int d{\stackrel{\_}{\rho }}^{\prime }{\stackrel{\stackrel{\_}{\_}}{G}}_P\left(\stackrel{\_}{r},{\stackrel{\_}{r}}^{\prime }\right)·{\stackrel{\_}{M}}_s\left({\stackrel{\_}{\rho }}^{\prime }\right)=-j\omega \varepsilon \int d{\stackrel{\_}{\rho }}^{\prime }{\stackrel{\stackrel{\_}{\_}}{G}}_P\\ \left(\stackrel{\_}{\rho }-{\stackrel{\_}{\rho }}_q,z,{\stackrel{\_}{\rho }}^{\prime }-{\stackrel{\_}{\rho }}_q,z^{\prime }\right)·\\ {\stackrel{\_}{M}}_s\left({\stackrel{\_}{\rho }}^{\prime }\right)\end{array}& \left(15\right)\end{array}$
Using the expression of the primary dyadic Green's function from Equations (7) and (8) in the above, the field coefficients incident on cylinder q are computed as

[0085] $\begin{array}{cc}{a}_{\ln }^{\mathrm{TM}\left(q\right)}=\frac{\eta \mathrm{j\omega }\varepsilon }{2d}\frac{{\left(-1\right)}^{n+1}}{k_{\rho l}^2}{f}_l\int d{\stackrel{\_}{\rho }}^{\prime }\times & \left(16\right)\\ {\stackrel{\_}{m}}_{-n}\left(k_{\rho l},k_{\mathrm{zl}},{\stackrel{\_}{\rho }}^{\prime }-{\stackrel{\_}{\rho }}_q\right)e^{jn\varphi \stackrel{\_}{{\rho }^{\prime }{\rho }_q}}·{\stackrel{\_}{M}}_s\left({\stackrel{\_}{\rho }}^{\prime }\right)& \mathrm{ }\\ =\frac{\eta \omega \varepsilon }{2d}\frac{{\left(-1\right)}^{n+1}}{k_{\rho l}^2}{f}_l\int d{\stackrel{\_}{\rho }}^{\prime }\times & \left(17\right)\\ {\stackrel{\_}{n}}_{-n}\left(k_{\rho l},k_{\mathrm{zl}},{\stackrel{\_}{\rho }}^{\prime }-{\stackrel{\_}{\rho }}_q\right)e^{jn\varphi \stackrel{\_}{{\rho }^{\prime }{\rho }_q}}·{\stackrel{\_}{M}}_s\left({\stackrel{\_}{\rho }}^{\prime }\right).& \end{array}$
After the Foldy-Lax multiple scattering equations are solved, the surface currents on the cylinder (p) are

[0086] $\begin{array}{cc}\begin{array}{c}{J}_s^{\left(p\right)}=\sum _{m,l}{w}_{\mathrm{lm}}^{\mathrm{TM}\left(p\right)}{\stackrel{\_}{J}}_s^{\mathrm{TM}}\left(k_{\rho l},k_{\mathrm{zl}},x\pm d/2\right)+\\ \sum _{m,l}{w}_{\mathrm{lm}}^{\mathrm{TE}\left(p\right)}{\stackrel{\_}{J}}_s^{\mathrm{TE}}\left(k_{\rho l},k_{\mathrm{zl}},z\pm d/2\right)\end{array}& \left(18\right)\end{array}$
where the modal surface currents are

[0087] $\begin{array}{cc}\begin{array}{c}{\stackrel{\_}{J}}_{\mathrm{sn}}^{\mathrm{TE}}\left(k_{\rho },k_z,z\right)=-\frac{2j/\eta \pi k_{\rho }a)e^{-jn\varphi }}{H_n^{{\left(2\right)}^{\prime }}\left(k_{\rho }a\right)}\times \\ [-z\frac{{\mathrm{nk}}_z}{\mathrm{ka}}\cos \left(k_{z}z\right)+\hat{\varphi }\frac{k_{\rho }^2}{k}j\sin \left(k_{z}z\right)]\end{array}& \left(19\right)\\ {\stackrel{\_}{J}}_{\mathrm{sn}}^{\mathrm{TM}}\left(k_{\rho },k_z,z\right)=\hat{z}\frac{\left(2/\eta \pi a\right)}{H_n^{\left(2\right)}\left(k_{\rho }a\right)}\left\{\cos k_{z}z\right\}{e}^{-jn\varphi }.& \left(20\right)\end{array}$
The Excitation of Magnetic Frill Currents

[0088] Consider magnetic frill currents at aperture j with {overscore (M)}_su at z′=d/2 and {overscore (M)}_sb at z′=−d/2:

[0089] $\begin{array}{cc}{\stackrel{\_}{M}}_{\mathrm{su}}\left({\stackrel{\_}{\rho }}^{\prime }\right)=-\frac{V_{\mathrm{ju}}}{{\stackrel{\_}{\rho }}^{\prime }-{\stackrel{\_}{\rho }}_j\ln \left(b/a\right)}{\hat{\varphi }}^{\prime }\frac{}{\rho {\rho }_j},\mathrm{for}a\le {\stackrel{\_}{\rho }}^{\prime }={\stackrel{\_}{\rho }}_j\le b& \left(21\right)\\ {\stackrel{\_}{M}}_{\mathrm{sb}}\left({\stackrel{\_}{\rho }}^{\prime }\right)=-\frac{V_{\mathrm{jb}}}{{\stackrel{\_}{\rho }}^{\prime }-{\stackrel{\_}{\rho }}_j\ln \left(b/a\right)}{\hat{\varphi }}^{\prime }\frac{}{\rho {\rho }_j},\mathrm{for}a\le {\stackrel{\_}{\rho }}^{\prime }={\stackrel{\_}{\rho }}_j\le b.& \left(22\right)\end{array}$
First consider the excitation terms of a_ln^TM(j) and a_ln^TE(j) for cylinder j.

[0090] Because {overscore (n)}₀(k_ρl, k_zl, {overscore (ρ)}′−{overscore (ρ)}_j)=0,

[0091] $\begin{array}{cc}{a}_{\ln }^{\mathrm{TM}\left(j\right)}=\frac{\mathrm{kj}}{d}\frac{{\left(-1\right)}^l}{k_{\rho l}^2}{f}_l\frac{\pi V_j}{\ln \left(b/a\right)}{\delta }_{\mathrm{n0}}\left[H_0^{\left(2\right)}\left(k_{\rho l}b\right)-H_0^{\left(2\right)}\left(k_{\rho l}a\right)\right]& \left(23\right)\end{array}$
a_ln^TE(j)=0  (24)

where V_j=V_ju for {overscore (M)}_su and V_j=V_jb for {overscore (M)}_sb.

[0092] Next, the excitation terms of a_ln^TM(q) and a_ln^TE(q) for cylinder q≠j are calculated. The translation addition theorem is used, since |{overscore (ρ)}′−{overscore (ρ)}_q|<|{overscore (ρ)}_j−{overscore (ρ)}_q|.

[0093] Thus, for q≠j,

[0094] $\begin{array}{cc}\begin{array}{c}{a}_{\ln }^{\mathrm{TM}\left(q\right)}=\frac{\mathrm{jk}}{2d}\frac{{\left(-1\right)}^{n+l}}{k_{\rho l}^2}{f}_l{H}_n^{\left(2\right)}\left(k_{\rho l}{\stackrel{\_}{\rho }}_q-{\stackrel{\_}{\rho }}_j\right)e^{jn\varphi \stackrel{\_}{{\rho }_q{\varrho }_j}}\times \\ \left[\frac{2\pi V_j}{\ln \left(b/a\right)}\left[J_0\left(k_{\rho l}b\right)-J_0\left(k_{\rho l}a\right)\right]\right]\end{array}& \left(25\right)\end{array}$
After the exciting field coefficients are calculated, the surface current density on cylinder p is calculated from the above relationship. Then the current on cylinder (p) is:

[0095] $\begin{array}{cc}\begin{array}{c}{I}^{\left(p\right)}={\int }_0^{2\pi }ad\varphi {\stackrel{\_}{J}}_s^{\left(p\right)}·\hat{z}\\ =2\pi a\sum _l{w}_{\mathrm{l0}}^{\mathrm{TM}\left(p\right)}\left(\hat{z}\right)k_{\rho l}\frac{\left(2/\pi k_{\rho l}a\right)}{H_0^{\left(2\right)}\left(k_{\rho l}a\right)}\frac{1}{\eta }\left\{\cos k_{\mathrm{zl}}\left(z\pm d/2\right)\right\}\end{array}& \left(26\right)\end{array}$
Next, the general case is extended, with N cylinders having voltages V_1u, V_2u, . . . , V_Nu, V_1b, V_2b, . . . and V_Nb. For sources at the upper aperture of z′=d/2, that is, V_1u, V_2u, . . . , V_Nu, the incident wave at cylinder q is:

[0096] $\begin{array}{cc}\begin{array}{c}{a}_{\ln }^{\mathrm{TM}\left(q\right)}=\frac{\mathrm{jk}}{2d}\frac{{\left(-1\right)}^{n+l}}{k_{\rho l}^2}{f}_l\frac{2\pi V_{\mathrm{qu}}}{\ln \left(b/a\right)}{\delta }_{\mathrm{n0}}\times \\ \left[H_0^{\left(2\right)}\left(k_{\rho l}b\right)-H_0^{\left(2\right)}\left(k_{\rho l}a\right)\right]+\sum _{j\ne q}^N\frac{\mathrm{jk}}{2d}\frac{{\left(-1\right)}^{n+l}}{k_{\rho l}^2}{f}_l{H}_n^{\left(2\right)}\\ \left(k_{\rho l}{\stackrel{\_}{\rho }}_q-{\stackrel{\_}{\rho }}_j\right)e^{jn\varphi \stackrel{\_}{{\rho }_q{\varrho }_j}}\times \\ \left[\frac{2\pi V_{\mathrm{ju}}}{\ln \left(b/a\right)}\left[J_0\left(k_{\rho l}b\right)-J_0\left(k_{\rho l}a\right)\right]\right].\end{array}& \left(27\right)\end{array}$
The Foldy-Lax equations are solved for w_ln^TM(q) and the currents are found at the cylinders, at z=+d/2. This current, I^uu, is the current at z=d/2 due to the source at z′=d/2:

[0097] $\begin{array}{cc}{I}^{\left(p\right)\mathrm{uu}}=\sum _l{w}_{\mathrm{l0}}^{\mathrm{TM}\left(p\right)}\frac{4}{\eta H_0^{\left(2\right)}\left(k_{\rho l}a\right)}\left(-1\right).& \left(28\right)\end{array}$
The currents are then found at the cylinders, at z=−d/2. This current, I^bu, is the current at −z=−d/2 due to the source at z′=d/2:

[0098] $\begin{array}{cc}{I}^{\left(p\right)\mathrm{bu}}=\sum _l{w}_{\mathrm{l0}}^{\mathrm{TM}\left(p\right)}\frac{4}{\eta H_0^{\left(2\right)}\left(k_{\rho l}a\right)}.& \left(29\right)\end{array}$
Next, consider the sources at the lower aperture of −z′=d/2, that is, V_1b, V_2b, . . . , and V_Nb. The Foldy-Lax equations are the same for this aperture.

[0099] Then the currents are found at the cylinders, at z=+d/2. This current, I^bb, is the current at the lower aperture z=−d/2 due to the sources of the lower aperture at z′=−d/2.

[0100] The total currents are
I^(p)u=I^(p)uu+I^(p)ub  (30)
I^(p)b=I^(p)bu+I^(p)bb.  (31)
Matrix Notation for the Interior Problem

[0101] Suppose that there are N vias, that l=L_max, and that multipoles are limited to n=+N_max. Note that different ls are uncoupled. Then the dimension of {overscore (w)}_l is M=N×(2N_max+1), when using a combined particle index p=1, 2, . . . , N and a multi-pole index n=−N_max, −N_max+1, . . . , 0, 1, N_max. A combined index is then α−(q, n)−(2N_max+1)×(q−1). Thus, α=1, 2, . . . , M. Let superscript T denote the transpose of a matrix:
{overscore (w)}_l^T=[w_l(−Nmax)⁽¹⁾w_l(−Nmax+1)⁽¹⁾. . . w_l(Nmax)⁽¹⁾w_l(−Nmax)⁽²⁾. . . w_l(Nmax)⁽²⁾. . . w_l(Nmax)^(N)].  (32)
Thus, the dimension of the Foldy-Lax equations is M, and it is independent of the number of modes. That is, it is independent of the number of basis functions in the z-direction, if modes are used. On the other hand, if N_z subsectional basis functions are used in the z-direction, then the dimension of the matrix equation is magnified N_z times. This decoupling drastically reduces the size of the matrix equation, particularly at high frequencies.

[0102] The Foldy-Lax matrix {double overscore (F)}_l is of dimension N(2N_max+1)×N(2N_max+1)=M×M:

[0103] $\begin{array}{cc}{\left(\left[{\stackrel{\_}{\stackrel{\_}{F}}}_l\right]\right)}_{\mathrm{qn},\mathrm{pm}}={\delta }_{\mathrm{nm}}{\delta }_{\mathrm{qp}}-\left(1-{\delta }_{\mathrm{pq}}\right)H_{n-m}^{\left(2\right)}\times \left(k_{\rho l}{\stackrel{\_}{\rho }}_p-{\stackrel{\_}{\rho }}_q\right)e^{j\left(n-m\right)\varphi {\stackrel{\_}{\rho }}_p{P}_q}{T}_m^{\left(N\right)}.& \left(33\right)\end{array}$
Define {double overscore (E)}_l to be of dimension M×N, so that

[0104] $\begin{array}{cc}{\left({\stackrel{\_}{\stackrel{\_}{E}}}_l\right)}_{\mathrm{qn},j}={\left(\left[{\stackrel{\_}{a}}_{\mathrm{lj}}\right]\right)}_{\mathrm{qn}}& \left(34\right)\end{array}$
where

[0105] $\begin{array}{cc}\begin{array}{c}{\left(\left[{\stackrel{\_}{a}}_{\mathrm{lj}}\right]\right)}_{\mathrm{qn}}=\frac{\mathrm{jk}}{2d}\frac{{\left(-1\right)}^l}{k_{\rho l}^2}{f}_l\frac{2\pi }{\ln \left(b/a\right)}{\delta }_{\mathrm{n0}}\times \\ \left[H_0^{\left(2\right)}\left(k_{\rho l}b\right)-H_0^{\left(2\right)}\left(k_{\rho l}a\right)\right],\text{if}q=j\\ =\frac{\mathrm{jk}}{2d}\frac{{\left(-1\right)}^{n+l}}{k_{\rho l}^2}{f}_l{H}_n^{\left(2\right)}\left(k_{\rho l}{\stackrel{\_}{\rho }}_q-{\stackrel{\_}{\rho }}_j\right)e^{jn\varphi \stackrel{\_}{{\rho }_q{\varrho }_j}}2\frac{\pi }{\ln \left(b/a\right)}\times \\ \left[J_0\left(k_{\rho l}b\right)-J_0\left(k_{\rho l}a\right)\right],\text{if}q\ne j.\end{array}& \left(35\right)\end{array}$
Let voltage vectors {overscore (V)}^u be of dimension N, with ({overscore (V)}^u)_j=V_ju.

[0106] In matrix form, the Foldy-Lax equation, for sources at z′=d/2, is:

[0107] $\begin{array}{cc}{\stackrel{\_}{\stackrel{\_}{F}}}_l{\stackrel{\_}{w}}_l^{\left(u\right)}={\stackrel{\_}{\stackrel{\_}{E}}}_l{\stackrel{\_}{V}}^u.& \left(36\right)\end{array}$
After the Foldy-Lax equation is solved, the currents can be computed. Let {double overscore (P)}_l be of dimension N×N:

[0108] $\begin{array}{cc}{\left({\stackrel{\_}{\stackrel{\_}{P}}}_l\right)}_{p,j}={\left({\stackrel{\_}{\stackrel{\_}{F}}}_l^{-1}{\stackrel{\_}{\stackrel{\_}{E}}}_l\right)}_{\mathrm{p0},j}.& \left(37\right)\end{array}$
The p0,j index is due to the fact that only the n=0th harmonic contributes to the current at the terminal, as indicated by Equation 26. Next, define current vectors of dimension N,Ī^uu, corresponding to the current at the upper aperture, at z=d/2, due to the sources at z′=d/2. Also, define current vectors of dimension N,Ī^bu, which is the current at the lower aperture, at z=−d/2, due to the sources at z′=d/2. Then,

[0109] $\begin{array}{cc}{\stackrel{\_}{I}}^{\mathrm{uu}}={\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{uu}}{\stackrel{\_}{V}}^u& \left(38\right)\\ {\stackrel{\_}{I}}^{\mathrm{bu}}={\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{bu}}{\stackrel{\_}{V}}^u& \left(39\right)\end{array}$
where {double overscore (Y)}^uu and {double overscore (Y)}^bu are of dimension N×N:

[0110] $\begin{array}{cc}{\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{uu}}=\sum _l{B}_l{\stackrel{\_}{\stackrel{\_}{P}}}_l& \left(40\right)\\ {\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{bu}}=\sum _l{D}_l{\stackrel{\_}{\stackrel{\_}{P}}}_l\text{with}& \left(41\right)\\ B_l=\frac{4{\left(-1\right)}^l}{\eta H_0^{\left(2\right)}\left(k_{\rho l}a\right)}& \left(42\right)\\ D_l=\frac{4}{\eta H_0^{\left(2\right)}\left(k_{\rho l}a\right)}.& \left(43\right)\end{array}$
Similarly, define current vectors of dimension N,Ī^ub, Which is the current at the upper aperture, at z=d/2 due to the sources at z′=d/2. Also, define current vectors of dimension N, Ī^bb, which is the current at the lower aperture, at z=−d/2 due to the sources at z′=−d/2.

[0111] Also, let voltage vectors {overscore (V)}^b be of dimension N. Then,

[0112] $\begin{array}{cc}{\stackrel{\_}{I}}^{\mathrm{ub}}={\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{ub}}{\stackrel{\_}{V}}^b& \left(44\right)\\ {\stackrel{\_}{I}}^{\mathrm{bb}}={\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{bb}}{\stackrel{\_}{V}}^b.& \left(45\right)\end{array}$
Because of reflection symmetry of the waveguide about z=0, the symmetry relation of the admittance matrix elements is {double overscore (Y)}^ub={double overscore (Y)}^bu and {double overscore (Y)}^bb={double overscore (Y)}^uu.

[0113] After the interior problem is solved, the solution is represented by

[0114] $\begin{array}{cc}\left[\begin{array}{c}{\stackrel{\_}{I}}^u\\ -{\stackrel{\_}{I}}^b\end{array}\right]=-\stackrel{\_}{\stackrel{\_}{Y}}\left[\begin{array}{c}-{\stackrel{\_}{V}}^u\\ {\stackrel{\_}{V}}^b\end{array}\right]& \left(46\right)\end{array}$
where {double overscore (Y)} is a 2N×2N matrix:

[0115] $\begin{array}{cc}\stackrel{\_}{\stackrel{\_}{Y}}=\left[\begin{array}{cc}{\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{uu}}& -{\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{ub}}\\ -{\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{bu}}& {\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{bb}}\end{array}\right]=\left[\begin{array}{cc}{\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{uu}}& -{\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{ub}}\\ -{\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{ub}}& {\stackrel{\_}{\stackrel{\_}{Y}}}^{\mathrm{uu}}\end{array}\right].& \left(47\right)\end{array}$
Solving the Exterior Problem

[0116] As discussed above, once the admittance matrix is determined by employing the Foldy-Lax equations as described in detail above, those of ordinary skill will recognize that an equivalent circuit can readily be determined. Optionally, that equivalent circuit can be used as input in SPICE solver applications, as noted above.

[0117] The next step in the present invention is to solve the exterior problem, as indicated in block 33 of FIG. 5. The solver for the exterior problem uses MoM and the matrix pencil technique to extract the S parameters of the exterior structures. FIG. 18 shows the exterior structure for the case of a wire. The exterior problem is decomposed into a short circuit problem and a wire antenna problem. The short circuit problem relates the positive and negative waves with the short circuit current at the via port; the wire antenna problem relates the positive and negative waves with the voltage at the via port. The combined problem gives the relation of the port voltage and current as follows:
B₁=Γ_SC1A₁+T_ant1V₁  (48)
I₁=I_SC1A₁+Y_ant1V₁
Combining the Interior & Exterior Problems to Obtain the Overall Scattering Matrix

[0118] FIG. 19 schematically shows the interior and exterior problem connected. The interior and exterior problems are combined by equating the port voltages and currents of the two problems. FIG. 19 includes power/ground plane 18a, power/ground plane 18b, vias 12, 12a and 12b, apertures 16, and traces 14. Note that as shown in FIG. 19, the coupling effects between vias 12a and 12b are being modeled. Note also that each aperture corresponds to a port, as discussed above. For each port, there are four unknowns: A, B, V, and I. Two equations derive from the exterior problem, and one equation from the interior problem. The fourth equation can derive from the connection between this port and other ports. If this port is not connected with other ports, then it is an input port or a matched output port. In either case, there are only three unknowns: B, V, and I for this port. By solving a set of matrix equations, the relations between the incident wave A_i and outgoing wave B_i can be obtained, yielding the scattering matrix of the whole structure. FIG. 20 graphically illustrates the coupling effect between two active vias in the presence of 256 passive vias, obtained by combining the exterior and interior problems in accord with the present invention.

[0119] A more exhaustive mathematical description of solving the exterior problem, in accord with the present invention, will now be provided. For a single via 12 in FIG. 4, the exterior problem for z>d/2 has a magnetic current opposite that of the interior problem, but with the same voltage

[0120] $\begin{array}{cc}{\stackrel{\_}{M}}_s^{\left(e\right)}\left({\stackrel{\_}{\rho }}^{\prime }\right)=\frac{V_{1b}}{{\stackrel{\_}{\rho }}^{\prime }-{\stackrel{\_}{\rho }}_1\ln \left(b/a\right)}{\hat{\varphi }}^{\prime }\frac{}{\rho {\rho }_1},\text{for}a\le {\stackrel{\_}{\rho }}^{\prime }-{\stackrel{\_}{\rho }}_1\le b.& \left(49\right)\end{array}$
Defining the current I_1u and I_1b as flowing out of the layer produces

[0121] $\begin{array}{cc}\left[\begin{array}{c}{B}_{1u}\\ I_{1u}\end{array}\right]=\left[\begin{array}{cc}{\Gamma }_{\mathrm{sc},1u}& T_{\mathrm{ant},1u}\\ I_{\mathrm{sc},1u}& Y_{\mathrm{ant},1u}\end{array}\right]\left[\begin{array}{c}{A}_{1u}\\ -V_{1u}\end{array}\right].& \left(50\right)\end{array}$

[0122] For the exterior problem of the lower port, a magnetic current is provided that is opposite the interior problem. To maintain the same matrix form of Equation 50, the current will be flowing out of the waveguide slab, and is equal to −I_1b:

[0123] $\begin{array}{cc}{\stackrel{\_}{M}}_s^{\left(e\right)}\left({\stackrel{\_}{\rho }}^{\prime }\right)=\frac{V_{1b}}{{\stackrel{\_}{\rho }}^{\prime }-{\stackrel{\_}{\rho }}_1\ln \left(b/a\right)}{\hat{\varphi }}^{\prime }\frac{}{\rho {\rho }_1}\text{for}a\le {\stackrel{\_}{\rho }}^{\prime }-{\stackrel{\_}{\rho }}_1\le b& \left(51\right)\\ \left[\begin{array}{c}{B}_{1b}\\ I_{1b}\end{array}\right]=\left[\begin{array}{cc}{\Gamma }_{\mathrm{sc},1b}& T_{\mathrm{ant},1b}\\ I_{\mathrm{sc},1b}& Y_{\mathrm{ant},1b}\end{array}\right]\left[\begin{array}{c}{A}_{1b}\\ V_{1b}\end{array}\right].& \left(52\right)\end{array}$

[0124] The definitions are extended to N upper and lower ports. Define N×N diagonal matrices with {double overscore (Γ)}_sc^u, {double overscore (Γ)}_sc^b, {double overscore (T)}_ant^u, {double overscore (T)}_ant^b, {double overscore (I)}_sc^u, {double overscore (I)}_sc^b, {double overscore (Y)}_ant^u, and N×1 column matrices Ā^u, Ā^b, {overscore (B)}^u, {overscore (B)}^b.

[0125] Solving the combined problem of the exterior and interior problem gives a 2N×2N scattering matrix of {double overscore (S)}:

[0126] $\begin{array}{cc}\left[\begin{array}{c}{\stackrel{\_}{B}}^u\\ {\stackrel{\_}{B}}^b\end{array}\right]=\stackrel{\_}{\stackrel{\_}{S}}\left[\begin{array}{c}{\stackrel{\_}{A}}^u\\ {\stackrel{\_}{A}}^b\end{array}\right]& \left(53\right)\\ \stackrel{\_}{\stackrel{\_}{S}}=\left[\begin{array}{cc}{\stackrel{\_}{\stackrel{\_}{S}}}_{\mathrm{uu}}& {\stackrel{\_}{\stackrel{\_}{S}}}_{\mathrm{ub}}\\ {\stackrel{\_}{\stackrel{\_}{S}}}_{\mathrm{bu}}& {\stackrel{\_}{\stackrel{\_}{S}}}_{\mathrm{bb}}\end{array}\right].& \left(54\right)\end{array}$
Simplified Exterior Problem with a Given Characteristic Admittance

[0127] Consider the simplified case of identical traces for the exterior problem with a given characteristic admittance Y₀ of the traces:

[0128] $\begin{array}{cc}\left[\begin{array}{cc}{\stackrel{\_}{\stackrel{\_}{\Gamma }}}_{\mathrm{sc}}^u& {\stackrel{\_}{\stackrel{\_}{\Gamma }}}_{\mathrm{ant}}^u\\ {\stackrel{\_}{\stackrel{\_}{I}}}_{\mathrm{sc}}^u& {\stackrel{\_}{\stackrel{\_}{Y}}}_{\mathrm{ant}}^u\end{array}\right]=\left[\begin{array}{cc}{\stackrel{\_}{\stackrel{\_}{\Gamma }}}_{\mathrm{sc}}^b& {\stackrel{\_}{\stackrel{\_}{T}}}_{\mathrm{ant}}^b\\ {\stackrel{\_}{\stackrel{\_}{I}}}_{\mathrm{sc}}^b& {\stackrel{\_}{\stackrel{\_}{Y}}}_{\mathrm{ant}}^b\end{array}\right]=\left[\begin{array}{cc}-{\stackrel{\_}{\stackrel{\_}{1}}}_N& {\stackrel{\_}{\stackrel{\_}{1}}}_N\\ -2Y_0{\stackrel{\_}{\stackrel{\_}{1}}}_N& Y_0{\stackrel{\_}{\stackrel{\_}{1}}}_N\end{array}\right].& \left(55\right)\end{array}$
Then the scattering matrix simplifies to:

[0129] $\begin{array}{cc}\stackrel{\_}{\stackrel{\_}{S}}={\left\{Y_0{\stackrel{\_}{\stackrel{\_}{1}}}_{2N}+\stackrel{\_}{\stackrel{\_}{Y}}\right\}}^{-1}\left\{Y_0{\stackrel{\_}{\stackrel{\_}{1}}}_{2N}-\stackrel{\_}{\stackrel{\_}{Y}}\right\}& \left(56\right)\end{array}$
where {double overscore (Y)} is the admittance matrix of the interior problem as given by Equation 47, and {double overscore (1)}_2N is the unit matrix of dimension 2N.

[0130] For the case of two vias, the S-matrix is:

[0131] $\begin{array}{cc}\stackrel{\_}{\stackrel{\_}{S}}=\left[\begin{array}{cccc}{S}_{1u1u}& S_{1u2u}& S_{1u1b}& S_{1u2b}\\ S_{2u1u}& S_{2u2u}& S_{2u1b}& S_{2u2b}\\ S_{1b1u}& S_{1b2u}& S_{1b1b}& S_{1b2u}\\ S_{2b1u}& S_{2b2u}& S_{2b1b}& S_{2b2b}\end{array}\right]& \left(57\right)\end{array}$

[0132] For the case of 256 vias, the S-matrix is of dimension 512×512, with two ports for each via.

Small Via Radius and Small Layer Thickness Approximations

[0133] For the case of a small via radius, keep only T₀^TM and discard all other T_n^TM. Then only the m=0 harmonics of the vias are coupled to each other, yielding the Foldy-Lax equation:

[0134] $\begin{array}{cc}{w}_{\mathrm{l0}}^{\mathrm{TM}\left(q\right)}=a_{\mathrm{l0}}^{\mathrm{TM}\left(q\right)}+\sum _{\begin{array}{c}p=1\\ p\ne q\end{array}}^N{H}_0^{\left(2\right)}\left(k_{\rho l}{\stackrel{\_}{\rho }}_p-{\stackrel{\_}{\rho }}_q\right)T_0^{\left(N\right)}{w}_{l0}^{\mathrm{TM}\left(p\right)}.& \left(58\right)\end{array}$
The dimension of the Foldy-Lax equation then becomes only N.