Efficient 3D High-Frequency Impedance Extraction ... - Semantic Scholar
Efficient 3D High-frequency Impedance Extraction for General Interconnects and Inductors Above a Layered Substrate Navin Srivastava
Roberto Suaya
Kaustav Banerjee
Mentor Graphics, Wilsonville, OR, USA navin [email protected]
Mentor Graphics, Grenoble, France roberto [email protected]
University of California, Santa Barbara, CA, USA [email protected]
Abstract—We present an efficient and highly accurate approach to highfrequency impedance extraction for VLSI interconnects and intentional on-chip inductors. The approach is based on a three-dimensional (3D) loop formalism that uses discrete complex images approximations applied to a quasi-magnetostatic treatment of the vector potential, resulting in closed-form expressions for the impedance matrix of current filaments in the presence of a multi-layer substrate. Populating the impedance (Z) matrix for 3D configurations of finite transverse dimensions (including nonManhattan wires and inductors) is computationally inexpensive, and includes substrate eddy current effects that become quantitatively important in the frequency regime beyond 20 GHz which is imminent at the 45 nm technology node onwards. The accuracy, as exemplified by the magnitude of inductor impedance |Z|, is within 5% of a full-wave electromagnetic field solver for frequencies up to 100 GHz, with an order of magnitude lower computation cost. The proposed method represents a core technology for incorporation into system level extraction of analog systems consisting of multiple inductors and nearby interconnects, for CMOS on-chip circuits in the nanometer era.
I. I NTRODUCTION As VLSI technology scales beyond the 45 nm node [1], faster transistors drive signals on interconnects with significant frequency components approaching 100 GHz. At the same time, unprecedented levels of functional integration enable large digital circuits, as well as sensitive analog circuits operating in the mm-wave frequency regime [1], to coexist on the same platform. While the frequency-dependent skin and proximity effects are already dealt with during system-level extraction of global on-chip interconnects in the 10-20 GHz range [2], hitherto unimportant electromagnetic phenomena are manifest beyond 20 GHz. At these frequencies, the magnetic fields generated by transient currents penetrating the conductive Silicon substrate result in eddy currents that have a significant impact on on-chip inductors as well as on timingcritical global interconnects. Accurate high-frequency characterization of inductors, including electromagnetic (EM) interactions with multiple nearby inductors and interconnects as often encountered in wireless applications, dictates that they be considered in toto including the substrate, as opposed to being analyzed as individual components in isolation. While computationally efficient analytical approaches have been proposed for impedance computation of 2D Manhattan-only interconnect structures [3], they are insufficient for on-chip inductors and 3D interconnects. Existing alternatives that are applicable to 3D geometries demand numerical quadratures for multi-dimension integrals such as those described in Section II. Accurate computation of the impedance matrix for multiple inductors and interconnects above a multi-layer substrate overwhelms EM field solvers. S-parameter measurements that are currently used for this purpose at 1-10 GHz severely limit design space exploration, and are quite challenging in the next frequency decade due to difficulties with accurate and reliable deembedding calibration. The demand is for computational techniques that correctly handle the substrate effect, while permitting the analysis of structures with complexity well beyond that of individual components. These
978-3-9810801-6-2/DATE10 © 2010 EDAA
considerations guide the search beyond finite element methods, such as HFSS [4], that can be used only for validation purposes. Full-wave field solvers based on boundary element techniques, such as Fastmaxwell [5], have many of the desired properties but lack the performance required for handling large structures. High-performance full-wave field solvers expand the two dimensional Fourier transform of the substrate Green’s function as a linear combination of complex exponentials, in an approximation called the Discrete Complex Images Method (DCIM) [6]. In FastMaxwell [5], DCIM is used in combination with higher order basis functions in a numerical approach to compute the full-wave substrate Green’s function, and applied to three-dimensional interconnects above a single-layer substrate. Quasi-static solvers that provide much better performance are commonly used at lower frequencies. The classical work of FastHenry [2] is the reference in magneto-quasi-static (MQS) applications, but it only allows for conductors and metallic ground planes, and cannot account for a large substrate with possibly high conductivity and permittivity. Our goal is to develop an analytical quasistatic approach capable of handling complex interconnect structures above a multi-layer substrate. To this end, we compute the impedance matrix of a configuration of wires using the Green’s function for a three-dimensional quasi-static magnetic dipole above a multi-layered substrate, and use DCIM to allow for analytical quadratures. The use of a magnetic dipole (closed current loop) [7] instead of an open current source guarantees the correct long distance behavior of electromagnetic fields. Since our computations are in the near-field regime (distances less than few mms) they are insensitive to surface waves that could make the DCIM inaccurate [8]. Furthermore, with the help of the superposition principle to represent finite three-dimensional interconnect loops as a continuous distribution of point-like magnetic dipoles, we obtain the self and mutual impedance of 3D interconnects above a multi-layered substrate in closed form. Our impedance computation methodology is described in Section III and implementation details in Section IV. In Section V, we compare our results with a full-wave field solver [4], after reintroducing the electric field generated by the displacement currents for simple structures where such comparisons are feasible. Besides the computational efficiency resulting from the closed form expressions derived here, our method seamlessly merges into existing loop impedance approaches based on the filaments decomposition of wires above the substrate, that work well at lower frequencies [9]. Hence, the work in this paper is a significant step towards realizing a system-level impedance extraction methodology covering a wide frequency spectrum from nearly D.C. to over 100 GHz, with reasonable computation cost. II. M AGNETIC DIPOLE VECTOR POTENTIAL G REEN ’ S FUNCTION For a general Manhattan interconnect, the collection of a signal wire and its return paths can be partitioned along their length to form bundles [7]. All wires in a bundle run parallel and have the same length. Each bundle consists of multiple current loops, each formed by a pair of
(x1,y2,z1)
z
Forward current on metal Mn
I ˆy x
y
(x2,y2,z2)
& p
Loop Area
0 where each Gs,R z,y corresponds to a single line current (monopole) on the right hand side of (3). To derive expressions for the Green’s function corresponding to the above magnetic dipole, it is convenient to first consider each current term in isolation. Using the 2D Fourier transform (in x and y), we obtain the monopole Green’s function as:
I ˆy Return current on metal Mn+2 ĭ
(x2,y1,z2)
(x1,y1,z1)
Fig. 1. Current loop formed by two parallel current filaments carrying current in opposite directions constitutes a planar magnetic dipole source.
parallel opposing current filaments. A current loop is the smallest atomic unit in our description of interconnect loop impedance. The physical equivalent of a planar current loop, as shown in Fig. 1, is a magnetic dipole [10]. Since on-chip conductors are confined to discrete metal layers in the x-y plane (neglecting vias which have small dimensions), the current loops they form are planar, although they may have arbitrary orientations (φ). For a planar loop carrying current I, its magnetic dipole moment p is perpendicular to the plane containing the loop, and its magnitude is given by [10]: | p| = I × LoopArea
(1)
Rotational symmetry allows us to choose the y-axis along the length of the wires forming the current loop. Hence, p has horizontal and vertical components: ˆ + pz zˆ = | p|[sin(φ)ˆ x + cos(φ)ˆ z] p = px x
(2)
0 Gs,R z,y (x, y, z, x , y , z ) =
∞ ∞ −∞ −∞
⎛
1 × (2π)2
⎞
ejkx (x−x )+jky (y−y ) ⎠ dkx dky ; ⎝ −μ × ejkz |z−z | + χN ejkz (z+z ) 2jkz kz2 = −kx2 − ky2
(6)
where kx and ky are the spectral variables corresponding to x and y, respectively, and χN can be found using the boundary conditions for the and magnetic field (H) at the substrate interface: magnetic flux (B)
R1 ) · zˆ
R0 − B (B
z=0
R0 − H R1 ) × zˆ
= 0; (H
z=0
= 0;
For a single layer substrate, we get: χ1 =
kz − kz 2 ; kz = −kx2 − ky2 − γ12 kz + kz
lim (ejkx a/2 − e−jkx a/2 )/a = lim (2j sin(kx a/2))/a = jkx
0 2 Gd,R z,y (x, y, z, x , y , z ) = ⎡ ⎤ a δ(x − x + )δ(z − z )− )δ(y − y (3) −2π ⎣ 2 ⎦ ; z, z > 0 lim a a→0 a δ(x − x − )δ(y − y )δ(z − z ) 2 Since there are no sources inside the substrate, the Green’s function in region Ri satisfies:
(2 − γi2 )GRi (x, y, z) = 0; z < 0
(4)
where = jωμ(σi + jω i ) for the i substrate layer . The dipole Green’s function in (3) can be expressed as: th
1
0 Gd,R z,y (x, y, z, x , y , z ) =
⎛ ⎞ a s,R0 −2π ⎝−Gz,y (x, y, z, x − 2 , y , z ) ⎠ lim a 0 a→0 μa + Gs,R ,y ,z ) z,y (x, y, z, x + 2
1 The
a→0
(9)
we get:
In 3D layered media, as in the case of a multi-layered substrate, the ¯ Green’s function for the vector potential is a second rank tensor, G, whose components are Gu,v , where u denotes the direction of the dipole source and v the direction of the resultant magnetic vector potential A. Our objective is to compute these quantities when the location of the source r , as well as that of the observation point r, are in the region R0 above the substrate. Consider two opposite currents I yˆ and −I yˆ in the x-y plane, centered at (x ,y ,z ) and separated by an infinitesimal distance a, lying above a multi-layered substrate. The two currents constitute a zˆ-directed magnetic dipole. The corresponding yˆ-directed vector potential Green’s 0 with I → ∞ and a → 0 such that the dipole moment function Gd,R z,y μ (Ia) = 1, satisfies: 2π
γi2
(8)
For an N -layer substrate, these boundary conditions must be applied at all interfaces between different layers [11]. Substituting (6) in (5), and using: a→0
A. Vector potential Green’s function in integral form
(7)
(5)
substrate layers have finite conductivity, and the permittivity of Si is large. Hence, γi2 cannot be neglected.
0 Gd,R z,y (x, y, z, x , y , z ) =
∞ ∞ −∞ −∞
⎛ kx
−μ × 2(2π)2
ejkx (x−x )+jky (y−y
⎝ kz
⎞
)
× ejkz |z−z | + χN ejkz (z+z
)
⎠ dkx dky
(10)
The above dipole Green’s function expression is the sum of a primary term corresponding to a magnetic dipole in free space, and a secondary term which is nearly identical in form (except for the coefficient χN ) arising due to the substrate: Fnet (r, r ) = Fprimary (r, r ) + χN Fsecondary (r, r )
(11)
Although the expressions for the primary and secondary field terms may vary depending on the source (such as electric or magnetic dipoles, or magnetic monopoles), the coefficient χN always appears unaltered [3], [6], [11]. B. Discrete complex images arising from the substrate Although analytical expressions are well-known for the primary field component in (10), analytical computation of the entire integral expression is hampered by the coefficient χN . Hence, the key to obtaining closed form expressions for the substrate Green’s function lies in finding a suitable approximation to χN that transforms the secondary field expression into an anaytically integrable form. An often favored approximation for χN - the DCIM [6], [12] - uses a linear combination of complex exponentials. The following modified DCIM extends several desirable properties to the approximation [3], including lower fitting time and smoothness across a wide range of frequencies: χN ≈
K k=1
βk e−αk (λ/γ1 ) ,
(12)
magnetic dipole source
magnetic dipole source
I
Indeed the same expression can be derived by using the definition of magnetic dipole moment [10] for a loop of area a lying in the x ˆ-ˆ y plane and carrying current I:
I z’
z=0
z’
Modified DCIM
ȕ 1I
K image dipoles
Multi-layer substrate
substrate
Į interface z’+ Ȗ1 1
z=0
1 p (r ) = 2
ȕ2I Į z’+ ȖK 1
ȕKI
where λ = kx2 + ky2 ∈ (0, ∞), γ1 = jωμ(σ1 + jω 1 ), and αk , βk are found by non-linear least squares fitting. The above approximation can be physically interpreted as the sum of K images of the source dipole representing the effect of the substrate, as shown in Fig. 2, and the net field in (11) takes the form: Fnet (r, r ) = Fprimary (r, r ) +
βk F˜primary (r, r , αk )
(13)
k=1
where F˜primary is the field due to an image of the source dipole and has a functional form similar to the term Fprimary . C. Analytical vector potential Green’s function in free space For a single elementary current in free space, the vector potential s,f ree is related to the corresponding Green’s function Gs,f ree (r, r ) A as:
s,f ree (r) = r )dr A Gs,f ree (r, r )J( (14) r ) is the current density at r (with components Jx along x where J( ˆ and Jy along yˆ, respectively). Similarly, for a dipole in free space, the vector d,f ree is written in terms of the dipole Green’s function as: potential A
d,f ree (r) = A
¯ d,f ree (r, r ) · p G (r )dr
(15)
¯ d,f ree (r, r ) is a where p (r ) is the dipole moment density at r and G second rank tensor. Now, for a single current in free space, the vector potential Green’s function is well-known [10]: 1 μ Gs,f ree (r, r ) = 4π |r − r | (16) |r − r | = (x − x )2 + (y − y )2 + (z − z )2 Noting that the vector potential Green’s function for a magnetic dipole consisting of two opposite currents I yˆ and −I yˆ, located at (x − a/2,y ,z ) and (x + a/2,y ,z ) separated by an infinitesimally small distance a, is given by the gradient between the Green’s functions corresponding to the two isolated currents, as shown in (5), we can use (16) to get an alternative representation of the dipole Green’s function in free space. Thus: ree (r, r ) = Gd,f z,y
⎛
μ 1⎜ ⎜ lim a→0 4π a ⎝
⎞
1 (x − x + a2 )2 + (y − y )2 + (z − z )2 1
−
(x −
x
−
a 2 ) 2
r )d3 r = zˆIa r × J(
(18)
and the Green’s function for a magnetic dipole [10]:
Fig. 2. Representation of the substrate as a series of images of the source magnetic dipole using the modified discrete complex images method [3].
K
+ (y −
y )2
+ (z −
⎟ ⎟ ⎠
ree = Gd,f z,y
μ (r − r ) × p 4π |r − r |3 | p|
Recognizing this simple form of the vector potential Green’s function (19) for on-chip interconnect currents in free space, we avoid the multiple dimension integrals in (10) and instead find analytical expressions for the dipole Green’s function in the presence of a multi-layered substrate in the following Section. III. I MPEDANCE E XTRACTION M ETHODOLOGY FOR 3D I NTERCONNECTS A BOVE A M ULTI - LAYER S UBSTRATE A. Analytical vector potential Green’s function for a 3D magnetic dipole in free space: The quasi-static magnetic vector potential at r due to a magnetic dipole located at a point r in free space is given by [10]: × (r − r ) r) = μ × p A( 4π |r − r |3
(x − x ) μ = 2 4π [(x − x ) + (y − y )2 + (z − z )2 ]3/2 ree ¯ d,f ree created by where Gd,f denotes the y-directed component of G z,y the z-directed quasi-static magnetic dipole.
(20)
¯ d,f ree (r, r )) gives the vector The vector potential Green’s function (G potential (A(r)) in free space, due to a unitary magnetic dipole (| p| = 1) ¯ d,f ree (r, r ), in terms of located at (r ). The relevant expressions for G ree (r, r ) (as explained before), are described below. its components Gd,f u,v Due to our choice of co-ordinates, the y-component of p is zero (2). Hence, ree ree ree = Gd,f = Gd,f =0 (21) Gd,f y,x y,y y,z ree is perpendicular Moreover, the cross product in (20) implies that Gd,f u,v to the source dipole moment p . In other words, the Green’s function due to an x-directed dipole source will have no component in the x-direction. Hence: ree ree ree = Gd,f = Gd,f =0 (22) Gd,f x,x y,y z,z
From the observations above, the vector potential Green’s function (in the x-y plane) at r = (x, y, z), due to a unit magnetic dipole source located at r = (x , y , z ) in free space, can be written as:
¯ d,f ree ¯ G
⎡
0
0 0 0
ree = ⎣ Gd,f x,y d,f ree Gx,z
⎤
ree Gd,f z,x ree ⎦ Gd,f z,y 0
(23)
Since we consider the interaction between wires that lie in the x-y plane ree along (vertical conductors or vias are ignored), the components Gd,f u,z d,f ree - the only non-zero component along zˆ are not relevant. Hence, Gx,z zˆ - is ignored. Substituting p from (2) in (20), and with | p| = 1, we get: ree (x, y, z, x , y , z ) = Gd,f x,y
μ sin(φ)(z − z ) 4π |r − r |3
(24)
ree (x, y, z, x , y , z ) = Gd,f z,y
μ cos(φ)(x − x ) 4π |r − r |3
(25)
(17)
z )2
(19)
ree (x, y, z, x , y , z ) = − Gd,f z,x
where the angle φ is shown in Fig. 1.
μ cos(φ)(y − y ) 4π |r − r |3
(26)
B. Analytical vector potential Green’s function for a 3D magnetic dipole above a multi-layer substrate: For the vector potential Green’s function in the presence of a multilayer substrate, we invoke the modified discrete complex images method [3] explained in Section II, to represent the substrate as a series of images of the source dipole (see Fig. 2). The procedure for computing accurate values for the parameters (αk ,βk ) in (12) is explained in the following section. The vector potential due to the combined effect of K +1 dipoles (source and its K images) is given by the superposition of the vector potential due to each dipole considered in isolation. Hence, the vector potential Green’s function at (x,y,z) due to a source dipole located at (x ,y ,z ) in the presence of a multi-layer substrate is given by: ¯ d,sub = G ¯ d,f ree + G
K
¯ d,img G k
(27)
r) along yˆ in (29), we get: Taking the components of A(
y4 M
d,f ree
C. Impedance for a finite current loop: The Green’s function shown above gives the vector potential due to a unit magnitude dipole source located at a point r = (x , y , z ). The vector potential at r = (x, y, z) due to a finite size current loop is given by the superposition of infinitesimally small point sources in the area occupied by the loop, each occupying an area 12 (dx × dy ) and carrying current I. The co-ordinates of the source (x ,y ,z ) span the rectangular (shaded) area shown in Fig. 1, whose extremities are defined by (x1 ,y1 ,z1 ), (x1 ,y2 ,z1 ), (x2 ,y1 ,z2 ) and (x2 ,y2 ,z2 ). Expressing the z-coordinate of the source in terms of x and φ: zx = z1 + (x − x1 )tan(φ), the vector potential at r = (x, y, z) due to the finite current loop in free space is given by (15), yielding:
ree d,f ree ree ˆAd,f yˆ Ad,f Gx,y (x, y, z) + AGz,y (x, y, z) + x Gz,x (x, y, z)
(29)
ree where each Ad,f Gu,v is given by:
y2 x2 =
ree Gd,f (x, y, z, x , y , zx )Idx dy u,v
d,f ree MG x,y
d,f ree MG z,y
r) · dl/I A(
M = ψ/I =
(31)
l
where dl is the length vector for an infinitesimal element of the victim conductor, r is the position vector for this element, and l is the contour along the conductor length. For a victim conductor oriented along yˆ extending from (x3 ,y3 ,z3 ) to (x3 ,y4 ,z3 ), (31) becomes:
y4 M
d,f ree
d,f ree (x3 , y, z3 ) · dˆ A y /I
= y3
μ = 4π
d,f ree MG tan2 (φ) z,y
μ = 4π
where, s=
3
(32)
+q
2
y =y2 x =x2
cos (φ)s − cos (φ)q
y =y1
x =x1
y =y2 x =x2 y =y1 x =x 1
√ m2 + v 2 − v v | + m2 + v 2 ; ln| √ 2 m2 + v 2 + v
⎛
zk cos(φ)tanh
⎜
q = sin(φ) ⎜ ⎝
+vtan−1
−1
√
t m2 + v 2
v t √ zk cos(φ) m2 + v 2
y=y4 (34) y=y3
y=y4
(35) y=y3
(36)
⎞ ⎟ ⎟ ⎠;
u = x3 − x ; v = y − y ; zk = (z3 − z1 ) − (x3 − x1 )tan(φ); t = usec(φ) + zk sin(φ); m2 = t2 + zk2 cos2 (φ) ¯ d,f ree in (15) must be In the presence of a multi-layer substrate, G ¯ d,sub from (27), which, when combined with (28), yields: replaced by G
y2 x2 d,sub (x, y, z) = A
+
K
¯ d,f ree (x, y, z, x , y , z ) · p G (r )dr x
y1 x1 y2 x2
βk
(37) ¯ d,f ree (x, y, z, x , y , z k ) · p G (r )dr x
y1 x1
= −(z1 + αk /γ1 ) + (x − x1 )tan(φ) is the z-coordinate where of the k image dipole. Since the integrands have the same form as those in free space, the subsequent integrals for the mutual impedance M d,sub in the presence of a multi-layer substrate can also be computed analytically: zxk th
d,f ree d,f ree + MG + M d,sub = MG x,y z,y
For the particular choice of co-ordinates used here (currents along yˆ), the ree term Ad,f Gz,x = 0. As expected, the resultant magnetic vector potential is directed along yˆ, just as the currents in the source loop. The mutual inductance between the source current loop and a victim conductor is given by:
(33)
All the integrals shown above can be evaluated in closed form:
(30)
y1 x1
dy
d,f ree d,f ree + MG = MG x,y z,y
k=1
d,f ree (x, y, z) = A
ree Ad,f Gu,v (x, y, z)
I y3
k=1
¯ d,f ree are given in (23)-(26). The expression for The expressions for G d,img ¯ Gk is easily obtained by relocating the source from (x ,y ,z ) to (x ,y ,−(z + αk /γ1 )), and multiplying by the linear coefficient βk : ¯ d,img (x, y, z, x , y , z ) = β G ¯ d,f ree (x, y, z, x , y , −z − αk ) (28) G k k γ1
=
ree d,f ree Ad,f Gx,y (x3 , y, z3 ) + AGz,y (x3 , y, z3 )
K
d,img d,img βk (MG + MG ) x,y z,y
(38)
k=1 d,img d,img and MG are given by the same expressions as those where MG x,y z,y d,f ree d,f ree for MGx,y and MGz,y , respectively, with the following modification to (36): (39) zk = (z3 + z1 + αk /γ1 ) − (x3 − x1 )tan(φ)
The equations above are sufficient to compute the mutual impedance of Manhattan interconnects by orienting the co-ordinate axes such that relevant conductors are parallel to yˆ. When computing the self inductance of a loop, the victim conductor coincides with the source and the center-to-center distance between the source and destination points (u = x3 − x in (36)) is then replaced by the geometric mean distance (GMD) [13] of a conductor with respect to itself (for rectangular crosssection GMD = elog(w+t)−3/2 [13]). The self impedance is then given by Zself = Rself + jωMself , where Rself = ρL( ws 1ts + ws 1ts ) 1 1 2 2 is the static resistance of the signal line (s1 ) and return path (s2 ) of the loop (ws1 , ws2 are widths and ts1 , ts2 are the thicknesses of the conductors).
Starting freq = fStart; Randomly assign starting values:
Besides Manhattan interconnects which are always inclined at right angles, we are also interested in computing the impedance of inductors, which may comprise conductor segments inclined at arbitrary angle θ. In this case, we choose a co-ordinate system that aligns the source current loop with yˆ. For the victim conductor extending from (x3 , y3 , z3 ) to (x4 , y4 , z3 ), the x-coordinate can be expressed in terms of its y-coordinate as: xy = x3 + (y − y3 )tan(θ). The mutual impedance is then given by: Mθd,f ree
μ μ = cos(φ)T1 + sin(φ)T2 ; 4π 4π
Run VP Get best fit
{O j , F j ( freq )}
Dk , Ek
D k (0,1)
Next freq = current freq +/- fStep; Next starting values = current best fits
START
NO
Reached end frequency? YES
Current freq = freq with min residual; Current best fits = fits at current freq
(40)
NO
Are all residuals < 1E-2 ?
YES
STOP
where,
−sinh−1
T1 =
y−y sxz
y3
y4
tan−1
T2 =
x =x2 y =y2 dy; x =x1
(xy − x )(y − y ) (z3 − zx )sxyz
sxyz =
(41)
y =y1
x =x2
y3
sxz =
Fig. 3. Algorithm flowchart for computing complex images approximation using Variable Projection method.
B. Impedance computation for realistic dimension interconnects
y =y2 dy;
x =x1
y =y1
(xy − x )2 + (z3 − zx )2 ;
(xy − x )2 + (y − y )2 + (z3 − zx )2 ;
xy = x3 + (y − y3 )tan(θ) For conductor loops in the x ˆ − yˆ plane (φ = 0), the integral in (40) R is computed analytically (with a symbolic integration tool Maple), although the resulting expression is too long to be included here. IV. I MPLEMENTATION OF I MPEDANCE E XTRACTION M ETHODOLOGY FOR VLSI I NTERCONNECTS
At frequencies and transverse dimensions where current distribution inside the conductors is non-uniform due to skin and proximity effects, each conductor is discretized into multiple thin filaments with piecewise constant current density for each filament [2]. The expressions for self and mutual impedance shown in Section III can then be used for such filaments. Each ith filament in the conductor carrying the forward current forms a loop with each j th filament in the conductor carrying the return current, resulting in L loops. The following linear system of equations is then solved to compute the net current I flowing through the two-conductor bundle:
⎡
R1 + jωM11 ⎢ .. ⎣ . jωML1
... .. . ...
⎤⎡
A. Modified discrete complex images for multi-layer substrates
⎡
⎤
1 .. ⎥ , . ⎦ 1
ηl = 1, ηl ∈ C
(42)
l
where Rl + jωMll is the self impedance of the lth loop, Mlp is the mutual impedance between the lth and pth loops, and ηl is the fraction of the total current I flowing through the lth loop. Finally, the impedance of the bundle is given by Z = I −1 . V. R ESULTS Fig. 4 shows the impedance of a square inductor configuration (combination of two orthogonal loops, Loop − x and Loop − y, shaded with different colors) lying above a 3-layer substrate. The input impedance is Cxy Loop-x
200
Loop-y
Rx/2 Lx/2
Ry/2
Ry/2
Rx/2 Lx/2
Ly/2
Cx
180
Ly/2
160
Cy
Re(Zin) HFSS Re(Zin) Proposed method + FC Im(Zin) HFSS Im(Zin) Proposed method + FC
140
Spice circuit to get Zin from loop impedance x
Loop-x
Loop-y
y 100
The advantage of the discrete complex images method in computing the substrate Green’s function is evident from the convenient analytical expressions for the mutual impedance between conductor currents that have been derived in Section III. The non-linear fit shown in (12) has a substantially different accuracy requirement in three dimensions visa-vis the two dimensional counterpart used in [3]. The 2D Fourier transform (10) needed to recover the co-ordinate space representation, or equivalently the 1D Fourier Bessel representation [11], is plagued with sign oscillations demanding higher accuracy in the non-linear least square fitting algorithm. We employ the Variable Projection (VP) algorithm [14], [15] to determine the complex images (12) in terms of the best fit parameters (αk , βk ) ∈ C. The VP algorithm has the advantage of allowing non-uniform sampling points in λ ∈ (0, ∞) to accurately capture rapid oscillations in χN . Moreover, random starting values are needed only for the non-linear parameters (αk ), since the linear coefficients (βk ) are uniquely determined once a solution to the non-linear parameter least square problem is found. We are interested in complex image representations of the substrate over a wide range of frequencies (20-100 GHz). We start at one end of the frequency spectrum (say 100 GHz) using random values between 0 and 1 as starting guesses for αk . It turns out that χN is a smooth function of ω, permitting us to use the best fit values for αk obtained at one frequency as a good starting guess for the parameters at an adjacent frequency. This strategy of ”continuation” is applied progressively over the entire frequency range of interest. If the desired level of accuracy is not achieved at all frequency points, this process is repeated by choosing as starting values the parameters at the frequency point with minimum residual error. This procedure is summarised in Figure 3.
⎤
jωM1L η1 I ⎥ ⎢ .. ⎥ ⎢ .. ⎦⎣ . ⎦ = ⎣ . RL + jωMLL ηL I
1E-3 Ohm-cm Thick = 1 µm 1E3 Ohm-cm Thick = 1 µm
1E0 Ohm-cm Thick = 500 µm
120 100 80 60 40 20
100
Length scale µm
Metal width/thick = 1.0 µm Ht above substrate = 1.5 µm
Zin (Ohms)
y4
Substrate profile
0
20
30
40
50
60
70
80
90
100
Frequency (GHz)
Fig. 4. Input impedance Zin for a square inductor above a 3-layer substrate, as a function frequency, in comparison with HFSS [4]. Zin is computed using R where the R, L parameters are computed the circuit shown above in HSpice, using our method and capacitances are obtained from FastCap [16]. Interconnect and substrate profiles are shown alongside.
R computed by simulating the circuit shown alongside using HSpice. Note that orthogonal loops have no mutual impedance. Comparison with the full-wave field solver HFSS [4] shows excellent agreement in the dominant imaginary part of Z. Fig. 5 shows the impedance (Z) of an octagonal inductor. In this case the octagon is composed of four loops: one each oriented along x ˆ and yˆ, and two inclined at 45◦ angles (non-Manhattan). While the mutual inductance between orthogonal loops is zero, that between inclined loops is computed as per (40). The mutual inductance is not shown in the circuit diagram to avoid clutter. The imaginary part of the impedance agrees well with HFSS. For both inductors, the maximum error in |Z| is less than 5%, with some discrepancy in the real part of Z beyond 60 GHz due to corner effects in capacitance and the possible small manifestations of full-wave effects. Fig. 6 shows the mutual impedance (|Z|) between two identical square inductors lying adjacent to each other above the same 3-layer substrate. The self impedance of one inductor is shown to compare the relative magnitude of the mutual impedance. Fig. 7 compares the impedance computation time for the two coupled inductors over a 3-layer substrate. The time taken for the impedance computation is shown in Fig. 7. In comparison with the full-wave
C14 C13 C12
Loop 1
C24 C23
Loop 2
C1
C34
Loop 3
C2
Loop 4
C3
C4
Spice circuit to get Zin from loop impedance (mutual inductances included but not shown here) 350 2 Lo op
4
45
250
Z (Ohms)
Loop1
Loop 1
°
70 µm
HFSS Proposed method + FC
300
op Lo
70
µm
Loop 3
R EFERENCES
150 Real(Z)
Lo op 2
4 op Lo
50
Metal width/thick = 1.0 µm Ht above substrate = 1.5 µm
0 20
30
40
50
60
70
80
Frequency (GHz)
90
100
Fig. 5. Input impedance Zin for an octagonal inductor above a 3-layer substrate, as a function frequency, in comparison with HFSS [4]. Zin is computed using R Interconnect profile is shown alongside. the circuit shown above in HSpice. Substrate profile is same as in Fig. 4. 100
VI. C ONCLUSIONS We have shown the computational opportunities offered by suitably generated discrete complex image expressions for the quasi-static magnetic vector potential Green’s function towards analytical impedance computation of general 3D VLSI interconnects, including non-Manhattan ones. We demonstrated the efficiency in the impedance extraction for spiral inductor structures at millimeter wavelengths, previously the realm of electromagnetic field solvers which use computationally expensive numerical methods. This work extends the reach of high-frequency VLSI interconnect impedance computation up to the 100 GHz frequency regime, by incorporating the effect of substrate eddy currents that have a non-negligible impact (up to 20%) in this frequency regime. Validation against a full-wave field solver shows that the dominant imaginary part of impedance is accurate over the relevant frequency range (20-100 GHz) where substrate effects are important, while some deviation is observed in the significantly smaller real part above 60 GHz. In terms of computation cost, our method shows an order of magnitude improvement over the full-wave field solver at each frequency point.
Imag(Z)
200
100 Loop3
field solver, our method gives an order of magnitude improvement in runtime for impedance computation at one frequency point. For our method, the time taken for the two steps - populating the impedance matrix in (42) and solving the linear system - are shown separately. A direct method is used here for solving the linear system. For systemlevel applications involving much larger systems, alternative accelerated solvers can be used, although this is beyond the scope of this paper where the fundamentals of the core algorithm are described. The following computation costs have not been included in the comparison with the full-wave field solver: the runtime for calculating the complex images corresponding to the substrate profile once per technology generation, the runtime for capacitance extraction (frequency independent) once per R to compute Zin for the interconnect geometry, and HSpicesimulation comparison.
[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]
200 180 100
|Z| (Ohms)
140
100
10
[14]
160
120 100 80 60
Proposed method HFSS
40 20
Length scale µm Metal width/thick = 1.0 µm Ht above substrate = 1.5 µm
[15] [16]
Zself
http://www.itrs.net/reports.html. ITRS, 2007. M. Kamon et al., IEEE TMTT, vol. 42(9), pp. 1750–1758, 1994. N. Srivastava et al., IEEE TCAD, vol. 28(7), pp. 1047–1060, 2009. User’s Guide: HFSSv10.0. Ansoft Corporation, 2005. T. Moselhy et al., DATE, 2007, pp. 1–6. X. Hu, Ph.D. dissertation, MIT, 2006. R. Escovar et al., IEEE TCAD, vol. 24(5), pp. 783–793, 2005. D. Fang et al., IEE Proc. H, vol. 135(5), pp. 297–303, 1988. R. Escovar, Ph.D. dissertation, UJF, Grenoble, France, 2006. J. Jackson, Classical Electrodynamics. New York, U.S.A.: Wiley, 1975. W. Chew, Waves and Fields in Inhomogeneous Media. New York: IEEE Press, 1995. M. Aksun, IEEE Trans. on MTT, vol. 44(5), pp. 651–658, 1996. F. Grover, Inductance Calculations: Working Formulas and Tables. New York, U.S.A.: Dover, 1962. G. Golub and V. Pereyra, SIAM J. Num. Anal., vol. 10(2), pp. 413–432, 1973. ——, Inverse Problems, vol. 19, pp. R1–R26, 2003. K. Nabors and J. White, IEEE TCAD, vol. 10(11), pp. 1447–1459, 1991.
0 20
Zmutual 30
40
50
60
70
Frequency (GHz)
80
90
100
Fig. 6. Mutual impedance |Z| between two square inductors above a 3-layer substrate, as a function frequency, in comparison with HFSS [4]. Self impedance is indicated for comparison of magnitude. Interconnect profile is shown alongside. Substrate profile is same as in Fig. 4.
MQS: R/L Time (s) Proposed Method Populating Solving linear Total [Z] system
21
75
97
Fullwave HFSS Time (s)
1230
Fig. 7. Runtime for inductor impedance computation for the structure in Fig. 6, at 80 GHz.
Roberto Suaya
Kaustav Banerjee
Mentor Graphics, Wilsonville, OR, USA navin [email protected]
Mentor Graphics, Grenoble, France roberto [email protected]
University of California, Santa Barbara, CA, USA [email protected]
Abstract—We present an efficient and highly accurate approach to highfrequency impedance extraction for VLSI interconnects and intentional on-chip inductors. The approach is based on a three-dimensional (3D) loop formalism that uses discrete complex images approximations applied to a quasi-magnetostatic treatment of the vector potential, resulting in closed-form expressions for the impedance matrix of current filaments in the presence of a multi-layer substrate. Populating the impedance (Z) matrix for 3D configurations of finite transverse dimensions (including nonManhattan wires and inductors) is computationally inexpensive, and includes substrate eddy current effects that become quantitatively important in the frequency regime beyond 20 GHz which is imminent at the 45 nm technology node onwards. The accuracy, as exemplified by the magnitude of inductor impedance |Z|, is within 5% of a full-wave electromagnetic field solver for frequencies up to 100 GHz, with an order of magnitude lower computation cost. The proposed method represents a core technology for incorporation into system level extraction of analog systems consisting of multiple inductors and nearby interconnects, for CMOS on-chip circuits in the nanometer era.
I. I NTRODUCTION As VLSI technology scales beyond the 45 nm node [1], faster transistors drive signals on interconnects with significant frequency components approaching 100 GHz. At the same time, unprecedented levels of functional integration enable large digital circuits, as well as sensitive analog circuits operating in the mm-wave frequency regime [1], to coexist on the same platform. While the frequency-dependent skin and proximity effects are already dealt with during system-level extraction of global on-chip interconnects in the 10-20 GHz range [2], hitherto unimportant electromagnetic phenomena are manifest beyond 20 GHz. At these frequencies, the magnetic fields generated by transient currents penetrating the conductive Silicon substrate result in eddy currents that have a significant impact on on-chip inductors as well as on timingcritical global interconnects. Accurate high-frequency characterization of inductors, including electromagnetic (EM) interactions with multiple nearby inductors and interconnects as often encountered in wireless applications, dictates that they be considered in toto including the substrate, as opposed to being analyzed as individual components in isolation. While computationally efficient analytical approaches have been proposed for impedance computation of 2D Manhattan-only interconnect structures [3], they are insufficient for on-chip inductors and 3D interconnects. Existing alternatives that are applicable to 3D geometries demand numerical quadratures for multi-dimension integrals such as those described in Section II. Accurate computation of the impedance matrix for multiple inductors and interconnects above a multi-layer substrate overwhelms EM field solvers. S-parameter measurements that are currently used for this purpose at 1-10 GHz severely limit design space exploration, and are quite challenging in the next frequency decade due to difficulties with accurate and reliable deembedding calibration. The demand is for computational techniques that correctly handle the substrate effect, while permitting the analysis of structures with complexity well beyond that of individual components. These
978-3-9810801-6-2/DATE10 © 2010 EDAA
considerations guide the search beyond finite element methods, such as HFSS [4], that can be used only for validation purposes. Full-wave field solvers based on boundary element techniques, such as Fastmaxwell [5], have many of the desired properties but lack the performance required for handling large structures. High-performance full-wave field solvers expand the two dimensional Fourier transform of the substrate Green’s function as a linear combination of complex exponentials, in an approximation called the Discrete Complex Images Method (DCIM) [6]. In FastMaxwell [5], DCIM is used in combination with higher order basis functions in a numerical approach to compute the full-wave substrate Green’s function, and applied to three-dimensional interconnects above a single-layer substrate. Quasi-static solvers that provide much better performance are commonly used at lower frequencies. The classical work of FastHenry [2] is the reference in magneto-quasi-static (MQS) applications, but it only allows for conductors and metallic ground planes, and cannot account for a large substrate with possibly high conductivity and permittivity. Our goal is to develop an analytical quasistatic approach capable of handling complex interconnect structures above a multi-layer substrate. To this end, we compute the impedance matrix of a configuration of wires using the Green’s function for a three-dimensional quasi-static magnetic dipole above a multi-layered substrate, and use DCIM to allow for analytical quadratures. The use of a magnetic dipole (closed current loop) [7] instead of an open current source guarantees the correct long distance behavior of electromagnetic fields. Since our computations are in the near-field regime (distances less than few mms) they are insensitive to surface waves that could make the DCIM inaccurate [8]. Furthermore, with the help of the superposition principle to represent finite three-dimensional interconnect loops as a continuous distribution of point-like magnetic dipoles, we obtain the self and mutual impedance of 3D interconnects above a multi-layered substrate in closed form. Our impedance computation methodology is described in Section III and implementation details in Section IV. In Section V, we compare our results with a full-wave field solver [4], after reintroducing the electric field generated by the displacement currents for simple structures where such comparisons are feasible. Besides the computational efficiency resulting from the closed form expressions derived here, our method seamlessly merges into existing loop impedance approaches based on the filaments decomposition of wires above the substrate, that work well at lower frequencies [9]. Hence, the work in this paper is a significant step towards realizing a system-level impedance extraction methodology covering a wide frequency spectrum from nearly D.C. to over 100 GHz, with reasonable computation cost. II. M AGNETIC DIPOLE VECTOR POTENTIAL G REEN ’ S FUNCTION For a general Manhattan interconnect, the collection of a signal wire and its return paths can be partitioned along their length to form bundles [7]. All wires in a bundle run parallel and have the same length. Each bundle consists of multiple current loops, each formed by a pair of
(x1,y2,z1)
z
Forward current on metal Mn
I ˆy x
y
(x2,y2,z2)
& p
Loop Area
0 where each Gs,R z,y corresponds to a single line current (monopole) on the right hand side of (3). To derive expressions for the Green’s function corresponding to the above magnetic dipole, it is convenient to first consider each current term in isolation. Using the 2D Fourier transform (in x and y), we obtain the monopole Green’s function as:
I ˆy Return current on metal Mn+2 ĭ
(x2,y1,z2)
(x1,y1,z1)
Fig. 1. Current loop formed by two parallel current filaments carrying current in opposite directions constitutes a planar magnetic dipole source.
parallel opposing current filaments. A current loop is the smallest atomic unit in our description of interconnect loop impedance. The physical equivalent of a planar current loop, as shown in Fig. 1, is a magnetic dipole [10]. Since on-chip conductors are confined to discrete metal layers in the x-y plane (neglecting vias which have small dimensions), the current loops they form are planar, although they may have arbitrary orientations (φ). For a planar loop carrying current I, its magnetic dipole moment p is perpendicular to the plane containing the loop, and its magnitude is given by [10]: | p| = I × LoopArea
(1)
Rotational symmetry allows us to choose the y-axis along the length of the wires forming the current loop. Hence, p has horizontal and vertical components: ˆ + pz zˆ = | p|[sin(φ)ˆ x + cos(φ)ˆ z] p = px x
(2)
0 Gs,R z,y (x, y, z, x , y , z ) =
∞ ∞ −∞ −∞
⎛
1 × (2π)2
⎞
ejkx (x−x )+jky (y−y ) ⎠ dkx dky ; ⎝ −μ × ejkz |z−z | + χN ejkz (z+z ) 2jkz kz2 = −kx2 − ky2
(6)
where kx and ky are the spectral variables corresponding to x and y, respectively, and χN can be found using the boundary conditions for the and magnetic field (H) at the substrate interface: magnetic flux (B)
R1 ) · zˆ
R0 − B (B
z=0
R0 − H R1 ) × zˆ
= 0; (H
z=0
= 0;
For a single layer substrate, we get: χ1 =
kz − kz 2 ; kz = −kx2 − ky2 − γ12 kz + kz
lim (ejkx a/2 − e−jkx a/2 )/a = lim (2j sin(kx a/2))/a = jkx
0 2 Gd,R z,y (x, y, z, x , y , z ) = ⎡ ⎤ a δ(x − x + )δ(z − z )− )δ(y − y (3) −2π ⎣ 2 ⎦ ; z, z > 0 lim a a→0 a δ(x − x − )δ(y − y )δ(z − z ) 2 Since there are no sources inside the substrate, the Green’s function in region Ri satisfies:
(2 − γi2 )GRi (x, y, z) = 0; z < 0
(4)
where = jωμ(σi + jω i ) for the i substrate layer . The dipole Green’s function in (3) can be expressed as: th
1
0 Gd,R z,y (x, y, z, x , y , z ) =
⎛ ⎞ a s,R0 −2π ⎝−Gz,y (x, y, z, x − 2 , y , z ) ⎠ lim a 0 a→0 μa + Gs,R ,y ,z ) z,y (x, y, z, x + 2
1 The
a→0
(9)
we get:
In 3D layered media, as in the case of a multi-layered substrate, the ¯ Green’s function for the vector potential is a second rank tensor, G, whose components are Gu,v , where u denotes the direction of the dipole source and v the direction of the resultant magnetic vector potential A. Our objective is to compute these quantities when the location of the source r , as well as that of the observation point r, are in the region R0 above the substrate. Consider two opposite currents I yˆ and −I yˆ in the x-y plane, centered at (x ,y ,z ) and separated by an infinitesimal distance a, lying above a multi-layered substrate. The two currents constitute a zˆ-directed magnetic dipole. The corresponding yˆ-directed vector potential Green’s 0 with I → ∞ and a → 0 such that the dipole moment function Gd,R z,y μ (Ia) = 1, satisfies: 2π
γi2
(8)
For an N -layer substrate, these boundary conditions must be applied at all interfaces between different layers [11]. Substituting (6) in (5), and using: a→0
A. Vector potential Green’s function in integral form
(7)
(5)
substrate layers have finite conductivity, and the permittivity of Si is large. Hence, γi2 cannot be neglected.
0 Gd,R z,y (x, y, z, x , y , z ) =
∞ ∞ −∞ −∞
⎛ kx
−μ × 2(2π)2
ejkx (x−x )+jky (y−y
⎝ kz
⎞
)
× ejkz |z−z | + χN ejkz (z+z
)
⎠ dkx dky
(10)
The above dipole Green’s function expression is the sum of a primary term corresponding to a magnetic dipole in free space, and a secondary term which is nearly identical in form (except for the coefficient χN ) arising due to the substrate: Fnet (r, r ) = Fprimary (r, r ) + χN Fsecondary (r, r )
(11)
Although the expressions for the primary and secondary field terms may vary depending on the source (such as electric or magnetic dipoles, or magnetic monopoles), the coefficient χN always appears unaltered [3], [6], [11]. B. Discrete complex images arising from the substrate Although analytical expressions are well-known for the primary field component in (10), analytical computation of the entire integral expression is hampered by the coefficient χN . Hence, the key to obtaining closed form expressions for the substrate Green’s function lies in finding a suitable approximation to χN that transforms the secondary field expression into an anaytically integrable form. An often favored approximation for χN - the DCIM [6], [12] - uses a linear combination of complex exponentials. The following modified DCIM extends several desirable properties to the approximation [3], including lower fitting time and smoothness across a wide range of frequencies: χN ≈
K k=1
βk e−αk (λ/γ1 ) ,
(12)
magnetic dipole source
magnetic dipole source
I
Indeed the same expression can be derived by using the definition of magnetic dipole moment [10] for a loop of area a lying in the x ˆ-ˆ y plane and carrying current I:
I z’
z=0
z’
Modified DCIM
ȕ 1I
K image dipoles
Multi-layer substrate
substrate
Į interface z’+ Ȗ1 1
z=0
1 p (r ) = 2
ȕ2I Į z’+ ȖK 1
ȕKI
where λ = kx2 + ky2 ∈ (0, ∞), γ1 = jωμ(σ1 + jω 1 ), and αk , βk are found by non-linear least squares fitting. The above approximation can be physically interpreted as the sum of K images of the source dipole representing the effect of the substrate, as shown in Fig. 2, and the net field in (11) takes the form: Fnet (r, r ) = Fprimary (r, r ) +
βk F˜primary (r, r , αk )
(13)
k=1
where F˜primary is the field due to an image of the source dipole and has a functional form similar to the term Fprimary . C. Analytical vector potential Green’s function in free space For a single elementary current in free space, the vector potential s,f ree is related to the corresponding Green’s function Gs,f ree (r, r ) A as:
s,f ree (r) = r )dr A Gs,f ree (r, r )J( (14) r ) is the current density at r (with components Jx along x where J( ˆ and Jy along yˆ, respectively). Similarly, for a dipole in free space, the vector d,f ree is written in terms of the dipole Green’s function as: potential A
d,f ree (r) = A
¯ d,f ree (r, r ) · p G (r )dr
(15)
¯ d,f ree (r, r ) is a where p (r ) is the dipole moment density at r and G second rank tensor. Now, for a single current in free space, the vector potential Green’s function is well-known [10]: 1 μ Gs,f ree (r, r ) = 4π |r − r | (16) |r − r | = (x − x )2 + (y − y )2 + (z − z )2 Noting that the vector potential Green’s function for a magnetic dipole consisting of two opposite currents I yˆ and −I yˆ, located at (x − a/2,y ,z ) and (x + a/2,y ,z ) separated by an infinitesimally small distance a, is given by the gradient between the Green’s functions corresponding to the two isolated currents, as shown in (5), we can use (16) to get an alternative representation of the dipole Green’s function in free space. Thus: ree (r, r ) = Gd,f z,y
⎛
μ 1⎜ ⎜ lim a→0 4π a ⎝
⎞
1 (x − x + a2 )2 + (y − y )2 + (z − z )2 1
−
(x −
x
−
a 2 ) 2
r )d3 r = zˆIa r × J(
(18)
and the Green’s function for a magnetic dipole [10]:
Fig. 2. Representation of the substrate as a series of images of the source magnetic dipole using the modified discrete complex images method [3].
K
+ (y −
y )2
+ (z −
⎟ ⎟ ⎠
ree = Gd,f z,y
μ (r − r ) × p 4π |r − r |3 | p|
Recognizing this simple form of the vector potential Green’s function (19) for on-chip interconnect currents in free space, we avoid the multiple dimension integrals in (10) and instead find analytical expressions for the dipole Green’s function in the presence of a multi-layered substrate in the following Section. III. I MPEDANCE E XTRACTION M ETHODOLOGY FOR 3D I NTERCONNECTS A BOVE A M ULTI - LAYER S UBSTRATE A. Analytical vector potential Green’s function for a 3D magnetic dipole in free space: The quasi-static magnetic vector potential at r due to a magnetic dipole located at a point r in free space is given by [10]: × (r − r ) r) = μ × p A( 4π |r − r |3
(x − x ) μ = 2 4π [(x − x ) + (y − y )2 + (z − z )2 ]3/2 ree ¯ d,f ree created by where Gd,f denotes the y-directed component of G z,y the z-directed quasi-static magnetic dipole.
(20)
¯ d,f ree (r, r )) gives the vector The vector potential Green’s function (G potential (A(r)) in free space, due to a unitary magnetic dipole (| p| = 1) ¯ d,f ree (r, r ), in terms of located at (r ). The relevant expressions for G ree (r, r ) (as explained before), are described below. its components Gd,f u,v Due to our choice of co-ordinates, the y-component of p is zero (2). Hence, ree ree ree = Gd,f = Gd,f =0 (21) Gd,f y,x y,y y,z ree is perpendicular Moreover, the cross product in (20) implies that Gd,f u,v to the source dipole moment p . In other words, the Green’s function due to an x-directed dipole source will have no component in the x-direction. Hence: ree ree ree = Gd,f = Gd,f =0 (22) Gd,f x,x y,y z,z
From the observations above, the vector potential Green’s function (in the x-y plane) at r = (x, y, z), due to a unit magnetic dipole source located at r = (x , y , z ) in free space, can be written as:
¯ d,f ree ¯ G
⎡
0
0 0 0
ree = ⎣ Gd,f x,y d,f ree Gx,z
⎤
ree Gd,f z,x ree ⎦ Gd,f z,y 0
(23)
Since we consider the interaction between wires that lie in the x-y plane ree along (vertical conductors or vias are ignored), the components Gd,f u,z d,f ree - the only non-zero component along zˆ are not relevant. Hence, Gx,z zˆ - is ignored. Substituting p from (2) in (20), and with | p| = 1, we get: ree (x, y, z, x , y , z ) = Gd,f x,y
μ sin(φ)(z − z ) 4π |r − r |3
(24)
ree (x, y, z, x , y , z ) = Gd,f z,y
μ cos(φ)(x − x ) 4π |r − r |3
(25)
(17)
z )2
(19)
ree (x, y, z, x , y , z ) = − Gd,f z,x
where the angle φ is shown in Fig. 1.
μ cos(φ)(y − y ) 4π |r − r |3
(26)
B. Analytical vector potential Green’s function for a 3D magnetic dipole above a multi-layer substrate: For the vector potential Green’s function in the presence of a multilayer substrate, we invoke the modified discrete complex images method [3] explained in Section II, to represent the substrate as a series of images of the source dipole (see Fig. 2). The procedure for computing accurate values for the parameters (αk ,βk ) in (12) is explained in the following section. The vector potential due to the combined effect of K +1 dipoles (source and its K images) is given by the superposition of the vector potential due to each dipole considered in isolation. Hence, the vector potential Green’s function at (x,y,z) due to a source dipole located at (x ,y ,z ) in the presence of a multi-layer substrate is given by: ¯ d,sub = G ¯ d,f ree + G
K
¯ d,img G k
(27)
r) along yˆ in (29), we get: Taking the components of A(
y4 M
d,f ree
C. Impedance for a finite current loop: The Green’s function shown above gives the vector potential due to a unit magnitude dipole source located at a point r = (x , y , z ). The vector potential at r = (x, y, z) due to a finite size current loop is given by the superposition of infinitesimally small point sources in the area occupied by the loop, each occupying an area 12 (dx × dy ) and carrying current I. The co-ordinates of the source (x ,y ,z ) span the rectangular (shaded) area shown in Fig. 1, whose extremities are defined by (x1 ,y1 ,z1 ), (x1 ,y2 ,z1 ), (x2 ,y1 ,z2 ) and (x2 ,y2 ,z2 ). Expressing the z-coordinate of the source in terms of x and φ: zx = z1 + (x − x1 )tan(φ), the vector potential at r = (x, y, z) due to the finite current loop in free space is given by (15), yielding:
ree d,f ree ree ˆAd,f yˆ Ad,f Gx,y (x, y, z) + AGz,y (x, y, z) + x Gz,x (x, y, z)
(29)
ree where each Ad,f Gu,v is given by:
y2 x2 =
ree Gd,f (x, y, z, x , y , zx )Idx dy u,v
d,f ree MG x,y
d,f ree MG z,y
r) · dl/I A(
M = ψ/I =
(31)
l
where dl is the length vector for an infinitesimal element of the victim conductor, r is the position vector for this element, and l is the contour along the conductor length. For a victim conductor oriented along yˆ extending from (x3 ,y3 ,z3 ) to (x3 ,y4 ,z3 ), (31) becomes:
y4 M
d,f ree
d,f ree (x3 , y, z3 ) · dˆ A y /I
= y3
μ = 4π
d,f ree MG tan2 (φ) z,y
μ = 4π
where, s=
3
(32)
+q
2
y =y2 x =x2
cos (φ)s − cos (φ)q
y =y1
x =x1
y =y2 x =x2 y =y1 x =x 1
√ m2 + v 2 − v v | + m2 + v 2 ; ln| √ 2 m2 + v 2 + v
⎛
zk cos(φ)tanh
⎜
q = sin(φ) ⎜ ⎝
+vtan−1
−1
√
t m2 + v 2
v t √ zk cos(φ) m2 + v 2
y=y4 (34) y=y3
y=y4
(35) y=y3
(36)
⎞ ⎟ ⎟ ⎠;
u = x3 − x ; v = y − y ; zk = (z3 − z1 ) − (x3 − x1 )tan(φ); t = usec(φ) + zk sin(φ); m2 = t2 + zk2 cos2 (φ) ¯ d,f ree in (15) must be In the presence of a multi-layer substrate, G ¯ d,sub from (27), which, when combined with (28), yields: replaced by G
y2 x2 d,sub (x, y, z) = A
+
K
¯ d,f ree (x, y, z, x , y , z ) · p G (r )dr x
y1 x1 y2 x2
βk
(37) ¯ d,f ree (x, y, z, x , y , z k ) · p G (r )dr x
y1 x1
= −(z1 + αk /γ1 ) + (x − x1 )tan(φ) is the z-coordinate where of the k image dipole. Since the integrands have the same form as those in free space, the subsequent integrals for the mutual impedance M d,sub in the presence of a multi-layer substrate can also be computed analytically: zxk th
d,f ree d,f ree + MG + M d,sub = MG x,y z,y
For the particular choice of co-ordinates used here (currents along yˆ), the ree term Ad,f Gz,x = 0. As expected, the resultant magnetic vector potential is directed along yˆ, just as the currents in the source loop. The mutual inductance between the source current loop and a victim conductor is given by:
(33)
All the integrals shown above can be evaluated in closed form:
(30)
y1 x1
dy
d,f ree d,f ree + MG = MG x,y z,y
k=1
d,f ree (x, y, z) = A
ree Ad,f Gu,v (x, y, z)
I y3
k=1
¯ d,f ree are given in (23)-(26). The expression for The expressions for G d,img ¯ Gk is easily obtained by relocating the source from (x ,y ,z ) to (x ,y ,−(z + αk /γ1 )), and multiplying by the linear coefficient βk : ¯ d,img (x, y, z, x , y , z ) = β G ¯ d,f ree (x, y, z, x , y , −z − αk ) (28) G k k γ1
=
ree d,f ree Ad,f Gx,y (x3 , y, z3 ) + AGz,y (x3 , y, z3 )
K
d,img d,img βk (MG + MG ) x,y z,y
(38)
k=1 d,img d,img and MG are given by the same expressions as those where MG x,y z,y d,f ree d,f ree for MGx,y and MGz,y , respectively, with the following modification to (36): (39) zk = (z3 + z1 + αk /γ1 ) − (x3 − x1 )tan(φ)
The equations above are sufficient to compute the mutual impedance of Manhattan interconnects by orienting the co-ordinate axes such that relevant conductors are parallel to yˆ. When computing the self inductance of a loop, the victim conductor coincides with the source and the center-to-center distance between the source and destination points (u = x3 − x in (36)) is then replaced by the geometric mean distance (GMD) [13] of a conductor with respect to itself (for rectangular crosssection GMD = elog(w+t)−3/2 [13]). The self impedance is then given by Zself = Rself + jωMself , where Rself = ρL( ws 1ts + ws 1ts ) 1 1 2 2 is the static resistance of the signal line (s1 ) and return path (s2 ) of the loop (ws1 , ws2 are widths and ts1 , ts2 are the thicknesses of the conductors).
Starting freq = fStart; Randomly assign starting values:
Besides Manhattan interconnects which are always inclined at right angles, we are also interested in computing the impedance of inductors, which may comprise conductor segments inclined at arbitrary angle θ. In this case, we choose a co-ordinate system that aligns the source current loop with yˆ. For the victim conductor extending from (x3 , y3 , z3 ) to (x4 , y4 , z3 ), the x-coordinate can be expressed in terms of its y-coordinate as: xy = x3 + (y − y3 )tan(θ). The mutual impedance is then given by: Mθd,f ree
μ μ = cos(φ)T1 + sin(φ)T2 ; 4π 4π
Run VP Get best fit
{O j , F j ( freq )}
Dk , Ek
D k (0,1)
Next freq = current freq +/- fStep; Next starting values = current best fits
START
NO
Reached end frequency? YES
Current freq = freq with min residual; Current best fits = fits at current freq
(40)
NO
Are all residuals < 1E-2 ?
YES
STOP
where,
−sinh−1
T1 =
y−y sxz
y3
y4
tan−1
T2 =
x =x2 y =y2 dy; x =x1
(xy − x )(y − y ) (z3 − zx )sxyz
sxyz =
(41)
y =y1
x =x2
y3
sxz =
Fig. 3. Algorithm flowchart for computing complex images approximation using Variable Projection method.
B. Impedance computation for realistic dimension interconnects
y =y2 dy;
x =x1
y =y1
(xy − x )2 + (z3 − zx )2 ;
(xy − x )2 + (y − y )2 + (z3 − zx )2 ;
xy = x3 + (y − y3 )tan(θ) For conductor loops in the x ˆ − yˆ plane (φ = 0), the integral in (40) R is computed analytically (with a symbolic integration tool Maple), although the resulting expression is too long to be included here. IV. I MPLEMENTATION OF I MPEDANCE E XTRACTION M ETHODOLOGY FOR VLSI I NTERCONNECTS
At frequencies and transverse dimensions where current distribution inside the conductors is non-uniform due to skin and proximity effects, each conductor is discretized into multiple thin filaments with piecewise constant current density for each filament [2]. The expressions for self and mutual impedance shown in Section III can then be used for such filaments. Each ith filament in the conductor carrying the forward current forms a loop with each j th filament in the conductor carrying the return current, resulting in L loops. The following linear system of equations is then solved to compute the net current I flowing through the two-conductor bundle:
⎡
R1 + jωM11 ⎢ .. ⎣ . jωML1
... .. . ...
⎤⎡
A. Modified discrete complex images for multi-layer substrates
⎡
⎤
1 .. ⎥ , . ⎦ 1
ηl = 1, ηl ∈ C
(42)
l
where Rl + jωMll is the self impedance of the lth loop, Mlp is the mutual impedance between the lth and pth loops, and ηl is the fraction of the total current I flowing through the lth loop. Finally, the impedance of the bundle is given by Z = I −1 . V. R ESULTS Fig. 4 shows the impedance of a square inductor configuration (combination of two orthogonal loops, Loop − x and Loop − y, shaded with different colors) lying above a 3-layer substrate. The input impedance is Cxy Loop-x
200
Loop-y
Rx/2 Lx/2
Ry/2
Ry/2
Rx/2 Lx/2
Ly/2
Cx
180
Ly/2
160
Cy
Re(Zin) HFSS Re(Zin) Proposed method + FC Im(Zin) HFSS Im(Zin) Proposed method + FC
140
Spice circuit to get Zin from loop impedance x
Loop-x
Loop-y
y 100
The advantage of the discrete complex images method in computing the substrate Green’s function is evident from the convenient analytical expressions for the mutual impedance between conductor currents that have been derived in Section III. The non-linear fit shown in (12) has a substantially different accuracy requirement in three dimensions visa-vis the two dimensional counterpart used in [3]. The 2D Fourier transform (10) needed to recover the co-ordinate space representation, or equivalently the 1D Fourier Bessel representation [11], is plagued with sign oscillations demanding higher accuracy in the non-linear least square fitting algorithm. We employ the Variable Projection (VP) algorithm [14], [15] to determine the complex images (12) in terms of the best fit parameters (αk , βk ) ∈ C. The VP algorithm has the advantage of allowing non-uniform sampling points in λ ∈ (0, ∞) to accurately capture rapid oscillations in χN . Moreover, random starting values are needed only for the non-linear parameters (αk ), since the linear coefficients (βk ) are uniquely determined once a solution to the non-linear parameter least square problem is found. We are interested in complex image representations of the substrate over a wide range of frequencies (20-100 GHz). We start at one end of the frequency spectrum (say 100 GHz) using random values between 0 and 1 as starting guesses for αk . It turns out that χN is a smooth function of ω, permitting us to use the best fit values for αk obtained at one frequency as a good starting guess for the parameters at an adjacent frequency. This strategy of ”continuation” is applied progressively over the entire frequency range of interest. If the desired level of accuracy is not achieved at all frequency points, this process is repeated by choosing as starting values the parameters at the frequency point with minimum residual error. This procedure is summarised in Figure 3.
⎤
jωM1L η1 I ⎥ ⎢ .. ⎥ ⎢ .. ⎦⎣ . ⎦ = ⎣ . RL + jωMLL ηL I
1E-3 Ohm-cm Thick = 1 µm 1E3 Ohm-cm Thick = 1 µm
1E0 Ohm-cm Thick = 500 µm
120 100 80 60 40 20
100
Length scale µm
Metal width/thick = 1.0 µm Ht above substrate = 1.5 µm
Zin (Ohms)
y4
Substrate profile
0
20
30
40
50
60
70
80
90
100
Frequency (GHz)
Fig. 4. Input impedance Zin for a square inductor above a 3-layer substrate, as a function frequency, in comparison with HFSS [4]. Zin is computed using R where the R, L parameters are computed the circuit shown above in HSpice, using our method and capacitances are obtained from FastCap [16]. Interconnect and substrate profiles are shown alongside.
R computed by simulating the circuit shown alongside using HSpice. Note that orthogonal loops have no mutual impedance. Comparison with the full-wave field solver HFSS [4] shows excellent agreement in the dominant imaginary part of Z. Fig. 5 shows the impedance (Z) of an octagonal inductor. In this case the octagon is composed of four loops: one each oriented along x ˆ and yˆ, and two inclined at 45◦ angles (non-Manhattan). While the mutual inductance between orthogonal loops is zero, that between inclined loops is computed as per (40). The mutual inductance is not shown in the circuit diagram to avoid clutter. The imaginary part of the impedance agrees well with HFSS. For both inductors, the maximum error in |Z| is less than 5%, with some discrepancy in the real part of Z beyond 60 GHz due to corner effects in capacitance and the possible small manifestations of full-wave effects. Fig. 6 shows the mutual impedance (|Z|) between two identical square inductors lying adjacent to each other above the same 3-layer substrate. The self impedance of one inductor is shown to compare the relative magnitude of the mutual impedance. Fig. 7 compares the impedance computation time for the two coupled inductors over a 3-layer substrate. The time taken for the impedance computation is shown in Fig. 7. In comparison with the full-wave
C14 C13 C12
Loop 1
C24 C23
Loop 2
C1
C34
Loop 3
C2
Loop 4
C3
C4
Spice circuit to get Zin from loop impedance (mutual inductances included but not shown here) 350 2 Lo op
4
45
250
Z (Ohms)
Loop1
Loop 1
°
70 µm
HFSS Proposed method + FC
300
op Lo
70
µm
Loop 3
R EFERENCES
150 Real(Z)
Lo op 2
4 op Lo
50
Metal width/thick = 1.0 µm Ht above substrate = 1.5 µm
0 20
30
40
50
60
70
80
Frequency (GHz)
90
100
Fig. 5. Input impedance Zin for an octagonal inductor above a 3-layer substrate, as a function frequency, in comparison with HFSS [4]. Zin is computed using R Interconnect profile is shown alongside. the circuit shown above in HSpice. Substrate profile is same as in Fig. 4. 100
VI. C ONCLUSIONS We have shown the computational opportunities offered by suitably generated discrete complex image expressions for the quasi-static magnetic vector potential Green’s function towards analytical impedance computation of general 3D VLSI interconnects, including non-Manhattan ones. We demonstrated the efficiency in the impedance extraction for spiral inductor structures at millimeter wavelengths, previously the realm of electromagnetic field solvers which use computationally expensive numerical methods. This work extends the reach of high-frequency VLSI interconnect impedance computation up to the 100 GHz frequency regime, by incorporating the effect of substrate eddy currents that have a non-negligible impact (up to 20%) in this frequency regime. Validation against a full-wave field solver shows that the dominant imaginary part of impedance is accurate over the relevant frequency range (20-100 GHz) where substrate effects are important, while some deviation is observed in the significantly smaller real part above 60 GHz. In terms of computation cost, our method shows an order of magnitude improvement over the full-wave field solver at each frequency point.
Imag(Z)
200
100 Loop3
field solver, our method gives an order of magnitude improvement in runtime for impedance computation at one frequency point. For our method, the time taken for the two steps - populating the impedance matrix in (42) and solving the linear system - are shown separately. A direct method is used here for solving the linear system. For systemlevel applications involving much larger systems, alternative accelerated solvers can be used, although this is beyond the scope of this paper where the fundamentals of the core algorithm are described. The following computation costs have not been included in the comparison with the full-wave field solver: the runtime for calculating the complex images corresponding to the substrate profile once per technology generation, the runtime for capacitance extraction (frequency independent) once per R to compute Zin for the interconnect geometry, and HSpicesimulation comparison.
[1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]
200 180 100
|Z| (Ohms)
140
100
10
[14]
160
120 100 80 60
Proposed method HFSS
40 20
Length scale µm Metal width/thick = 1.0 µm Ht above substrate = 1.5 µm
[15] [16]
Zself
http://www.itrs.net/reports.html. ITRS, 2007. M. Kamon et al., IEEE TMTT, vol. 42(9), pp. 1750–1758, 1994. N. Srivastava et al., IEEE TCAD, vol. 28(7), pp. 1047–1060, 2009. User’s Guide: HFSSv10.0. Ansoft Corporation, 2005. T. Moselhy et al., DATE, 2007, pp. 1–6. X. Hu, Ph.D. dissertation, MIT, 2006. R. Escovar et al., IEEE TCAD, vol. 24(5), pp. 783–793, 2005. D. Fang et al., IEE Proc. H, vol. 135(5), pp. 297–303, 1988. R. Escovar, Ph.D. dissertation, UJF, Grenoble, France, 2006. J. Jackson, Classical Electrodynamics. New York, U.S.A.: Wiley, 1975. W. Chew, Waves and Fields in Inhomogeneous Media. New York: IEEE Press, 1995. M. Aksun, IEEE Trans. on MTT, vol. 44(5), pp. 651–658, 1996. F. Grover, Inductance Calculations: Working Formulas and Tables. New York, U.S.A.: Dover, 1962. G. Golub and V. Pereyra, SIAM J. Num. Anal., vol. 10(2), pp. 413–432, 1973. ——, Inverse Problems, vol. 19, pp. R1–R26, 2003. K. Nabors and J. White, IEEE TCAD, vol. 10(11), pp. 1447–1459, 1991.
0 20
Zmutual 30
40
50
60
70
Frequency (GHz)
80
90
100
Fig. 6. Mutual impedance |Z| between two square inductors above a 3-layer substrate, as a function frequency, in comparison with HFSS [4]. Self impedance is indicated for comparison of magnitude. Interconnect profile is shown alongside. Substrate profile is same as in Fig. 4.
MQS: R/L Time (s) Proposed Method Populating Solving linear Total [Z] system
21
75
97
Fullwave HFSS Time (s)
1230
Fig. 7. Runtime for inductor impedance computation for the structure in Fig. 6, at 80 GHz.
