Time- and Frequency-Domain Evaluation of Stochastic Parameters on Signal Lines

This paper focuses on the derivation of enhanced transmission-line models allowing to describe, in time and frequency domain, a realistic interconnect with the inclusion of external uncertainties, like process variations or routing and layout uncertainties. The proposed method, which is based on the expansion of the well-known telegraph equations in terms of orthogonal polynomials, turns out to be accurate and more efficient than alternative solutions like Monte Carlo method in determining the transmission-line response sensitivity to parameters variability. Moreover, an implementation into standard circuit analysis tools such as SPICE is possible. Two application examples based on PCB structures of common use in commercial packages conclude the paper.


Introduction
Nowadays, the numerical simulation of interconnect structures is a fundamental step in the design phase due to the urging necessity to perform right-the-first-time designs. In this regard, several tools are available, although they are usually deterministic and this represents a strong limitation whenever manufacturing tolerances or uncertainties on design parameters cannot be neglected. In this framework, the stochastic analysis is a tool that is extremely useful in the early design phase for the prediction of the system performance and for setting realistic design margins. Relevant examples are provided by the process-induced variability or routing decisions that unavoidably impact on the performance of PCB lines.
The typical resource allowing to collect quantitative information on the statistical behavior of the circuit response is based on the application of the brute-force Monte Carlo (MC) method. Such an approach, however, is computationally expensive, and this fact prevents us from their application to the analysis of complex realistic structures.
Recently, an effective solution that overcomes the previous limitation has been proposed. This methodology is based on the polynomial chaos (PC) theory and on the representation of the stochastic solution of a dynamical circuit in terms of orthogonal polynomials [1]. PC technique enjoys applications in several domains of Physics, including computational electromagnetics [2]- [7]. We limit ourselves to mention that the authors of this contribution proposed an extension of PC theory to distributed structures described by transmission-line equations [8]. Furthermore, a SPICE implementation of this methodology and its extension to lumped circuit elements has been recently provided [9]. This paper extends the PC theory to the stochastic simulation of a realistic interconnected structure consisting of a cascade connection of distributed multiconductor transmission lines and lumped multiport circuits. Both frequency-and time-domain results are presented, being the latter obtained through either harmonic superposition or SPICE simulation. In addition, an extension in the formulation allowing to account for variability in transmission line lengths will be presented. The feasibility and strength of the proposed approach are then further validated by means of the stochastic simulation of two realistic structures, for which the effects of uncertainties in the system parameters are predicted and analyzed.

Polynomial Chaos Primer
This section provides a brief overview of the PC method. The idea underlying this technique is the spectral expansion of a stochastic function (intended as a given function of a random variable) in terms of a truncated series of orthogonal polynomials. Within this framework, any function H, carrying the effects of variability, can be approximated by means of the following truncated series [10] where {ϕ k } are suitable orthogonal polynomials expressed in terms of the random variable ξ. The above expression is defined by the class of the orthogonal bases, by the number of terms P + 1 (limited to the range from 2 to 20 for practical applications, depending also on the number of random variables considered) and by the expansion coefficients H k . The choice of the orthogonal basis relies on the distribution of the random variables. Uncertainties arising from fabrication tolerances turn out to be properly characterized in terms of Gaussian variability. Hence, the most appropriate orthogonal functions for the expansion (1) are the Hermite polynomials, the first three being ϕ 0 = 1, ϕ 1 = ξ and ϕ 2 = ξ 2 − 1, where ξ is the standard normal random variable with zero mean and unit variance. It is relevant to remark that any random parameter in the system, e.g., the substrate permittivity ε r , can be related to ξ as follows whereε r and σ ε are the mean value and standard deviation, respectively. The orthogonality property of Hermite polynomials is expressed by where δ kj is the Kronecker delta and < ·, · > denotes the inner product in the Hilbert space of the variable ξ with Gaussian weighting function, i.e., When the variability is more properly characterized in terms of uniform distribution, Legendre polynomials ϕ 0 = 1, ϕ 1 = ξ, ϕ 2 = 3 2 ξ 2 − 1 2 , etc., can be used. In this case, ξ in uniform in the range [−1, 1] and the inner product is defined as With the above definitions, the expansion coefficients H k of (1) are computed via the projection of H onto the orthogonal components ϕ k . It is worth noting that relation (1), which is a known nonlinear function of the random variable ξ, can be used to predict the probability density function (PDF) of H(ξ) via numerical simulation or analytical formulae [11].
The basic results of PC theory outlined above can be extended to the case of multiple independent random variables. However, for the sake of brevity, the formal development (consisting in the application of orthogonality relations to build higher dimensional polynomials as the product combination of polynomials in one variable) is omitted here and readers are referred to [8].

