Stress assessment for device performance in three-dimensional IC: linked package-scale/die-scale/ feature-scale simulation flow

. Potential challenges with managing mechanical stress and the consequent effects on device performance for advanced three-dimensional (3-D) IC technologies are outlined. The growing need in a simulation-based design verification flow capable of analyzing a design of 3-D IC stacks and detecting across-die out-of-spec variations in MOSFET electrical characteristics caused by the die thinning and stacking-induced mechanical stress is addressed. The development of a multiscale simulation methodology for managing mechanical stresses during a sequence of designs of 3-D IC dies, stacks, and packages is focused. A set of physics-based compact models for a multiscale simulation is proposed to assess the mechanical stress across the device layers in silicon chips stacked and packaged with the 2.5D interposer-based, and true 3-D through silicon via-based technology. A simulation flow is developed for the hot-spot checking in different types of devices/circuits such as digital, analog, analog matching, memory, IO, characterized by different sensitivities to the stress-induced mobility variations. A calibration technique based on fitting to measured electrical characteristics of the test-chip devices is presented. The limited characterization or measurement capabilities for 3-D IC stacks and a strict “ good die ” requirement make this type of analysis critical in order to achieve an acceptable level of functional and parametric yield. © The Authors


Introduction
][3] Three-dimensional (3-D) IC stacking technologies (including 2.5D interposer-based approaches), employing thinned wafer/through silicon via (TSV) structures, are novel solutions that result in reduced floor space, higher bandwidth, and reduced energy consumption.To enable 3-D chip stacking, new components were adopted by interconnect technology: through-silicon vias to provide connectivity between the back-end-of-line (BEoL) and back-chip redistribution layer (BRDL) metallization of some of the stacked dies, and solder or metal bumps and pillars for interconnecting the neighboring dies.The processing of high-density TSV structures through thinned dies and subsequent 3-D stacking is a promising technological alternative to the traditional two-dimensional (2-D) lithography/etch scaling.However, several issues have to be addressed to guarantee the needed product performance and reliability.Among them are management of mechanical stress and high Joule heating as well as reliability issues such as electromigration and stress migration. 4Many of the process steps employed by 3-D IC technology act as stress sources that can affect the chip performance and reliability.These are TSV etching and filling, wafer/die thinning, wafer bumping, high-temperature solder reflow, chip stacking, etc. 3,5 Hence, it is important to have a capability to accurately assess the stress generated during 3-D IC stacking.
The traditional chip-package interaction (CPI) concerns are related to reliability issues caused by high peeling and shear stresses (hundreds of MPa), which are able to cause cracking, delamination, etc., resulting in shorts or opens.These stresses represent a serious risk for chips' durability; nevertheless, the existing measurement and screening procedures enable the failure modes detection, allowing to eliminate the detected threats.The failure risk can be minimized by using appropriate package assembly materials, which can decrease substantially the CPI stress.Analysis of the reliability issues is out of scope of this paper.
The stress simulator described here reflects the needs of designers and manufacturers of 3-D IC chip packages to control deviations of design-in MOSFET parameters that cause degradation of the chip performance.Even relatively low values of CPI stress (less than hundred MPa) can cause parametric failures of circuits due to change of charge carrier mobility in transistor's channel regions.Stress gradients across regions where circuits are located can worsen the chip performance even if the stress level is not high.Due to the lack of measurement and screening procedures for detecting mobility shift in transistors, the developed simulator becomes an additional metrology tool for 3-D IC package inspection.
It is a very challenging task to get an entire picture of the modification of the stress distribution across device layers caused by 3-D IC technology.Different scales of the stress distributions must be taken into account: die scale (i.e., several millimeters) global stress variations which are generated due to thin die bending; bump and TSV scale (i.e., 10 to 100 μm) local stress variations, which are induced by the package component assembly; device scale (10 to 100 nm) variations due to transistor layer nonuniformity.Traditional methods such as finite-difference analysis and finite-element analysis (FEA) cannot be employed for a simulation of the transistor channel stress distribution across a die due to the size of a model, which can easily reach hundreds of millions degrees of freedom.The established FEA-based stress simulators have typically been used for addressing the traditional CPI effect where a silicon chip was modeled as a homogeneous isotropic piece of silicon.Details of chip structures (layer information, layout, etc.) have not been considered, and the problem with calculating the transistor-to-transistor intrachannel stress variation and consequent variation in transistor electrical characteristics has not been addressed yet. 6In order to be able to consider these effects, the stress simulation methodology should be capable to resolve scales of the order of a transistor size (approximately nanometers) and, to account for all major internal (layout-induced) and external (e.g., packaging) stress sources affecting a particular device.Compact modelbased approaches for the layout-induced stress effects, which typically employ the empirical modeling, 7,8 cannot take into account the package-induced variations in transistor characteristics.Because of the lack of physics-based foundation, this kind of modeling cannot provide a link to the physics-based package-scale simulation in order to include a CPI-induced stress loading.A look-up table methodology is not practical due to a large size of each device's surrounding layout area (radius of up to 5 μm) that should be accounted for a correct stress prediction.It will require the generation of an enormous amount of local layout configurations around a gate in order to obtain a proper representation.Therefore, the development of compact physics-based models is the only possible solution to achieve the ultimate goal of predicting and simulating transistor-to-transistor stress variation across a device layer.These compact models are based on analytical treatment of the elasticity problems: approximate solutions of the basic differential equations are used to generate a set of simple algebraic equations allowing fast full-chip analysis of stress values in transistor channels.Resolved stress variation should be converted further into the variation of the device electrical characteristics.
This paper describes the recently developed physicsbased simulator that predicts variation in transistor's electrical characteristics, caused by CPI.The simulator provides an interface between layout formats (GDS II, OASIS) and FEA-based package-scale models.

