Abstract—This paper introduces an improved electrostatic discharge (ESD) system-level transient simulation modeling method and discusses its validation using IEC 61000-4-2 ESD pulses on a real-world product. The system model is composed of high current and broadband (up to 3 GHz) models of R, L, C, ferrite beads, diodes, and integrated circuit IO pins. A complex return path model is the key to correctly model the system’s response to the IEC excitation. The model includes energy-limited time-dependent IC damage models. A power–time integral method is introduced to accurately determine if a junction would experience thermal runaway under an arbitrary injection waveform. The proposed method does not require knowledge of the junction’s microscopic geometry, material information, defect location, or melting temperature.
Index Terms—Common mode, electromagnetic compatibility (EMC), electrostatic discharge (ESD), human metal model (HMM), IEC 61000-4-2, system efficient ESD design (SEED), transmission line pulser (TLP).
Recent studies have shown that system-level electrostatic discharge (ESD) simulation can serve as a powerful tool for analyzing ESD performance . The simulation enables the design of reliable protection on the first attempt and avoids the need for repeated design optimization tests.
The concept of ESD simulation has been promoted as an option in system-level ESD efficient design (SEED) . SEED emphasizes the analysis of the interaction between the quasistatic I–V curve of a vulnerable pin and the pin’s external protection. Gossner et al. applied SEED for analyzing an IO pin’s response to ESD for different on-board protection solutions . Monnereau et al. extended the modeling framework by adding trace and package models, and validated their method with an inverter circuit under a 100 ns transmission line pulser (TLP) excitation . Li et al. previously published a hard error analysis of a cellphone’s keyboard illumination circuit based on a 35 ns TLP source . Orr et. al used similar method to characterize IC pins .
Although the SEED simulation offers greatly improved system-level ESD design, some issues remain unresolved. A TLP-excited system simulation may not substitute IEC/HMM  excited analysis in certain cases. TLP-based simulation results may be valid when the damage is caused by the IEC’s second peak (residue portion), which has a long duration and can be mimicked by a TLP pulse. It does not reflect the consequences of the first few nanoseconds of an IEC excitation. In addition, as an IEC waveform passes through a complicated system, the resulting injection on a vulnerable part could be in an arbitrary shape and, thereby, break a fixed empirical HMM-IEC relation . A TLP source is not suitable for modeling soft error, near-field coupling or signal integrity (SI) problems caused by an ESD injection. Due to the above reasons, many researches have moved forward to perform ESD transient simulation under IEC/HMM excitations –.
It is difficult to convert a TLP-based simulation into an IEC setup directly by substituting the TLP model with an ESD gun model, based on a real product system. Compared to a TLPbased model, an IEC source-based setup requires more sophisticated modeling on the current return path in order to achieve an accurate circuit response under ESD tests. Furthermore, intensive use of flex-printed circuits (FPCs) for connecting multiple PCBs creates complex return paths. Among the recent publications that researched system-level IEC simulation, some showed less accurate results compared to measurement, especially at the very first nanosecond, e.g., . Some demonstrated excellent modeling results, but the investigated problems were only at the PCB level rather than the real product level due to the lack of complex return path structures, e.g.,  .
In addition to modeling the PCB-based and IC internal ESD protection structures, a failure criterion is needed. Using only a TLP-derived constant failure, current threshold  may be insufficient if this threshold is only surpassed for a few nanoseconds. This will be the case if the initial peak of the ESD current surpasses the threshold but the second peak remains below it. As Notermans et al. concluded in , “For a real system, dynamic failure must be taken into account as well.” Particularly, it will be shown later in this paper that a complex network could introduce an oscillatory current waveform inside the system, thereby making a constant current threshold inapplicable. The dynamic failure according to junction overheat was investigated by Wunsch and Bell , who characterized the failure model with the tested pulse–power relationships. Later, several researchers such as Yiquan Cao, applied the thermal failure model in ESD scenarios .
In the study presented in this paper, we modeled a cell phone circuit in realistic IEC testing scenarios. The state of the art of this paper includes the following four parts. First, typical components (R, L, C, ferrite beads, and semiconductor devices) under high-current and high-frequency excitations are modeled. Second, a detailed model of the complex return path inside the phone is presented. Finally, a time-dependent destruction model and power–time integral method is introduced to accurately determine if a junction would suffer thermal damage under an arbitrary injection waveform. The checking algorithm is an extension to Taska’s work .
The remainder of this paper is organized as follows. Section II describes the product under investigation. The test systems and methods for creating the model are introduced in Section III. The component models are shown in Section IV. Section IV presents the semiconductor’s failure model and discusses the development of the thermal runaway criterion of a junction under an arbitrary waveform. Section VI mainly discusses the ESD gun model and common-mode path modeling. Section VII shows the validation of the system-level model and the model’s application for hard error analysis.
II. SYSTEM UNDER INVESTIGATION
A vulnerable keypad backlight LED circuit in a smart phone, as shown in Fig. 1, was investigated. The driver IC controlled the LED’s brightness by varying the IO pin’s state. All component information will be kept confidential because of intellectual property constraints. ESD tests indicated that the LED was a sensitive zapping point. During product-level tests, airmode discharge sometimes struck through the aperture between the plastic buttons that covered the LED and coupled into the illumination circuits.
At first glance, the circuit’s behavior under ESD appeared somewhat complex for the following reasons: 1) L–Cpairs could cause resonance; 2) ferrite beads and capacitors may saturate or show nonlinear behavior under high-current injection; and 3) the keyboard PCB was connected to the main PCB through an FPC, which introduced a complex return path for the ESD current.
III. MODELING METHODOLOGY
A component model was created based on an RF model and a device model obtained under high current, as shown in Fig. 2. This combination ensures sufficient accuracy under IEC 61000- 4-2 or HMM excitations. Based on the 0.7–1 ns rise time and the response of nonlinear elements, a modeling bandwidth of 3 GHz was selected. Z-parameters were used to obtain the RF model.
The high-current I–V curves were extracted using a 15∼40 ns adjustable TLP pulse (see Fig. 3), which is long enough to extract a stabilized result yet avoid damaging the DUT. To control the parasitics of the test setup, inductances were minimized, e.g., a circular arrangement of five 10-Ω resistors was used to create a broadband 2-Ω current measurement shunt.
IV. COMPONENT MODELS
A. Semiconductor Devices
Similar measurements were used to determine the VI behavior of LEDs, Zener diodes, and IC pins. The only difference was that the IC was powered to ensure the same operating conditions as those encountered during system-level testing.
The Zener diode’s transient I–V curve appears in Fig. 4, as does a behavioral model developed by fitting this curve. Diode 11 defined the I–V characteristics of the Zener diode under negative pulses applied to its cathode; diode 10 and the switch (actually, a voltage controlled resistor) determined the positive I–V characteristics. Diode 9 was used as a unidirectional switch to separate the positive and negative pulse injections.
The capacitance of the Zener diode was measured using a vector-network analyzer (VNA). Due to its large value of 25 pF, it was determined that the diode would carry most of the current during the first nanoseconds of the ESD pulse.
The LED model (see Fig. 5) is based on a similar concept. It has two parts: the factory-provided SPICE model for nominal current conditions, and two voltage-controlled resistances to mimic the high-current I–V behavior. The factory model already included the capacitance, so no external RF model is needed here.
The IO pin on the driver IC was modeled as a three-terminal device. First, a TLP was used to obtain the power clamp of the Vcc network (Diode 3 in Fig. 6). Then, the high-side (DIODE1) and low-side (DIODE2) protection diodes of the IO pin were measured by applying positive and negative pulses to the IO pin, respectively. Finally, using a VNA, the values of the linear components (C, R, and C33) were derived. The 300-pF-power rail capacitor is a combination of junction, gate, and metallization capacitance. The system contains a large 2-μF on-board capacitor placed in parallel.
B. Capacitors, Ferrite Beads, and Inductors
The voltage across a capacitor may lead to sparking, capacitor breakdown, and a recoverable change in the capacitance value , . Fig. 7 shows the voltage and current of a 10-V-rated 10-nF X7R capacitor that was excited with a 15 ns 3 kV TLP. Although the charge current was constant, a nonlinear voltage increase occurred. This indicates that the capacitance decreased as the voltage increased. The capacitance variation over time, or C(t), can be calculated from the measured voltage and charging current waveform
The C–V behavior was approximated by an arc–tangent function (2) to account for this C–V behavior, although other researchers have shown that quartic functions can work equally well 
where A–D tune the model, as shown in Fig. 8. For this specific capacitor, the best match was achieved at A = 18, B = 2.2, C = 2.8, and D = 7.
Using the arc–tangent function, together with equivalentseries-resistance and equivalent-series-inductance obtained from measured Z-parameters, a complete capacitor model can be created in Agilent’s Advanced Design System .
Not all capacitors behave nonlinearly under ESD. The low dielectric constant of NP0 ceramic will show little or no nonlinearity; however, the low capacitance values achievable with small-package NP0 capacitors may spark over.
Similar to capacitors, ferrites may exhibit saturation or other nonlinear behavior under high-current conditions. The nonlinear inductance can be approximated using the following:
In certain cases, the additional high-frequency noise on the measured I(t) may cause dI/dt to change significantly, thereby interfering with the calculated L(t). To calculate the L(t), one could either perform low-pass filtering on the tested raw data, or use
V (t)dt I (t) to calculate.
The inductance–current relationship can be modeled by a nonlinear arc–tangent function, as used in capacitor modeling. Here, we used an alternative method, a quartic equation, for modeling
where I stands for the current flow through the nonlinear inductor; L0 is the initial/nominal inductance; and Lsat represents the saturated inductance. A = 2 and B = 1 for the specific ferrite we tested. Fig. 9 shows the modeled curve of a ferrite with an equivalent 60-nH inductance that can be saturated to 20 nH.
The complete model of the ferrite appears in Fig. 10. Besides the nonlinear inductance model (SDD1P), other linear models can express the effect of the capacitance and loss following Yu’s topology . These linear parts usually can be found in a device’s datasheet and can be checked by measuring the S-parameters. This model does not take hysteresis into account because the ferrite bead uses soft magnetic materials that exhibit no relevant hysteresis .
V. DYNAMIC DESTRUCTION THRESHOLD MODELING
A. Failure Power Models
To determine if a specific ESD will damage a device, its robustness threshold must be known. As discussed previously, a simple current threshold may not be sufficient; a dynamic threshold will better predict complex waveforms, such as an HMM discharge. Using a TLP with a varying pulse width, the damage threshold function (see Fig. 11) was created. The TLP current decreased as the pulse length increased, indicating that the device was energy limited.
Semiconductor devices under electrical over stress (EOS) have many microscopic failure mechanisms, e.g., surface breakdown around a junction and internal body breakdown through a junction. However, as Wunsch noted in , most failure mechanisms are linked primarily to the junction temperature. The widely used junction thermal model was developed by Wunsch and Bell , and later, Taska . Their thermal analysis yielded the failure power (P) per unit junction area (A) as a function of the rectangular pulse width (tp ):
where K1 , K2 , and K are design-specific parameters that relate to the junction material and conductivities. The resulting curve of (5) appears in Fig. 12.
The parameters K1 , K2 , and K may not always be derived explicitly from junction design because in many applications, the material information and junction geometries are not known. They can be determined, however, by fitting the measured curves, as shown in Figs. 13 and 14.
B. Failure Criteria
To determine device failure under time-varying waveform P(τ ) based on the knowledge of the TLP tested failure power/time relationship P0 (t), one can identify whether or not any portion in P(τ ) injected the same amount of energy as a certain destructive rectangular pulse.
This idea can be derived from heat transfer equation 
where T is the junction temperature, ρ is the density, Cp is the specific heat capacity, D is the thermal diffusivity, and q(t) is the heating rate per unit volume. The Green’s function, or the solution to this function, is
The Green’s function is known as the impulse response in both the time and spatial domains. As an injection source P(r
,τ ) heating a defect volume Δ, the temperature at an observation location r (the vulnerable point) at time t can be written as 
where T0 is the initial ambient temperature.
A rectangular pulse with an amplitude of P0 and a duration of tf can damage a semiconductor junction because the failure point temperature reaches the failure temperature Tc
If an arbitrary injection profile that starts at an arbitrary time τ0 can also generate the same amount of heat within a duration of tf
This arbitrary waveform can be considered destructive. Therefore, the heat contribution of this arbitrary waveform to its equivalent rectangular pulse can be related as
The rectangular pulse failure power P0 is a function of duration tf (the failure power–time model in Section V-A), so the failure criterion is written as
Note that the power–time integral must be performed in an assumed failure time span tf ; otherwise, the integral of heat transfer function G cannot be eliminated. This is intuitive; if the injected arbitrary wave’s energy reaches P0 (tf ) ∗ tf over a longer span than tf , the junction temperature may still be lower than Tc because more heat has dissipated.
an Tc because more heat has dissipated. Equation (12) allows a devices’ thermal failure to be evaluated without knowing its material, geometry, failure location, or melting temperature. Only its tested failure model P0 (tf ) and simulated time-varying power profile P(τ ) are needed. Equation (12) can be implemented with the following algorithm:
Equation (12) can be simplified further if τ0 = 0, or, if the highest power portion always occurs at the beginning of an injection (usually the case for an ESD event). The criterion, therefore, is simplified as
The interception point of the left and right sides of (13) stands for the failure time and destructive injection energy (but not the energy that heats the defect region).
Section VII contains examples of applying the failure criterion.
VI. SYSTEM-LEVEL SETUP, ESD GUN MODEL, AND COMMON-MODE MODELING
A. System-Level Test Setup and Modeling
A contact-mode discharge on the DUT setup is shown in Fig. 15. A cellphone’s battery charging cord, filtered with a ferrite, was connected to the cell phone’s USB port as part of the return path. The cord’s shielding at the other end was shorted to a large metal plane.
In such a test setup, the ESD current return path (commonmode path) and the ESD generator should be modeled in order to correctly calculate the ESD current within the circuitry under investigation.
B. System-Level Grounding Model
For the system test setup shown in Fig. 15, the connection between the cell phone’s ground (metal frame) and the main ground plate can be modeled as shown in Fig. 16. The transmission lines TL1 and TL2 modeled the IO and Vcc nets on the double-sided flex circuit, respectively. The characteristic impedance was measured as 45Ω with a TDR. This impedance can also be calculated from the flex’s 2-D cross-sectional geometry. TL1 and TL2 were not referenced to the same metal; instead, their left sides were connected to the keyboard PCB’s local ground, and their right sides were shorted to the main PCB’s reference plane.
The transmission line TL3 modeled the flex’s ground metal relative to the cellphone’s body frame metal. The characteristic impedance of this common-mode path was measured as 120Ω.
C. ESD Gun Model
An ESD generator, TESEQ NSG 438 , was used in this project. Its equivalent circuit model appears in Fig. 17, which was developed based on Wang’s topology .
VII. SYSTEM-LEVEL SIMULATION RESULTS
A. System Model Validation
The system model was constructed by inserting all of the circuit models developed as described in Section IV, as well as the ESD gun model, into the system scheme shown in Fig. 16.
Fig. 18 shows one of our most challenging validation setups used for checking the model’s credibility and the robustness of the modeling methods. A Tektronix CT-6 probe was inserted in front of the IO pin to measure the ESD current flowing into the IC. To allow the current probe to be placed, an 8-mm-long wire was soldered in-series to the IO pin. This wire introduced an additional 4-nH inductance. The simulated current conformed to measurements reasonably well (see Fig. 19). The difference between simulated and measured results can be quantitatively described with the feature selective validation technique – .
B. Application of the System Model for ESD Hard-Error Analysis
1) Transient Current Flows Into the LED: One objective was to determine the conditions under which the LED would suffer damage. Calculating the destruction criteria (12) on the simulated power profile and the LED’s failure model, respectively, showed that under +14 kV, the LED would be damaged (see Fig. 20), which agreed with our tested result. The checking algorithm also showed that under 15-kV injection, the damage would occur within the first 5 ns; under 14 kV, the damage occurred at 32 ns. Fig. 21 shows the result of the simplified checking algorithm (13).
2) Thermal Failure of the Driver IC: Another objective was to analyze the conditions under which the driver IC could survive without any external protection (same setup as shown in Fig. 18, but with a 10-nF nonlinear capacitor in parallel to the LED to avoid LED destruction).
By applying (12), we were able to predict that the driver IC could survive under 15 kV but would not withstand a 16-kV injection (see Fig. 22). This prediction also agreed with our tested results.
The transient response of a real cell phone product under IEC 61000-4-2 excitation was modeled. The proposed method features both high voltage/current and high speed (up to 3 GHz) modeling of typical components, including R, L, C, ferrite, diodes, and IC pins, as well as a complex return path model. The simulation result resembled the tested waveform at both the first and second peaks of the IEC excitation.
The time-dependent destruction threshold of a semiconductor device can be obtained from the tested Wunsch–Bell model with rectangular waveforms. This model accounts for thermal-related junction failures, which have been proven to be the primary cause of a semiconductor junction’s failure mode under EOS.
To determine the device failure under an arbitrary waveform based on the knowledge of the TLP tested failure power/time relationship, one can identify whether or not any portion in the arbitrary waveform P(τ ) injected the same amount of energy as a certain destructive rectangular pulse. Our proposed checking algorithm (12) and its simplified version (13) can be applied for IEC excitation scenarios. For other injection profiles in which the power peak does not occur at the very beginning of the whole waveform, (13) cannot be applied. It must be noted that the checking algorithm is based on a thermal failure model. In rare cases when a component is vulnerable to voltage breakdown, one needs to compare the simulated voltage profile to the TLP The proposed model is very suitable for both pre- and postdesign analysis due to its high computational efficiency. An engineer can quickly understand the holes in a design as long as off-the-shelf circuit models and failure threshold models can be provided readily by the device vendor. Besides using an extracted equivalent circuit model, one may model the return path more precisely with its geometry at the cost of computing time, as what have been done in . A component can always be characterized with automated TLP and VNA measurements. characterized voltage threshold to predict if the device would damage. However, when a component is not available, one could model it from its geometry and material . In addition to failure analysis, the system model also can be used to analyze ESD-induced interference in SI problems, with an additional coupling path model.