Stochastic Model for Distributed Lines
Under the quasi-TEM assumption (i.e., the electric and magnetic fields are nearly orthogonal to the propagation direction), which is valid for electrically-small cross-sections and weak conductor losses, Maxwell's equations reduce to two-dimensional problems. The wave propagation along distributed transmission lines is equivalently governed by the following differential equations in the Laplace domain [12]: also known as telegrapher's equations. The above equations are derived from Faraday's law and where c i is the contour between the ith conductor and the reference, while contour c ′ i surrounds the ith conductor just off its surface.
Moreover, Z = R + sL and Y = G + sC are the p.u.l. impedance and admittance matrices, depending on the geometrical and material properties of the structure. The entries of the p.u.l. inductance matrix L, conductance matrix G and of the capacitance matrix C are related to the fields external to the conductors and are determined as the static field solutions in the transverse plane assuming perfect conductors. The entries of the p.u.l. resistance matrix R are instead governed by the fields interior to the imperfect condcutors. Using classical circuit theory, (6) is equivalent, from an electrical point of view, to the RLGC model in Fig. 1, with (capital letters denoting the entries of the p.u.l. matrices in (6)). The solution of (6) is given by the combination of its chain parameter matrix (CPM), that writes where expm denotes the matrix exponential, with the boundary conditions defined by the port electrical relations at the line terminations. The CPM relates voltages and currents at the two extremities z = 0 and z = L, i.e. ] .
(13) By defining for notation convenience X(z, s) = [V(z, s), I(z, s)] T , (13) can be rewritten as When the problem becomes stochastic, in order to account for the uncertainties affecting the guiding structure, we must consider the p.u.l. parameters as random quantities, with entries depending on the random variable ξ. In turn, (6) becomes a stochastic differential equation, leading to randomly-varying voltages and currents along the line. Therefore, they also depend on ξ.
The expansion (1) of the p.u.l parameters and of the unknown voltage and current variables in terms of Hermite polynomials, yields a modified version of (6), whose second row becomes where a second-order expansion (i.e., P = 2) is assumed; the expansion coefficients of electrical variables and of p.u.l. parameters are readily identifiable in the above equation.
Projection of (15) and of the companion relation arising from the first row of (6) on the first three Hermite polynomials leads to the following augmented system, where the random variable ξ does not appear explicitly, due to the integral projection form given in (4): .
T collect the different coefficients of the polynomial chaos expansion of the voltage and current variables. The new p.u.l. matrixỸ turns out to bẽ and a similar relation holds for matrixZ. It is worth noting that (16) is analogous to (6) and plays the role of the set of equations of a multiconductor transmission line with a number of conductors that is (P + 1) times larger than those of the original line. It should be also remarked that the increment of the equation number is not detrimental for the method, since for small values of P (as typically occurs in practice), the additional overhead in handling the augmented equations is much less than the time required to run a large number of MC simulations.
As far as the solution of the stochastic problem is concerned, the augmented equation (16) is used in place of (6), as well as the corresponding CPM, that becomes The augmented CPM relates the coefficients of the voltage and current variables (i.e.,X(z, s) = [Ṽ(z, s),Ĩ(z, s)] T ) at the line extremities. Readers are referred to [8] for additional details on the formulation for multiple random variables.