Simulation Methodology
The proposed simulation methodology/flow, which will be discussed in detail below, results in the development of a new type of design verification tools.The tool is capable of analyzing any 3-D IC die stack design with regard to out-of-spec variations in device electrical characteristics caused by mechanical stress generated by warpages of the stacked dies, by solder bumps pressure generated in a course of die stacking, as well as by mismatch of thermomechanical properties of TSV and silicon die bulk.The target of the developed simulation flow is the calculation of across-die distributions of the stress normal components fσ x ; σ y ; σ z g (i.e., diagonal components of stress tensor) inside transistor channels and its conversion into the stress-modified electrical characteristics of the transistors.
This simulation flow, shown schematically in Fig. 1, can be described as a sequence of the simulation steps performed with different simulation tools: 1. FEA-based package-scale simulation of the stack deformation, which is originated by the package assembly.Output: fields of displacements on the surfaces of all stacked dies.2. Compact model for bump-induced displacements.
Output: fields of total displacements, caused by die warpage and bump-die contact mechanics, on the die surfaces.3. Compact models for nonuniform, coordinate-dependent mechanical properties of BEoL and BRDL interconnects and silicon bulk with TSVs. 4. FEA-based die-scale simulation of the strains generated in the die by the surface displacements calculated at the step 2; material properties are described by matrices created at the step 3. Output: CPI-induced strain and stress distributions across the device layer. 5. Compact model for TSV-induced stresses.6. Compact model-based simulation of the transistor-totransistor stress variation across the device layer, and the conversion of the final stress into the mobility multiplier (transistor instant parameter MULU0).This stress variation is caused by a layout-induced relaxation of the stress simulated at the previous step.Output: components of CPI and TSV induced stresses, and instant parameters for each analysed transistor; the annotated SPICE netlist.
Both FEA simulations and compact modelling are performed assuming an elastic behavior of all materials involved in package assembly.The possibilities of plastic deformation of solder bumps and viscoelastic behaviour of interlayer dielectric and underfill are not considered in the presented simulation flow.Viscoelastic properties of underfill layer above glass transition temperature T g are not modeled assuming that the major input into stress caused by the bumps contact interactions with the dies (thin multilayered plates) was developed at temperatures below T g .Within the accepted linear theory of elasticity, the superposition of stresses induced by various assembly components can be used to obtain the total stress distribution.According to the flow schematics (Fig. 1), the stress created by TSVs is considered separately and is assumed to be independent of wafer deformation.Obviously, TSVs cannot be resolved in package scale FEA simulations (due to their small sizes), which implies that the possible effect of TSVs on assemblyinduced wafer deformation is neglected.It should be mentioned that this effect is captured at the FEA die-scale simulation step by introducing the nonuniform, coordinatedependent mechanical properties of the silicon die bulk layer governed by TSV locations.
The described simulations methodology has resulted in the development of a flow for assessment of stress-induced variations in the carrier (electrons and holes) mobility in MOSFET channels, which, when it is plugged into the standard design flow, can be used for the design hot-spots screening and for accurate simulation of a variety of chip characteristics such as performance, leakage, power, etc.

