Fast Transient Analysis of Power Distribution
Network Modeled by Unstructured Meshes byUsing Locally Implicit Latency Insertion Method
Shingo Okada #1 , Hiroki Kurobe # , Tadatoshi Sekine ∗2 , Hideki Asai #∗3#

Graduate School of Engineering, Dept. of
Systems Eng., Shizuoka University
3-5-1 Johoku, Naka-ku, Hamamatsu-shi, 432-8561, Japan1

okada@tzasai7.sys.eng.shizuoka.ac.jp∗

Dept. of Systems Eng., Shizuoka University
3-5-1 Johoku, Naka-ku, Hamamatsu-shi, 432-8561, Japan2

sekine@tzasai7.sys.eng.shizuoka.ac.jp
3hideasai@sys.eng.shizuoka.ac.jp

Abstract—This paper proposes a locally implicit latency insertion method (LILIM), which is a suitable method for the fast
simulation of an arbitrary shaped power distribution network(PDN) modeled by triangular meshes. First, an efﬁcient modelingmethod based on triangular mesh is reviewed and we refer tothe limitation of the LIM for the meshed PDN analysis. Next,
in order to overcome the problem, we formulate the LILIM bycombining the efﬁcient modeling and the locally implicit schemes.
Finally, the numerical results show that the LILIM is applicableand efﬁcient for the simulation of the PDN analysis.
Index Terms—Delaunay triangular mesh, latency insertion
method, locally implicit ﬁnite difference method, power distribution network, time-domain analysis, Voronoi tessellation.

i j2

i jn j

i j1

B

R jk

L jk

E jk
-

Gj

Cj H j
(a)

vj

i jk

+

vk
(b)

I. I NTRODUCTION

Fig. 1. Network topologies required for LIM. The arrows correspond todirections of the branch currents. (a) Node topology. (b) Branch topology.

Recent semiconductor packaging technologies have enabled
integrated circuits (ICs) to be assembled within an extremely
small area. To provide power to each system, an irregularshaped power distribution networks are arranged in chips,
packages, and boards. However, the simultaneous switching
noise (SSN) induced by increased power density and shorter
signal edges directly degrades the power/signal integrity(PI/SI) of the recent low-voltage PDN. In order to guarantee
PI/SI of the system composed of the chips, packages, andboards, the PDN must be modeled and simulated efﬁciently
and accurately. One of the useful approaches for modelingand simulation is the ﬁnite different method using square
meshes [1]. This method can simulate efﬁciently the PDN withapertures. However, modeling of the PDN using square meshes
suffers from the staircase approximation error when the PDN
contains small features such as via holes, small apertures,
plane gaps and so on. Another approch for the PDN modelingis one using triangular meshes instead of the square ones [2],[3]. This type of modeling method usually uses the Delaunary
triangulation, which genarates the Delaunay mesh and Voronoi
tessellation, which are mutually orthogonal to each other. This
property of the mesh geometry is important in the sense thatelectric and magnetic ﬁelds are orthogonal to each other. Thesemodeling methods achieve successfully to reduce the total
number of unknown variables without degrading the accuracy
for the modeling of the ﬁne geometric structure compared withthat of the square meshes.
On the other hand, the latency insertion method (LIM) has
been proposed as one of the fast circuit simulation techniques[4]. The LIM is suitable to simulate an equivalent circuit
extracted by the modeling methods with triangular meshes
since the equivalent circuit satisﬁes the required topologies,
where each node and branch contain the grounded capacitance
C and the serial inductance L as shown in Fig. 1, respectively.
However, the LIM has the constraint of time step size to satisfy
the numerical stability condition similar to the FDTD method
since the method is based on the explicit leapfrog scheme [4],[5], [6]. Concretely, the maximum time step size is limited by
the minimum values of C and L in the computational domain.
The values of C and L in the equivalent circuit are proportional
to the size of each triangular mesh. Thus, the maximum time
step size of the LIM is extremely dependent on these values

b=4

b=5
da5
da6

la5

la4

da1

la1

La6

b=2

(a)

La4
b=3
a

20

b=6
La1

da2

b=1

A

10

35

La5

da4

Aa la3 b=3
a
da3
la2

b=6 la6

b=4

b=5

b=1

La3
La2

x
C (Unit: mm)

y

z

50

10 20
50

Ca

b=2

25

B

(b)

Fig. 2. Voronoi diagram and triangular mesh for conductor media. (a) Double
circle, dash lines, and solid lines indicate the Voronoi point, the Voronoi
region, and triangular mesh. (b) Equivalent circuit, where grounded resistance
G and series resistance R are not shown, of triangular mesh, where vertices
and edges of triangular mesh are corresponding to nodes and branches.

100
Fig. 3.

Top view of an example PDN modeled by triangular meshes.

0

500

0

since the vicinities of ﬁne structures are modeled by exploiting
the tiny size of triangular meshes. As a consequence, a large
amount of computational cost is required even if the LIM,
which is much faster than matrix-based SPICE-like simulators,
is employed.
In this paper, we propose the locally implicit (LI) LIM
to overcome the inherent problem of the LIM. The LILIM
employs a locally implicit scheme instead of the globally
explicit leapfrog scheme used in the original LIM The similar
approaches have been applied to the FDTD method [7], [8].
These approaches apply the implicit scheme to one direction inwhich the space is discretized into smaller cells than those in
the other directions. By contrast, in our proposed LILIM, theimplicit scheme is applied to the local area which is composedof small values of reactance elements.
The remainder part of the paper is organized as follows.
In Section II, the modeling method by using the triangularmeshes is reviewed. Section III describes the formulations
and features of the proposed LILIM. Section IV demonstrates
accuracy and efﬁciency of the proposed method by showing
some numerical examples, and conclusions are given in Section V.
II. PDN M ODELING BY U SING T RIANGULAR M ESHESIn modeling of the arbitrary shaped PDN, square meshes
must be globally tiny even if the ﬁne structure is contained in
the local region. This fact leads that the number of unknown
variables increases drastically. Additionally, the square meshscheme suffers from the staircase approximation on the edgeof the geometric structure. In order to cope with the above
inherent problems, the efﬁcient modeling methods based on
triangular meshes have been proposed [2], [3]. These methods
successfully reduce the total number of unknowns without
degrading the accuracy since the triangular mesh scheme is
ﬂexible to capture the complex structures by resizing the meshsizes locally. These triangular meshes are lead by Delaunay
Triangularization and Voronoi Tessellation as shown in Fig.
2(a). In Fig. 2(a), the gray hexagonal part indicates the area
Aa of Voronoi region surrounding Voronoi point a, and dab

1000

1500

Block circuit 1
Block circuit 2
Block circuit 4

500

Block circuit 5
Block circuit 81000 Block circuit 3
Block circuit 6

1500

Block circuit 7

Block circuit 9

0

1000

Fig. 4.

Coefﬁcient matrix structure of {v} in the LILIM.

and Iab indicate the length of the Voronoi edge and the lengthof the edge in triangular mesh, respectively. The equivalentcircuit is extracted from the object meshed triangularly, as
shown in Fig. 2(b), where the Voronoi point and the side oftriangular mesh correspond to the node and branch topologies.
And the circuit is composed of the grounded capacitances andinductances of which the values are derived by the following
manner:Aa
dab
Ca = ε , Lab = µh
,
(1)
h
lab
where ε, µ, and h indicate permittivity, permeability, anddielectric thickness. It is conﬁrmed that the equivalent circuit
extracted from the triangular meshes has the topologies, inwhich every node and branch contains the grounded capacitance C and the serial inductance L. Therefore, the equivalentcircuit is suitable to be solved by using the basic LIM. For
example (1) means that the value of Ca is proportional to
Aa , and the arbitrary shaped PDN model is constructed by
the various values of capacitance according to the mesh sizes.
Thus, the maximum time step size of the basic LIM is strictly

TABLE I
CPU TIME

1.01
LIM
LILIM (Cth = 0.5 pF)
Voltage[V]MethodLIM
LILIM

CPU Time [s]59.71916.406

Speed-up
1.00
3.64

TABLE IIMEMORY UTILIZATION

Method
LIM
LILIM

1

0

2
Time[ns]Fig. 5.

