Afef Bouchaala

Lionel Courau

Olivier Bonnaud

Philippe Galy

W-element RLGC matrices calculation for power distribution planes modeling using MCTL matrix method

Keywords: multilayered PDN, preditcive model, MCTL, W-Element, RLGC Matrices

I. Introduction

Modern electronic systems require a large number of integrated circuits (ICs), more Input/Output (I/O) connections, faster operating clock frequencies associated to low cost and high performance integration, with increased functionality while preserving reliability. In addition, the new generations of electronic products are involving more mixed signal because of the integration of digital, RF, optical and micro-electro-mechanical functions on a single chip or module, which pose tremendous challenges for designers. As a result, maintaining the Signal Integrity (SI) and Power Integrity (PI) for future systems is becoming one of the most important issues. In fact, the number of failures caused by SI/PI problems is on rise because existing methodologies and tools cannot address these issues successfully and at early design phases.

Advanced electronic packages and Printed Circuit Boards (PCB) have been developed in the recent decades to match these new challenges. They are essentially formed by Power Distribution Planes which play a very important role by serving as conduit for the transportation of current, providing charge to the switching circuits at high frequencies and support return currents for the signal lines referenced to them. Ideally, these planes should exhibit low impedance over a large frequency range of operation so that the transient currents induced by simultaneous switching of digital circuits do not lead to excessive noise propagation over the PDN [START_REF] Roy | Delay-Extraction-Based Waveform Relaxation Algorithm for Fast Transient Analysis of Power Distribution Networks[END_REF]. However with the increase in clock speed, the scaling of supply voltage and high switching speed of logic circuits, effects like ground bounce, EM interferences and Simultaneous Switching Noise (SSN) are arising in the PDNs can quickly lead to undesirable voltage fluctuations and propagation delays in chip, board and packaging levels. Hence, Power Distribution Planes are considered as critical area for Electromagnetic Compatibility (EMC) and for Power Integrity verification of high speed packages, and they need more and more investigations in order to analyze and predict unwanted noises [START_REF] Swaminathan | Power integrity Modeling and Design for Semiconductor and Systems[END_REF]. Using proper predictive models, these problems can be identified and eliminated and thus costly tests measurements and redesign phases can be avoided.

Typical PDNs are composed of metal planes stacked on top of each other separated by low-loss dielectrics. Since each layer formed by metal planes with the low-loss insulator can act as a cavity, the PDNs are highly resonant structures [START_REF] Lee | Analysis and suppression of SSN noise coupling between Power/Ground plane cavities through cutouts in multilayer packages and PCBs[END_REF]. To completely characterize such structures through time-domain analysis, a tremendous amount of time is required for a simulation. Hence, the frequency-domain analysis of package PDNs is more beneficial.

PDNs' modeling methods in the frequency domain can be roughly classified into three categories: i) numerical full-wave approaches; ii) analytical approaches; iii) hybrid methods. PDNs have traditionally been modeled using the numerical full-wave techniques based on resolving Maxwell's equations, such as the Finite-Element Method (FEM) [START_REF] Polucarpou | The Finite-Element method for moedling circuits and interconnects for electronic packaging[END_REF], and Finite Difference Time-Domain method (FDTD) [START_REF] Choi | Modeling of power sup-10mm sion ply noise in large chips using the circuit-based finite difference time-domain method[END_REF] to name a few. These full-wave tools are certainly able to handle all of the structures discussed previously with high accurate results, but they come at a major computational cost, and this can be prohibitive for efficient characterization of complex structures, especially when there are small features and discontinuities in the package. In addition, a complete layout is required in order to launch these tools. For these reasons, designers relegate their use to final verification, stage at which design iterations are expensive. Analytical methodologies have long been used for this purpose as well. The most popular ones are based on discretizing the planar structure into a mesh of squared unit cells and each unit cell can be represented either by lumped elements [START_REF] Kim | Modeling of irregular shaped power distribution planes using transmission matrix method[END_REF] or by four simple transmission lines [START_REF] Lee | Modeling and Analysis of Multichip module power supply planes[END_REF]. These methods do not require tremendous time and memory; however, they are limited to single pairs of Power/Ground planes. However, to avoid the computational cost of full-wave methods and the geometrical limitations of the analytical methods, a more efficient approach consists to combine both analytical and numerical techniques in one hybrid technique like in the cavity resonator method [START_REF] Kim | Modeling and measurement of inter level electromagnetic coupling and fringing effect in a hierarchical power distribution network using segmentation method with resonant cavity model[END_REF]. In these cases, the layout is decomposed into traces, vias, planes, and circuits. These elements are sent off to specifically tuned solvers optimized for these structures, and their results are integrated back together into comprehensive S-parameters. These techniques provide nearly full-wave accuracy, while at the same time they result in very large-scale problems to be handled in a reasonable amount of time. Moreover, a PDN layout is still needed, so it is difficult to get predictive electrical model.

A typical multilayered PDN contains stacked Power/Ground planes as shown in Fig. 1. Inspection of this structure reveals two structural characteristics. Firstly, the PDN is composed of multiple metal layers; therefore, it can be modeled by Multi-Conductor Transmission Lines (MCTLs). Secondly, it has a periodic geometry that can be analyzed by a meshing of MCTLs [START_REF] Paul | Analysis of Multiconductor Transmission Lines[END_REF]. MCTL can be modeled with conventional RLGC (Resistance, Inductance, conductance, and Capacitance) lumped elements [START_REF] Ruelhi | Progress in the methodologies for the electrical modeling of interconnects and electronic packages[END_REF]. However, the accuracy while using these elements is limited to scenarios where the rise time of the excitation is much higher than the propagation time over the planes. Moreover, by involving a SPICE tool, the modeling of high frequency effects like skin effect requires additional lumped elements, which can quickly lead to a prohibitive memory size and to excessive runtime costs. Alternative techniques to model MCTLs have also been proposed [START_REF] Yu | Passive multipoint moment matching model order reduction algorithm on multi port distributed interconnect networks[END_REF]. However, all the above models may still lead to large circuit matrices due to the introduction of internal nodes by representing each transmission line segment in the Modified Nodal Analysis (MNA) formulation used by SPICE. In this work we have used the HSPICE circuit simulator with which it is possible to model a multi-conductor lossy frequency dependent transmission line through the so-called W-Element model [START_REF]HSPICE User Guide: Signal Integrity, Version E-2010.12[END_REF].

The main contribution of this paper consists in providing a constructive model for efficient predictive analysis of multilayered PDN while having only some technological and geometrical information about the package. The proposed model is based on W-Element representation of calibrated unit MCTL which can be used for modeling different PDN sizes. Novel equations of RLGC Matrices are given and effective self and mutual inductance are proposed to efficiently describe the inductive interactions among coupled signal lines. This set of novel RLGC equations strengthens accuracy up to 50 GHz compared to full-wave solvers. This paper is organized as follows: section II describes in details our modeling methodology that will be called "MCTL matrix". In section III, numerical examples will be given to illustrate the validation and the efficiency of the proposed method in comparison with a full wave solver based on method of Moments, emphasizing the gain of CPU time and memory obtained by the MCTL Matrix method. Section IV concludes the paper with some prospects for future work.

II. MCTL Predictive Modeling Method

Advanced electronic packages are consisted of multiple layers in order to reduce the parasites of the PDN (e.g., to reduce the inductance of the planes), these layers can be allocated to power and ground in an alternating manner so that multiple plane pairs can exist in a package or board. Fig. 1 shows a typical Power/ Ground planes structure. This structure can be discretized into multilayered squared unit cells. For a good accuracy, the unit cell size must be about twenty times less than the wavelength of the highest frequency of interest, otherwise refection errors may occur. To avoid confusion, there should be no conflict regarding a common global reference terminal for the definition of voltages at the interconnection of the unit cells. The MCTL matrix model overcomes this practical problem by defining multilayered unit cell models that have the same ground reference such that interconnection of the unit cells becomes straightforward. Fig. 2 shows an equivalent circuit model for a sample unit cell including four planes, where the bottom plane is set as the common reference terminal.

A big misconception in multiple plane pair structures is considering that the planes are isolated from each other. For a single plane pair, the skin effect is assumed to be dominant and the solution is extended to multilayered structures by assuming zero coupling between the plane layers through the conductor. However, espe-their use to final verification, stage at which design iterations are expensive.

Analytical methodologies have long been used for this purpose as well. The most popular ones are based on discretizing the planar structure into a mesh of squared unit cells and each unit cell can be represented either by lumped elements [START_REF] Kim | Modeling of irregular shaped power distribution planes using transmission matrix method[END_REF] or by four simple transmission lines [START_REF] Lee | Modeling and Analysis of Multichip module power supply planes[END_REF]. These methods do not require tremendous time and memory; however they are limited to single pairs of Power/Ground planes.