Package-Scale FEA-Based Simulation
Assembly of 3-D IC stack involves the wafer-to-wafer, or die-to-wafer or die-to-die mounting.Mechanical and electrical integrities of the stacked dies are achieved by employment of solder balls/metal pillars and TSVs (Fig. 2).Large flip-chip (FC) bumps connect the tier1 die to the substrate, while the set of TSVs and micro-bumps provide the die-to-die connection.The assembly of 2.5D packages is similar but involves only a die-to-interposer mounting.The encapsulation of the dies is performed by employment of an underfill resin in the gaps between dies, and finally molding the whole stack.
Curing the package at high temperatures (∼150°C) and subsequent cooling down to room temperature create substantial thermal load, resulting in the bending (warpage) of the dies, due to mismatch of the thermomechanical parameters of different layers.This process generates large stresses in the dies, mainly due to large values of the coefficient of thermal expansion (CTE) of the underfill, which results in "shrinking" of silicon die in lateral directions.FEA-based simulation, which is a traditional method for analyzing a mechanical behavior of the package, is included in the developed flow.Another possible approach is an approximate analytical analysis of stresses in multilayered structures. 9,10A review of different types of analytical techniques available for resolving the stress-related problems in silicon technology was presented in Ref. 11.Although the analytical modeling can enable a fast computing of some of the stress components in 3-D IC stack, it is expected to be less accurate than the FEA analysis in the case when tiers of different lateral sizes are stacked in the package as shown in Fig. 2. Obviously, the FEA-based package simulation cannot resolve geometrically the effect of hundreds of FC-bumps and microbumps.Therefore, all layers in the stack are considered as the homogeneous thin plates.At this step, the thermomechanical characteristics of underfill layers are defined using the so-called rule of mixtures. 12This means that the bumps are assumed to be "smeared" throughout the entire layer.The mechanical characteristics of this "smeared" layer are obtained by the volume averaging of the corresponding properties of bumps and underfill in accordance with their volume fractions.
An example of the simulated distribution of the stress components in two-tier stack across tier1 is shown in Fig. 3.This stress, which is generated due to mismatch in thermomechanical properties of the constituent layers, is accompanied by the warpage of the dies.With the thicknesses of the layers in the package and values of material properties collected in Table 1, the values of lateral components of the stress exceed 100 MPa.Larger stress can be generated with decreasing die thickness.Evidently, stress distribution is nonuniform near the edges of the analysed (tier1) die.In addition, a strong nonuniformity is observed in the region under the periphery of upper (tier2) die, which is caused by vertical pressure created by tier2 die edges.