Uncertainty in the Interconnect Length
The above strategy can take into consideration any variation in the p.u.l. parameters, i.e., in the interconnect crosssection. Nonetheless, they can be adapted to account also for uncertainties in the interconnect length, whose inclusion is of paramount importance when the routing decisions are not a priori known [14].
From both Fig. 1 and (12), it is clear that in principle nothing changes if Z, Y and L are replaced by Z ′ = Z · L, Y ′ = Y · L and L ′ = 1 m, respectively. Yet, thanks to this definition, variations of the length L can be associated to Z ′ and Y ′ , that act as the p.u.l. parameters of an equivalent line of length 1 m.
As an example, if the random interconnect length can be expressed by means of the following series it is possible to study the equivalent line whose random p.u.l. parameters are given by which are analogous to (1). Therefore, the same procedure described above can be used. In the general case, i.e., when the p.u.l. parameters are random themselves and also expressed according to (1), the following relations can be used to determine the expansion coefficients of the equivalent p.u.l. parameters Z ′ and Y ′ : where I denotes the identity matrix, whileZ andỸ are the augmented parameters associated to the original line.

Stochastic Model for Lumped Blocks
We now suppose to have a lumped block described by transmission equations: where X 1,2 (s) = [V 1,2 (s), I 1,2 (s)] T and When we intend to include the variability of its parameters, we must consider the elements as random quantities and expand the above equations in terms of Hermite orthogonal polynomials, as already done for the case of distributed lines. By introducing a new random variable η, possibly different from the variable ξ that affects the distributed part, the second-order expansion of (22) yields (24) Again, projection onto the first three Hermite polynomials leads to the following deterministic augmented system T collects the coefficients of the PC expansion of the voltage and current variables. The four blocks ofT C turn out to bẽ with i, j = 1, 2. Also in this case, the augmented equations belong to the same class as the initial ones, and the system matrices have the same structure as in (17). The extension to multiple random variables is in principle obtainable as for the case of distributed structures.

Stochastic Simulation of Cascaded Structures
This section summarizes the procedure for the stochastic simulation of a complex cascaded structure. The proposed strategy is the following: (i) generate augmented stochastic models for the different parts composing the cascaded structure, by expanding their characteristics as described in Sections 3 and 4; (ii) simulate the entire structure in the frequency domain by suitably concatenating these models. If a time-domain calculation is desired, in principle this can be accomplished through harmonic superposition.

Incorporation of the Boundary Conditions
The simulation of a cascaded structure amounts to combining the port electrical relations of the two terminal elements defining sources and loads with the overall CPM, resulting from the cascaded connection of the intermediate elements.
According to the properties of the CPM, the overall matrix is given by the product of the matrices of the individual blocks. For the stochastic case, this leads tõ It should be noted that, in general, the generation of each extended model in (27) must be carried out including all the random variables in the problem, whether they affect all blocks or only a part of them, in order to have a consistent representation of all blocks for the concatenation.
When dealing with augmented models, also the port relations -although here no variability is assumed on the terminations -must be written in an extended form, in order to allow their combination with (27). Assuming to express the port excitations in terms of Thévenin equivalents, the port relations at the left-and righ-end sides become . .]) collect the port excitation voltages and the port impedances, respectively. These deterministic terms can be considered as a zero-order expansion, thus proportional to ϕ 0 . Expansion of (28) and the subsequent projection onto Hermite polynomials leads to where in this specific caseṼ S (s) = [V S (s), 0 . . . 0] T , while theZ S (s) are diagonal, due to the absence of variability on the termination networks. If needed, the above equations can be suitably modified to account for variability in the characteristics of the source elements. The solution of the stochastic system is achieved via the standard procedure used for combining the boundary conditions with the equations given by the CPM [12]. As an example, the left-and right-end currents can be computed as where whereas the voltages are obtained from (29).

Frequency-Domain Stochastic Analysis
Once the unknown voltages and currents are computed, the quantitative information on the spreading of circuit responses can be readily obtained from the analytical expression of the unknowns. As an example, the frequencydomain solution of the right-end ith terminal voltage V Ri , arising from (29) and (30) with a single random variable ξ and P = 2, is where the second numerical index denotes the expansion term. The above relation can be used for instance to compute the PDF of |V Ri (jω, ξ)|.

Time-Domain Stochastic Analysis
Several techniques can be adopted to apply PC modeling also to time-domain simulation, depending on the waveform and circuit characteristics.