However, to avoid the computational cost of full-wave methods and the geometrical limitations of the analytical methods, a more efficient approach consists to combine both analytical and numerical techniques in one hybrid technique like in the cavity resonator method [START_REF] Kim | Modeling and measurement of inter level electromagnetic coupling and fringing effect in a hierarchical power distribution network using segmentation method with resonant cavity model[END_REF]. In these cases, the layout is decomposed into traces, vias, planes, and circuits. These elements are sent off to specifically tuned solvers optimized for these structures, and their results are integrated back together into comprehensive S-parameters. These techniques provide nearly full-wave accuracy, while at the same time they result in very large-scale problems to be handled in a reasonable amount of time. Moreover a PDN layout is still needed, so it is difficult to get predictive electrical model.

A typical multilayered PDN contains stacked Power/Ground planes as shown in Fig. 1. Inspection of this structure reveals two structural characteristics. Firstly, the PDN is composed of multiple metal layers; therefore, it can be modeled by Multi-Conductor Transmission Lines (MCTLs). Secondly, it has a periodic geometry that can be analyzed by a meshing of MCTLs [START_REF] Paul | Analysis of Multiconductor Transmission Lines[END_REF]. MCTL can be modeled with conventional RLGC (Resistance, Inductance, conductance, and Capacitance) lumped elements [START_REF] Ruelhi | Progress in the methodologies for the electrical modeling of interconnects and electronic packages[END_REF]. However, the accuracy while using these elements is limited to scenarios where the rise time of the excitation is much higher than the propagation time over the planes. Moreover, by involving a SPICE tool, the modeling of high frequency effects like skin effect requires additional lumped elements, which can quickly lead to a prohibitive memory size and to excessive runtime costs. Alternative techniques to model MCTLs have also been proposed [START_REF] Yu | Passive multipoint moment matching model order reduction algorithm on multi port distributed interconnect networks[END_REF]. However, all the above models may still lead to large circuit matrices due to the introduction of internal nodes by representing each transmission line segment in the Modified Nodal Analysis (MNA) formulation used by SPICE. In this work we have used the HSPICE circuit simulator with which it is possible to model a multi-conductor lossy frequency dependent transmission line through the so-called W-Element model [START_REF]HSPICE User Guide: Signal Integrity, Version E-2010.12[END_REF].

The main contribution of this paper consists in providing a constructive model for efficient predictive analysis of multilayered PDN while having only some technological of novel RLGC equations strengthens accuracy up to 50 GHz compared to full-wave solvers. This paper is organized as follows: section II describes in details our modeling methodology that will be called "MCTL matrix". In section III, numerical examples will be given to illustrate the validation and the efficiency of the proposed method in comparison with a full wave solver based on method of Moments, emphasizing the gain of CPU time and memory obtained by the MCTL Matrix method. Section IV concludes the paper with some prospects for future work.

II. MCTL PREDICTIVE MODELING METHOD

Advanced electronic packages are consisted of multiple layers in order to reduce the parasites of the PDN (e.g., to reduce the inductance of the planes), these layers can be allocated to power and ground in an alternating manner so that multiple plane pairs can exist in a package or board. Fig. 1 shows a typical Power/Ground planes structure. This structure can be discretized into multilayered squared unit cells. For a good accuracy, the unit cell size must be about twenty times less than the wavelength of the highest frequency of interest, otherwise refection errors may occur. To avoid confusion, there should be no conflict regarding a common global reference terminal for the definition of voltages at the interconnection of the unit cells. The MCTL matrix model overcomes this practical problem by defining multilayered unit cell models that have the same ground reference such that interconnection of the unit cells becomes straightforward. Fig. 2 shows an equivalent circuit model for a sample unit cell including four planes, where the bottom plane is set as the common reference terminal. cially when the planes resonate, substantial coupling between the plane layers can occur through the magnetic fields penetrating the solid conductor. As a result, each plane assigns the plane below it as its local reference plane. In other words, considering a current on the bottom of the i th plane, the return path for the same current is considered to exist on the top of the (i+1) th plane below, regardless its nature. This can be proven in Fig. 3 which represents a 3D visualization of forward current on the top layer and the return current on the plane just below. This structure was simulated using FEM EMPro 3D simulator. It is clear that when exciting the top plane, we get nearly the same intensity for the return currents on the plane just below, despite the ground reference was assigned to the bottom plane.

MCTL can be represented by W-Element model in a SPICE simulator like HSPICE. W-Element has numerous advantages compared to the HSPICE U-Element. First, it does not create any spurious ringing as that is produced by the U-Element in the time domain. Moreover, at higher frequencies, some phenomena become increasingly predominant such as skin effect losses, dielectric losses, dispersion and radiation losses, especially when lossy dielectric materials like FR4 is used. The W-Elements are best to accurately take into account these frequency dependent parameters [START_REF] Yu | Passive multipoint moment matching model order reduction algorithm on multi port distributed interconnect networks[END_REF].

MCTL Matrix Modeling Flow

At this stage, it is important to note that our goal from this study is to provide a predictive model for Power/Ground planes without the need of a representative layout. In other words, only some key information about PDN geometry (planes 'length, width, and thickness, port locations, etc.) and stack-up technology (dielectric constant, loss tangent, etc.) are needed to construct a MCTL-based grid as described. As a consequence, we can get our package or PCB model as early as the specification phase, thing that other models cannot provide. The modeling flow is described in Fig. 4. Being based on the chip design constraints, we start by fixing the geometrical parameters (length, width, thickness, port locations, number of layers, shape of Power/Ground planes, etc.) and technological ones (dielectric constant, loss tangent, etc.). Given the PDN size and with the maximum frequency of interest, we can conclude the size of the unit cell, L u , and then we can deduce the total number of unit cells required to construct the PDN predictive model. The technological parameters of the PDN (conductor and dielectric properties) will determine the substrate to which the unit MCTL will be associated. Once these parameters are fixed, we define a MCTL line with L u as length, associated to the substrate. A W-Element Extraction can be then performed to provide RLGC lossy matrices, but thanks to our approach we can obtain these matrices with simple calculations, which will be developed in the next section.

In order to reduce the time for building the equivalent circuit for the entire PDN's Planes, which is time consuming for organizing millions of nodes and elements manually, a software program in Perl was developed to generate the Power/Ground grid model automatically. The program takes the geometry of the layout and the properties of the material as inputs. It is easily customizable in terms of the meshing size and the number of coupled layers of the PDN to meet the application requirements.

RLGC Matrices Calculation for the W-Element

The aim of this section is to build accurate per-unit-length (p.u.l) RLGC matrices of the W-Element, using frequencydependent equations. As mentioned previously, a W-Element is a lossy multi-conductor frequency-dependent transmission line, based on a novel state-of-the art simulation method; they are simulated very fast by SPICE simulators, and are accurate and robust enough especially when we need to model high frequency dependent loss. The W-Element is organized in terms of coupled conductors, and there is no limit on the number of coupled conductors. Each of these conductors has two terminals (one at each end of the conductor). A reference conductor is always assured, the number of conductors is therefore related to the number of W-Element nodes as n=2*(c+1), where c is the number of conductors (including the reference) and n is the number of nodes. Once the length of the unit cell is set, we can approximate the width of each line to half of the line length. Adding the permittivity and loss tangent of the dielectric, the p.u.l RLGC elements can be easily computed.

single plane pair, the skin effect is assumed to be dominant and the solution is extended to multilayered structures by assuming zero coupling between the plane layers through the conductor. However, especially when the planes resonate, substantial coupling between the plane layers can occur through the magnetic fields penetrating the solid conductor. As a result, each plane assigns the plane below it as its local reference plane. In other words, considering a current on the bottom of the i th plane, the return path for the same current is considered to exist on the top of the (i+1) th plane below, regardless its nature. This can be proven in Fig. 3 which represents a 3D visualization of forward current on the top layer and the return current on the plane just below. This structure was simulated using FEM EMPro 3D simulator. It is clear that when exciting the top plane, we get nearly the same intensity for the return currents on the plane just below, despite the ground reference was assigned to the bottom plane. MCTL can be represented by W-Element model in a SPICE simulator like HSPICE. W-Element has numerous advantages compared to the HSPICE U-Element. First, it does not create any spurious ringing as that is produced by the U-Element in the time domain. Moreover, at higher frequencies, some phenomena become increasingly predominant such as skin effect losses, dielectric losses, dispersion and radiation layers, shape of Power/Ground planes, etc.) and technological ones (dielectric constant, loss tangent, etc.). Given the PDN size and with the maximum frequency of interest, we can conclude the size of the unit cell, L u , and then we can deduce the total number of unit cells required to construct the PDN predictive model. The technological parameters of the PDN (conductor and dielectric properties) will determine the substrate to which the unit MCTL will be associated. Once these parameters are fixed, we define a MCTL line with L u as length, associated to the substrate. A W-Element Extraction can be then performed to provide RLGC lossy matrices, but thanks to our approach we can obtain these matrices with simple calculations, which will be developed in the next section. In order to reduce the time for building the equivalent circuit for the entire PDN's Planes, which is time consuming for organizing millions of nodes and elements manually, a software program in Perl was developed to generate the Power/Ground grid model automatically. The program takes conductor. However, especially when the planes resonate, substantial coupling between the plane layers can occur through the magnetic fields penetrating the solid conductor. As a result, each plane assigns the plane below it as its local reference plane. In other words, considering a current on the bottom of the i th plane, the return path for the same current is considered to exist on the top of the (i+1) th plane below, regardless its nature. This can be proven in Fig. 3 which represents a 3D visualization of forward current on the top layer and the return current on the plane just below. This structure was simulated using FEM EMPro 3D simulator. It is clear that when exciting the top plane, we get nearly the same intensity for the return currents on the plane just below, despite the ground reference was assigned to the bottom plane. MCTL can be represented by W-Element model in a SPICE simulator like HSPICE. W-Element has numerous advantages compared to the HSPICE U-Element. First, it does not create any spurious ringing as that is produced by the U-Element in the time domain. Moreover, at higher frequencies, some phenomena become increasingly predominant such as skin effect losses, dielectric losses, dispersion and radiation conclude the size of the unit cell, L u , and then we can deduce the total number of unit cells required to construct the PDN predictive model. The technological parameters of the PDN (conductor and dielectric properties) will determine the substrate to which the unit MCTL will be associated. Once these parameters are fixed, we define a MCTL line with L u as length, associated to the substrate. A W-Element Extraction can be then performed to provide RLGC lossy matrices, but thanks to our approach we can obtain these matrices with simple calculations, which will be developed in the next section. In order to reduce the time for building the equivalent circuit for the entire PDN's Planes, which is time consuming for organizing millions of nodes and elements manually, a software program in Perl was developed to generate the Power/Ground grid model automatically. The program takes