Compact Model for Bump-Induced Displacements
Although the FEA simulation allows us to obtain the detailed variation of the warpage-induced stresses (Fig. 3), it does not take into account the local deformation/stress generated by bumps.To capture this deformation, an analytical modelling linked with the package-scale FEA simulation was developed.Fields of warpage-related displacements u W ðx; yÞ ¼ fu W top ðx; yÞ; u W bot ðx; yÞg of the top and bottom surfaces of all layers calculated by FEA tool can be used for linking the package-scale simulations with the bump-scale model.These displacements allow calculating analytically the approximate local deformations generated around each bump.Here and below the displacement u implies 3-D vector u ¼ fu x ; u y ; u z g, which defines the stress components in the employed approximation of elastic materials.
Figure 4 shows the local deformation that is generated around the FC bump due to chip cooling from the processing temperature down to operating temperature.Bumps, consisting of intermetallic compounds (in the case of FC bumps) or copper (micro-bumps), are characterized by the substantially larger Young's modulus (E b ) than the underfill material (E u ), i.e., E b > E u , but by smaller CTE (i.e., α b < α u ).As a result, the placement of the bump into the underfill layer should modify the warpage-related displacement fields both in lateral (x; y) and vertical (z) directions 13 (Fig. 4).Due to a mismatch between the bump height shrinkage and the vertical deformation of the underfill Δu W z at the bump location, both caused by the cooling, the bump acts as an indenter deforming both the substrate and die subsurface regions in the case of FC bumps.A vertical load generated by the bump can be written as Here, H is the bump height, Δu W z ∕H represents warpagerelated vertical strain at the bump location, and ΔT describes the thermal load.Similarly, due to warpage-induced lateral displacements u W xðyÞ , each bump generates the following tangential load where D is the bump diameter, and u W xðyÞ ∕D represents warpage-related lateral strain in the bump region.It should be mentioned that the thermal mismatch of bump and underfill is able to reduce the tangential load provided by Eq. ( 2).Besides, an extra load is generated by the shear deformation of the bumps, which originates due to different CTE of the silicon die and of the substrate.This shear component of bump deformation is more essential for the bumps located near the die edges.Formulas (1) and ( 2) establish the link between the bump compact model-induced displacements and the warpage-related displacements calculated with the FEA model.
The additional deformation of the die surface u bump due to loading Eqs. ( 1) and ( 2) can be calculated by solving the corresponding problems of contact mechanics. 14,15Additional field of z-displacements at the die surface that is generated by a single bump is Here, r is the distance to the center of the bump, KðD∕2rÞ and EðD∕2rÞ are the complete elliptic integrals of the first and second kinds correspondingly, with modulus D∕2r, and E γ and ν γ are the Young's modulus and Poisson's factor of the material which is in contact with the bump.Taking into account a fast decay of the bump contribution to the total vertical displacements at large distance r ≫ D, we have introduced an effective "radius of bump interaction" r eff , which restricts the size of the area affected by the bump.This r eff estimates the bump model accuracy.Similarly to Eq. ( 3), the lateral displacements caused by loads (1) and ( 2) can be obtained by using the force-balance equation based analysis.
Finally, the total displacements at the die surfaces can be obtained by adding the bump-induced displacements u bump to warpage-induced one u W : Equation ( 4) shows that all bumps at r ¼ r i within the region jr − r i j ≤ r eff affect the displacement generated at this point.