Harmonic Superposition
The time-domain response can be readily obtained from the frequency-domain solution by considering a periodic input source and expressing it in terms of a truncated Fourier series where c 0 is the signal average over one period and c n is the complex Fourier coefficient for the n-th harmonic at frequency f n . The system being linear, its time-domain behaviour is in principle obtainable by the superposition of the analyses carried out for all signal harmonics. For the individual solution at frequency f n , the voltage sources appearing in (28), are replaced by their n-th harmonic components, i.e., E i (s n ) = E i (j2πf n ) = c i,1,2 . The timedomain expression of the output voltage v Ri (t, ξ) is then obtained as a linear superposition of harmonics: where v 0 is the DC term, obtained from a DC calculation, and the output coefficients V Ri (j2πf n , ξ) are computed according to (32). The linearity of Fourier decomposition assures that the PC structure is preserved also for the timedomain expression of the output: v where with k = 0, 1, 2.

SPICE Implementation
An alternative solution, also suitable for non-periodic inputs, is the implementation of the augmented equations (16) and (25) into standard analysis tools such as SPICE [9]. For distributed lossless lines, this can be easily achieved by means of the generalized Branin's equivalent for multinconductor transmission lines [12], which is suitable for handling also non-symmetric p.u.l. matrices, as typically occurs in the generated augmented systems (cfr. (17)). Readers are referred to [9] for the derivation of SPICE stochastic models for linear lumped elements. Of course, the created models are also suited to carry out frequency-domain simulations.
Besides new augmented circuit models, also new boundary conditions (i.e., terminations) must be defined in order to solve the overall resulting system (27). As shown in Sec. 5.1, when considering linear lumped terminations in the form of Thévenin or Norton equivalents and when no variability is attributed to them (as considered in this contribution), this simply amounts to replicating the original termination impedances on the additional conductors. On the other hand, no replication is required for independent sources.

Applications and numerical results
In this section, the proposed technique is applied to the analysis of two different PCB structures. Both frequencydomain and time-domain simulations are presented, carried out by means of both MATLAB and SPICE. The randomness is provided by the substrate parameters, i.e., h and ε r , that are considered to be the same for both the transmission lines, as well as by the lumped capacitance C. These parameters are considered as three independent Gaussian random variables with a relative standard deviation of 10%. The total number of terms P + 1 (corresponding also to the magnification of the size with respect to the original system) is given by where n is the number of random variables and p is the order of accuracy, that represents the maximum degree of the polynomials used for the expansion. Three augmented models are built for the distributed lines and the intermediate lumped section, which are indicated by their transmission matrices T TL1,2 and T C , respectively. The approximate relations given in [15] were used to numerically compute the PC expansion of the p.u.l. parameters of the coupled microstrips, whereas the expansion of the CPM for the lumped block is obtained analytically.
In the first example the structure is analyzed in frequency domain. Figure 4 shows the transfer functions between the voltage source and the two right-end terminations. The black thick lines represent the response of the structure for the nominal values of its parameters, while the thinner black lines indicate the limits of the ±3σ bound (σ being the standard deviation) determined from the results of the proposed technique. Finally, a qualitative set of 100  Magnitudes of |V d1 (jω)/E 1 (jω)| and |V d2 (jω)/E 1 (jω)| for the variability of substrate parameters and connector capacitance. Solid black thick line: deterministic response; solid black thin line: ±3σ limits of the third-order PC expansion; gray lines: a sample of responses obtained by means of the MC method (limited to 100 curves, for graph readability).
MC simulations is plotted using gray lines. Clearly, the parameter variations lead to a spread in the transfer function, that is well predicted by the estimated 3σ limits.
Often the knowledge of the standard deviation represents a limited information, since the quantitative information about how the values are distributed is missing. Nonetheless, from the analytical PC model we can also obtain the probability density function of the system responses. Figure 5 compares the PDFs of |V d2 /E 1 | computed at three different frequencies by means of the PC technique with those generated after 40,000 MC simulations. The frequencies correspond to the dashed vertical lines shown in Fig. 4.
The good agreement between the actual and the predicted PDFs and, in particular, the accuracy in reproducing the tails and the large variability of non-Gaussian shapes of the reference distributions, confirm the potential of the proposed method. For this example, it is also clear that a PC expansion with order 3 (i.e., P + 1 = 20) is already accurate enough to capture the dominant statistical information of the system response. Figure 6 shows the results for a time-domain analysis on the same structure. In this case the time-domain voltage source e 1 (t) is a trapezoidal wave with an amplitude of 1 V, rise and fall times of 100 ps, duty cycle 50% and total period of 5 ns. Harmonic superposition with N = 30 harmonics is considered for the calculations.  : v d1 (t) and v d2 (t). Solid black thick line: deterministic response; solid black thin line: ±3σ limits of the third-order PC expansion; gray lines: a sample of responses obtained by means of the MC method (limited to 100 curves, for graph readability). and frequency domain. For fairness, the table includes also the overhead introduced by PC due to the building of the augmented matrices. Nevertheless, the speed-up observed is over two orders of magnitude.