blocks. n denotes the time index, and ∆t is the time step size.
In (4), the node voltage and the branch currents in branch
block are arranged at the same time step. Thus, to update theunknown variables, simultaneous equations must be solved.
n+ 1
iam 2 given by the implicit difference of (3), namely (5),)
Lam n− 21∆t ( n+ 12
n+ 1
n+ 1
iam 2 =
iam +va
− vm 2 ,
(5)
βam
βam

4

Waveform results at the port C.

limited by the reactance element values as mentioned in [5].
III. LILIM FOR T HE M ESHED PDN
In order to simulate efﬁciently the arbitrary shaped PDN
modeled by triangular meshes, we propose and formulate
the LILIM by employing locally implicit scheme insteadof globally explicit leapfrog scheme. First, the Kirchhoff’s
current law (KCL) and Kirchhoff’s voltage law (KVL) areapplied to the node a and the branch ab in Fig. 2, and thedifferential equations are given by
a

dva
Ca
+ Ga va
dt

= −

NB
∑

iab ,

Memory utilization[ M byte]
43
145

(2)

b=1

diab
(3)
+ Rab iab = va − vb ,dt
where va , iab , and NBa indicate the voltage at the node a, thecurrent ﬂowing from the node a to the node b, and the total
number of branches connecting to the node a. Next, the termof the branch currents connecting to the node a in (2) is split
into two terms, which are accompanied with the capacitancelarger than the threshold capacitance Cth and less than Cth . In
this paper, the former can be regarded as branch topology, on
the other hand, the latter is deﬁned as branch block. Next, byapplying the leapfrog scheme [4], [5] and the implicit scheme
to the branch topology and the branch block respectively in(2), then, the following equation is derived;
( N
)
NS
L
∑
∑
Ca n− 12
∆t
n+ 21
n+ 21
n
va
−
iam −
iam
,
(4)
=
+
va
αa
αa
m=1
m=1
Lab

where αa = (Ca + ∆tGa ), NL + NS = NBa , and NL and
NS are the total numbers of the branch topologies and branch

is substituted into the terms of branch blocks in (4), whereβam = (Lam + ∆tRam ). Then, (4) is transformed to
(
)
)
NS
NS (
∆t ∑
∆t ∑
∆t
∆t n+ 12
n+ 12
1+
va
−
vm
αa m=1 βam
αa m=1 βam
[ N
)]
NS (
L
∑
∑
∆t
Ca n− 12
Lam n− 12
n
va
+
−
iam
. (6)
=
iam −
αa
αa
βam
m=1
m=1
In the left hand side of (6), there are plural unknown node
voltages at the (n + 12 )-th step. The node accompanied withn+ 1

the voltage variable vm 2 is connected to the node a throughthe branch block.
This is caused by applying the implicit scheme to thecurrents in the branch blocks. Thus, in our proposed LILIM,
we deﬁne the circuit in which the nodes are coupled with eachother as a block circuit. By deriving the node-based equations
in a block circuit, then, the updating formula of the block
circuit is given by
1(([I] + ∆t[K1 ]) + [K2 ][M ]) {v}n+ 21= [K3 ]{v}n− 2 + {b},(7)where {v}, {b}, [Kj=1,2,3 ], [M ], and [I] indicate the unknown
voltage vector in the block circuit, the known vector, thecoefﬁcient matrices, the incidence matrix, and the unit matrix,
respectively. The dimensions of the vectors and the matricesin (7) are deﬁned as{v}, {b} ∈ <Nblock ,[Kj ], [M ], [I] ∈ <Nblock ×Nblock ,where Nblock denotes the number of the nodes in one block
circuit. Equation (7) is solved by using a direct matrix solver,e.g. LU decomposition method. The computational cost is
relatively smaller than the cost of the SPICE-like algorithm
although the calculation cost increases as the size of the blockcircuit becomes larger, because the direct method is used only
for the local block circuits similarly to the block-LIM [9].
In the proposed LILIM, the voltage variables in each block
circuit are updated by (7). After that, the current variables inthe branch block are updated explicitly by using (5). On theother hand, the unknown variables in the domain except forthe block circuits are updated explicitly based on the following
basic LIM updating formulas
Nan+ 1
va 2 =

B
∑
Ca
∆t
n− 1
va 2 −
inab ,Ca + ∆tGa
Ca + ∆tGa

(8)

b=1

Lab − ∆tRab n∆t n+ 12
n+ 1
iab +
(va
− vb 2 ).
(9)
Lab
Lab
The LILIM has the advantage similar to that of the HIEFDTD method [7], [8], [10], which is locally free from the
CFL condition. The maximum time step size of the LILIM
does not depend on the value of the reactance elements inthe block circuit. Thus, the LILIM can perform the efﬁcient
circuit simulation of the the arbitrary shaped network withvarious latencies.
in+1
ab =

IV. N UMERICAL R ESULTS
In order to estimate the efﬁciency of the LILIM, an example
PDN has been simulated by the basic LIM and LILIM. The
shape and dimension of the PDN are shown in Fig. 3. The
height of the dielectrics between the conductor planes is h =0.6 mm, and its relative permittivity εr is 4.2. The conductor
plane is discretized using the triangular meshes with 1890
nodes. The value of Cth is derived empirically and equals
to 0.05 pF. In this case, 1890 nodes are classiﬁed into 165
node topologies and 9 block circuits, of which the numbersof the internal nodes are 1, 3, 16, 840, 1, 1, 860, 2, and1, respectively. Fig. 4 shows the coefﬁcient matrix of theunknown voltages generated by (7) and (8). In Fig. 4, eachcoefﬁcient matrix of the block circuit is factored using the
LU decomposition before the transient analysis. The maximum
time step sizes used in the LIM and LILIM are 0.0267 ps and3.15 ps, respectively. These values are evaluated by assuming
the minimum value of capacitance and the value of Cth are lessthan the upper bound mentioned in [5]. The constant voltage
source (= 1.0 V) is connected to the port A in Fig. 3. In
addition, the input current source is connected to the port B
in Fig. 3, and is the triangle wave of which the pulse value is0.005 A, delay time is 0.1 ns, rise and fall times are 0.1 ns
and 0.3 ns, and period is 0.5 ns. The simulation interval is [0ns, 5 ns], and we observed the voltage waveform at the port
C in Fig. 3.
The waveform results are illustrated in Fig. 5. It is conﬁrmedthat the waveform obtained from the LIM and LILIM agree
well with each other. CPU times and memmory utilization of
the LIM and LILIM are shown in Table I and II, respectively.
From Table I, we can seee that the LILIM is about 3.64
times faster than the LIM. As a result, the proposed method ismore effective than the conventional methods for the arbitrary
shaped PDN analysis with appropriate accracy.

V. C ONCLUSION
In this paper, an efﬁcient LIM-based algorithm for themeshed PDN analysis based on locally implicit scheme has
been proposed. The the time step size limitation in the LIM
was alleviated by exploiting the advantages of the locally
implicit-based algorithm. As a consequence, it has been veriﬁed that the LILIM can perform the fast circuit simulations
using the larger time step sizes which are prohibited in thebasic LIM analysis. The numerical results showed that the
LILIM was about 3.64 times faster than the basic LIM withoutloosing accuracy for the simulation of the the arbitrary shaped
PDN.
ACKNOWLEDGMENT
This work has been partially supported by The Semiconductor Technology Academic Research Center (STARC).
R EFERENCES[1] A. Engin, K. Bharath, and M. Swaminathan, “Multilayered ﬁnitedifference method (MFDM) for modeling of package and printed circuitboard planes,” Electromagnetic Compatibility, IEEE Transactions on,
vol. 49, no. 2, pp. 441 –447, may 2007.
[2] K.-B. Wu, G.-H. Shiue, W.-D. Guo, C.-M. Lin, and R.-B. Wu,“Delaunay–voronoi modeling of power–ground planes with source port
correction,” IEEE Trans. Adv. Packag., vol. 31, no. 2, pp. 303–310, May
2008.
[3] J. Y. Choi and M. Swaminathan, “An effective modeling method for
multi–scale and multilayered power/ground plane structures,” in Proc.
IEEE ECTC 2011, May 2011, pp. 477–483.
[4] J. E. Schutt-Ainé, “Latency insertion method (LIM) for the fast transient
simulation of large networks,” IEEE Trans. Circuits Syst. I, vol. 48, pp.
81–89, Jan. 2001.
[5] S. N. Lalgudi and M. Swaminathan, “Analytical stability condition of
the latency insertion method for nonuniform GLC circuits,” IEEE Trans.
Circuits Syst. II, vol. 55, no. 9, pp. 937–941, Sep. 2008.
[6] K. S. Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” IEEE Trans. Antennas
Propag., vol. AP-14, pp. 302–307, May 1966.
[7] J. Chen and J. Wang, “A three-dimensional semi–implicit FDTD scheme
for calculation of shielding effectiveness of enclosure with thin slots,”
IEEE Trans. Electromagn. Compat., vol. 49, no. 2, pp. 354–360, May
2007.
[8] M. Unno and H. Asai, “HIE-FDTD method for hybrid system withlumped elements and conductive media,” IEEE Microw. Wireless Compon. Lett., vol. 21, no. 9, pp. 453–455, Sep. 2011.
[9] T. Sekine and H. Asai, “Block latency insertion method (block-LIM)
for fast transient simulation of tightly coupled transmission lines,” IEEE
Trans. Electromagn. Compat., vol. 53, pp. 193–201, Feb. 2011.
[10] M. Unno, S. Aono, and H. Asai, “GPU-based massively parallel 3D HIE-FDTD method for high-speed electromagnetic ﬁeld simulation,”
IEEE Trans. Electromagn. Compat., vol. 54, no. 4, pp. 912 –921, Aug.
2012.

Power Distribution Network Modeling for 3-D ICs with TSV Arrays
Chi-Kai Shen

Yi-Chang Lu

Yih-Peng Chiou

Tai-Yu Cheng

Tzong-Lin Wu

Graduate Institute ofCommunication Eng.,National Taiwan
University,
Taipei 10617, Taiwan

Graduate Institute ofElectronics Eng.,National Taiwan
University,
Taipei 10617,
TaiwanGraduate Institute of
Photonics &Optoelectronics,
National Taiwan
University
Taipei 10617, Taiwan

Graduate Institute ofCommunication Eng.,National Taiwan
University,
Taipei 10617, Taiwan

Graduate Institute ofCommunication Eng.,National Taiwan
University,
Taipei 10617, Taiwan
wtl@cc.ee.ntu.edu.tw

Three-dimensional integrated circuit (3-D IC) provides a
possible way to follow the prediction of Moore’s Law andachieve design goals, which includes high transmission rate,
low power consumption, compact size, and high integration ofheterogeneous functionality. In power distribution networks(PDNs) of 3-D IC, simultaneous switching noise (SSN) willcause even more serious problems in power integrity (PI),signal integrity (SI) and electromagnetic interference (EMI)because noise coupling could occur in any direction, not juston one specific plane. To deal with these noise issues,
full-wave simulation is one of the methods being adopted toaccurately analyze PDNs in 3-D ICs. However, because timeand memory consumption are concerned, equivalent circuit
models are often used to approximate the results of full-wave
simulation instead.
As shown in Fig. 1(a), PDN modeling in 3-D IC can be
categorized into meshed-type PDN and via-type PDN
according to the shape. Meshed-type PDN comprises
power/ground grids on the chip and on the interposer, while
via-type PDN includes TSVs and bumps. For the mesh-type
PDN, as shown Fig. 1(b), the conventional equivalent circuits
for power/ground grids are series resistance and inductancealong a wire, and mutual capacitance as well as mutual
inductance between two wires. For the equivalent circuit of
via-type PDN, as shown in Fig. 1(c), the arrays of TSVs andbumps can be modeled similarly except for the conductivity Gof silicon substrate.
In previous researches, constructing equivalent circuits ofgrids and TSVs separately and then combining them together
is the typical modeling approach for 3-D IC PDNs. However,it assumes that there is no coupling between the power/groundgrids and TSVs. In fact, the silicon substrate with low
conductivity at low frequency has electrical coupling to themetals atop. This phenomenon becomes significant if the areaof the metal is very large, such as high-density power/groundgrids.
This coupling phenomenon can be demonstrated with asimple simulation, as shown in Fig. 2(a). A 2-D cross-section
structure with two conductors in ANSYS Q2D environment,
which is used to compare with the structure with two
conductors located above the silicon substrate, as shown in
Fig. 2(b). Since the silicon acts as a conductor at low
frequency range, we can substitute it with a PEC material. The
simulation results are shown in Table I. The numbers revealthat there is strong coupling between the metal and the siliconsurface, but the coupling between the metals is much weaker
than the case with no silicon surface Therefore, once we

(a)

(b)

(c)

Fig. 1. Structure and equivalent circuit model of PDN in 3-D IC.
(a) PDN in a 3-D IC (b) Equivalent circuit of mesh-type PDN.
(c) Equivalent circuit of via-type PDN.

combine the models of power/ground grids and the TSVs
together, the coupling between the metals in power/groundgrids can be neglected, while only the coupling capacitances
Cm_sub between the metals and the silicon substrate with TSVs
inside need to be considered. In our case, these capacitance
Cm_sub can be calculated approximately by the formula ofmicrostrip line capacitance by treating the silicon surface as
the ground.
The coupling node insertion method (CNIM) for array-type
TSVs is proposed here to create the additional nodes forconnecting the coupling capacitance Cm_sub between the silicon
substrate and the power/ground grids. In this method, we

(a)

(b)
≈
Fig. 4. Test structure for integration of mesh-type and via-type
PDN and its port position.
Fig. 2. Simulation for testing the effect of the silicon surface.
(a) Without silicon surface. (b) With silicon surface.
TABLE I
Table for summarizing the simulation results in
Fig. 2(a) and Fig. 2(b)Without Silicon Surface
With Silicon Surface(pF/m, Fig. 4(a))(pF/m, Fig. 4(b))

Cm_rod59.412

11.496

Cm_subNone

243.29(a)	
  Fig. 3. The 3-D scheme for deciding where the capacitance islocated between the metal and the silicon surface.

partition the Y elements of TSV model from the typical one totwo series elements with each one composed of the parallel
2CSi and 2GSi, as shown in Fig. 3. The node connected to theelements Cox and CSi is regarded as the silicon substrate node.
On the other hand, the nodes between the Rrod and Ls,rod on the
power/ground grids is regarded as the metal nodes. Therefore,
the two ends of the Cm,sub elements can be connected to thesenodes of the nearby TSV and power/ground grids respectively.
To demonstrate the full-PDN model and the effect of CNIM,
we construct a test structure of a passive silicon interposerwith metals in the redistribution layer (RDL) and connected toTSVs. As shown in Fig. 4, the test structure is composed of
16-by-16 power/ground grids on the RDL stacked on a 6-by-6
TSV array. In Fig. 4, the ports are assigned with one (Port 1)
located on the power/ground grids and three (Port 2 to 4)
located at the opposite end of the TSVs.
Two of numerical results for self-impedance (Z11) andtransfer impedance (Z24) are shown in Fig. 5. The results fromthe full-wave simulation and the equivalent models
with/without the coupling capacitances are shown in eachsub-figure. The results from the equivalent model are in good(b)

Fig. 5. Numerical results of self and transfer impedances of the test
structure in Fig. 4. (a) |Z11|. (b) |Z 24|.

consistency with full-wave simulation from DC to 40 GHz.
The results shown in Fig. 7 demonstrate that considering thecoupling effect using our CNIM model will greatly improve
the accuracy of simulation results. Using our model, thecalculated resonance frequency and magnitude of impedance
are both closer to those obtained from full-wave simulator.
Therefore, the proposed CNIM can offer an effective way
to revise the equivalent circuit model, obtain more accurate
results in simulating PDN in 3-D IC and also lower the timeand memory consumption compared with full-wave
simulation.

An Explicit and Unconditionally Stable Finite
Difference Scheme for the Fast Transient Analysis ofa Power Distribution Network
Norio Nishizaki#1, Tadatoshi Sekine*2, Hideki Asai*#3#

Graduate School of Engineering, Dept. of
Systems Eng., Shizuoka University
3-5-1 Johoku, Naka-ku, Hamamatsu-shi, 432-8561, Japan1nishizaki@tzasai7.sys.eng.shizuoka.ac.jp*

Dept. of Systems Eng., Shizuoka University
3-5-1 Johoku, Naka-ku, Hamamatsu-shi, 432-8561, Japan2

sekine@tzasai7.sys.eng.shizuoka.ac.jp
3hideasai@sys.eng.shizuoka.ac.jp

Abstract— Generally, an explicit difference scheme has thenumerical stability condition which constrains the time step size,
and this condition depends on the lowest reactance component. If
the time step size that does not fulfill the numerical stability
condition is used, unstable modes arise, and we fail the transient
analysis. In this paper, we propose an explicit andunconditionally stable finite difference scheme in which theconstraint of time step size is overcome by removing the unstablemodes. Numerical results show that the proposed method is quite
effective to perform the fast transient simulation of power
distribution networks with extremely small apertures.

I. INTRODUCTION
Significant advances of the semiconductor packaging
technology have provided a complicated design of high-speed
and high-density electronic circuits. Therefore, various effects
depending on the high-frequency characteristics of signals
cause unexpected behaviors of chips and packages on a
printed circuit board. In addition, as supply voltages for thesecomponents are becoming low, the simultaneous switching
noise (SSN) seriously affects the power integrity of a power
distribution network (PDN). Therefore, it becomes importantto verify the PDN properly in order to ensure the power
integrity of the advanced packaging system.
In order to simulate the SSN by performing the circuit
analysis, equivalent circuit models of the PDN are usuallyprovided by commercially-based extractors or hand-madeextraction techniques [1]. In general, if there exist fine objects
or apertures in the PDN, a large number or small meshes are
required for capturing the exact shapes of the fine structures.
As a result, these fine meshes induce a large network andextremely small reactance elements.
There exist several explicit finite difference schemes in thecircuit [2, 3] and electromagnetic simulations [4, 5] while they
are basically based on primary finite difference methods [6].
These explicit schemes usually have low computational costs
and can perform the fast simulations of large problems.
Especially, the latency insertion method (LIM) and extended

algorithms of it get considerable attention in the field of thefast transient analysis of a large-scale network [2, 3]. However,
the explicit scheme has the numerical stability condition
which strictly limits the maximum time step size used in thetransient analysis. In order to alleviate such a numericalstability problem, the recent researches on the finite difference
scheme pay much attention to unconditionally stable schemes[7, 8, 9]. Although the existing unconditionally stable schemes
are more efficient than the conventional explicit method, theschemes partially include time-consuming matrix calculations.
In this paper, we propose an fully-explicit andunconditionally stable difference scheme. In the proposed
method, the root cause of the instability of the explicit method
is specified by using the eigenvalue decomposition. Then, theinstability is eliminated by performing simple algebraic
operations.
The rest of the paper is organized as follows. Section II
discribes the explicit finite difference scheme in the circuit
simulation and its numerical stability condition followed by
the proposed explicit and unconditionally stable scheme. The
accuracy and CPU time of the proposed method are verified
by performing the transient simulations of the PDN withextremely-small apertures in Section III. We conclude thepaper in Section IV.
II. PROPOSED METHOD
We assume that the circuit to be analysed is composed of anumber of node and branch topologies shown in Fig. 1. The
node topology is composed of the parallel-connected capacitor
and current source, and the branch topology consists of theseries-connected inductor and voltage source. By using
Kirchhoff's current and voltage lows, the circuit equation
associated with the network can be written in the vectormatrix form as followsd
C v + Mi = h ,
(1)
dt

d
i − MT v = e ,(2)dtwhere C and L are the capacitance and inductance matrices, v
and i are the voltage and current vectors, h and e are thecurrent source and voltage source vectors, and M is theincidence matrix. Then, differentiating (1) with respect to timeand substituting (2) into (1) lead tod2dC 2 v + ML−1M T v = h − ML−1e .
(3)
dt
dt
Equation (3) can be regarded as the second order ordinary differential
equation and rewritten asΤd2
u + Su = j ,
dt 2

I i2

I i1

L

I ik

Vi

I i3

L ij
＋

Hi

Ci

Vi

－

Vj

I ij

(a)
(b)
Fig. 1. The network topologies constituting the lossless network. The arrows
correspond to the directions of the branch currents. (a) Node topology. (b)Branch topology.

(4)

2 mm

where

T = C,(5)(6)
S = ML−1M T ,(7)
u = v,d−1(8)
j = h − ML e ,dt
T is the N × N symmetric positive definite matrix, S is the
N × N positive semidefinite matrix, u is the unknown vector,
and j is the current excitation vector. The generalized
eigenvalue problem related to (4) is written in the following
form [10]Sφ = λ Tφ .
(9)In this case, N eigenvalues (λ1 , λ 2 , ..., λ N ) and corresponding
eigenvectors φ can be obtained by solving (9). By defining

2 mm

100 mm

100 mm(a)(b)
Fig. 2. The rectangular meshes of the PDN with extremely fine apertures.
The gray cells have rectangular shape, and the black ones represent theapertures. (a) The top-view. (b) The closeup of the central part.

the matrix Φ = [φ1 φ 2 ... φ N ] , which stores the N eigenvectors
in its colmuns, and premultiplying both sides of (4) by ΦT ,(4) can be rewritten asd2
y + Dy = ΦT j ,
(10)
dt 2
where u = Φy , Φ T TΦ = I [11], I is the unit matrix,
D = Φ T SΦ = diag{λ1 , λ2 , ..., λN } , diag{⋅ } denotes a diagonal

matrix, and λi and φi (i = 1, 2, ..., N ) are the i-th eigenvalue
and eigenvector, respectively. Additionally, y is the unknown
coefficient vector which is derived by projecting u on the
eigenspace Φ , and φi corresponds to each mode obtained by
the mode decomposition of (4) [10]. Applying the explicit
central-difference scheme to (10) and rearranging the equation
lead to the following updating formula of yy n+1 = 2y n − y n−1 − ∆t 2 Dy n + ∆t 2ΦT j ,(11)where n is the time index, and ∆t is the time step size. In thiscase, ∆t must fulfill the numerical stability condition [10],(12)∆t 2λ2i ≤ 4, i = 1, 2, ..., N .
BecauseE ij

λi is an angular resonance frequency associatedwith (4), if high-frequency components are contained in the
numerical solutions, ∆t becomes a very small value.
Therefore, the calculation cost increases dramatically, and theabove explicit method is not effective especially for thesimulation of the object which has extremely fine structures.

w

2 Lpucld

C puc

2 L puc

4

2 L puc

C puc
4

C puc

Dielectrics

Perfect electric conductor

4

2 L puc

C puc
4

(a)
(b)
Fig. 3. The unit-cell constituting the PDN. (a) The unit-cell and (b) its Pitype equivalent circuit.

In the proposed method, the explicit updating formula (11) is
improved as follows: First, in order to obtain the accurate
numerical solutions in the transient analysis, we make include
the component related to not less than ten times higher
frequency than the maximum one, f max , into the solutions. To
do so, the time step size ∆t has to be satisfied with1~∆t ≤ ∆ tmax =f max .
(13)
10~
Next, in the case that ∆t = ∆ tmax , the eigenvectors stored inΦ are divided into two groups, namely, Φ l and Φ h so that
Φ = [Φ l Φ h ] : eigenvalues associated with the eigenvectors inΦ h do not fulfil (12), whereas those associated with Φ l do.
The unknowns in y are also divided into two groups related toΦ l and Φ h so that y = [y Tl y Th ]T . In other words, wecategorize the stable modes as Φ l and the unstable modes asΦ h . Then, removing the unstable modes Φ h and redefining

III. NUMERICAL RESULTS
The PDN illustrated in Fig. 2 is analysed to verify theaccuracy and efficiency of the proposed method. The PDN is
composed of the parallel perfect electric conductor planes andFR4 between the planes. The size of the PDN is 100
mm × 100 mm, and the thickness of FR4 is 0.5 mm. Since 16
square apertures, each of which is 0.2 mm on a side, exist on
the central part of the PDN, the narrower grid spacing thanthat in the other part is used around the apertures. The relative
permittivity of FR4 is ε r = 4.0 . In this example, we assume
that the maximum frequency f max = 5.0 GHz , and the wave
propagation speed cr = 1.5 × 108 m/s in the media, and theminimum wavelength becomes 30 mm. To perform accurate
simulations, the maximum cell size is set to 2.0 mm, which isless than one-tenth of the wavelength. Since the apertures are
extremely small compared to the PDN, the low reactancecomponents are produced. As a result, the PDN is divided
using nonuniform rectangular meshes into 65 × 65 unit cells,
one of which is shown in Fig. 3(a). The unit cell is modeled as
the Pi-type equivalent circuit illustrated in Fig. 3(b). In ourcase, the per-unit-cell inductance Lpuc and capacitance Cpuc in
the Pi-model are calculated as follows:
dl
,
(14)
L puc = µ 0w
wl
C puc = ε 0ε r,
(15)dwhere µ 0 and ε 0 are the permeability and permittivity of thevacuum, l is the cell length, w is the cell width, and d is the
distance between the planes. The number of the nodes of theequivalent circuit is 4221, and therefore, 4221 eigenvalues arecalculated from S and T. Fig. 5 shows the time step sizes,
each of which is the maximum value satisfying (12) andrelated to the eigenvalue λi . According to (13) and f max , the~time step size ∆ tmax is calculated to 20 ps and shown as the
dot-dashed line in Fig. 5. Clearly, 3813 time step sizes are less~~than ∆ tmax , and therefore, the unstable modes arise if ∆ tmax isused in the existing explicit scheme. As a result, the existing
scheme has to use the minimum time step size ∆t1 = 0.9 ps ,which increases the computational cost of the transient
simulation. These very small time step sizes come from low
reactance components in the equivalent circuit related to theextremely fine structures in the computational domain, namely,
PDN. On the other hand, by removing the unstable modes,~~∆ tmax , which is about 20 times larger than ∆ t1 , can be used in
the proposed method. A current source is appended to thenode at the bottom-left corner of the PDN as an input

−4

10

−6

10
Time step size

Φ and y so that Φ = Φ l and y = y l lead to the stable andexplicit updating formulation in the form of (11). thenumerical solutions obtained from the proposed updating
formula are calculated stably without the direct matrix~
operations even if ∆ tmax is used.

−8

10

~
∆ tmax
−10

10

−12

10

0

2000

4000

Index of the eigenvalue
Fig. 4. Index of the eigenvalue versus time step size ∆ti .

excitation. The input has the triangular waveform of which therise and fall times are 1 ns, the amplitude is 0.01 A. We
perform the transient analysis by using the proposed method,
LIM, and HSPICE.
Figs. 5, 6, and 7 show the voltage waveforms at the topright corner obtained by setting the time step size to 0.9 ps or20 ps. In these figures, the waveform obtained using HSPICE
with the time step size of 1 ps is shown as Reference. As
illustrated in Fig. 5, the waveform of HSPICE changes withthe time step size except for the case of 1 ps, and the reference
waveform agrees with that obtained from the LIM. Therefore,
we regard the reference waveform as the most accurate one inthis verification. From Fig. 6, we can see that the waveform of
the proposed method with 0.9 ps time step size completely
agrees with the reference waveform. In addition, in Fig. 7, theproposed method can provide the numerically stable solutions
with appropriate accuracy even if 20 ps is used. Table I shows
the CPU time in the case that the simulation interval is from 0
ns to 100 ns and the speed-up ratio between HSPICE with~∆ tmax = 20 ps and the other methods. From Table I, theproposed method is about 1048 times faster than the LIM andabout 21462 times faster than HSPICE. As a result, it isconfirmed that the proposed method is as accurate as theconventional solvers and much more efficient than them in thesimulation of the PDN with the extremely-small apertures.
IV. CONCLUSION
In this paper, we have proposed the method which canovercome the time step size limitation of the existing explicit
difference scheme by removing the unstable modes. The
unstable modes were specified by using the eigenvalue
decomposition and defining the maximum frequency ofinterest. The numerical results showed that the proposed
method was comparable to the LIM and HSPICE from aviewpoint of the accuracy. In addition, it has been confirmedthat our approach was much faster than HSPICE and could
reduce the computational cost of the existing scheme
effectively.

0.03

0.02

0.02
Voltage [V]Voltage [V]

0.03

0.01

0
0

Reference

HSPICE ( Δt = 0.9 ps)

LIM

HSPICE ( Δt = 20 ps)

2

0.01
Proposed method
Reference004Time [ns]Fig. 5. Comparison of the voltage waveforms between HSPICE and
LIM.

2

4

Time [ns]Fig.7. Voltage waveforms in the case that Δt = 20 ps.

TABLE 1

0.03

CPU TIME

MethodVoltage [V]CPU time (ms)

Speed-up

0.9

LIM

20

N/A

N/A

HSPICE

0.981435300.04

20343390

1

0.02

0.01

∆t (ps)

Proposed

3140

109

20

16

21462

0.9

16765

20

Proposed method
Reference

0
0

2

[7]

4

Time [ns]Fig.6. Voltage waveforms in the case that Δt = 0.9 ps.

[8]

ACKNOWLEDGMENT
This work is partially supported by Japan Society for the
Promotion of Science (JSPS).

[9]REFERENCES
[1]

[2]

[3]

[4]

[5]

[6]

A. E. Engin, K. Bharath, and M. Swaminathan, “Multilayered
finitedifference method (MFDM) for modeling of package and printed
circuit board planes,” IEEE Trans. Electromagn. Compat., vol. 49, no.
2, pp. 441–447, May 2007.
J. E. Schutt-Aine, “Latency insertion method (LIM) for the fasttransient simulation of large networks,” IEEE Trans. Circuits Syst. I,
vol. 48, no. 1, pp. 81–89, Jan. 2001.
T. Sekine and H. Asai, “Block latency insertion method (block-LIM)
for fast transient simulation of tightly coupled transmission lines,”
IEEE Trans. Electromagn. Compat., vol. 53, no. 1, pp. 193–201, Feb.
2011.
A. Taflove and S. C. Hagness, Computational Electrodynamics: The
Finite-Difference Time-Domain Method, 3rd ed. Boston: Artech House
Inc., June 2005.
D. Jiao and J. M. Jin, “Finite element analysis in time domain,” in The
Finite Element Method in Electromagnetics,” New York: John Wiley& Sons, 2002, pp. 529-584..
J. Thomas, “Numerical partial differential equations: Finite difference
methods.” New York, NY: Springer, Sep. 1995.

[10]

[11]

E. L. Tan, "Fundamental schemes for efficient unconditionally stable
implicit finite-difference time-domain methods," IEEE Trans.
Antennas Propag., vol. 56, no. 1, pp. 170-177, Jan 2008.
H. Kurobe, T. Sekine, and H. Asai, "Alternating direction explicitlatency insertion method (ADE-LIM) for the fast transient simulationof transmission lines," IEEE Trans. Compon., Packag., Manuf.
Technol., vol. 2, no. 5, pp. 783-792, May 2012.
H. Kurobe, T. Sekine, and H. Asai, "Locally implicit LIM for thesimulation of PDN modeled by triangular meshes," IEEE Microw.
Wireless Compon. Lett., vol. 22, no. 6, pp. 291-293, Jun. 2012.
Q. He and D. Jiao, “An explicit time-domain finite-element methodthat is unconditionally stable,” in Proc. APSURSI 2011, Spokane, WA,
Jul. 2011, pp. 2976-2979.
Gilbert Strang, Linear Algebra and Its Applications, Academic Press,
Inc,. New York: 1976.

Fast Transient Analysis of Multiconductor
Transmission Lines Using Nodal Block
Relaxation (NBR) Method
Takahiro Takasaki†,* , Tadatoshi Sekine††, **, and Hideki Asai†, ††, ***†Department of Systems Engineering, Graduate School of Engineering, Shizuoka University††
Department of Systems Engineering, Shizuoka University
3-5-1, Johoku, Naka-ku, Hamamatsu-shi, 432-8561 Japan
E-mail: {*takasaki, **sekine}@tzasai7.sys.eng.shizuoka.ac.jp***
hideasai@sys.eng.shizuoka.ac.jp

Abstract— This paper describes a nodal block relaxation(NBR) method in which the block relaxation method is applied
to the nodal circuit equation restructured from the modified
nodal equation. First, the circuit equation formulated by the
RLCG-MNA method is transformed into the nodal equation
composed of only node voltage variables. As a result, the block
relaxation method becomes available in the numerical solutionsof tightly coupled multiconductor transmission lines. Finally,
numerical results show that the proposed method is efficient for
the fast simulation of multiconductor transmission lines.
Fig. 1.

I.

INTRODUCTION

In order to estimate signal delay, reflection, and cross-talkon tightly coupled multiconductor transmission lines (MTLs),
circuit-based modeling and simulation techniques areadequate in terms of accuracy and efficiency. Althoughgeneral-purpose SPICE-like simulators are applicable to sucha simulation, their algorithms with sparse matrix techniques
are not efficient in the transient analysis of a coupled network.
Recently, the block latency insertion method (block-LIM) has
been proposed for the fast simulation of the MTLs, but it has
a strict numerical stability condition, which may reduce theefficiency of the block-LIM significantly [1]. On the other
hand, the nodal analysis (NA) approach used in [2] is usefulfor the simulation of interconnects because it eliminates thecurrent variables so that only the node voltages are variables
in a circuit equation. However, the NA approach has not been
applied to the coupled network explicitly.
In this paper, we apply the NA approach to the equivalentcircuit of the MTLs along with the block-relaxation method.
The nodal block relaxation (NBR) method enables to exploit
several times larger time step size than that of the block-LIM
and works effectively to separate locally dense parts in the
MTLs. And its iterative solution converges rapidly due to thepositive definite coefficient matrix of the circuit equation
obtained by the NA formulation. The algorithm of the
proposed method is similar to the block-LIM, but it has
superior features related to the numerical stability andconvergence.
The remainder part of the paper is organized as follows. In
Section II, we begin with the RLCG-MNA formulation andgenerate the iterative updating formula of the proposed

The equivalent circuit of MTLs which is composed of 128
transmission lines.

method by combining the NA approach and block relaxation
method. Section III shows some numerical results, and theconclusion is given in Section IV.
II.

THE PROPOSED METHOD

The equivalent circuit of the MTLs is shown in Fig. 1. To
derive the circuit equation associated with the network in Fig.
1, we first apply the RLCG-MNA method [3] to the circuit.
The RLCG-MNA method can generate a low-dimensional
circuit equation by eliminating the node voltage variablecorresponding to the intermediate node in every resistorinductor series-connected branch as follows.

 G M TM   v  C 0  d  v  h ,R   i   0 L  dt  i   e (1)where v, i, h, and e are the node voltage, branch current,
independent current source, and independent voltage source
vectors, and R, L, C, and G are the resistance, inductance,
capacitance, and conductance matrices. The matrix M is theincidence matrix which represents the relation of connection
between the voltages and currents.
In the next step, applying the backward difference method
at the (n+1)-th step to (1) leads to 1t C  G

T
 M

  v n 1   1t C

L  R   i n 1   0M1
t0   v n  h n 1 
   
 , (2)L  i n   e n 1 

1
t

0

128

256

384

512

640

768

0

128

256

384

512

640

768

0

0

0

0

128

128

128

256

256

256

384

384

384

512

512

512

640

640

640

768

768

(a)
Fig. 2.

128

256

384

512

640

768

768

(b)

(c)

Structural representations of coefficient matrices.
(a) T in (5). (b) P in (6). (c) Q in (6).

where Δt is the time step size, and n denotes the time step. In
this case, because the coefficient matrix in the left-hand sideof (2) is not a positive definite matrix, it is not suitable for arelaxation-based numerical simulation.
In order to overcome this drawback, we follow the NA
approach: First, (2) is partitioned into two equations1
 1
 n 1n 1
Cv n  h n 1 ,
 C  G  v  Mi 

t

t



(3)

1
 1

 M T v n 1   L  R i n 1 
Li n  e n 1 .
t
 t


(4)

Next, eliminating in+1 from (3) by using (4) leads to
Tv n 1 

1
ˆ 1 Li n  s ,Cv n  MK
t
t

(5)

where
-1

T

1
ˆ e n 1 .
ˆ MT , Kˆ   1 L  R  , s  h n 1  MK
C  G  MK
t
 t


The coefficient matrix T in (5) is shown in Fig. 2(a), and we
can see that T is a block tridiagonal matrix. Furthermore, it
can be proven that T is the symmetric positive definite matrix,
and therefore, the convergence of the block relaxation method
is guaranteed. Subsequently, in order to apply the block
relaxation method to (5), T is split into the block lower
triangular matrix P and upper triangular matrix Q as shown in
Fig. 2(b) and (c). Then the iterative updating formula isobtained and written as
Pv n 1,( k 1)  Qvn 1,( k ) 1
ˆ 1 Li n  s ,Cv n  MK
t
t

(6)

where k is the iteration counter. As described above, because
T generated by the NA approach is the symmetric positive
definite matrix, T ensures the convergence of the iterative
solutions. After updating the node voltages, the unknown
current variables are updated by using (4).
III.

NUMERICAL RESULTS

We perform the transient analyses of the circuit shown in Fig.
1 by using the proposed method and block-LIM. We set theFig. 3.

Waveform results of vobs.

TABLE ICOMPARISON OF CPU TIME
Methodblock-LIM (Δtmax)
Proposed method (Δtmax, Nk=1)
Proposed method (5Δtmax, Nk=1)
Proposed method (5Δtmax, Nk=2)

CPU time (s)
38.44
45.29
9.11
11.44

Speedup
1.00
0.85
4.21
3.36

element values as follows: R = 100 Ω, L = 1.0 nH, C = 0.01
pF, G = 0.02 S, Lm = 0.03 nH, and Cm = 1.0 fF, where Lm andCm are the mutual inductance and mutual capacitance.
Furthermore, that we use L/100 and Lm/100 as the inductanceand mutual inductance within the dotted rectangle in Fig. 1.
The input current sources Hin1 and Hin2 are trapezoidal pulses.
These sources are the same trapezoidal pulse of which theinitial value is 0 A, pulse value is 0.02 A, rise and fall times
are 0.2 ns, pulse width is 1.0 ns and period is 2.2 ns, except
for the delay times, which are 0.5 ns and 0.6 ns, respectively.
The number of the iterations used in the proposed method isfixed and equal to Nk. The simulation interval is from 0 ns to
4.5 ns.
The waveform results of vobs are illustrated in Fig. (3),where Δtmax is the maximum time step size used in the blockLIM. In Fig. 3, it is confirmed that the waveform results ofthe proposed method are numerically stable and agree withthat of the block-LIM. The CPU times and speed-up ratios tothe block-LIM are shown in Table I. In Table I, it isconfirmed that the proposed method with 5Δtmax and Nk = 2 isabout 3.36 times faster than the block-LIM without losing
accuracy. Note that the proposed method is substantially fast
because the block-LIM itself is more than ten times fasterthan HSPICE.
IV. CONCLUSION
In this paper, we have proposed the efficient circuit
simulation technique, NBR, based on the NA approach andblock relaxation method. In the proposed method, the lowdimensional circuit equation was generated by combining the
RLCG-MNA method and the NA approach. Finally, NBR
method was applied to solve the resultant equation.
Numerical results showed that the proposed method wasabout 3.36 times faster than the block-LIM with appropriate
accuracy.
REFERENCES[1] T. Sekine and H. Asai, “Block latency insertion method (blockLIM) for fast transient simulation of tightly coupledtransmission lines,” IEEE Trans. Electromagn. Compat., vol 53,
no. 1, pp. 193–201, Feb. 2011.
[2] T. H. Chen, C. Luk, and C. C.-P. Chen, “INDUCTWISE:
Inductance-wise interconnect simulator and extractor,” IEEE
Trans. Comput.-Aided Des. Integr. Circuits Syst., vol. 22, no. 7,
pp. 884–894, Jul. 2003.
[3] Y. Tanji, T. Watanabe, H. Kubota, and H. Asai, “Large scale rlc
circuit analysis using RLCG-MNA formulation,” in Proc.
Date ’06, Munich, Germany, May 2006, pp. 45–46.

A Novel Broadband Bandstop Filter Design
Tsui-Wei Weng and Tzong-Lin Wu
Department of Electrical Engineering
National Taiwan University
10617 Taipei, Taiwan
r00942011@ntu.edu.tw, wtl@cc.ee.ntu.edu.tw
Abstract—A novel wide stopband and compact bandstop filter(BSF) is proposed and realized based on defected groundstructure (DGS) and meandered signal line. A condition fortransmission zeros is derived. It can be observed that theequivalent circuit model reasonably agrees with the full-wave
simulation and measurement results. Also, it is found that the
proposed filter has high rejection over 25 dB in its stopband. The
filter has a compact size of 0.21 λg × 0.19 λg, where λg is the
wavelength of the center frequency of the stopband, and its -10
dB fractional bandwidth (FBW) is up to 90%.

Ground
St

b

Patch
p1

INTRODUCTION

Defected ground structure (DGS) has been widely used in
microwave and millimeter wave circuits recently for size
reduction. It is realized by etching defects on the ground plane.
The function of DGS is that the defect can disturb returncurrents, which results in additional inductance andcapacitance to create resonance at some frequency. The
resonance property provides a band-rejection characteristic,
which can be directly used in designing BSFs.
Several approaches [1]-[9] have been proposed for BSF
design based on DGS. Conventional DGS has the shape ofdumbbell [1] or spiral [2]. Later on, various shapes have beenbrought out and discussed, such as UV slots [3], double-H
shaped [4], modified dumbbell shaped [5], T shaped [6], E
shaped [7], and modified-E shaped [8]. However, the common
problem of [1]-[8] is that only one transmission zero (TZ) is
generated by a unit cell, which in turn results in narrow
stopband bandwidth of BSF. To achieve a wider stopband, one
method is cascading several unit cells, like three cascaded Vslots DGS in [3], which is at the expense of occupying larger
surface area. On the other hand, double equilateral U shaped
DGS in [9] can generate two TZs. Nevertheless, the rejection
between two TZ is not high; therefore, if -10 dB rejection isrequired, the unit cell in [9] still cannot create a wide stopband.
In this paper, a novel BSF is proposed by using a DGS andmeandered signal line. The proposed BSF can generate morethan one TZ to achieve wide stopband. In addition, thecorresponding equivalent circuit model is constructed andanalyzed, and the TZs condition is derived and discussed.
Finally, a test sample is designed and fabricated. The surface
area of the proposed BSF is 7.6 mm × 6.9 mm, which is small.
It is only 0.04 λg2, in terms of the wavelength propagating in

Bridgep2

Keywords- Bandstop filter (BSF); defected ground structure(DGS).

I.

G

b

G

G
G

(a)
W
Signalh
Substrate
Groundt(b)
Figure 1. The proposed BSF: (a) Top view (b) Side view

the substrate. Besides, it has a wide stopband, which covers
from 2.24 GHz to 5.93 GHz below -10 dB, and possesses alarge FBW equal to 90% by measurement.
II.

STRUCTURE AND DESIGN CONCEPTThe proposed BSF is implemented in a two-layered
structure, and its top view and side view are shown in Figs.

Figure 2. Test sample of the proposed BSF

displays the return loss and the insertion loss of measurement
results and simulation results (full-wave and equivalent circuit
model), respectively. Good agreement between measurement
results and full-wave simulation is observed. The cutoff
frequency is defined by -10 dB, and a wide stopband from
2.24 GHz to 5.93 GHz for measurement can be observed inFig. 3.

0|S11|

S-para

-10
-20

|S21|

-30
Measurement
Full-wave
Model

-40

REFERENCES

-50
0

1

2
3
4
Frequency (GHz)5

6

Figure 3. Return loss and insertion loss of the proposed BSF.

1(a) and (b). As can be shown in Fig. 1(a), the filter
comprises a meandered signal trace on the top layer, and aground plane with DGS on the bottom layer. The DGS iscomposed of two parts, one is Bridge, and the other is Patch.
Bridge is a connection of Ground and Patch, and the width ofbridge is denoted as b. On the other hand, Patch is a floating
ground with slots surrounded, and p1, p2 are used to describe
the size of Patch, and G is the gap between Patch and Ground.
In addition, total length of the meandered signal line, line
width, substrate height and metal thickness are denoted as St,
W, h, and t, respectively.
As we know, multiple TZs can produce a wide stopband;
therefore, it is essential to know how TZs are introduced by
the proposed filter. The TZ condition is derived as follows:2
1  4 2 LDGS CDGS f Z 2 ( LDGS  Lm )sin EM  f Z . (1)
ZM
III.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

SIMULATION AND MEASUREMENTA BSF test sample is fabricated with FR4 substrate andoccupies an area of 7.6 mm × 6.9 mm, as shown in Fig. 2. The
dielectric constant is 4.3, and the loss tangent is 0.02. Fig. 3[9]D. Ahn, J. Park, C. Kim, J. Kim, Y. Qian, and T. Itoh, “A design of the
low-pass filter using the novel microstrip defected ground structure,”
IEEE Trans. Microw. Theory Tech., vol. 49, no. 1, pp. 86–93, Jan. 2001.
C. Kim, J. Lim, S. Nam, K. Kang, J. Park, G. Kim, and D. Ahn, “Theequivalent circuit modeling of defected ground structure with spiralshape,” in IEEE MTT-S Int. Microw. Symp. Dig., Jun. 2002, vol. 3, pp.
2125–2128.
D. J. Woo, T. K. Lee, J. W. Lee, C. S. Pyo, W. K. Choi, "Novel U-slot
and V-slot DGSs for bandstop filter with improved Q factor," IEEE
Transactions on Microwave Theory and Techniques, vol. 54, no.6, pp.
2840-2847, June 2006.
M. K. Mandal and S. S. , “A novel defected ground structure for planar
circuits,” IEEE Microw. Wireless Compon. Lett., vol. 16, no. 2, pp. 93–
95, Feb. 2006.
A. M. E. Safwat, F. Podevin, P. Ferrari, and A. Vilcot, “Tunablebandstop defected ground structure resonator using reconfigurable
dumbbell-shaped coplanar waveguide,” IEEE Trans. Microw. Theory
Tech.,vol. 54, no. 9, pp. 3559–3564, Sep. 2006.
X. Wang, B. Wang, H. Zhang, and K. J. Chen, “A tunable bandstop
resonator based on a compact slotted ground structure,” IEEE
Trans.Microw. Theory Tech., vol. 55, no. 9, pp. 1912–1917, Sep. 2007.
S.-Y. Huang and Y.-H. Lee, "A compact E-shaped patterned groundstructure and its applications to tunable bandstop resonator, " IEEE
Trans. on MTT, Vol. 57, No. 3, 657-666, 2009.
S.U. Rehman, A.F. Sheta, and M. Alkanhal, " Compact bandstop filter
using defected ground structure (DGS), " 2011 Saudi International
Electronics, Communications and Photonics Conference (SIECPC), pp.
1-4, 2011.
S.-W. Ting, K.-W. Tam, and R. P. Martins," Miniaturized microstrip
lowpass filter with wide stopband using double equilateral U-shaped
defected ground structure, " IEEE Microw. Wireless Compon. Lett., vol.
16, no. 5, pp. 240-242, May. 2006.

A Solution and Its Estimation Method for Slotcrossing Signals to Reduce ISI-increased Crosstalk
Yu-Jen Chang1, Chiu-Chih Chou, Hao-Hsiang Chuang, and Tzong-Lin Wu3
Department of Electrical Engineering and Graduate Institute of Communication Engineering
National Taiwan University, Taipei, 10617, Taiwan
1r00942014@ntu.edu.tw3
wtl@cc.ee.ntu.edu.tw
Abstract—A branched reflector technique for the suppression ofslot-induced crosstalk especially including the effect of intersymbol interference (ISI) is newly proposed. A fast algorithm forthe performance evaluation of the proposed technique is also
provided. To fulfil this algorithm, a flow chart is suggested toquickly estimate the slot-induced crosstalk for arbitrary system
excitation with any operating speed. Further, the estimatedresults obtained are compared with the measured ones. Goodagreement between them confirms the validity of the proposed
algorithm. Also, the effectiveness of the proposed technique iswell demonstrated based on these obtained results.

I. INTRODUCTION
In order to satisfy the demand for high performance, lowcost and compact electronic products nowadays, more andmore highly integrated multi-chip modules (MCM), systemin-a-packages (SIP), and high-density printed circuit boards(PCB) are employed in design. For such, the isolation of
power/ground noise and the separation of different voltage
supplies related to different chips become critical. Splitting
power/ground planes has been verified to be an effective andlow-cost solution to these problems [1], [2]. However, highspeed signal lines sometimes unavoidably cross the slotted
power/ground planes under high-density routing. Manysignal-integrity issues such as reflection, crosstalk, andradiated emission thus arise [3]-[5].
Many solutions have been proposed to reduce the slotinduced crosstalk stated above. The most commonly usedmethod is to add stitching capacitors [5], which can provide areturn current path at high frequency range and still maintainthe isolation between the split power/ground planes at thedirect current (DC) level. However, these additional
capacitors increase the fabrication cost. Furthermore, theeffective series inductance (ESL) of the component-type
capacitor invalidates the capacitor at high frequency range andthe potential LC resonance may result in longer duration ofcrosstalk.
Recently, a branched reflector technique for slot-induced
crosstalk suppression has been proposed [6]. Its physicalmechanism and design method have been investigated only
for the step input source; nevertheless, the crosstalk may behighly increased due to the inter-symbol interference (ISI). In
this paper, the design method of branched reflector technique
is refined and able to reduce the slot-induced crosstalkespecially increased by the ISI effect, for which the pseudo
random bit sequence (PRBS) is considered as the input source.

To estimate this performance, a new fast algorithm evolved
from the previous one [7] and its flow chart are also proposed.
This new algorithm can take the ISI effect into account andrapidly assess the crosstalk for arbitrary system excitation
with any operating speed.
II. PROPOSED TECHNIQUE AND ESTIMATION METHOD
A. Proposed Solution and Its Equivalent Circuit Model
Fig. 1(a) shows the physical structure of the solution
board with a branched reflector. On the top layer, two
identical microstrip lines with length 2∙L2 and width Wm areprinted. The distance between these two lines is 2∙Ls2. The
bottom layer is a ground plane. A slot line with width Ws isetched in the middle of this plane and both ends of the slot
line are kept a distance L1 to the plane edges. The spacing
between each end of the slot line and its adjacent microstrip
line is Ls1. An open stub with width Wbr and length Lbr is
printed on the top layer and connected to the ground plane bya via in the center of the two lines. This open stub is called
branched reflector hereafter. The branched reflector canhighly decrease the slot-induced crosstalk if Lbr and Wbr areproperly designed [6].
Fig. 1(b) shows the equivalent circuit model of thestructure shown in Fig. 1(a). In this model, the characteristic
impedance of each microstrip line is represented by Zm. For
simplicity, the frequency-dependant characteristic impedanceof the slot line is assumed to be a constant Zs while thecharacteristic impedance of the branched reflector is Zbr. From
this equivalent circuit model, it can be easily observed that theslot-induced near-end and far-end crosstalk are identical inmagnitude. Therefore, only the near-end crosstalk will beconsidered in the following sections.
B. Brief Introduction for Algorithm
The first step of the proposed algorithm is to obtain thestep response of the considered system. This response can beused to rapidly predict the maximum ISI-increased crosstalk
incorporating with the following steps of this algorithm. Oncethe step response is obtained, the maximum crosstalk voltage(Vmax) for the PRBS input with a specific bit-rate can becalculated as [7](1)Vmax  max{Vmax, n },
n

B’

A’

V
L1

1st UI
T

Ls1

V1,n

m-th UI
T
Vm,n

Wm
Ls2
Wbr

Group nLs2 LbrWm(n-1)T/N

Ls1
Ws

L2

B

A

t

Fig. 2. A typical step response.

L1

L2(n-1)T/N

metal
substrate
via

A-A’

B-B’
(a)

L2

Zm
Vin

Ls1

Zs

L2

Zm

Zm
Ls2

Zs

Ls2

Zs

Zm

Lbr
Zbr

Zm

Vnear

Zm

Zm

L2

L2
Microstrip line
Slot line
Branched reflector

Ls1

Zs

Vfar

Zm

(b)
Fig. 1. The proposed structure: (a) its top and side views, and (b) its
corresponding equivalent circuit model.

where Vmax,n denotes the maximum voltage of the group n (seeFig. 2) and can be expressed as
Vmax, n   M i , n   m j , n .
(2)
i

j

In Fig. 2, Mi,n denotes the ith local maximum voltage in the setof group n, while mj,n denotes the jth local minimum voltage inthe same set. The duration of each unit interval (UI) is denoted
as T and the number of the sampling points in each UI isrepresented by N.
Fig. 3 shows the flow chart of the proposed algorithm. In
this chart, three steps surrounded by a dashed line aresuggested by the previous research [7]. According to this chart,many steps are essential to fast and successful estimation of
the ISI-increased crosstalk for arbitrary system excitation withany operating speed. First, the step response of the crosstalkfor the system can be obtained by circuit simulation or evenmeasurement. Second, the UI size is determined by theoperating bit-rate of the system. Then, (1) and (2) can beapplied to get the maximum crosstalk with respect to the given

Fig. 3. The flow chart of the proposed algorithm.

step response. Third, if there is any other bit rate in theinterested range, the iterating procedure proceeds with the UI
size changed. It is worth mentioning that by using thisalgorithm, the influence of ISI-increased crosstalk withdifferent operating speed can be well evaluated without
running the time-consuming simulation for step response
again. Furthermore, it provides a straightforward assessment
for efficiently analyzing the proposed technique.
III. PREDICTION OF WORST ISI-INCREASED CROSSTALKIn this section, the effectiveness of the proposed techniqueand its evaluation following the algorithm as described in Part
B of Section II will be well verified. Then, some important
properties related to this technique will be revealed.
A. Numerical Result
Table I summarizes the physical dimension of the proposed
structure (see Fig. 1(a)), which would be viewed as a solution
board for demonstration. The reference board without thebranched reflector has the same dimension as the solution

TABLE I
THE VALUES OF STRUCTURE PARAMETERS
Variable
L1

Length (mm)Variable
Lbr

Length (mm)

25

Ls1

25

Wm

1.5

Ls 2

75

Ws

1

L2

50

Wbr

1.5

20

The reduction ratio is positive in almost all the bit-rate range
except in the range from 0.75 to 0.9 Gbps, which means thatthe branched reflector can perform well except in this narrow
range. However, if we only observe the step response shownin Fig. 4(a), it seems that the reduction ratio will always be
about 50%. This indicates that the prediction only based on
the step response is not sufficient. In fact, the peak at 0.82
Gbps is primarily caused by the half-wavelength resonance ofthe slot. Nonetheless, the length of the branched reflector, Lbr,
can slightly affect the resonant frequency. The peak can be
moved to lower/higher frequency by increasing/decreasing thevalue of Lbr. It is worth noting that the increase/decrease of Lbrshould be kept small because the change will degrade theimprovement of Vmax at other bit rates.
In Fig. 5(a), the frequency response of Vmax can be divided
into three characteristic regions by two frequency points, fL
and fH. These frequencies can be approximated as
1
fL  ,
(4)
TL

1
.
(5)
TH
In (4), TL is defined by the duration time from the peak valueof the step response to 5% of its peak value. In this case, thevalue of TL is about 10 ns and the resultant fL is 0.1 Gbps.
Meanwhile, TH in (5) is defined by the time difference
between the first transition V1,s (V1,r) and the second transition
V2,s (V2,r) of the step response. In this case, the value of TH isabout 0.28 ns and the resultant fH is 3.5 Gbps.
The frequency response of Vmax is almost constant in thetwo regions below fL and above fH, while it varies in the region
between fL and fH. When the operating speed is below fL, theripple of voltage waveform will decay to an insignificant
value before the transition of next bit. Therefore, the ISI
between the adjacent bits is negligible. Indeed, the value of
Vmax in this range is the same as that of the step response. Fig.
5(b) shows Vmax,r and Vmax,s, respectively, for the reference andsolution boards at the region below fL. Comparing Vmax,r (Vmax,s)in this figure with V1,r (V1,s) observed from the step response
in Fig. 4 (b), good agreement can be seen, which indicates theminor influence of ISI effect.
On the other hand, Vmax is almost frequency-independentand shows the largest value when the bit rate of PRBS is
higher than fH. The reason can be well explained by (1), whichis the predicting formula of Vmax. When the bit rate is fastenough, all the local maxima and minima of the step response
will be used to calculate Vmax in (1). Therefore, the value of
Vmax for small UI (high bit rate) is constant and also is thelargest value among all frequency ranges.
fH 

(a)

(b)
Fig. 4. The near-end crosstalk: (a) the whole waveform, and (b) the enlarged
partial waveform.

board. The equivalent circuit model in Fig. 1(b) would be used
to obtain the required step response. The characteristic
impedances for Zm, Zs, and Zbr are 50 Ω, 90 Ω, and 50 Ω,
respectively. A step voltage source (Vs) with a magnitude of
0.5 V is connected to the input end of the upper microstrip
line and the circuit simulator [8] is used for the time-domain
simulation. Fig. 4(a) shows the step response of the near-end
crosstalk (Vnear) for the solution and reference boards, and Fig.
4(b) illustrates the zoom-in view of Fig. 4 (a). Then, these two
step responses could be applied to the estimation algorithm.
Fig. 5 shows the maximum crosstalk (Vmax) with ISI effect,
which is predicted by the newly proposed algorithm.
Comparing the result of the solution board (Vmax,s) with that of
the reference board (Vmax,r), the reduction ratio (R) can be
defined asV Vmax, s(3)
R  max, r
100 (%).
Vmax, rB. Experimental Validation
In order to demonstrate the effectiveness of the proposed
technique, the solution and reference boards with the same
physical dimension as summarized in Table I are fabricatedusing typical PCB process. Each board has a thickness of 0.8mm and a relative dielectric constant of 4.4. The experimental
results are measured by a time-domain reflectometer (TDR,

(a)
(a)

(b)
(b)
Fig. 5. Maximum crosstalk with the ISI effect: (a) the entire frequency range,and (b) the enlarged low-frequency region.

Tektronix CSA 8000B). The PRBS with a magnitude of 0.5 V
and a bit rate of 1.2 Gbps generated by a pattern generator(Anritsu MP1763) having a 50-Ω source resistance is fed
intothe input end of the microstrip line. All the other ends ofthe two microstrip lines are terminated by a 50-Ω resistor. The
measured results of the near-end crosstalk for the reference
and solution boards are illustrated in Figs. 6(a) and 6(b),respectively. Both of them are displayed in the same way withthe same data-rate of 1.2 Gbps. As shown in these figures, themaximum near-end crosstalk Vmax for the solution board is
48.75 mV and that for the reference board is 82 mV, whichcorresponds to a reduction ratio (R) of about 41%. Meanwhile,
the reduction ratio estimated by the proposed algorithm asshown in Fig. 5(a) is about 40%, which is in good agreementwith the measurement results.
IV. CONCLUSION
Since the slot-crossing signal lines are inevitably used inmodern package and PCB design, a solution to the slotinduced crosstalk suppression is in urgent need. The branchedreflector technique is newly proposed and well demonstrated
to suppress the slot-induced crosstalk including the ISI effect.
Moreover, the fast prediction algorithm and its flow chart areprovided for evaluating the performance of this technique.
Through the experiment, the measurement results obtained
also show good agreement with the algorithm prediction.
Consequently, this new algorithm can be applied to rapidly

Fig. 6. Measured near-end crosstalk with the PRBS input: (a) the referenceboard, and (b) the solution board.

evaluate the ISI-increased slot-induced crosstalk for arbitrary
system excitation with any operating speed.
REFERENCES
[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

C. Wei, F. Jun, R. Yong, S. Hao, J. L. Drewniak, and R. E. DuBroff,“DC power bus noise isolation with power-plane segmentation,” IEEE
Trans. Electromagn. Compat., vol. 45, no. 2, pp. 436–443, May 2003.
T.-L. Wu, H.-H. Chuang, and T.-K. Wang, “Overview of power
integrity solutions on package and PCB: decoupling and EBG
isolation,” IEEE Trans. Electromag. Compat., vol. 52, no.2, pp. 346356, May 2010.
H.-J. Liaw and H. Merkelo, “Signal integrity issues at split ground andpower planes,” IEEE Electron. Compon. Technol. Conf., Orlando, FL,
May 1996, pp. 752-755.
H.-H. Chuang and T.-L. Wu, “A novel ground resonator technique toreduce common-mode radiation on slot-crossing differential signals,”IEEE Microw. Wireless Compon. Lett., vol. 20, no. 12, pp. 660-662,
Dec. 2010.
F. Xiao, Y. Nakada, K. Murano, and Y. Kami, “Crosstalk analysis
model for traces crossing split ground plane and its reduction by
stitching capacitor,” IEICE Trans. Electron., vol. J89-C, no. 11, pp.
885-893, 2006.
H.-H. Chuang, C.-C. Chou, Y.-J. Chang, and T.-L. Wu, “A branchedreflector technique to reduce crosstalk between slot-crossing signallines,” IEEE Microw. Wireless Compon. Lett., submitted for
publication.
Y.-S Cheng and R.-B Wu, “Direct eye diagram optimization for
arbitrary transmission lines using FIR filter,” IEEE Trans. Compon.,
Packag. Manuf. Technol., vol. 1, no. 8, pp. 1250-1258, Aug. 2011.
Agilent Technol., ADS, ver. 2009.

Circuit/Electromagnetic Hybrid Simulation ofElectrostatic Discharge in Contact Discharge Mode
Tsuyoshi Takada

Tadatoshi Sekine

Hideki Asai

Dept. of Systems Eng.,
Shizuoka University
3-5-1 Johoku, Naka-ku,
Hamamatsu-shi, 432-8561 Japan
takada@tzasai7.sys.eng.shizuoka.ac.jp

Dept. of Information Science and Tech.,
Graduate School of Science and Tech.,
Shizuoka University
3-5-1 Johoku, Naka-ku,
Hamamatsu-shi, 432-8561 Japan
sekine@tzasai7.sys.eng.shizuoka.ac.jp

Dept. of Systems Eng.,
Shizuoka University
3-5-1 Johoku, Naka-ku,
Hamamatsu-shi, 432-8561 Japanhideasai@sys.eng.shizuoka.ac.jp

Abstract—In this paper, we describe a hybrid simulation
technique for electrostatic discharge (ESD) events in contact
discharge mode. The proposed technique is based on circuit,electrostatic field, and full-wave simulations. Because the hybrid
simulation can take a low computational cost relative to existing
full-wave techniques, it can reduce the total CPU time and theamount of memory consumption.

I. I NTRODUCTION
Electrostatic discharge (ESD) causes failures and damage
to electronic equipments because the current produced by ahuman-to-metal (HTM) discharge has high-frequency components due to the fast rise time of the current waveform. In
order to ensure immunity-aware electromagnetic compatibility(EMC) design, it is necessary to simulate the actual ESD
events along the IEC 61000-4-2 [1]. Although some papers describe about the numerical modeling and simulation techniques
for the ESD [2]–[5], it is still one of challenging issues in the
EMC community to deal with the physical phenomena of thedischarge and include them for the numerical simulation.
The three-dimensional electromagnetic simulation provides
full-wave solutions as spatially-distributed field components
generated by discharging [4]–[7]. However, the computationalcost, the CPU time and memory requirement of the electromagnetic analysis tend to be large compared to the circuitbased simulation. On the other hand, some circuit modelsof ESD generators have been developed [5] and [8]. Since
the circuit model includes the elements which represent only
major effects during the discharging process, its computational
cost is much lower than that of the full-wave simulation.
In this paper, we describe a hybrid technique combining
the circuit and full-wave simulations for the ESD events inthe contact discharge mode. The proposed technique receives
benefits of the low computational cost of the circuit simulation
and, at the same time, can obtain the electromagnetic field
distribution from the full-wave simulation.
II. E XISTING FDTD S OLVER FOR THE ESD E VENTSThe existing full-wave simulation techniques applied tothe immunity test use either a commercially-based full-wave
simulator or an FDTD solver [4]–[7]. In this paper, we follow

the FDTD simulation proposed in [6] with some modifications.
In [6], a Noiseken ESD generator and a test setup according to
the IEC standard are modeled in a computational domain. The
domain is discretized into the Yee grid, and the whole domain
to be simulated is shown in Fig. 1(a). In addition, the detailed
structure of the ESD generator is also illustrated in Fig. 1(b).
The ESD generator includes the capacitor unit, resistor, andelectrode. The capacitor unit and electrode are aluminium andthe resistor has the conductivity of 2.41 S/m. In addition to
the actual objects, we insert some lumped elements into threeparts, namely, charging, switching, and target parts as shownin Fig. 1(b). Furthermore, free space surrounding the domain
is truncated by the absorbing boundary which is placed at thearea immediately external to the domain.
The transient simulation performed in [6] is separated into
two phases, namely, charging and discharging phases. In boththe phases, a time-dependent voltage source Vc is inserted into
the charging part, and its value is defined as
)
 (
−t

(0 ≤ t ≤ τ0 )V0 1−e τVc (t) =, (1)() 1∫ t−τ0 V0 1−e τ −ic (ξ)dξ (τ0 < t)Cs τ0
where τ is the time constant of charging, τ0 is the stop timeof charging, V0 is the high voltage supply, Cs is the valueof the storage capacitor, and ic is the current flowing throughVc . In fact, the first equation in (1) is used for the charging
process, and the second one relates to the discharging process.
Additionally, the target, into which the current is injected, ismodeled as the load resistor RL and inserted into the target
part shown in Fig. 1(b). As for the switching part, the capacitor
Csw is used in both the phases whereas the inductor Lsw
is inserted only in the discharging phase. In the simulations
performed in Section IV, we set the parameters so that τ = 30
ns, τ0 = 5τ , V0 = 5 kV, Csw = 0.8 pF, Lsw = 70 nH, and
RL = 2 Ω.
III. H YBRID S IMULATION FOR THE ESD E VENTS
The proposed technique is separated into three processes,
and their outlines are as follows:Capacitor unit

Ground strap

CsStandard

Cp

330

Resistor
170

Rd

308

254

Electrode
54

32

Cg

60

448

ig

Lsw Switch

Ct

Lt

it

Tip

RL Target

60
Fig. 2.

The proposed equivalent circuit model of the ESD generator.

Charging part40

TABLE IPARAMETERS IN THE CIRCUIT MODEL

Vc
Category

Switching partParameter

Value

Storage cap. (pF)

Cs

150

Discharge res. (Ω)

Rd

330

Cap. (pF)

Cp

0.9

Cap. (pF)

Csw

0.8

Ind. (nH)

Lsw

70

Tip

Ind. (nH)Lt

11

Tip-GndCap. (pF)

Ct

7.3

Gun-GndCap. (pF)Cg

4.3

Gnd strap

Ind. (nH)

Lg

600

TargetRes. (Ω)

RL

2

Excitation

Voltage source

Vs–

Standard
Parasitics

Csw

88 94

Lsw
Switch8
Electrode
Target part
RL

50

Csw

Gun-Gnd Tip-Gnd

40

Resistor

Lg

108[unit: mm]

Ground plane(a)Capacitor unit

Vs

Gnd strap Parasitics

Element type

[unit: mm](b)Fig. 1. The dimensions and structures of the computational domain including
the ESD generator, ground strap, target, and large metal plate. The absorbing
boundary exists around the domain. (a) The dimensions of the whole system.
(b) The detailed dimension and structure of the ESD generator.

1) The estimation of input currents: we perform the transient simulation of the equivalent circuit model of an
ESD generator to obtain the current injected into a target
and one flowing through a ground strap. These currents
are used as inputs for the discharging phase.
2) The calculation of initial electric field distribution: given
voltages on the conductors of a storage capacitor, potentials distributed over a computational domain arecalculated by solving the Laplace’s equation. The initialvalues of the electric fields used in the next process arecalculated from the potentials.
3) Full-wave simulation for the discharging phase: theelectromagnetic fields generated by the discharge current flowing through the target are calculated using the
FDTD-based modeling and simulation.
In fact, the first and second processes are preparation for thethird process, and therefore, we do not care the order of thefirst two processes. The above three processes are describedin detail below.
A. Circuit Simulation for the Estimation of the Input Currents
The conventional circuit models of the ESD generator in[4], [5], [8], [9] mainly include a storage capacitor Cs = 150
pF and a discharge resister Rd = 330 Ω according to the
IEC standard. Although the values of circuit parameters varydepending on a product type of an ESD generator, the existingB. Electrostatic Field Simulation for the Calculation of Initial
Electric Field Distribution
After the ESD generator is charged, it is expected that
an electric field is generated around the generator due tothe charge distributed over the conductor surfaces of thestorage capacitor. Therefore, we need to know electric field
distribution before starting the transient simulation for thedischarging phase. In order to obtain the distribution after
charging, we begin with the Laplace’s equation in threedimensional space∂ 2 φ(x, y, z) ∂ 2 φ(x, y, z) ∂ 2 φ(x, y, z)
+
+
= 0,
∂x2
∂y 2
∂z 2