Effect of Nonuniform Interconnect, BRDL and
Silicon Bulk Mechanical Characteristics on the CPI-Induced Stress Simulated with the Die-Scale FEA Model Strain and stress distribution across a device layer, which is a layer located in the silicon interior just nanometers below  the interface with the BEoL interconnect, should be calculated with a FEA tool by implementing the field of displacements (4) as the boundary conditions (BCs) for the die faces.At this step, each die is considered as a multilayer stack consisting of a silicon layer, BEoL interconnect, and BRDL layers.These layers are characterized by the nonuniform spatial distributions of the elastic and thermal properties determined by their layouts.An assessment of additional stress variation caused by the nonuniformity of the mechanical properties of these layers is implemented in the developed simulation flow.
A calculation methodology of the effective Young's modulus, Poisson's ratio, and CTE as functions of metal density ρ M in each metal layer has been developed based on the theory of mechanical properties of anisotropic composite materials. 12This methodology requires a division of all considered composite layers into a number of bins.Average values of thermomechanical properties of the composite material are calculated for each bin on the basis of the average metal density ρ M inside the bin.The size of the bin should be chosen based on required simulation accuracy: the finer partitioning provides more accurate results at the expense of the run time.Depending on routing direction of the metal wires, the Young's modulus of i'th bin at j'th layer should be calculated using one of the following formulas: 12 where the symbols k and ⊥ denote directions parallel and perpendicular to routing direction correspondingly; E M and E D are Young's modulus of metal and interlayer dielectric.The similar approach is used for calculating Poisson's ratio and CTE of all composite layers.FEA model is used for the die-scale simulation of the stress components σ ¼ fσ x ; σ y ; σ z g.In this simulation, the BCs are represented by previously calculated die face displacements, and the material properties are represented by the calculated position-dependent mechanical properties of the composite layers: σðx; y; z ∈ devicelayerÞ ¼ f½u surf ðx; yÞ; Eðx; yÞ; νðx; yÞ; αðx; yÞ: (6) Figure 5 demonstrates the distributions of stress components across the device layer of tier1 of the stack shown in Fig. 3(a).The studied region is located far from the die edges; therefore, the obtained periodical distribution of the stress components reflects a periodical pattern in bump locations.Across die variation of the BEoL mechanical properties results in the stress distribution blurring.The colormaps, showing stress values in Fig. 5, were obtained using the material properties presented in Table 1 and the properties of the bumps, interlayer dielectric (ILD), and copper presented in Table 2.

Compact Model for TSV-Induced Stress
TSV fabrication generates strain in the surrounding silicon: a thermal load ΔT generated by cooling the chip down from copper anneal temperature to room temperature creates a thermal mismatch strain ε th ¼ ðα Cu − α Si ÞΔT due to large CTE difference between copper and silicon.The distribution of radial and circumferential stress components around the TSV, at a distance r > D TSV ∕2, where D TSV is the TSV diameter, can be approximated by Lame formula: 14 Due to axial symmetry, the stress components of Eq. ( 7) in polar coordinates are independent of azimuthal angle θ.Equation ( 7) is valid when the TSV is placed far enough from the die periphery (at a distance much larger than D TSV ), where the boundary effects can be neglected.It should be mentioned that for more accurate description of TSV-induced stress distribution, the deformation of silicon surface due to prevailing vertical shrinking of the TSV in comparison with silicon should be considered. 16However, it can be shown that for the device layer, which is located very close to the silicon/BEoL interface, this surface effect results in just Table 2 Thermomechanical properties of bumps and back-end-ofline (BEoL) interconnect materials.a renormalization of the parameter ε th .In the discussed simulation flow, the value of this parameter should be determined by calibration.This should validate the employment of the Eq. ( 7) as an appropriate approximation for the effect of TSV on device characteristics.Vertical component of stress σ z is essential only in the immediate vicinity of TSV and can be neglected since the devices are not allowed to be placed inside the so-called "keep-out-zones" around TSVs, where zone size is prescribed by the design rules.
The radial and circumferential components [Eq.( 7)] of TSV-induced stress can be transformed to the Cartesian components σ x ; σ y by This transformation provides These stress components should be determined in the same grid points, which were used for calculation of CPI stress components [Eq.( 5)].Finally, assuming that CPI and TSV stresses are independent and contribute additively, the total stress distribution is determined as