a. DC and AC Resistance Matrices

The R matrices are composed of R 0 and R s parameters. The R 0 parameter is the DC resistance of the transmission line. In the R 0 matrix, the diagonal elements represent the resistance of each conductor and the off-diagonal elements are null.

R 0 can be calculated using the equation ( 1)

Where l u is the length of the MCTL, σ is the conductivity of the conductor, and e is the thickness of the conductor which is in reality the thickness of the Power/Ground planes.

The R s parameter is the skin effect resistance of the transmission line. In the R s matrix the diagonal elements represent the AC resis- tance of each conductor in isolation from the other. The off diagonal elements are null. The AC series resistance is scaled with the square root of the frequency.

R s can be computed using the equation ( 2)

Note that the R 0 and R s resistances equations in the circuit formu- lation used in the method of the reference [START_REF] Lee | Modeling and Analysis of Multichip module power supply planes[END_REF] to model a simple planes pair, contains a factor two to account of the resistance of the both planes. This factor is omitted since we calculate the resistance of each layer separately.

b. Capacitance Matrix

The C parameter represents the self and the coupling capacitance between the conductors. Basically, each diagonal term in Max-wellian form is the sum of the entire row, so for example the C 11maxwellian term is actually equal to C 10 +C 12 , where C 10 is the capacitance of the conductor 1 to the global ground. The off-diagonal terms are equal to the minus of the mutual elements, so for example C 21 -maxwellian is equal to -C 21 .

C ij can be computed as:

Where ε r is the relative permittivity of the dielectric material, ε0 represents the permittivity of the free-space, and d is the separation between conductors i and j.

c. Conductance Matrix

The G parameter is the dielectric-loss conductance per-unitlength. The G matrix has the same general form as the C Matrix G ij can be expressed as:

Where tanδ represents the loss tangent of the dielectric material and C ij the capacitance p.u.l.

d. Partial and Mutual Inductance

The most confusing, subtle, and important parameter in highspeed packaging and interconnect design is the inductance. It plays a key role in the origin of SSN, and in crosstalk between transmission line structures. This part is devoted for a new comprehensive theory of p.u.l inductance in Inductance is calculated by using the partial concept which was firstly introduced by Dr. Clayton Paul and Dr. Ruehli [14] [15].

Traditionally, inductance elements in the matrix can be calculated as the capacitance matrix. While the capacitive coupling between non-adjacent conductors can often be ignored, mutual inductive coupling is a long range issue and cannot be ignored in non-adjacent conductors. So, for the L matrix calculation the diagonal elements are actually the self-partial inductance minus the mutual inductance from other conductors. Coupling with the rest of the lines in the mesh can be ignored since the length line L u is much larger than the separation between conductors of the same MCTL.

The off-diagonal parameters on each column are the mutual inductance between the conductor and the remaining conductors of the MCTL and which are shown in Fig. 6.

approximate the width of each line to half of the line length.

Adding the permittivity and loss tangent of the dielectric, the p.u.l RLGC elements can be easily computed. The R matrices are composed of R 0 and R s parameters. The R 0 parameter is the DC resistance of the transmission line. In the R 0 matrix, the diagonal elements represent the resistance of each conductor and the off-diagonal elements are null.

𝑅𝑅 0 = [ 𝑅𝑅 011 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 𝑅𝑅 0𝑖𝑖𝑖𝑖 ]
R 0 can be calculated using the equation ( 1)

R 0ii = ( 𝑙𝑙 𝑢𝑢 𝜎𝜎.𝑙𝑙 𝑢𝑢 .𝑒𝑒 ) /𝑙𝑙 𝑢𝑢 (1)
simple planes pair, contains a factor two to account of the resistance of the both planes. This factor is omitted since we calculate the resistance of each layer separately.

b. Capacitance Matrix

The C parameter represents the self and the coupling capacitance between the conductors. Basically, each diagonal term in Maxwellian form is the sum of the entire row, so for example the C 11 -maxwellian term is actually equal to C 10 +C 12 , where C 10 is the capacitance of the conductor 1 to the global ground. The off-diagonal terms are equal to the minus of the mutual elements, so for example C 21 -maxwellian is equal to -C 21 . C ij can be computed as:

𝑐𝑐 𝑖𝑖𝑖𝑖 = ( 𝜀𝜀 𝑟𝑟 .𝜀𝜀 0 .𝑙𝑙 𝑢𝑢 .𝑙𝑙 𝑢𝑢 𝑑𝑑 ) 𝑙𝑙 𝑢𝑢 ⁄ (3)
Where ε r is the relative permittivity of the dielectric material, 𝜀𝜀 0 represents the permittivity of the free-space, and d is the separation between conductors i and j. p.u.l RLGC elements can be easily computed. The R matrices are composed of R 0 and R s parameters. The R 0 parameter is the DC resistance of the transmission line. In the R 0 matrix, the diagonal elements represent the resistance of each conductor and the off-diagonal elements are null.

𝑅𝑅 0 = [ 𝑅𝑅 011 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 𝑅𝑅 0𝑖𝑖𝑖𝑖 ]
R 0 can be calculated using the equation ( 1)

R 0ii = ( 𝑙𝑙 𝑢𝑢 𝜎𝜎.𝑙𝑙 𝑢𝑢 .𝑒𝑒 ) /𝑙𝑙 𝑢𝑢 (1)
calculate the resistance of each layer separately.

b. Capacitance Matrix

The C parameter represents the self and the coupling capacitance between the conductors. Basically, each diagonal term in Maxwellian form is the sum of the entire row, so for example the C 11 -maxwellian term is actually equal to C 10 +C 12 , where C 10 is the capacitance of the conductor 1 to the global ground. The off-diagonal terms are equal to the minus of the mutual elements, so for example C 21 -maxwellian is equal to -C 21 . C ij can be computed as:

𝑐𝑐 𝑖𝑖𝑖𝑖 = ( 𝜀𝜀 𝑟𝑟 .𝜀𝜀 0 .𝑙𝑙 𝑢𝑢 .𝑙𝑙 𝑢𝑢 𝑑𝑑 ) 𝑙𝑙 𝑢𝑢 ⁄ (3)
Where ε r is the relative permittivity of the dielectric material, 𝜀𝜀 0 represents the permittivity of the free-space, and d is the separation between conductors i and j.

p.u.l RLGC elements can be easily computed. The R matrices are composed of R 0 and R s parameters. The R 0 parameter is the DC resistance of the transmission line. In the R 0 matrix, the diagonal elements represent the resistance of each conductor and the off-diagonal elements are null.

𝑅𝑅 0 = [ 𝑅𝑅 011 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 𝑅𝑅 0𝑖𝑖𝑖𝑖 ]
R 0 can be calculated using the equation ( 1)