(2)

where φ is the scalar potential. In this paper, we solve (2)numerically by means of the finite difference method (FDM)[10]. Then, we get a set of simultaneous equations. To solve
the simultaneous equations, we adopt a conventional iterativesolver, the successive over relaxation (SOR) method [11]. The
grid used in the FDM is compatible with the Yee grid used in
the FDTD method in the discharge phase.
C. Full-Wave Simulation for the Discharging Phase
In the third process, the full-wave simulation is performed
using those pre-calculated values. Actually, we adopt the
FDTD method and simulate the reduced computational domainshown in Fig. 3 instead of the domain defined in Fig. 1(a).
The reduced domain is derived by removing the upper part of
Fig. 1(a). Before starting the transient simulation, we set the

Ground strap

Electrode

32

54
448

60

18

108
[unit: mm]

it
RL

ig

Fig. 3. The dimensions and structures of the reduced computational domainused for the full-wave simulation of the proposed technique.

20
Discharge Current [A]models deal with some common components; the materialof the generator; parasitic coupling; ground strap; and target
objects. Especially, because the return path of the discharge
current and the parasitic elements directly affect the waveformof the discharge current, we have to include them in theequivalent circuit model.
We synthesized the circuit model of the ESD generator in
the previous section and it is illustrated in Fig. 2. And therelated circuit parameters are shown in Table I. The circuit
is composed of the elements associated with the RC circuit
used in the IEC standard, parasitic elements, materials of thegenerator body, discharge tip, and return path. In this synthesis,
the return path is represented by two elements. The first one
is the capacitor Cg , which affects the current value and therising time before the first peak. The other is the inductor Lg ,which corresponds to the ground strap and affects the current
value and the appearance time of the second peak. To performthe transient simulation by using a SPICE-like simulator, it isassumed that the network is driven by the independent voltage
source Vs which is given by a step function. This means thatthe real discharging process of the charge is replaced by the
charging process of the storage capacitor [9].
The currents are calcurated by the transient simulation ofFig. 2. We can obtain the current it flowing through the tip
inductor Lt and the current ig flowing through the ground strap
inductor Lg specified in Fig. 2. We use these currents as the
inputs for the full-wave simulation in the proposed method.