Structure #2
Layout and routing uncertainties may significantly affect signal integrity in interconnects at the package and board level as well [14]. Let us consider the geometry in Fig. 8, that represents the wiring between two components on the same substrate. The cross-section of this microstrip interconnect is shown in Fig. 9. Several sources of uncertainty may be present. First, routing decisions as well as the displacement between the two components impact on the in-  terconnect length. For high-speed applications, the general trend is to minimize the path in order to avoid degradations. Therefore, a variation in the range 5 ± 0.5 cm is a reasonable assumption. In the example of Fig. 8, a second source of uncertainty is provided by the pad pitch for the two components, that in fact defines the trace separation. In this application, it is supposed to vary in the range 2.5 ± 1 mm. Finally, a fine control of the substrate properties is sometimes prohibitive, depending on the process technology. For instance, variations of the permittivity in the range 4.5 ± 0.25 are likely to be expected. These three sources of variation represent an example of the most influential that can be expected in a problem like the one sketched in Fig.  8. All the remaining parameters are considered to assume fixed values, i.e., w = 1 mm, h = 0.2 mm and t k = 0.03 mm.
In order to validate the flexibility of the proposed approach, a SPICE polynomial chaos model for the interconnect in Fig. 8 is created that accounts for the aforementioned variations, which are all considered as uniform. As such, Legendre polynomials represent the optimal choice  for the basis {ϕ k }. First, a transient simulation is performed. The center trace is driven at the near-end side with a current source in shunt connection with a 50-ohm resistance. The source waveform is a Gaussian pulse of peak amplitude of 0.1 A, occurring at about t = 320 ps, and 3 dB width of about 60 ps. All the other terminations are RC loads with a 50-ohm resistance and a capacitance of 10 pF. Figure 10 shows the mean value and the standard deviation for the signal propagated at the far-end side v F E (t), computed with HSPICE both from the PC model and by means of 10,000 simulations carried out with the available MC feature. Figure 11 shows the PDF of the far-end crosstalk v F EXT (t) computed at t = 650 ps instead. Of the two PC curves, the one marked (ii) is computed by increasing the expansion order from 2 to 3 w.r.t. the one marked (i), thus showing that the number of terms can be effectively increased to achieve a better accuracy. A frequency sweep is also performed for the structure under consideration. In this case, the transient current source is replaced by a phasor with amplitude 1 A and phase zero. All the other parameters are left unchanged. The interconnect system is simulated up to 5 GHz. Mean value and ±3σ estimations are reported in Fig. 12 for the far-end crosstalk. A set of MC responses is included as well, thus confirming that the 3σ bound provides a valuable information about the spread of the response. Figure 13 shows a PDF of the signal propagated on the driven trace instead.
Finally, Tab. 2 collects the computational times and the memory required by the simulations. From these results it is clear that PC allows to obtain the same information as MC, with comparable accuracy but a great efficiency improvement.

Conclusions
This paper addresses the uncertainty-aware analysis for the design of realistic board-level interconnections. The pro-   posed approach is based on the expansion of the governing equations onto a basis of orthogonal polynomials. The result is a modified model that inherently takes possible random variations of the interconnect parameters into account. Nonetheless, the generated stochastic models can be implemented and solved into standard circuit analysis tools such as SPICE. The method offers comparable accuracy and improved efficiency with respect to alternate methods like Monte Carlo simulations. Two realistic application examples involving high-speed PCB links are used to demonstrate the feasibility and strength of the approach.