R 0ii = ( 𝑙𝑙 𝑢𝑢 𝜎𝜎.𝑙𝑙 𝑢𝑢 .𝑒𝑒 ) /𝑙𝑙 𝑢𝑢 (1)
calculate the resistance of each layer separately.

b. Capacitance Matrix

The C parameter represents the self and the coupling capacitance between the conductors. Basically, each diagonal term in Maxwellian form is the sum of the entire row, so for example the C 11 -maxwellian term is actually equal to C 10 +C 12 , where C 10 is the capacitance of the conductor 1 to the global ground. The off-diagonal terms are equal to the minus of the mutual elements, so for example C 21 -maxwellian is equal to -C 21 . C ij can be computed as:

𝑐𝑐 𝑖𝑖𝑖𝑖 = ( 𝜀𝜀 𝑟𝑟 .𝜀𝜀 0 .𝑙𝑙 𝑢𝑢 .𝑙𝑙 𝑢𝑢 𝑑𝑑 ) 𝑙𝑙 𝑢𝑢 ⁄ (3)
Where ε r is the relative permittivity of the dielectric material, 𝜀𝜀 0 represents the permittivity of the free-space, and d is the separation between conductors i and j.

Where l u is the length of the MCTL, σ is the conductivity of the conductor, and e is the thickness of the conductor which is in reality the thickness of the Power/Ground planes. The R s parameter is the skin effect resistance of the transmission line. In the R s matrix the diagonal elements represent the AC resistance of each conductor in isolation from the other. The off diagonal elements are null. The AC series resistance is scaled with the square root of the frequency.

𝑅𝑅 𝑠𝑠 = [ 𝑅𝑅 𝑠𝑠11 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 𝑅𝑅 𝑠𝑠𝑖𝑖𝑖𝑖 ]
R s can be computed using the equation ( 2)

R si = (√ 𝜋𝜋.𝑓𝑓.𝜇𝜇 0 0.6321.𝜎𝜎.𝑙𝑙 𝑢𝑢 2 ) 𝑙𝑙 𝑢𝑢 ⁄ (2)
Note that the R 0 and R s resistances equations in the circuit formulation used in the method of the reference [START_REF] Lee | Modeling and Analysis of Multichip module power supply planes[END_REF] to model a simple planes pair, contains a factor two to account of the resistance of the both planes. This factor is omitted since we calculate the resistance of each layer separately.

b. Capacitance Matrix

The C parameter represents the self and the coupling capacitance between the conductors. Basically, each diagonal term in Maxwellian form is the sum of the entire row, so for example the C 11 -maxwellian term is actually equal to C 10 +C [START_REF]HSPICE User Guide: Signal Integrity, Version E-2010.12[END_REF] , where C 10 is the capacitance of the conductor 1 to the global ground. The off-diagonal terms are equal to the minus of the Where l u is the length of the MCTL, σ is the conductivity of the conductor, and e is the thickness of the conductor which is in reality the thickness of the Power/Ground planes. The R s parameter is the skin effect resistance of the transmission line. In the R s matrix the diagonal elements represent the AC resistance of each conductor in isolation from the other. The off diagonal elements are null. The AC series resistance is scaled with the square root of the frequency.

𝑅𝑅 𝑠𝑠 = [ 𝑅𝑅 𝑠𝑠11 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 𝑅𝑅 𝑠𝑠𝑖𝑖𝑖𝑖 ]
R s can be computed using the equation ( 2)

R si = (√ 𝜋𝜋.𝑓𝑓.𝜇𝜇 0 0.6321.𝜎𝜎.𝑙𝑙 𝑢𝑢 2 ) 𝑙𝑙 𝑢𝑢 ⁄ (2)
Note that the R 0 and R s resistances equations in the circuit formulation used in the method of the reference [START_REF] Lee | Modeling and Analysis of Multichip module power supply planes[END_REF] to model a simple planes pair, contains a factor two to account of the resistance of the both planes. This factor is omitted since we calculate the resistance of each layer separately.

b. Capacitance Matrix

The C parameter represents the self and the coupling capacitance between the conductors. Basically, each diagonal term in Maxwellian form is the sum of the entire row, so for example the C 11 -maxwellian term is actually equal to C 10 +C [START_REF]HSPICE User Guide: Signal Integrity, Version E-2010.12[END_REF] , where C 10 is the capacitance of the conductor 1 to the global The R matrices are composed of R 0 and R s parameters. The R 0 parameter is the DC resistance of the transmission line. In the R 0 matrix, the diagonal elements represent the resistance of each conductor and the off-diagonal elements are null.

𝑅𝑅 0 = [ 𝑅𝑅 011 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 𝑅𝑅 0𝑖𝑖𝑖𝑖 ]
R 0 can be calculated using the equation ( 1)

R 0ii = ( 𝑙𝑙 𝑢𝑢 𝜎𝜎.𝑙𝑙 𝑢𝑢 .𝑒𝑒 ) /𝑙𝑙 𝑢𝑢 (1)
The C parameter represents the self and the coupling capacitance between the conductors. Basically, each diagonal term in Maxwellian form is the sum of the entire row, so for example the C 11 -maxwellian term is actually equal to C 10 +C [START_REF]HSPICE User Guide: Signal Integrity, Version E-2010.12[END_REF] , where C 10 is the capacitance of the conductor 1 to the global ground. The off-diagonal terms are equal to the minus of the mutual elements, so for example C 21 -maxwellian is equal to -C 21 .

Where ε r is the relative permittivity of the dielectric material, 𝜀𝜀 0 represents the permittivity of the free-space, and d is the separation between conductors i and j. The R matrices are composed of R 0 and R s parameters. The R 0 parameter is the DC resistance of the transmission line. In the R 0 matrix, the diagonal elements represent the resistance of each conductor and the off-diagonal elements are null.

𝑅𝑅 0 = [ 𝑅𝑅 011 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 𝑅𝑅 0𝑖𝑖𝑖𝑖 ]
R 0 can be calculated using the equation ( 1)

R 0ii = ( 𝑙𝑙 𝑢𝑢 𝜎𝜎.𝑙𝑙 𝑢𝑢 .𝑒𝑒 ) /𝑙𝑙 𝑢𝑢 (1)
The C parameter represents the self and the coupling capacitance between the conductors. Basically, each diagonal term in Maxwellian form is the sum of the entire row, so for example the C 11 -maxwellian term is actually equal to C 10 +C [START_REF]HSPICE User Guide: Signal Integrity, Version E-2010.12[END_REF] , where C 10 is the capacitance of the conductor 1 to the global ground. The off-diagonal terms are equal to the minus of the mutual elements, so for example C 21 -maxwellian is equal to -C 21 . G ij can be expressed as:

𝐺𝐺 𝑖𝑖𝑖𝑖 = 𝜋𝜋. 𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡. 𝑐𝑐 𝑖𝑖𝑖𝑖 (4)
Where tanδ represents the loss tangent of the dielectric material and C ij the capacitance p.u.l.

d. Partial and Mutual Inductance

The most confusing, subtle, and important parameter in high-speed packaging and interconnect design is the inductance. It plays a key role in the origin of SSN, and in crosstalk between transmission line structures. This part is devoted for a new comprehensive theory of p.u.l inductance in RLGC Matrix. Inductance is calculated by using the partial concept which was firstly introduced by Dr. Clayton and Dr. Ruehli [START_REF] Ruehli | Inductance calculations in c complex integrated circuit environment[END_REF] [START_REF] Paul | Modeling electromagnetic interference properties of printed circuit boards[END_REF]. Traditionally, inductance elements in the matrix can be calculated as the capacitance matrix. While the capacitive coupling between non-adjacent conductors can often be ignored, mutual inductive coupling is a long range issue and cannot be ignored in non-adjacent conductors. So, for the L matrix calculation the diagonal elements are actually the self-partial inductance minus the mutual inductance from other conductors. Coupling with the rest of the lines in the mesh can be ignored since the length line L u is much larger than the separation between conductors of the same MCTL. The off-diagonal parameters on each column are the mutual inductance between the conductor and the remaining conductors of the MCTL and which are shown in Fig. 6 G ij can be expressed as:

𝐺𝐺 𝑖𝑖𝑖𝑖 = 𝜋𝜋. 𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡. 𝑐𝑐 𝑖𝑖𝑖𝑖 (4)
Where tanδ represents the loss tangent of the dielectric material and C ij the capacitance p.u.l.

d. Partial and Mutual Inductance