10
Proposed method
Existing method
Mathematical model
IEC standard

00

Fig. 4.
target.

30
Time [ns]60

The waveform results of the discharge current flowing through theinitial values of the electric field obtained from the electrostatic
field simulation in the second process. Then, the transient
responses are calculated step by step along the time marchingprocedure of the FDTD method. In this case, excitations areexactly the same as the currents it and ig obtained fromthe circuit simulation in the first process. These currents areinserted as independent current sources as illustrated in Fig.
3. Clearly, because the computational domain in our method
is smaller than that of the existing FDTD solver, it is expectedthat the proposed method can reduce the computational cost
significantly.
IV. N UMERICAL R ESULTS AND D ISCUSSION
We perform transient simulations using the existing FDTD
solver and the proposed hybrid technique. In the FDTD-based
simulation of both the proposed and existing methods, theobjects and space to be analyzed are discretized using a cubic
cell of 2 mm on a side. We also use a perfectly matched layer(PML) which is 16 cells thick as the absorbing boundary to
truncate the computational domain.
The waveform results of the discharge currents are shownin Fig. 4(a). For comparison, a waveform calculated by amathematical model described in [4] is also plotted as areference one in these figures. Each of the values from the
proposed method is nealy equal to that from the existing one,and all the values are within the tolerance defined in the
IEC standard. From the above results, it is confirmed that the
waveform of the discharge current obtained by the proposed
method is as accurate as the one by the existing full-wave