Simulation of the Transistor Intrachannel Stress Components
The generated across-die distributions of stress components make it possible to calculate the averaged stress components inside transistor channels by using a simple interpolation procedure.However, at a device scale, the composite nature of the device layer becomes important.In fact, this layer represents a sequence of silicon islands separated by the shallow trench isolation (STI) regions filled by deposited silicon oxide.Keeping in mind a drastic difference in the mechanical characteristics of silicon and silicon oxide, as well as the fact that both the Si islands (i.e., the regions where the transistors are located) and STI regions are characterized by wide dispersion in sizes and shapes, we can conclude that strain distribution should depend on the local layout configurations of the device layer.FEA-based die-scale simulations of the CPI stress described in Sec.2.3 cannot resolve the billions of shapes existing in devices layout.Therefore, in the next step of the multiscale simulation flow, a compact model-based calculation of the device-scale stress variations should be performed.This compact model considers the strain/stress distribution, calculated with Eq. ( 9), as an intial distribution that should be changed (relaxed) in accordance with the layout geometries.To calculate this stress redistribution, each transistor channel and the neighboring layout are portitioned on a set of "cut-layers" (both in x-and y-directions).Each cut-layer consists of the device layer and the silicon bulk, as shown in Fig. 6(a).The device layer consists of the sequence of segments representing slices of silicon islands and STIs separating silicon islands.Each i'th segment in the device "composite" layer is characterized by its length L i , Young's modulus E i and Poisson's ratio ν i .The stresses given by Eq. ( 9) are considered as the "initial" stresses in each segment of the cutline.The difference in elastic properties of the neighbouring segments results in redistribution of this stress: each segment edge experiences additional lateral displacement u 0 i due to the action of the forces i .These displacements can be obtained from the solution of the force balance equation that takes into account the interaction between adjacent segments and the traction with the silicon bulk.Initial stress redistribution can be described as generation of an additional stress σ 0 , which can be obtained from solution of the corresponding force balance equation.For example, for each cutlayer directed along x-axis, the force balance equation for calculating the stress component σ 0 x is the following where τ xz is the shear stress component.Assuming that the vertical displacements everywhere in the plane are small in comparison with the lateral displacements, we can employ the following representation of stresses as functions of displacements It reduces the problem to the solution of the following equation for the lateral displacement μ 0 x : (a) (b) As it was shown in Refs.17 and 18, this equation can be further reduced to the system of linear equations for the lateral displacements at each node of the considered cut-line: Here, E i and E iþ1 are the Young's modulus of i'th and (i þ 1)'th segments; E 0 0 i and E 0 i are the known functions of the materials properties and segment geometries; ε 0 i is the initial strain in i'th node.These equations clearly demonstrate that the coupling between the global stress load (TSVs, packaging, bumps, etc.) and the resulted layoutinduced stress distribution is introduced through the initial strain distribution calculated with the combined FEA-compact model.Solution of the system of Eq. ( 13) provides the distribution of displacements along all considered cutlines, which after standard transformation [Eqs.( 14) and ( 15)] generate the distribution of stress components everywhere in the cut-lines, and particularly inside the transistor channels, Fig. 6(b).