The most confusing, subtle, and important parameter in high-speed packaging and interconnect design is the inductance. It plays a key role in the origin of SSN, and in crosstalk between transmission line structures. This part is devoted for a new comprehensive theory of p.u.l inductance in RLGC Matrix. Inductance is calculated by using the partial concept which was firstly introduced by Dr. Clayton and Dr. Ruehli [START_REF] Ruehli | Inductance calculations in c complex integrated circuit environment[END_REF] [START_REF] Paul | Modeling electromagnetic interference properties of printed circuit boards[END_REF]. Traditionally, inductance elements in the matrix can be calculated as the capacitance matrix. While the capacitive coupling between non-adjacent conductors can often be ignored, mutual inductive coupling is a long range issue and cannot be ignored in non-adjacent conductors. So, for the L matrix calculation the diagonal elements are actually the self-partial inductance minus the mutual inductance from other conductors. Coupling with the rest of the lines in the mesh can be ignored since the length line L u is much larger than the separation between conductors of the same MCTL. The off-diagonal parameters on each column are the mutual inductance between the conductor and the remaining conductors of the MCTL and which are shown in Fig. 6 The general form of the L inductance is presented as follows:

Traditionally, the mutual inductance is calculated as the product µ 0 .d, but this is only valid for important dielectric thickness. Based on the work of Dr. Clayton Paul on parallel conductors of rectangular cross section, and since the dielectric separation is much smaller than the length of the conductors, which is the case for the MCTL, the mutual-partial inductance p.u.l can be given by: Where d is the dielectric separation between two conductors. On the other hand, the self-partial-inductance p.u.l is given by: Where w and t are respectively the width and the thickness of the conductor.

e. The Mesh Schema

The mesh of the model containing all the RLGC elements is presented in the Fig. 7. The interconnection of the different MCTLs or the RLGC elements with each other does not require any special consideration and can be done in a straightforward manner. Note that Fig. 7 contains two different MCTLs elements: Internal elements for all the internal paths of the grid, whereas the external elements are only for the edges. The inductance and resistances of the external elements are adjusted by a factor of two, to take into account the fact that they have no neighboring planes. The elements on the boundary are assigned with a half of the capacitance of an internal cell.

III. Numerical Results

In this section we present scattering parameters results of different shapes and sizes of multilayered Power/Ground planes with different stack-up definition in order to demonstrate the validity and the efficiency of the proposed MCTL based model for PDNs Power/Ground planes. Four test structures were considered, the three first ones are simple two-layered structures and the last one is an extension to a third dimension PDN. Results have been compared to Keysight Momentum EM Solver results. All the simulations were performed on an Intel Xeon workstation with a 3GHz CPU and 3.25 GBytes of RAM.

Two Layered Power/Ground Planes a. Thin Dielectric Layer

We consider the structure shown in Fig. 8. It consists of 3mm by 3mm of Power/Ground plane pair. Metal layers are made of copper (σ = 5.87 x 10 7 S/m) of thickness 20μm, and separated by an FR4 with permittivity ε r = 4.8, Tanδ = 0.002, and thickness d = 30μm. Port1 and Port2 are located on the middle of each side. The mesh density was fixed to 30 MCTLs/row and so L u = 0.1 mm and W u = 0.05mm.

RLGC Matrices calculated at 10GHz are given by:

The modeled results of S 11 and S 12 (magnitude only) between Port1 and Port2 are illustrated in Fig. 9. Good agreement with Momentum full-wave solver is evident up to 40 GHz with some deviation in detail above that frequency. Potential sources of high frequency error include first order causes like the size of the unit cell (length of the MCTL), second order sources like dispersion, dielectric loss, and edge radiation.

inductance between the conductor and the remaining conductors of the MCTL and which are shown in Fig. 6. Traditionally, the mutual inductance is calculated as the product 𝜇𝜇 0 . 𝑑𝑑, but this is only valid for important dielectric Fig. 7: MCTL method mesh schema.

III. NUMERICAL RESULTS

In this section we present scattering parameters results of different shapes and sizes of multilayered Power/Ground planes with different stack-up definition in order to demonstrate the validity and the efficiency of the proposed MCTL based model for PDNs Power/Ground planes. Four test structures were considered, the three first ones are simple two-layered structures and the last one is an extension to a inductance between the conductor and the remaining conductors of the MCTL and which are shown in Fig. 6. Traditionally, the mutual inductance is calculated as the product 𝜇𝜇 0 . 𝑑𝑑, but this is only valid for important dielectric Based on the work of Dr. Clayton on parallel conductors of rectangular cross section, and since the dielectric separation is much smaller than the length of the conductors, which is the case for the MCTL, the mutual-partial inductance p.u.l can be given by:

𝑀𝑀 𝑖𝑖𝑖𝑖 ≈ 𝜇𝜇 2𝜋𝜋 [ln ( 2𝑙𝑙 𝑑𝑑 ) -1] for 𝑑𝑑 𝑙𝑙 ≪ 1 ( 5
)
Where d is the dielectric separation between two conductors. On the other hand, the self-partial-inductance p.u.l is given by:

𝐿𝐿 𝑖𝑖𝑖𝑖 ≈ 𝜇𝜇 2𝜋𝜋 [ln ( 8𝑙𝑙 𝑤𝑤 ) -1] for 𝑤𝑤 𝑙𝑙 ≪ 1 and 𝑤𝑤 𝑖𝑖 ≪ 1 (6)
Where w and t are respectively the width and the thickness of the conductor.

e. The Mesh Schema

The mesh of the model containing all the RLGC elements is presented in the Fig. 7. The interconnection of the different MCTLs or the RLGC elements with each other does not require any special consideration and can be done in a straightforward manner. Note that Fig. 7 contains two different MCTLs elements: Internal elements for all the internal paths of the grid, whereas the external elements are only for the edges. The inductance and resistances of the external elements are adjusted by a factor of two, to take into account the fact that they have no neighboring planes. The elements on the boundary are assigned with a half of the capacitance of an internal cell. Based on the work of Dr. Clayton on parallel conductors of rectangular cross section, and since the dielectric separation is much smaller than the length of the conductors, which is the case for the MCTL, the mutual-partial inductance p.u.l can be given by:

𝑀𝑀 𝑖𝑖𝑖𝑖 ≈ 𝜇𝜇 2𝜋𝜋 [ln ( 2𝑙𝑙 𝑑𝑑 ) -1] for 𝑑𝑑 𝑙𝑙 ≪ 1 ( 5
)
Where d is the dielectric separation between two conductors. On the other hand, the self-partial-inductance p.u.l is given by:

𝐿𝐿 𝑖𝑖𝑖𝑖 ≈ 𝜇𝜇 2𝜋𝜋 [ln ( 8𝑙𝑙 𝑤𝑤 ) -1] for 𝑤𝑤 𝑙𝑙 ≪ 1 and 𝑤𝑤 𝑖𝑖 ≪ 1 (6)
Where w and t are respectively the width and the thickness of the conductor.

e. The Mesh Schema

The mesh of the model containing all the RLGC elements is presented in the Fig. 7. The interconnection of the different MCTLs or the RLGC elements with each other does not require any special consideration and can be done in a straightforward manner. Note that Fig. 7 contains two different MCTLs elements: Internal elements for all the internal paths of the grid, whereas the external elements are only for the edges. The inductance and resistances of the external elements are adjusted by a factor of two, to take into account the fact that they have no neighboring planes. The elements on the boundary are assigned with a half of the capacitance of an internal cell. Based on the work of Dr. Clayton on parallel conductors of rectangular cross section, and since the dielectric separation is much smaller than the length of the conductors, which is the case for the MCTL, the mutual-partial inductance p.u.l can be given by:

𝑀𝑀 𝑖𝑖𝑖𝑖 ≈ 𝜇𝜇 2𝜋𝜋 [ln ( 2𝑙𝑙 𝑑𝑑 ) -1] for 𝑑𝑑 𝑙𝑙 ≪ 1 ( 5
)
Where d is the dielectric separation between two conductors. On the other hand, the self-partial-inductance p.u.l is given by:

𝐿𝐿 𝑖𝑖𝑖𝑖 ≈ 𝜇𝜇 2𝜋𝜋 [ln ( 8𝑙𝑙 𝑤𝑤 ) -1] for 𝑤𝑤 𝑙𝑙 ≪ 1 and 𝑤𝑤 𝑖𝑖 ≪ 1 (6)
Where w and t are respectively the width and the thickness of the conductor.

e. The Mesh Schema