Magnetic field [A/m]

0

−100

Proposed method
Existing method

−200
0

30
Time [ns]60Electric field [kV/m](a)confirmed that our technique can perform about 4.4 times
faster simulation than the existing FDTD solver.
Although the proposed method is more efficient and can
produce the adequate waveform of the discharge current, theaccuracy of the field components is not good enough. The
difference between the results from the existing and the hybrid
solvers is caused by the simple hybridation in the proposed
technique. Therefore, we have to improve the ways to separate
the system to be analyzed into subdomains and combine theresults obtained from the separated subdomains. This consideration is necessary for the hybrid simulation to provide moreaccurate results of electromagnetic field distribution generated
by the ESD events.
ACKNOWLEDGEMENT0

We would like to express the special thanks to SESAME
Technology Inc. for the technical advice.

−50

R EFERENCES
Proposed method
Existing method

−100
0

30
Time [ns]60(b)Fig. 5. The waveform results of field intensities. (a) Magnetic field intensity.
(b) Electric field intensity.

solver.
The electromagnetic fields are also calculated and we observe those at the obsevation points which are placed at the
distance of 50 mm from the discharge point. The observed
electric field and magnetic field intensities are illustrated inFig. 5. From Fig. 5(a), it is confirmed that the waveformof the magnetic field intensity from the proposed method isapproximately similar to that from the existing method. In
contrast, there is a difference between the electric field intensities shown in Fig. 5(b) after the first peak. We presume that
one of the reasons for the difference between the electric field
intensities is because the current path in our computational
domain does not include the ground strap, and thereby weignore the electromagnetic effect from the ground strap in the
reduced domain.
From the numerical result, the proposed method can reduce
the whole computational time and is totally about 4.4 times
faster than the existing solver. The amount of the memory
used during the simulation of our solver is 515 296 KB andthat of the existing one is 1 758 832 KB. This means that the
proposed method can reduce the required memory capacity by
70 %.
V. C ONCLUSIONS AND F UTURE W ORKS
We have proposed the hybrid simulation technique for thesimulation of the ESD events in the contact discharge mode.
The proposed method was described separately divided into
the three processes. From the numerical results, we have[1] IEC 61000-4-2 Ed. 2.0 Electromagnetic Compatibility (EMC)—Part 4-2:
Testing and Measurement Techniques—Electrostatic Discharge Immunity
Test, Dec. 2008.
[2] R. Chundru, D. Pommerenke, K. Wang, T. V. Doren, F. P. Centola, andJ. S. Huang, “Characterization of human metal ESD reference discharge
event and correlation of generator parameters to failure levels–part I:
Reference event,” IEEE Trans. Electromagn. Compat., vol. 46, no. 4,
pp. 498–504, Nov. 2004.
[3] ——, “Characterization of human metal ESD reference discharge eventand correlation of generator parameters to failure levels–part II: Correlation of generator parameters to failure levels,” IEEE Trans. Electromagn.
Compat., vol. 46, no. 4, pp. 505–511, Nov. 2004.
[4] K. Wang, D. Pommerenke, R. Chundru, T. V. Doren, J. L. Drewniak,
and A. Shashindranath, “Numerical modeling of electrostatic discharge
generators,” IEEE Trans. Electromagn. Compat., vol. 45, no. 2, pp. 258–
271, May 2003.
[5] S. Caniggia and F. Maradei, “Circuit and numerical modeling of electrostatic discharge generators,” IEEE Trans. Ind. Appl., vol. 42, no. 6,
pp. 1350–1357, Nov. 2006.
[6] O. Fujiwara, X. Zhang, and Y. Yamanaka, “FDTD simulation of electrostatic discharge current by testing,” IEICE Trans. Commun (Japanese
Edition), vol. J86-B, no. 11, pp. 2390–2396, Nov. 2003.
[7] C. Qing, J. Koo, A. Nandy, D. Pommerenke, J. S. Lee, and B. S. Seol,“Advanced full wave ESD generator model for system level coupling
simulation,” in Proc. IEEE EMC Symposium 2008, Detroit, MI, Aug.
2008, pp. 1–6.
[8] H. Tanaka, O. Fujiwara, and Y. Yamanaka, “A circuit approach tosimulate discharge, current injected in contact with an ESD-gun,” inProc. IEEE EMC Symposium 2002, Minneapolis, MN, May 2002, pp.
1–6.
[9] J. Koo, Q. Cai, G. Muchaidze, A. Martwick, K. Wang, and D. Pommerenke, “Frequency-domain measurement method for the analysis of
ESD generators and coupling,” IEEE Trans. Electromagn. Compat.,vol. 49, no. 3, pp. 504–511, Aug. 2007.
[10] M. N. O. Sadiku, Numerical Techniques in Electromagnetics with
MATLAB, 3rd ed. Boca Raton, FL: CRC Press, Apr. 2009.
[11] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed. Philadelphia, PA: SIAM, Apr. 2003.

