2.1.

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.

Fig. 2

Schematics of 3-D chip stack.

Curing the package at high temperatures ($\sim 150°\mathrm{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,10 A 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.¹² This 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.

Fig. 3

(a) The architecture of the analyzed die; (b) distribution of global stress components in tier1 across the subsurface region indicated in (a) by the dashed line.

Table 1

Thicknesses and thermomechanical properties of materials in the 3-D package used in FEA simulations.

2.2.

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 ${\mathit{u}}^W(x,y)=\{{\mathit{u}}_{\mathrm{top}}^W(x,y),{\mathit{u}}_{\text{bot}}^W(x,y)\}$ 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 $\mathit{u}=\{u_x,u_y,u_z\}$, 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., ${\alpha }_b<{\alpha }_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¹³ (Fig. 4). Due to a mismatch between the bump height shrinkage and the vertical deformation of the underfill $\mathrm{\Delta }{u}_z^W$ 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

Fig. 4

Bump induced deformation of silicon.

Here, $H$ is the bump height, $\mathrm{\Delta }{u}_z^W/H$ represents warpage-related vertical strain at the bump location, and $\mathrm{\Delta }T$ describes the thermal load. Similarly, due to warpage-induced lateral displacements $u_{x(y)}^W$, each bump generates the following tangential load

The additional deformation of the die surface $u^{\mathrm{bump}}$ due to loading Eqs. (1) and (2) can be calculated by solving the corresponding problems of contact mechanics.^14,15 Additional 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_{\gamma }$ and ${\nu }_{\gamma }$ 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\gg D$, we have introduced an effective “radius of bump interaction” $r_{\mathrm{eff}}$, which restricts the size of the area affected by the bump. This $r_{\mathrm{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 ${\mathit{u}}^{\mathrm{bump}}$ to warpage-induced one ${\mathit{u}}^W$:

Equation (4) shows that all bumps at $\mathit{r}={\mathit{r}}_i$ within the region $|\mathit{r}-{\mathit{r}}_i|\le r_{\mathrm{eff}}$ affect the displacement generated at this point.

2.3.

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 ${\rho }_M$ in each metal layer has been developed based on the theory of mechanical properties of anisotropic composite materials.¹² This 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 ${\rho }_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:¹²

FEA model is used for the die-scale simulation of the stress components $\sigma =\{{\sigma }_x,{\sigma }_y,{\sigma }_z\}$. 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:

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.

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\times 1300{\mu \text{m}}^2$.

Table 2

Thermomechanical properties of bumps and back-end-of-line (BEoL) interconnect materials.

2.4.

Compact Model for TSV-Induced Stress

TSV fabrication generates strain in the surrounding silicon: a thermal load $\mathrm{\Delta }T$ generated by cooling the chip down from copper anneal temperature to room temperature creates a thermal mismatch strain ${\epsilon }_{\mathrm{th}}=({\alpha }_{\mathrm{Cu}}-{\alpha }_{\mathrm{Si}})\mathrm{\Delta }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_{\mathrm{TSV}}/2$, where $D_{\mathrm{TSV}}$ is the TSV diameter, can be approximated by Lame formula:¹⁴

Due to axial symmetry, the stress components of Eq. (7) in polar coordinates are independent of azimuthal angle $\theta$. Equation (7) is valid when the TSV is placed far enough from the die periphery (at a distance much larger than $D_{\mathrm{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.¹⁶ However, 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 a renormalization of the parameter ${\epsilon }_{\mathrm{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 ${\sigma }_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 ${\sigma }_x,{\sigma }_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

2.5.

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 ${\nu }_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_i^{\prime }$ due to the action of the forces $F_i\sim (1-E_{\mathrm{STI}}/E_{\mathrm{Si}}){\sigma }_i^{\mathrm{init}}$. 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 $\sigma \text{'}$, which can be obtained from solution of the corresponding force balance equation. For example, for each cut-layer directed along $x$-axis, the force balance equation for calculating the stress component ${\sigma }_x^{\prime }$ is the following

Fig. 6

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

It reduces the problem to the solution of the following equation for the lateral displacement ${\mu }_x^{\prime }$:

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_i^{\prime \prime }$ and $E_i^{\text{'}}$ are the known functions of the materials properties and segment geometries; ${\epsilon }_i^0$ 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 layout-induced 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 cut-lines, 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).

2.6.

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 well-determined piezoresistance coefficients:¹⁹

Here, $\mathrm{\Delta }U=U(\sigma )-U0$, $U(\sigma )$ is the low-field mobility in the stressed silicon, $U0$ is the mobility at the zero stress condition;²⁰ the signs and values of piezoresistance coefficients ${\pi }_x,{\pi }_y,{\pi }_z$ are known to be dependent on crystallographic orientation of silicon surface and on the transistor channel orientation.^19,21 Correspondingly, 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 $\langle 110\rangle$ channel direction. The magnitude of the mobility change depends on the pitch of bump pattern and on device location relative to bumps.

Fig. 7

Mobility variation in nMOS and pMOS with $\langle 110\rangle$ channel orientation caused by bump pattern.

SPICE netlist annotation with the calculated instant parameters MULU0 makes it possible to calculate stress-modified 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-${\mathrm{Si}}_{1-x}{\mathrm{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.