The mesh of the model containing all the RLGC elements is presented in the Fig. 7. The interconnection of the different MCTLs or the RLGC elements with each other does not require any special consideration and can be done in a straightforward manner. Note that Fig. 7 contains two different MCTLs elements: Internal elements for all the internal paths of the grid, whereas the external elements are only for the edges. The inductance and resistances of the external elements are adjusted by a factor of two, to take into account the fact that they have no neighboring planes. The elements on the boundary are assigned with a half of the capacitance of an internal cell. 𝐺𝐺 𝑖𝑖𝑖𝑖𝑖𝑖 (𝑆𝑆/𝑚𝑚) = [ 8.92𝑒𝑒 -13 -8.92𝑒𝑒 -13 -8.92𝑒𝑒 -13 8.92𝑒𝑒 -13 ]; 𝐺𝐺 𝑒𝑒𝑒𝑒𝑖𝑖 (𝑆𝑆/𝑚𝑚) = [ 4.39𝑒𝑒 -13 -4.39𝑒𝑒 -13 -4.39𝑒𝑒 -13 4.39𝑒𝑒 -13 ] The modeled results of S 11 and S 12 (magnitude only) between Port1 and Port2 are illustrated in Fig. 9. Good agreement with Momentum full-wave solver is evident up to 40 GHz with some deviation in detail above that frequency. Potential sources of high frequency error include first order causes like the size of the unit cell (length of the MCTL), second order sources like dispersion, dielectric loss, and edge radiation. The modeled results of S 11 and S 12 (magnitude only) between Port1 and Port2 are illustrated in Fig. 9. Good agreement with Momentum full-wave solver is evident up to 40 GHz with some deviation in detail above that frequency. Potential sources of high frequency error include first order causes like the size of the unit cell (length of the MCTL), second order sources like dispersion, dielectric loss, and edge radiation. As mentioned in section I, MCTL Matrix method enables memory and computation time savings. The proposed model substantially reduces the CPU run time that requires only 10.89 seconds on the Intel Xeon machine to simulate this Power/ Ground plane pair, for which Momentum full-wave solver would have taken 32 min for a mesh density of 50 cells per wavelength. Our approach can reduce also by about minimum six times the memory amount required compared to Momentum. Indeed, it uses only 191,02Mbytes to simulate this Power/ Ground Planes pair, however, Momentum took about 1156 Mbytes. The detailed comparison is summarized in Table I. Note that these results have been improved comparing to our last publication [START_REF] Bouchaala | Efficient method for predictive package Power/Ground planes modeling using Multi-Conductor Transmission Lines[END_REF] thanks to the exact calculations of the RLGC matrix which allowed rapid convergence. This huge difference between the resources usage comes from the long process for creating and simulating structures in a full-wave solver: first of all, engineer must draw the layout or imports it from another design simulator. Substrate characteristics must be defined layer by layer if we do not have a ready one. User must also assign the ports definition. Once the design is ready, Momentum calculates the Green's functions that characterize the substrate and generate circuit mesh for a specified frequency range. To study the efficiency of the method to model different stack-up, we consider the case of dielectric thickness of 100μm and the permittivity was fixed to 5. Conductor thickness is assumed to be 50μm.

RLGC Matrices are given by:

The reflection and transmission parameters are given by Fig. 10.

We can easily see that we keep the same accuracy of the model as for thin dielectric. So our model is totally independent of the technological parameters and can be used for any substrate definition. Another advantage is that the simulation time is invariant since the number of layers remains the same. MCTL Matrix has used the same CPU and memory resources as for the thin structure to simulate the thick dielectric.

c. Invariant Unit Cell For Different Power/Ground Planes Sizes

Another advantage of the MCTL Matrix method is that all Power/ Ground planes dimensions can be simulated by using the same unit cell topology and sizing which makes our method easy customizable in terms of the design requirements. For instance, let's consider the stack-up of the previous example while changing the dimensions to 10mm by 10mm. This can be translated by a simple modification of the number of horizontal and vertical elements in the Perl software. The structure was simulated up to 40GHz and the numbers of horizontal and vertical elements are assumed to be 100 MCTL on each side, having the same length and RLGC matrices definition as the previous example. From Fig. 11, good agreement is evident up to 40GHz with some deviation in details above that frequency. As mentioned in section I, MCTL Matrix method enables memory and computation time savings. The proposed model substantially reduces the CPU run time that requires only 10.89 seconds on the Intel Xeon machine to simulate this Power/Ground plane pair, for which Momentum full-wave solver would have taken 32 min for a mesh density of 50 cells per wavelength. Our approach can reduce also by about minimum six times the memory amount required compared to Momentum. Indeed, it uses only 191,02Mbytes to simulate this Power/Ground Planes pair, however, Momentum took about 1156 Mbytes. The detailed comparison is summarized in Table I. Note that these results have been improved comparing to our last publication [START_REF] Bouchaala | Efficient method for predictive package Power/Ground planes modeling using Multi-Conductor Transmission Lines[END_REF] thanks to the exact calculations of the RLGC matrix which allowed rapid convergence. This huge difference between the resources usage comes from the long process for creating and simulating structures in a full-wave solver: first of all, engineer must draw the layout or imports it from another design simulator. Substrate characteristics must be defined layer by layer if we do not have a ready one. User must also assign the ports definition. Once the design is ready, Momentum calculates the Green's functions that characterize the substrate and generate circuit mesh for a specified frequency range. To study the efficiency of the method to model different stack-up, we consider the case of dielectric thickness of 100μm and the permittivity was fixed to 5. Conductor thickness is assumed to be 50μm.