Output Signal Synthesis of Nonlinear Digital Channel
Using Superposition of Multiple Edge Responses
Chiu-Chih Chou, Tzong-Lin Wu
Graduate Institute of Communication Engineering
National Taiwan University
Taipei, Republic of ChinaI.

INTRODUCTION

The pseudo-random bit sequence (PRBS) simulation takes
more and more time due to the increasing circuit complexity.
Several methods using superposition technique have then been
proposed to obtain the output waveform quickly. Thesemethods began with expressing an arbitrary input data b[n] as alinear combination of elementary functions (e.g., unit step) as:

x (t )   (b[n ]  b[n  1]) u(t  nUI )
n 0







  b[n ] (u(t  nUI )  u(t  ( n  1)UI )),n 0where x(t) is the input waveform, UI the unit interval, and u(t)the unit step function, assuming b[n]=0 for n<0. The next stepis to analytically express the output as a function of b[n].
Casper et al. [2], starting from the second line of (1), built theoutput using single bit response (SBR), or ‘010’ response, asy (t )   b[n]  (t  nUI ), 




n 0

where y(t) is the output waveform and φ(t) the SBR. Shi et al.
[3], starting from the first line of (1), constructed the output
waveform using one rising transition response, or ‘01’ response,and one falling transition response, or ‘10’ response, asy (t )   | b[n]  b[n  1] | s b[ n ] (t  nUI ), 



n 0

where sb[n] is the ‘01’ or ‘10’ response if b[n]=1 or 0. Ren et al.
[1] constructed the output waveform using multiple edge
responses (MER) asy(t )   | b[n]  b[n  1] | S b[ nm]...b[ n1] (t  nUI ),  n0

where S(t) is the particular rising or falling edge response whenthe m preceding bits are b[n-m] to b[n-1]. We term m the orderof the method. From (2~4), we can see the superposition
principle: the output signal can be constructed using only a few
simple responses.

Limitation exists, though. A completely linear system isassumed by (2), and (3) as well except that the ‘01’ and ‘10’
responses are allowed to be different. Only (4), theoretically, ispromised to work for arbitrary nonlinear system, as long as theorder m is chosen high enough.
II.

BRIEF INTRODUCTION TO THE ALGORITHM

A. Superposition of Multiple Edge Response
Definition: An edge response is defined as the output
response when the input undergoes a ‘01’ or ’10’ transition
with some particular preceding bit pattern. For a linear system,the ‘10’ response is exactly the inverse of the ‘01’ response,
and all ‘01’ or ‘10’ responses are always the same no matter
what happened previously. For a nonlinear system, neither of
the above properties is true. For example, the ‘01’ responses(the underlined portions) to the input ‘001’ and ‘101’ are in fact
different because a ‘10’ transition just before the target ‘01’
transition in the second input will alter the system greatly, e.g.
by causing a voltage ripples in the PDN, as illustrated in Fig. 1.
For such systems, errors occur if one uses (2) or (3) to build theoutput. The MER method cleverly solves this problem by using
more than two edge responses to build the output. Thus, thenonlinear effects such as the voltage ripples in PDN due toearlier transitions can be stored in these edge responses.
Example: one can choose to use 4 edge responses, as
defined in Table I, to construct the output waveform. In thistable, yabc(t) denotes the output signal to the input bit pattern‘abc’. These 4 responses together correspond to a 2nd-order
MER since all possible input bit patterns, up to 2 bits before thetarget ‘01’ or ‘10’ transition, i.e., the last transition, areconsidered. As an example, to construct the output due to input‘00101101’ using 2nd order MER, according to (4), we havey(t)=r 001(t)+f 010(t−UI)+r 101(t−2UI)+f 110(t−4UI)+r 101(t−5UI).
In general, an m’th order MER stores a total of 2m edge
responses. Edge responses up to order 5 are listed in Table II,
where r abc01(t) = yabc01(t) – yabc00(t).
The reason that the MER method can include nonlinear effects
is that: Nonlinearity decays to zero in the end. For example,
voltage ripples on the PDN due to earlier transitions cannot lastforever. The PDN voltage will saturate in the end. It can be
expected that rising transitions r 101(t) and r 001(t) are different,
but r10000001(t) and r001(t) shall be quite similar. Thus by storingenough amounts of edge responses, nonlinear effects can bewell captured. This is the basic idea of MER.

Voltage (V)

1.5

1.0

'001'

0.5

'101' − '100'
0.0
0

0.1

0.2

Time (ns)0.3

Figure 1. Different ‘01’ transitions.

DEFINITION OF 2ND ORDER BP RESPONSESTABLE I.
Name

Definition

Name

Definitionr 001(t)y 001(t)f 110(t)y 110(t)r 101(t)y 101(t) − y 100(t)f 010(t)y 010(t) − y 011(t)Figure 2. The test circuit consists of 4 CMOS inverters as an output driverand 2 T-model RLC power distribution networks.
1.8
1.6

Order

Rising Transitions a

Falling Transitions

1

r 01

f 10

2

101

010

r

f

Voltage (V)BP RESPONSES UP TO 5TH ORDERTABLE II.

1.4
1.2
1.0

r01

0.8

r101

0.6

r0101

0.4

r101010.2

3

r 1001, r 0101

f 0110, f 1010

4

r 10001, r 01001, r 01101, r 10101

f 01110, f 10110, f 10010, f 01010

5

r 100001, r 010001, r 011001, r 101001
r 011101, r 101101, r 100101, r 010101

f 011110, f 101110, f 100110, f 010110
f 100010, f 010010, f 011010, f 101010

0.0
0

0.1

Time (ns)

0.2

0.3

Figure 3. Four BP responses of the test circuit. The difference between
subsequent waveforms decreases when the order increases.