Stress-Induced Variation in Device Electrical Characteristics
The piezoresistive effect describes change in the electrical resistivity of a semiconductor or metal when mechanical strain is applied.In the case of semiconductors, this resistivity change is caused by the strain-induced modification of their band structure and results in mobility changes of the conductive electrons and holes.Stress-induced modification of the charge carrier mobility is calculated based on welldetermined piezoresistance coefficients: 19 Here, ΔU ¼ UðσÞ − U0, UðσÞ is the low-field mobility in the stressed silicon, U0 is the mobility at the zero stress condition; 20 the signs and values of piezoresistance coefficients π x ; π y ; π z are known to be dependent on crystallographic orientation of silicon surface and on the transistor channel orientation. 19,21Correspondingly, for the SPICE instance parameter MULU0, we can write Figure 7 demonstrates the simulated mobility distributions in NMOS and PMOS transistors caused by an FC bump array for the case of h110i channel direction.The magnitude of the mobility change depends on the pitch of bump pattern and on device location relative to bumps.SPICE netlist annotation with the calculated instant parameters MULU0 makes it possible to calculate stressmodified electrical characteristics of any device in the analyzed design by running SPICE simulator.It should be mentioned that the target of the developed simulation flow is the assessment of the effect of stresses associated with the 3-D IC technology on chip performance.These stresses are the unintentional and unwanted addition to the stresses already existing in the device layer generated by intentional stress sources, such as stressed layers (CESLs), source-drain epi-Si 1−x Ge x , stress memorization, etc., and unintentional ones (STI, residual, etc.).The total stress-induced variation of the transistor characteristics is caused by a combination of the layout-induced and package (CPI/TSV)-induced stresses.An instant parameter MULU0_layout describing the effect of the layout-induced stresses on the mobility can be calculated with the foundry calibrated SPICE models.A combined effect of these two types of stresses on the device characteristics can be obtained from SPICE simulations by employing the netlist annotated with MULU0 that represents the product of MULU0_layout and MULU0_3D.3 Outline of the Calibration and Validation Procedures For calibration of the parameters involved into the developed compact models, the measured MOSFETs drain current should be used.Stress-induced variations of drain current are assumed to be caused by mobility variation, which is related to stress by means of relation. 16A representative set of transistors should be used for calibration, which means that the chosen transistors should be located on different distances from the die edges, corners, TSVs, and bumps.The drain current I d can be measured and simulated either in linear or in saturation regimes.As it was mentioned in the previous section, the 3-D stacking induced stress is not the only cause of I d variations.Layout-proximity effects (e.g., stress induced by built-in stressors, process/litho variations, wall-proximity effect, etc.) are able to introduce substantial variations of the device characteristics and should be treated properly.Therefore, the difference of the measured current I meas  16): where the coefficients π I i are related to piezoresistance coefficients π i by means of special transformation matrix that can be obtained by SPICE simulations. 7Then, the calibration engine tunes all the model parameters by minimizing the target function: The calibration results demonstrated below were obtained with the measurements done on a single-chip package where a 100-μm-thin silicon chip was mounted on the package substrate using FC bumps.Since this package did not contain TSVs, only the die warpage-and bump pressure-related effects were considered.The thicknesses and properties of the involved materials are listed in Tables 1 and 2.
Figure 8 demonstrates the schematics of the structures used for the model calibration.They contain FC bump arrays located: Fig. 8(a) shows near the die edge and Fig. 8(b) shows near the die corner.MOSFET devices, located at different distances from FC bumps, were used for calibration.All these devices represent identical p-type transistors, having 1 × 1 μm channel sizes.The identity of transistor sizes and the identical neighboring layouts have allowed us to exclude layout-induced variations in mobility, and thus to avoid the issues related to the accuracy of SPICE models of layout-proximity effects.The calibration results are demonstrated in Fig. 9.The relative changes of the drain current [i.e., %ΔI d ¼ ðΔI d ∕I d Þ • 100%] for the selected p-type transistors are in the range of ~(3 to 7)% for the die-edge region, and ∼ð−3 to 3Þ% for the die corner.These changes are caused by the demonstrated lateral stresses (having average values ∼200 MPa).The effect of vertical stress component is much smaller due to small value of the corresponding piezoresistance coefficient: for the considered h110i orientation of the channels the values of the piezoresistance coefficients are π x ¼ 0.718; The smaller values of the current variation in die corner region are due to higher values of y-component of stress tensor, which reduces the current and for some transistors can surpass the positive effect of the x-component.The obtained correlation between the measured and simulated values clearly demonstrates the calibration capability of the developed compact-model-based simulation flow.Fit of the model predictions to the measurements for all transistors, demonstrated in Fig. 9(c), is characterized by mean-square deviation of ∼0.5% and maximal deviation of ∼1% (in terms of relative change of the drain current).This corresponds to the stress prediction accuracy of ∼15 MPa.
4 Hot-Spot Analysis Developed model, once calibrated with the test-chip data from a particular foundry and a technology node, can be employed for analyzing any design of the chips as long as the foundry manufacturing process and technology node are unchanged.It addresses the transistor-to-transistor variation of electrical characteristics caused by 3-D IC technology-induced stresses.However, full-chip detailed analysis of each transistor performance can be expensive in terms of time and CPU resources.So, it seems reasonable to have a capability of fast hot-spot analysis, or in other words, a capability for fast detecting the devices or circuits for which the CPI/TSV-induced stress impact can destroy the nets functionalities.Potential impact of the most stressed regions can be realized by means of the hot-spot analysis, consisting of the following steps: 1.For each circuit in the design (IO, digital, analog, etc.), the corresponding failing criteria are specified, i.e., the maximal acceptable current variation ðΔI d ∕I d Þ threshold is provided as an input parameter.2. A coarse grid covering the die is constructed.The bin size of the grid is comparable with the characteristic size of the analyzed circuits.3.For each bin, the maximal value of the CPI/TSV stress is obtained using the calculated stress distribution [Eq.( 9)].  as a hot spot and can be visualized by means of color maps as displayed in Fig. 10. 5.A detailed analysis of stress impact on the device/circuit characteristics should be performed for this bin as it was described in Secs.2.5 and 2.6.Figure 11(a) depicts visualization of detailed analysis results.For the highlighted transistor, the calculated mobility variation MULU0 is shown, as well as the stress components induced by TSV, CPI (denoted as "packaging"), and the total stress (sigma).Figure 11(b) demonstrates the obtained results when the methodology of Secs.2.5 and 2.6 was employed for the analysis of the layout-induced stress effect on the device performance: the developed capabilities have allowed to