RLGC Matrices are given by: 𝐿𝐿 𝑖𝑖𝑖𝑖𝑖𝑖 (𝐻𝐻/𝑚𝑚) = [ 6.88𝑒𝑒 -07 9.96𝑒𝑒 -09 9.96𝑒𝑒 -09 6.88𝑒𝑒 -07 ] ; 𝐿𝐿 𝑒𝑒𝑒𝑒𝑖𝑖 (𝐻𝐻/𝑚𝑚) = [ 13.76𝑒𝑒 -07 1.9𝑒𝑒 -08 𝐺𝐺 𝑖𝑖𝑖𝑖𝑖𝑖 (𝑆𝑆/𝑚𝑚) = [ 2.77𝑒𝑒 -13 -2.77𝑒𝑒 -13 -2.77𝑒𝑒 -13 2.77𝑒𝑒 -13 ];𝐺𝐺 𝑒𝑒𝑒𝑒𝑖𝑖 (𝑆𝑆/𝑚𝑚) = [ 1.38𝑒𝑒 -13 -1.38𝑒𝑒 -13 -1.38𝑒𝑒 -13 1.38𝑒𝑒 -13 ]

The reflection and transmission parameters are given by Fig. 10. We can easily see that we keep the same accuracy of the model as for thin dielectric. So our model is totally independent of the technological parameters and can be used for any substrate definition. Another advantage is that the simulation time is invariant since the number of layers remains the same. MCTL Matrix has used the same CPU and memory resources as for the thin structure to simulate the thick dielectric.

( meshing density of 50 cells/wavelength. For the memory resources, only 305.67 Mbytes are needed by our approach, which corresponds to about 8 times lower than Momentum that used 2544 Mbytes. So the complexity has increased sharply with Momentum, while MCTL Matrix has kept nearly the same complexity, without forgetting the required changes and adjustments on the layout structure which can be arduous for complex packages. In addition, a full substrate re-analysis is needed before performing the electromagnetic meshing. Fortunately, changing structure size in the MCTL Matrix can be performed by a simple update in the Perl software.

Extension to Multilayered Power/Ground Planes

The MCTL matrix method can be extended to a third dimension without any limit of the number of layers. In a single plane pair, the fundamental noise coupling occurs in the horizontal direction but also in the vertical direction. The vertical noise generation is emphasized in the multilayered packages and need to be carefully modeled.

Let's consider the three-layered structure depicted in Fig. 12. It consists of three solid planes of 20μm for the thickness, separated by 30μm of dielectric with permittivity equal to 5 and Tanδ = 0.002.

the same. MCTL Matrix has used the same CPU and memory resources as for the thin structure to simulate the thick dielectric.

(a)

agreement is evident up to 40GHz with some deviation details above that frequency. The deviation is important high frequency comparing to the previous test structure. Th can be explained by higher edge radiation since the boundari are larger. From computation time cost point of view, MCT matrix method took 34.66 seconds to simulate this structu however, Momentum simulated it in more than one hour wi a meshing density of 50 cells/wavelength. For the memo resources, only 305.67 Mbytes are needed by our approac which corresponds to about 8 times lower than Momentu that used 2544 Mbytes. So the complexity has increas sharply with Momentum, while MCTL Matrix has kept near the same complexity, without forgetting the required chang and adjustments on the layout structure which can be arduo for complex packages. In addition a full substrate re-analy is needed before performing the electromagnetic meshin Fortunately, changing structure size in the MCTL Matrix c be performed by a simple update in the Perl software. Another advantage of the MCTL Matrix method is that all Power/Ground planes dimensions can be simulated by using the same unit cell topology and sizing which makes our method easy customizable in terms of the design requirements. For instance, let's consider the stack-up of the previous example while changing the dimensions to 10mm by 10mm. This can be translated by a simple modification of the number of horizontal and vertical elements in the Perl software. The structure was simulated up to 40GHz and the numbers of horizontal and vertical elements are assumed to be 100 MCTL on each side, having the same length and RLGC matrices definition as the previous example. From Fig. 11, good agreement is evident up to 40GHz with some deviation in details above that frequency. The deviation is important at high frequency comparing to the previous test structure. This can be explained by higher edge radiation since the boundaries are larger. From computation time cost point of view, MCTL matrix method took 34.66 seconds to simulate this structure, however, Momentum simulated it in more than one hour with a meshing density of 50 cells/wavelength. For the memory resources, only 305.67 Mbytes are needed by our approach,

Extension to Multilayered Power/Ground Planes

The MCTL matrix method can be extended to a third dimension without any limit of the number of layers. In a single plane pair, the fundamental noise coupling occurs in the horizontal direction but also in the vertical direction. The vertical noise generation is emphasized in the multilayered packages and need to be carefully modeled. Let's consider the three-layered structure depicted in Fig. 12. It consists of three solid planes of 20μm for the thickness, separated by 30μm of dielectric with permittivity equal to 5 and Tanδ = 0.002. The p.u.l RLGC Matrices has been calculated at 10GHz as following.

Reflection and transmission scattering parameters were simulated up to 50 GHz. Fig. 13 shows that there is an excellent agreement up to 40 GHz regarding the reflection and the transmission parameters obtained from Momentum and the proposed method. Some deviation is evident beyond 40GHz as expected. For the computation cost, one can expect that the CPU time will increase with the layer number. Nonetheless, the simulation took only about 11 seconds to simulate this three-layered structure, which corresponds nearly to the same time taken to simulate a two-layered Power/ Ground Planes. However, Momentum simulator has spent more than one hour to analyze the structure. Nearly constant CPU time has been verified to simulate a five-layered structure. This constitutes a main advantage of our approach; the designer can avoid the tremendous computing time that a full-wave simulator would take to extend the analysis to higher number of layers.

IV. Conclusion

This paper described the physical principles and implementation of the MCTL Matrix method. The proposed method has been successfully implemented for a variety of planes structures and geometries and has been shown to yield significant savings in memory and run time costs without sacrificing accuracy. The simple implementation and the easy customization of the method over a full-wave simulator make it a promising approach, but the main contribution of the method remains in the fact that the electrical performance of the PDN in a high speed electronic product can be predicted before the design is started. In other words, Power Integrity can be checked at the pre-layout design stage. Consequently, the electronic product can be designed with this methodology to satisfy the Power Integrity specifications which lead to a shorter design cycle and rapid time-to-market. Our approach can be used as well for modeling on chip Power/Ground grids since they have already the shapes and the dimensions as the MCTLs. The p.u.l RLGC Matrices has been calculated at 10GHz as following.

𝐿𝐿 𝑖𝑖𝑖𝑖𝑖𝑖 (𝐻𝐻/𝑚𝑚) = [ 2.68𝑒𝑒 -07 1.79𝑒𝑒 -07 1.79𝑒𝑒 -07 3.5𝑒𝑒 -07 0.40𝑒𝑒 -07 1.79𝑒𝑒 -07 0.40𝑒𝑒 -07 1.79𝑒𝑒 -07 2.68𝑒𝑒 -07 ] 𝐿𝐿 𝑒𝑒𝑒𝑒𝑖𝑖 (𝐻𝐻/𝑚𝑚) = [ 5.36𝑒𝑒 -07 3.58𝑒𝑒 -07 3.58𝑒𝑒 -07 7𝑒𝑒 -07 0.82𝑒𝑒 -07 3.58𝑒𝑒 -07 0.82𝑒𝑒 -07 3.58𝑒𝑒 -07 5.36𝑒𝑒 -2.3𝑒𝑒 -13 -2.3𝑒𝑒 - 13 1.85𝑒𝑒 -12 -9.24𝑒𝑒 -13 -9.24𝑒𝑒 -13 -9.24𝑒𝑒 -13 -9.24𝑒𝑒 - 13 1.38𝑒𝑒 -12 ] 𝐺𝐺 𝑒𝑒𝑒𝑒𝑖𝑖 (𝑆𝑆/𝑚𝑚) = [ 6.9𝑒𝑒 -13 1.15𝑒𝑒 -13 -1.15𝑒𝑒 -13 4.59𝑒𝑒 -14 -4.62𝑒𝑒 -13 -4.62𝑒𝑒 -13 -4.62𝑒𝑒 -13 -4.62𝑒𝑒 -13 6.9𝑒𝑒 -13 ] main advantage of our approach; the designer can avoid the tremendous computing time that a full-wave simulator would take to extend the analysis to higher number of layers. the e to s shor can sinc MCT [START_REF] Roy | Delay-Extraction-Based Waveform Relaxation Algorithm for Fast Transient Analysis of Power Distribution Networks[END_REF] [2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

Fig. 1 :

 1 Fig. 1: Typical multilayered Power/Ground Planes. This architecture is the base of the proposed simulation.

Fig. 1 :

 1 Fig.1: Typical multilayered Power/Ground Planes. This architecture is the base of the proposed simulation.

Fig. 2 :Fig. 2 :

 22 Fig. 2: Multilayered Power/Ground planes discretization into MCTL unit cells.A big misconception in multiple plane pair structures is

Fig. 3 :

 3 Fig. 3: FEM EMPro current distribution intensity.

Fig. 4 :

 4 Fig. 4: MCTL Matrix method modeling flow. It includes MCTL calibration, W-element extraction, HSPICE simulation.

Fig. 3 :

 3 Fig. 3: FEM EMPro current distribution intensity.

Fig. 3 :

 3 Fig. 3: FEM EMPro current distribution intensity.

Fig. 4 :

 4 Fig. 4: MCTL Matrix method modeling flow. It includes MCTL calibration, W-element extraction, HSPICE simulation.

Fig. 4 :

 4 Fig. 4: MCTL Matrix method modeling flow. It includes MCTL calibration, W-element extraction, HSPICE simulation.

Fig. 5 :

 5 Fig. 5: Multi-conductor transmission line segment of length Lu.The length Lu is supposed much higher than the wavelength of the high frequency signal

Fig. 5 :

 5 Fig. 5: Multi-conductor transmission line segment of length L u .The length L u is supposed much higher than the wavelength of the high frequency signal.

Fig. 5 :

 5 Fig. 5: Multi-conductor transmission line segment of length Lu.The length Lu is supposed much higher than the wavelength of the high frequency signal

Fig. 5 :

 5 Fig. 5: Multi-conductor transmission line segment of length Lu.The length Lu is supposed much higher than the wavelength of the high frequency signal

Fig. 5 :

 5 Fig. 5: Multi-conductor transmission line segment of length Lu.The length Lu is supposed much higher than the wavelength of the high frequency signal

Fig. 5 :

 5 Fig. 5: Multi-conductor transmission line segment of length Lu.The length Lu is supposed much higher than the wavelength of the high frequency signal

Fig. 6 :

 6 Fig. 6: Mutual near and far mutual inductance.The general form of the L inductance is presented as follows:

Fig. 6 :

 6 Fig. 6: Mutual near and far mutual inductance.The general form of the L inductance is presented as follows:

Fig. 7 :

 7 Fig. 7: MCTL method mesh schema.

Fig. 7 :

 7 Fig. 7: MCTL method mesh schema.

Fig. 7 :

 7 Fig. 7: MCTL method mesh schema.

Fig. 7 :Fig. 7 :

 77 Fig.7: MCTL method mesh schema.

Fig. 6 :

 6 Fig. 6: Near and far mutual inductances supposed much higher than the wavelength of the high frequency signal

Fig. 8 :

 8 Fig. 8: 3mm by 3mm Power/Ground Plane pair. RLGC Matrices calculated at 10GHz are given by: 𝐿𝐿 𝑖𝑖𝑖𝑖𝑖𝑖 (𝐻𝐻/𝑚𝑚) = [ 3.51𝑒𝑒 -07 1.79𝑒𝑒 -07 1.79𝑒𝑒 -07 3.51𝑒𝑒 -07 ] ; 𝐿𝐿 𝑒𝑒𝑒𝑒𝑖𝑖 (𝐻𝐻/𝑚𝑚) = [ 7.02𝑒𝑒 -07 3.58𝑒𝑒 -07 3.58𝑒𝑒 -07 7.02𝑒𝑒 -07 ] 𝐶𝐶 𝑖𝑖𝑖𝑖 (𝐹𝐹/𝑚𝑚) = [ 1.41𝑒𝑒 -10 -1.41𝑒𝑒 -10 -1.41𝑒𝑒 -10 1.41𝑒𝑒 -10 ]; 𝐶𝐶 𝑒𝑒𝑒𝑒 (𝐹𝐹/𝑚𝑚) = [ 0.7𝑒𝑒 -10 -0.7𝑒𝑒 -10 -0.7𝑒𝑒 -10 0.7𝑒𝑒 -10 ] 𝑅𝑅0 𝑖𝑖𝑖𝑖𝑖𝑖 (Ω/𝑚𝑚) = [ 17.24 0 0 17.24 ] ; 𝑅𝑅0 𝑒𝑒𝑒𝑒𝑖𝑖 (Ω/𝑚𝑚) = [ 34.45 0 0 34.45 ] 𝑅𝑅𝑅𝑅 𝑖𝑖𝑖𝑖𝑖𝑖 (Ω/𝑚𝑚) = [ 3.28 0 0 3.28 ] ; 𝑅𝑅𝑅𝑅 𝑒𝑒𝑒𝑒𝑖𝑖 (Ω/𝑚𝑚) = [ 6.56 0 0 6.56 ]

Fig. 8 :

 8 Fig.8: 3mm by 3mm Power/Ground Plane pair.

Fig. 8 :

 8 Fig. 8: 3mm by 3mm Power/Ground Plane pair. RLGC Matrices calculated at 10GHz are given by: 𝐿𝐿 𝑖𝑖𝑖𝑖𝑖𝑖 (𝐻𝐻/𝑚𝑚) = [ 3.51𝑒𝑒 -07 1.79𝑒𝑒 -07 1.79𝑒𝑒 -07 3.51𝑒𝑒 -07 ] ; 𝐿𝐿 𝑒𝑒𝑒𝑒𝑖𝑖 (𝐻𝐻/𝑚𝑚) = [ 7.02𝑒𝑒 -07 3.58𝑒𝑒 -07 3.58𝑒𝑒 -07 7.02𝑒𝑒 -07 ] 𝐶𝐶 𝑖𝑖𝑖𝑖 (𝐹𝐹/𝑚𝑚) = [ 1.41𝑒𝑒 -10 -1.41𝑒𝑒 -10 -1.41𝑒𝑒 -10 1.41𝑒𝑒 -10 ]; 𝐶𝐶 𝑒𝑒𝑒𝑒 (𝐹𝐹/𝑚𝑚) = [ 0.7𝑒𝑒 -10 -0.7𝑒𝑒 -10 -0.7𝑒𝑒 -10 0.7𝑒𝑒 -10 ] 𝑅𝑅0 𝑖𝑖𝑖𝑖𝑖𝑖 (Ω/𝑚𝑚) = [ 17.24 0 0 17.24 ] ; 𝑅𝑅0 𝑒𝑒𝑒𝑒𝑖𝑖 (Ω/𝑚𝑚) = [ 34.45 0 0 34.45 ] 𝑅𝑅𝑅𝑅 𝑖𝑖𝑖𝑖𝑖𝑖 (Ω/𝑚𝑚) = [ 3.28 0 0 3.28 ] ; 𝑅𝑅𝑅𝑅 𝑒𝑒𝑒𝑒𝑖𝑖 (Ω/𝑚𝑚) = [ 6.56 0 0 6.56 ]𝐺𝐺 𝑖𝑖𝑖𝑖𝑖𝑖 (𝑆𝑆/𝑚𝑚) = [ 8.92𝑒𝑒 -13 -8.92𝑒𝑒 -13 -8.92𝑒𝑒 -13 8.92𝑒𝑒-13 ]; 𝐺𝐺 𝑒𝑒𝑒𝑒𝑖𝑖 (𝑆𝑆/𝑚𝑚) = [ 4.39𝑒𝑒 -13 -4.39𝑒𝑒 -13 -4.39𝑒𝑒 -13 4.39𝑒𝑒-13 ]

Fig. 9 :

 9 Fig. 9: Scattering parameters of thin dielectric layer (a) Reflection parameters; (b) transmission parameters

FREQUENCYFig. 9 :

 9 Fig.9: Scattering parameters of thin dielectric layer (a) Reflection parameters; (b) Transmission parameters.

FREQUENCYFig. 10 :

 10 Fig. 10: Scattering parameters of thick dielectric layer reflection parameters; (b) transmission parameters. One can notice the good agreement between both techniques. c. Invariant Unit Cell for Different Power/Ground Planes Sizes

FREQUENCYFig. 10 :Fig. 11 :

 1011 Fig.10: Scattering parameters of thick dielectric layer (a) Reflection parameters ; (b) Transmission parameters. One can notice the good agreement between both techniques.

Fig. 11 :Fig. 12 :FREQUENCYFig. 12 :

 111212 Fig.11: Simulated and modeled S-parameters for 10mm by 10mm Power/ Ground planes (a) Reflection parameters; (b) Transmission parameters.

Fig. 13 :

 13 Fig. 13: Simulated and modeled S-parameters for three-layered Power/Ground planes (a) Reflection parameters; (b) Transmission parameters. The results are very close but the computation time is much lower.

Fig. 13 :

 13 Fig.13: Simulated and modeled S-parameters for three-layered Power/Ground planes (a) Reflection parameters; (b) Transmission parameters. The results are very close but the computation time is much lower.

  Where ε r is the relative permittivity of the dielectric material, 𝜀𝜀 0 represents the permittivity of the free-space, and d is the separation between conductors i and j.

	𝐶𝐶 = [	𝐶𝐶 10 + 𝐶𝐶 12 + ⋯ + 𝐶𝐶 1𝑖𝑖 -𝐶𝐶 12 ⋮ -𝐶𝐶 1𝑖𝑖	-𝐶𝐶 12 𝐶𝐶 20 + 𝐶𝐶 12 + ⋯ + 𝐶𝐶 2𝑖𝑖 ⋮ …	… -𝐶𝐶 1𝑖𝑖 -𝐶𝐶 2𝑖𝑖 ⋮ … -𝐶𝐶 𝑖𝑖𝑖𝑖
			𝑐𝑐 𝑖𝑖𝑖𝑖 =	( 𝜀𝜀 𝑟𝑟 .𝜀𝜀 0 .𝑙𝑙 𝑢𝑢 .𝑙𝑙 𝑢𝑢 𝑑𝑑	) ⁄ 𝑙𝑙 𝑢𝑢	(3)
	c. Conductance Matrix
	The G parameter is the dielectric-loss conductance per-unit-length. The G matrix has the same general form as the C Matrix
	𝐺𝐺 = [	𝐺𝐺 10 + 𝐺𝐺 12 + ⋯ + 𝐺𝐺 1𝑖𝑖 -𝐺𝐺 12 ⋮ -𝐺𝐺 1𝑖𝑖		-𝐺𝐺 12 𝐺𝐺 20 + 𝐺𝐺 12 + ⋯ + 𝐺𝐺 2𝑖𝑖 ⋮ …	… -𝐺𝐺 1𝑖𝑖 -𝐺𝐺 2𝑖𝑖 … -𝐺𝐺 𝑖𝑖𝑖𝑖 ⋮	]

]

C ij can be computed as:

  .

	c. Conductance Matrix The G parameter is the dielectric-loss conductance per-unit-length. The G matrix has the same general form as the C Matrix	thick cond diele cond mutu
	𝐺𝐺 = [	𝐺𝐺 10 + 𝐺𝐺 12 + ⋯ + 𝐺𝐺 1𝑖𝑖 -𝐺𝐺 12 ⋮ -𝐺𝐺 1𝑖𝑖	-𝐺𝐺 12 𝐺𝐺 20 + 𝐺𝐺 12 + ⋯ + 𝐺𝐺 2𝑖𝑖 ⋮ …	… -𝐺𝐺 1𝑖𝑖 -𝐺𝐺 2𝑖𝑖 … -𝐺𝐺 𝑖𝑖𝑖𝑖 ⋮	]	Whe On t
						𝐿𝐿 𝑖𝑖𝑖𝑖
						Whe the c
						is pr MCT requ strai MCT of th edge are a that boun inter

  .

	thic cond diel cond mut
	Whe On t
	𝐿𝐿 𝑖𝑖
	Whe the c
	is p MC requ strai MC of t edge are that bou inte

Table I CPU Runtime and Memory Cost Comparison.

 I

	Method	Momentum	MCTL Matrix
	CPU runtime	32 min	10.89s
	Memory cost	1156 Mbytes	191.02 Mbytes
	b. Thick Dielectric Layer

  The deviation is important at high frequency comparing to the previous test structure. This can be explained by higher edge radiation since the boundaries are larger. From computation time cost point of view, MCTL matrix method took 34.66 seconds to simulate this structure, however, Momentum simulated it in more than one hour with a

	to ns Hz
	of are nd 02, the 30
	13 3 ]
	ly) od to cy. der L), ge

Table I CPU

 I Runtime and Memory Cost Comparison

	Method CPU runtime Memory cost	Momentum 32 min 1156 Mbytes	MCTL Matrix 10.89s 191.02 Mbytes
	b. Thick Dielectric Layer