a. All ‘(t)’s are omitted for simplicity.

III.

TABLE III.

TEST RESULTS

A. Test Circuit
The test circuit is a 4-stage CMOS inverter with power/ground nets connected to PDN before to the ideal power/ground as shown in Fig. 2. The PDN is an RLC T model andthe values are typical of a package level PDN (see, for example,[4]). Board level PDN is omitted deliberately to make the PDN
worse. The input data rate is 10 Gb/s with rise/fall time 10 ps.
B. Difference between Edge Responses
The edge responses r01, r101, r0101, and r10101 shown in Fig. 3
illustrate that the differences between subsequent waveforms
become smaller when the order increases, which follows thebasic idea of MER. A quantitative comparison for risingtransitions up to order 5 is listed in Table III, where thedifference DRxxx,yyy is defined as the absolute area between
rising transitions rxxx, and ryyy as

DR0101,10101=6.2

DR xxx, yyy   r xxx  r yyy dt. 

DR0101,100101=4.8

DR01,101= 24.5a
DR101,01101=5.8

DR01101,101101=5.4

DR101,011101=5.6
DR1001,01001=8.0

DR01001,101001=5.7

DR01,1001=10.1
DR1001,011001=5.1
DR01,10001=5.0

DR10001,010001=6.3

DR01,100001=4.6
a. The units are all ps×V.

REFERENCES
[1]


[3]

Table III is organized in a manner that each row itself makes
up a comparable series. For example, in the first row, DR01,101=24.5 quantifies the impact of the ‘10’ transition 1 bit earlier,and DR101,0101 = 8.7 indicates that the influence of the ‘01’
transition 2 bits earlier is smaller. DR0101,10101 and DR10101,010101
are even smaller, which conforms to the basic idea of MER.
Most of the rows follow this trend.

DR10101,010101=6.0

DR101,0101=8.7[2]COMPARISON OF RISING TRANSITIONS[4]

J. Ren and D. Oh, “Multiple edge responses for fast and accurate systemsimulations,” IEEE Trans. Adv. Packag., vol.31, no. 4, pp. 741-748, Nov.
2008.
B.K. Casper, M. Haycock, and R. Mooney, “An accurate and efficient
analysis method for multi-Gb/s chip-to-chip signaling schemes,” inIEEE Symp. VLSI Circuits, Jun. 2002, pp. 54-57.
R. Shi, W. Yu, Y. Zhu, E.S. Kuh, and C.K. Cheng, “Efficient andaccurate eye diagram prediction for high speed signaling,” in Proc.
IEEE Int. Conf. Computer Aided Design, San Jose, CA, Nov 2008, pp.
655-661
K. Iokibe, Y. Yano, and Y. Toyota, “Insertion of parallel RL circuits
into power distribution network for simultaneous switching current
reduction and power integrity,” in Proc. IEEE Asia-Pacific Internation
Symposium on Electromagnetic Compatibility, Sentosa, Singapore, May
2012, pp. 417-420

Macromodeling and Circuit Simulation of High-Speed Interconnects
Based on Vector Fitting and Equivalent Circuit Synthesis
Daisuke Honda†,1 , Tadatoshi Sekine‡,2 , and Hideki Asai†,‡,3†

Dept. of Systems Eng., Graduate School of Eng., Shizuoka University‡
Dept. of Systems Eng., Shizuoka University
3–5–1 Johoku, Naka-ku, Hamamatsu-shi, 432-8561, Japan
Phone/FAX:+81-53-478-1237/+81-53-478-1269
E-mail: { 1 dicek, 2 sekine}@tzasai7.sys.eng.shizuoka.ac.jp, 3 hideasai@sys.eng.shizuoka.ac.jp

Abstract

nents interconnected with each other, the S parameters arewidely used because they are useful to represent the relaThis paper describes a way to approximate the frequency re- tions among systems from a viewpoint of the characteristic
sponse of interconnects by using a rational function and syn- impedance.
thesize a SPICE-compatible equivalent circuit from the funcIn This paper, we describe a way to include the discretetion. We introduce the vector fitting (VF) to derive the ra- data of the frequency response to circuit simulations [4–6].
tional function in the pole-residue form from the dataset of The way that we introduce consists of three procedures: a
the frequency response. Then, the state equation which has rational function approximation of the frequency response;
only a real number is realized using the Jordan canonical form state-space realization; and a SPICE-compatible equivalentand proper similarity transformation. The SPICE-compatible circuit synthesis.
circuit is synthesized from the state equation with dependent
The remainder part of the paper is organized as follows. In
sources in addition to RC lumped elements. We provide the Section II, we introduce the vector fitting (VF) to derive thenumerical results for an example system which represents a rational function approximation of the frequency response.
simplified transmission line and discuss about the availability Then, in Section III, the state equation is derived from theof the VF and equivalent circuit realization technique.
rational function. Subsequently, The equivalent circuit associated with the state equation is obtained in Section IV. The
availability of the above techniques is demonstrated by show1. Introduction
ing some numerical results in Section V. Finally, we conclude
the paper in Section VI.
One of the problems in the recent advanced integration
technologies is that the space to which electronic circuits aremounted is limited to a small area. In this situation, inter- 2. Rational Approximation Using Vector Fitting (VF)ferences among them become significant, especially, those
We adopt VF to derive the rational function in the poleinside mobile and in-vehicle electronic devices. Therefore,
residueform from the dataset of the frequency response beit becomes a challenging issue to ensure signal integrity (SI)
cause
the
function obtained by VF is accurate and alwaysof the electronic circuits because of the signal delay, reflection, skin effect, and cross talk caused by the high-frequency stable [7, 8]. In VF, it is assumed that a frequency responsecomponents of the signal [1]. In order to lead SI-aware de- function F (s) is approximated as the rational function in thesign, it is necessary to consider not only the electrical be- pole-residue form
havior of the circuit but also the electromagnetic characterisN∑
rntics of surrounding components. For this purpose, it becomes+d(1)
F (s) ≈
more and more important to develop the modeling and sims − pnn=1
ulation techniques for high-speed interconnects in the nextgeneration circuit design [2, 3].
where N is the order of the rational function, s = jω, j is
In general, the characteristics of the high-speed intercon- the imaginary unit, ω is the angular frequency, and d is thenects can be represented by frequency responses, such as the real constant while the poles pn and residues rn are either
Z, Y, and S-parameters, which tend to be provided in the form real or complex conjugate pairs. The VF algorithm separates
of discrete data from measurements or full-wave simulations. the calculations of the poles, residues, and d in (1) into two
In particular, in the case that we deal with electronic compo- processes: the poles are obtained in the first process, and theresidues and d are calculated in the second one. To identify
the poles, we first define an unknown weighting function σ(s)and a weighted function σ(s)F (s) as rational functions asfollows:σ(s)σ(s)F (s)

=

=

N
∑

r̃n
+1
s
−
p̃n
n=1
N
∑

rn
+d
s
−
p̃n
n=1

(2)

(3)