Conclusions
The paper describes the developed multiscale simulation methodology and flow for the assessment of the stressinduced performance variation in 3-D IC chips.The proposed approach allows the linkage between the packagescale FEA formats and the chip layout formats by means of the developed physics-based compact models.The accuracy of the stress prediction is enhanced by developing and implementing the compact models for the effective thermomechanical properties of the composite BEoL interconnect and BRDL layers.Finally, the compact model of the device-scale stress relaxation allows taking into account every individual transistor and its neighboring layout content.Depending on the design stage where the analysis is made, the simulation results can be used for hot-spot analysis, or for back annotating the SPICE netlist with instant parameter MULU0 allowing more accurate circuit simulations can be done.

Fig. 5
Fig. 5 An example of the die-scale simulation results: distributions of stress components across the device layer of tier1 die.The obtained patterns are due to the effect of FC bumps of 80-μm diameter.The size of the demonstrated region is 1900 × 1300 μm 2 .

Fig. 6
Fig. 6 Schematics of the cut-line: (a) vertical slice, (b) top view demonstrating device partitioning into cutlines.Here, light polygons are silicon islands, narrow lines are the poly gates, small squares are contacts, and dark background is STI.

Fig. 7
Fig. 7 Mobility variation in nMOS and pMOS with h110i channel orientation caused by bump pattern.

d,
and the current calculated with the foundry calibrated SPICE model I SPICE d should be used for calibration: ΔI meas d ¼ I meas d − I SPICE d .The predicted values of the current variations ΔI simul d caused by the 3-D stress in CMOS channels are related to the simulated stress components by the expression similar to Eq. (

15 Fig. 8
Fig. 8 Schematics of the test structures used for model calibration: die edge (a), and die corner (b) regions.

Fig. 9 Fig. 10
Fig. 9 Model calibration results: 3-D stress components and the induced drain current variations for p-type transistors, located: (a) in the die edge region; (b) in the die corner region; (c) overall fit of the model prediction to measurements.

Table 1
Thicknesses and thermomechanical properties of materials in the 3-D package used in FEA simulations.
4. Maximal current variation ðΔI d ∕I d Þ max is determined for each bin.If this variation exceeds the threshold, i.e., ðΔI d ∕I d Þ max > ðΔI d ∕I d Þ threshold , for any type of circuits presented in this bin, then the bin is recognized