where p̃n and r̃n are the poles and residues of σ(s). Note thatthe poles of the weighted function are the same as those ofσ(s), and its residues are the same as those in (1). Substituting (2) into (3) and rearranging the equation lead to
)
(N
N
∑
∑ r̃n
rn
+d−
F (s) = F (s).
(4)
s − p̃ns − p̃nn=1
n=1

where pr and rr are the real poles and residues, pc and rc are
the complex poles and residues, ∗ denotes the complex conjugate, Nr is the number of the real poles, Nc is the number ofthe complex conjugate pairs, and N = Nr + 2Nc . Assuming
each element of the transfer function matrix
H(s) = C(sI − A)−1 B + D(8)has the form of (7), the associated state equation which has
M state variables and P ports can be written as follows:{ẋ = Ax + Bu
y = Cx + Du

(9)where x is the state variable vector, u and y are the input andoutput vectors, A is M × M diagonal matrix containing thepoles, B is M × P matrix containing 0 or 1, C is P × M
matrix containing the residues, D is P × P matrix containing the real constants. More specifically, (9) is written as theGiven p̃n as known initial poles and F (s) at K frequency
Jordan canonical form
points, we get K linear equations in the form of (4). As sembling those equations leads to the overdetermined linear   ẋ   A
0
0
xr
Br

r
r

 
problemẋc,1  =  0 Ac 0  xc,1  +  Bc  u



ẋc,2
0
0  A∗c  xc,2B∗cΦz = f
(5)
(10)


[
] xr

∗



where Φ is K × (2N + 1) matrix, z is the unknown vector  y = Cr Cc Cc
xc,1 + Duof which the size is 2N + 1, and f is the vector of which the xc,2size is K. Since the unknown vector is composed of rn , d,and r̃n , we can obtain r̃n by solving (5). In fact, the poles pn where the submatrices and vectors with the subscript r areof (1) are equal to the zeros of σ(s), and these zeros can be related to the real number, and those with the subscript c arefound by calculating the eigenvalues of a matrix
related to the complex number. In order to make (10) have
only real numbers, we perform a similarity transformation.
Q = Â − b̂ĉ
(6)
The transformation is introduced by replacing x by J−1 x̃ andwhere Â is the diagonal matrix of which the elements are premultiplying J by the first equation of (10), where J is ap̃n , b̂ is the column vector of ones, and ĉ is the row vector transformation matrix
containing r̃n . The eigenvalues of (6) provide a new set ofI 00poles, and the new poles are used as p̃n in (4) to construct the
I 
J= 0 I
(11)
overdetermined problem (5) again. By iterating solving (5)0jI−jI
and calculating the eigenvalues of (6), the new set of poles
converges to pn . Once pn is obtained, the accurate values of
the residues rn and d can be calculated by solving the other and I is the unit matrix. As a result, we get the real state
overdetermined problem, which is derived from (1) in a sim- equation composed only of real numbers
 
 


ilar manner to the way to derive (5) from (4).
x̃˙ r
Ar
0
0
x̃r



  x̃˙ c,1  =  0
Re(Ac ) Im(Ac )  x̃c,1 
3. State-Space Realization
 x̃˙
0
−Re(A
x̃c,2

c ) Re(Ac )
c,2






B
r
The rational function (1) obtained by VF can be rewritten(12)+  2Re(Bc )  u
in detail as follows:
−2Im(B
)
c




Nr

∑

rr,n
x̃r

[
]

F (s) =

y = Cr Re(Cc ) Im(Cc )  x̃c,1  + Dus − pr,n


n=1
x̃c,2
)
Nc (
∗
∑
rc,n
rc,n
+
+
+ d (7)s − pc,ns − p∗c,n
where x̃ = [x̃Tr x̃Tc,1 x̃Tc,2 ]T .
n=1

Port 1 i1

Z0

Z0

v1

V = 2b1 Z 0
v + Z 0i1
V= 1
2 Z0

b1

i2 Port 2v2

V = 2b2 Z 0
V=

v2 + Z 0i2
2 Z0

b2

a1

Port 1

1Ω

1 nH

1Ω

1 nH

0.1 nF

Vin

Port 20.01 pF

a2

c11 x1 c12 x2 d12 a2
1
R=
d11

d 21a1 c22 x2 c21 x1
R=

1
d 22

x1

x2
C =1

a12 x2 b11a1 b12 a2

C =1

b22 a2 b21a1 a21 x1

1
R=
a11

1
R=
a22

Fig. 1: The equivalent circuit realization of the S-parameters
4. Equivalent Circuit Synthesis

where ap and bp are the incident and reflected waves at theport p, and they are related to the port voltage vp and current
ip as
vp + Z0 ip
√
,
2 Z0

bp =

of which the transconductances are a12 , b11 , and b12 , respectively, are connected in parallel to each other. Fig. 1 shows
the schematic of the circuit associated with (13) and (14). Ina similar manner, we can apply the above synthesis technique
to arbitrary state space equations.
5. Numerical Results

It has been reported that the circuit synthesized using only
passive RLC elements may be nonpassive because the values of the elements may become negative [9, 10]. Therefore,
we adopt another technique, which uses control sources in
addition to RC elements and can always synthesize the network that is able to be analyzed by using SPICE-like simulators [4, 11]. To illustrate this technique, we use the simple
state equation of which the transfer function represents the Sparameter. The example system has two ports and two statevariables as [ ] [][ ] [][ ]ẋ1
a11 a12 x1
b
b
a1


=
+ 11 12

ẋ
a
a
x
b
b
2
21
22
2
21
22
[ ] [
][ ] [] [ a2 ] (13)
b
c
c
x
d
d
a1

1
11
12
1
11
12

=
+

b2
c21 c22 x2
d21 d22a2

ap =Fig. 2: The example circuit

vp − Z0 ip
√
2 Z0

(14)

As an example application of the above approaches, wedeal with the circuit illustrated in Fig. 2. In Fig. 2, the partsurrounded by the dotted line represents the equivalent circuit
model of a simple transmission line. We regard this part as thetwo-port network to which the voltage source Vin is connected
at Port 1 and the capacitor of 0.01 pF is connected at Port 2.
The values of the other elements are shown in Fig. 2.
First, the S-parameters of the two-port network and theirrational approximation obtained by VF are plotted in Fig. 3.
In this case, we provide three rational functions of which theorders N are 2, 3, and 6, respectively. As shown in Fig. 3, theS-parameters calculated from the functions of which the orders are 3 and 6 completely agree with the analytical solution
in the given frequency range while the results from the lower
order function is slightly different from the others.
Next, we perform the transient simulations using LTspicefor the original circuit in Fig. 2 and the equivalent circuits
synthesized from the rational functions. For this purpose, thewaveform of Vin is set to the trapezoidal pulse of which therise and fall times are 0.5 ns and the pulse width is 10 ns.
The transient responses observed at Port 2 are plotted in Fig.
4. As we can see in Fig. 4, the waveform from the function
with N = 3 agrees well with one from the original circuit. In
contrast, the transient response from the function with N = 6
does not completely agree with one from the original circuit
though the rational function is accurate in the frequency domain. This may be because the higher order function has anerror at the high frequency range which does not shown in
Fig. 3. The waveform from the lower order function also does
not agree with the original one as we expected.

where Z0 is the reference impedance. In order to solve the ordinary differential equations (ODEs) (13) along with (14) byusing the SPICE-like simulator, we assume that each ODE
in (13) fulfills Kirchhoff’s current law (KCL), and the equations in (14) are derived from Kirchhoff’s voltage law (KVL).
This can be achieved by regarding the state variable and ap as
the voltages, the coefficient related to the state variable as resistance or the transconductance of a voltage control current
source (VCCS), and the coefficient related to the derivativeof the state variable as capacitance. For instance, by assuming the first equation of (13) is the KCL equation related to 6. Conclusion
the node with the voltage x1 , we can synthesize the circuit in
In this paper, we have introduce one of the ways to include
which the capacitor of 1 F, resistor of 1/a11 Ω, and VCCSs,

References

1
S11,S22[1] C. R. Paul, Analysis of Multiconductor Transmission
Lines, 2nd ed. Hoboken, NJ: Wiley-IEEE Press, Oct.
2007.

Magnitude

AnalyticalFitted (N = 2)
Fitted (N = 3)
Fitted (N = 6)[2] R. Achar and M. S. Nakhla, “Simulation of high-speed
interconnects,” in Proc. IEEE, vol. 89, May 2001, pp.
693–728.

0
S12,S210

1[3] A. E. Ruehli and A. C. Cagellaris, “Progress in themethodologies for the electrical modeling of interconnects and electronic packages,” in Proc. IEEE, vol. 89,
May 2001, pp. 740–771.

2

Frequency [GHz]Fig. 3: The S-parameters of the two-port network in Fig. 2

Voltage[V]

5

0

0

[5] S. Caniggia and F. Maradei, “Circuit and numericalmodeling of electrostatic discharge generators,” IEEE
Trans. Ind. Appl., vol. 42, no. 6, pp. 1350–1357, Nov.
2006.

Original
Synthesized (N = 2)
Synthesized (N = 3)
Synthesized (N = 6)
10Time [ns]Fig. 4: The transient response at Port 2[4] E.-P. Li, E.-X. Liu, L.-W. Li, and M.-S. Leong, “A
coupled efficient and systematic full-wave time-domain
macromodeling and circuit simulation method for signal
integrity analysis of high-speed interconnects,” IEEE
Trans. Adv. Packag., vol. 27, no. 1, pp. 213–223, Feb.
2004.

20

[6] T. Sekine, H. Asai, and J. S. Lee, “Unified circuit modeling technique for the simulation of electrostatic discharge (ESD) injected by an ESD generator,” in IEEE
EMC Symposium 2012, Pittsburgh, PA, Aug. 2012, pp.
340–345.
[7] B. Gustavsen, “Rational approximation of frequency
domain responses by vector fitting,” IEEE Trans. Power
Del., vol. 14, no. 3, pp. 1052–1061, Jul. 1999.

[8] ——, The Vector Fitting Web Site. [Online]. Available:
http://www.energy.sintef.no/Produkt/VECTFIT/
the discrete data of the frequency response to the circuit simulation. The VF algorithm was adopted to derive the ratio[9] A. Scarlatti and C. L. Holloway, “An equivalent
nal function approximation in the pole-residue form, and thetransmission-line model containing dispersion for highSPICE-compatible equivalent circuit was obtained from the
speed digital lines–with an FDTD implementation,”state equation associated with the rational function. From theIEEE Trans. Electromagn. Compat., vol. 43, no. 4, pp.
numerical results, it has been confirmed that the S-parameters504–514, Nov. 2001.
of the simple transmission line can be approximated accurately by using VF. Moreover, the equivalent circuit derived [10] T.-L. Wu, C.-C. Kuo, H.-C. Chang, and J.-S. Hsieh,“A novel systematic approach for equivalent model exusing the circuit synthesis technique could be applied to thetraction of embedded high-speed interconnects in timetransient simulation performed by the conventional SPICEdomain,” IEEE Trans. Electromagn. Compat., vol. 45,
based simulator. We saw that the order of the rational funcno. 3, pp. 493–501, Aug. 2003.
tion should not be too low or too high in order to produce theaccurate results in the time domain.
[11] R. Neumayer, A. Stelzer, F. Haslinger, and R. Weigel,“On the synthesis of equivalent-circuit models for multiports characterized by frequency-dependent parameAcknowledgment
ters,” IEEE Trans. Microw. Theory Tech., vol. 50, no. 12,
This work is partially supported by Japan Society for the Propp. 2789–2796, Dec. 2002.
motion of Science (JSPS).

A Compact and Broadband Forward-Wave Directional Coupler with
Arbitrary Coupling Levels
Tong-Hong Lin, Tzong-Lin Wu
Department of Electrical Engineering
Graduate Institute of Communication Engineering
National Taiwan University
Taipei 10617, Taiwan
wtl@cc.ee.ntu.edu.tw
Abstract — A novel broadband and compact forward-wave
directional coupler is proposed. The broadband impedance
matching of the even- and odd-mode of the proposed coupler canbe achieved, and the phase difference between the even- and oddmode can keep approximately at a constant level during a broadfrequency range. The equivalent models are constructed to
predict the propagation properties based on the Bloch-Floquet
theorem. This three-layer structure is fabricated on the Duroid
substrate. The coupling level is -3.8 ± 0.5 dB from 1.97 GHz to
2.85 GHz and the corresponding fractional bandwidth is 37%.
Besides, the size is only 0.4 λg. Compared to the previous works,
the proposed structure can be designed with wider operational
bandwidth and a smaller size.
Index Terms – forward-wave directional coupler. broadband.
(a)I. INTRODUCTION
Forward-wave directional coupler (FWDC) is an importantcomponent in communication system and is often used for
power dividing and power combining. It is often applied in
the circuit design for performance enhancement, such as filter,
amplifiers, mixers, and modulators [1]. The coupling level ofthe forward coupler is decided by the phase difference (PD)between the even and odd mode while the input impedancesof the two modes are matched perfectly [2]. Therefore themain difficulties in realizing a broadband FWDC are: (i)achieving broadband input impedance matching of the evenand odd mode and (ii) keeping the PDs between the even andodd mode at a constant value during a broad frequency range.
Several couplers have been proposed, such as the periodic
shunt stubs along the coupled-line [3], [4], the microstrip
coupled-line with pattern ground structure [2], [5]. They allcontribute to the enhancement of PDs, but they are only usedin narrow bandwidth applications because of theirmonotonously increased PDs. Asymmetrical coupler isanother method to enhance the operational bandwidth, but theeven- and odd-mode decomposition can not be used toanalyze it [6], [7]. Therefore the analysis of asymmetricalcoupler becomes more complicated.
In this paper, a broadband and symmetrical forward
coupler is proposed. The even- and odd-mode Bloch
impedances of the proposed FWDC can be designed to thevalues which are very close to 50 Ω. Therefore, thebroadband impedance matching of the two modes can beachieved. The PDs between two modes also keeps(b)Fig. 1. (a) Unit cell of the proposed forward-wave directional coupler. (b)The front view.

approximately at a constant level within a broad frequency
range. Due to these two characteristics, a FWDC can be
designed with a broad bandwidth.
II. DESIGN CONCEPT AND THEORETICAL MODELINGThe proposed FWDC is designed based on the concept of
periodic structures, and the configuration of the proposed unit
cell is shown in Fig. 1. The middle layer in Fig. 1(a) istransparent to some degree so that the bottom layer withlighter color can be seen. As shown in Fig. 1, it is a threelayer structure and the middle layer serves as the reference
plane. There are a transmission line and a stub at the center of
the transmission line on the top layer. Besides, the pattern on
the bottom layer is the same as the top layer, and they areconnected with each other by a via at the end of the stubs. The
width and length of the transmission line can be represented
as w and l, respectively. The width of the shunt stub is t andthe length is s. The thickness between each layer isrepresented as h and the diameter of the connecting via is d.

TABLE IGEOMETRICAL DIMENSIONS AND THE CORRESPONDING PARAMETERS
Geometricaldimensions (mm)l = 18e = 3.4
w = 1.8
f = 4.2
t = 0.16
h = 0.813
s = 9.7
d = 0.28

Parameters of theequivalent models (at 2.5 GHz)
ZM = 50 Ω
ZStub = 133 Ω
θM = 45oθS = 45o

ZM θ M

ZStub

ZM θM

θS
(a)
(a)

ZM θ M

ZStub

ZM θM

θS

(b)
Fig. 2. Equivalent models of a unit cell (a) Odd mode. (b) Even mode.

The rectangular hole which etches on the middle layer toensure that the connecting via did not short to the middle
layer is e in width and f in length. The geometrical dimensionsof a unit cell are summarized in Table I.
Since the proposed structure is a symmetric structure, theeven- and odd- mode decomposition can be used to synthesize
it. The perfect electric conductor (PEC) plane and the perfectmagnetic conductor (PMC) plane are placed on the plane ofsymmetry, the middle layer, for cases of the odd mode and theeven mode, respectively. The equivalent models of bothmodes are shown in Fig. 2. As shown in Fig. 2(a), the oddmode structure can be modeled as two sections of microstrip
lines with characteristic impedances ZM and electrical lengths
θM and a short stub at the center. The stub is shorted to thereference plane because a PEC plane is placed on the middle
layer. The corresponding characteristic impedance andelectrical length of the short stub are ZStub and θS, respectively.
The even-mode equivalent model is shown in Fig. 2(b). It issimilar to the equivalent model of the odd mode except forthe stub at the center. The stub is now an open stub because a
PMC plane is placed on the middle layer. As shown in Fig.
2(b), the parameters using to represent the characteristic
impedances and electrical lengths are the same as the odd
mode. The extracted parameters are also summarized in TableI.
The even- and odd- mode electrical lengths and Bloch(b)
Fig. 3. Simulation results for the (a) even- and odd-mode electrical lengths
and (b) even- and odd-mode Bloch impedances.

impedances of a unit cell are important for designing aperiodic structure and they can be derived from the equivalentmodels by the Bloch-Floquet theorem [1]. The results are
given as



r
2



r
2



odd  cos1 cos  2 M   z sin  2 M  cot  S even  cos 1 cos  2 M   z sin  2 M  tan  S Z odd Z even 

jZ M sin  M  2cos  M  rz sin  M cot  S A B
jZ M sin  M  2cos  M  rz sin  M tan  S CDA    rz 2  cot  S   1 sin 2 (2 M )2

B  rz sin  2M  cos(2M ) cot SC    rz 2  tan  S   1 sin 2 (2 M )D  rz sin  2M  cos(2M ) tan  S
2

(4a)
(4b)
(4c)
(4d)
(4e)
(4f)
(4g)
(4h)

Fig. 4. Photographs of the fabricated 3-dB FWDC. (εr = 3.55)

range. The phase of the odd mode keeps at 0 degree from 1
GHz to 1.3 GHz for the same reason and starts to increase
from the cut-off frequency 1.3 GHz. It also increases fasternear the cut-off frequency and then turn to increases linearly
from 1.7 GHz to 4 GHz. The coupling level of a FWDC isbased on the PD of the two modes, so the PD is also plotted inFig. 3(a). As shown in Fig. 3(a), the even- and odd- mode
phases are parallel to each other from 1.7 GHz to 3.3 GHz.
Therefore, the PD between the even and odd modes keeps at
22.5o within this frequency range. As demonstrated in Fig.
3(b), the even-mode Bloch impedance decreases graduallyfrom 45 Ω to 0 Ω and there is a flat region from 1 GHz to 3.1
GHz where the impedance keeps approximately at 45 Ω. The
odd-mode impedance also decreases gradually from 100 Ω to55 Ω and keeps approximately at 57 Ω within a flat region
from 1.9 GHz to 4 GHz. Because of these flat regions, theproposed structure can achieve good broadband matching.
Recalling that two conditions to achieve a broadband FWDC
are broadband constant PDs between the even and odd modes
and broadband impedance matching. As demonstrated in Fig.
3, both conditions are achieved. As a result, a broadband
FWDC can be expected.

(a)III. MEASURED RESULTS(b)
Fig. 5. Measurement and full-wave simulation of the proposed 3-dB FWDC(a) S parameters and (b) phase imbalance.

where rz is the impedance ratio and equals to ZM/ZStub. The
frequency dependence of  odd ,  even , Z odd , Z even according to(4a)-(4d) are shown in Fig. 3. The full-wave simulation is also
included for validation and they show good agreement. The
values of the parameters of the equivalent models and thegeometric dimensions used in the full-wave simulation can befound in Table I. As shown in Fig. 3(a), the phase of the evenmode increases linearly with the frequency from 1 GHz to
3.3GHz and then increases faster near the cut-off frequency
3.7 GHz. The phase of the even mode keeps at 180 degreefrom 3.7 GHz to 4 GHz because the even-mode Bloch
impedance becomes pure imaginary and corresponds toevanescent rather than propagating within this frequency

A test sample of the proposed broadband 3-dB FWDC withthe dimension shown in Table I is fabricated using Ro4003.
The relative permittivity and loss tangent are 3.55 and 0.0027,
respectively. The photographs and the port names are shownin Fig. 4. As shown in Fig. 4, the test sample is miniaturized
by meandering the stubs and the total length is 28.51 mm (0.4λg) excluding the additional structure for measurement. The
measured and full-wave simulated scattering parameters areshown in Fig. 5(a), and they show good agreement. As shownin Fig. 5(a), the return loss and isolation are below -15 dB
from 1.94 GHz to 3.28 GHz. Therefore, broadband matchingis achieved. Moreover, the measured S21 and S31 are -3.8 ± 0.5
dB from 1.97 GHz to 2.85 GHz. The measured and simulated
phase imbalance between port2 and port 3 are demonstrated
in Fig. 5(b). They show good agreement except for the higher
frequencies. This is because the glue used to connect between
layers is FR4, the effect of loss become significant at higher
frequencies. As shown in Fig. 5(b), the measured phaseimbalance stays within -90o ± 2o from 1.78 GHz to 2.92 GHz.
IV. CONCLUSION
A novel broadband and compact FWDC is proposed. The
broadband impedance matching of both the even- and oddmode can be achieved. Besides, the phase difference between
the even and odd modes can keep at a constant level within abroad frequency range. Therefore, a broadband FWDC can be
designed. The measured S21 and S31 are -3.8 ± 0.5 dB from1.97 GHz to 2.85 GHz and the corresponding fractional
bandwidth is 37%.

REFERENCES[1] D. M. Pozar, Microwave Engineering, New York: Wiley, 2004,
pp. 333-364.
[2] S. K. Hsu, J. C. Yen, and T. L. Wu, “A novel compact forwardwave directional coupler design using periodical patterned
ground structure (Periodical style—Accepted for publication),”IEEE Trans. Microw. Theory Tech., vol. 59, no. 5, pp. 12491257, May. 2011.
[3] T. Fujii, and I. Ohta, “Size-reduction of coupled-microstrip 3
dB forward couplers by loading with periodic shunt capacitive
stubs,” in IEEE MTT-S Int.Microw.Symp. Dig., Jun 2005,
pp.1235-1238.
[4] J. H. Park and Y. Lee, “Improved capacitive loading method
for miniaturization of 0 dB forward-wave directional couplers(Periodical style—Accepted for publication),” IEEE Microw.
Wireless Compon. Lett., to be published.
[5] C. C. Chang, Y. Qian, T. Itoh, ”Enhanced forward coupling
phenomena between microstrip lines on periodically patterned
ground plane,” in IEEE MTT-S Int.Microw.Symp. Dig., May
2001, pp. 2039-2042.
[6] S. Uysal, C. W. Turner, and J. Watkins, “Nonuniform
transmission line codirectional couplers for hybrid MIMIC andsuperconductive applications,” IEEE Trans. Microw. Theory
Tech., vol. 42, no. 3, pp. 407-414, Mar. 1994.
[7] S. Uysal, and J. Watkins, ”Novel microstrip multifunction
directional couplers and filters for microwave and millimeterwave applications,” IEEE Trans. Microw. Theory Tech., vol. 39,
no. 6, pp. 977–985, Jun. 1991.
[8] C. Caloz and T. Itoh, Electromagnetic Metamaterials :Transmission Line Theory and Microwave Applications. New
York: Wiley, 2006, pp. 124-127.
[9] W. R. Eisenstadt and Y. Eo, “S-parameter-based ICinterconnect transmission line characterization,” IEEE Trans.
Compon., Hybrids, Manuf. Technol., vol. 15, no. 4, pp. 483490, Aug. 1992.
[10] R. Mongia, I. Bahl, P. Bhartia, and J. Hong, RF andMicrowave Coupled-Line Circuits, Boston / London : Artech
House, 2007.