1.
Introduction
By providing additional vertical interconnection resources, through-silicon via technology provides an attractive route to “More than Moore” [1], resulting in the smaller footprint, higher integration density, lower latency and power consumption[2]. Though TSV is a promising way to realize 3-D IC and has many benefits over traditional 2-D IC, the reliability of TSV is much more serious since TSV structure is more complex and more sensitive to the local multi-physics coupling. Especially for the thermal-stress coupling, the coefficient of thermal expansion (CTE) mismatch induced by high temperature difference will generate residual stress in TSV and silicon substrate, resulting in the TSV structure cracking and even the degradation of transistor mobility[3, 4].
To investigate the effect of thermal-stress on the reliability of 3-D IC, the primary work is configuring out the stress distribution in TSV and substrate as accurate as possible. Experimental testing, such as Micro-Raman[5], is the most accurate way to approach the real situation, but test-chips and optical instruments need high cost. Though the finite element analysis (FEA) simulation is a good choice with the lower cost, FEA is computationally expensive and infeasible for full-chip scale stress analysis. Considering the cost and efficiency, most previous works focus on the analytical modeling of TSV thermal-stress. A uniaxial thermal-stress model of TSV is proposed in Ref. [6], but it has a calculation error because the real TSV stress is biaxial. The biaxial stress model is proposed under axisymmetric assumption[7] and then is refined considering the traction–free condition in Ref. [8]. A semi-analytical stress model of TSV is proposed and applied in the calculation of multi-TSVs thermal-stress at a given point on the substrate[9]. The above works make a great contribution on exploring TSV thermal-stress modeling, but their models are all established under the isotropic assumption, which has the computation error in the stress distribution of silicon substrate. Actually, silicon is an anisotropic material[10]. The stress distribution in substrate should be calculated based on the anisotropic elastic theory.
To correct the calculation error induced by isotropic assumption, we propose an anisotropic thermal-stress model for TSV. The modeling principle and related assumptions are introduced firstly. Then, the derivation of the stress model is presented. Finally, we use Synopsys TCAD-S Interconnect to verify the accuracy of the proposed model.
2.
Modeling principles and assumptions
2.1
Plane strain and plane stress assumptions
To simplify the modeling process and get the analytical solution easier, a 3-D modeling issue is always reduced to be a 2-D issue. As a 3-D structure, the stress distribution on the surface of the substrate is mostly concerned because the device layer and top insulation layer suffer from stress a lot. The transistor performance drift and insulation cracking (TSV extrusion) induced by TSV stress are two critical reliability problems[11, 12]. To model the stress distribution on the substrate surface, it is feasible to reduce TSV thermal-stress modeling from a 3-D problem to a 2-D problem.
In elastic theory, two classic assumptions are used to solve the 2-D plane problem: plane strain assumption and plane stress assumption. Their specific definitions are given as follows:
Plane strain assumption: when the dimension of an object along the z-axis is much longer than its cross section along orthogonal axes directions, the displacements and stresses are independent of the z-axis. Accordingly, the strain εiz = 0 (i = x, y, z), but the normal stress along the z-axis σzz is non-zero.
Plane stress assumption: when the thickness of an object along the z-axis is much smaller than the other dimensions, the stress and displacement are independent of the z-axis. Accordingly, the stress σzz = 0, but the normal strain along the z-axis εzz is non-zero.
Apparently, TSV and the sidewall insulation layer are long bodies which satisfy the plane strain assumption, while the device layer satisfies the plane stress assumption.
2.2
Modeling method for TSV thermal-stress
Based on the plane strain and plane stress assumptions, the modeling of TSV thermal-stress can be separated into two steps, as shown in Fig. 1.
onerror="this.onerror=null;this.src='http://www.jos.ac.cn/fileBDTXB/journal/article/jos/2018/2/PIC/17050014-1.jpg'" class="figure_img" id="Figure1"/>
Download
Larger image
PowerPoint slide
Figure1.
The solving steps of anisotropic thermal-stress of TSV.
The first step only considers TSV and insulation layer. The radius of TSV is a and the thickness of insulation layer is b-a. The barrier layer and seed layer are neglected because their thickness is much smaller compared with the radius of TSV, which almost has no contribution to the total stress in the substrate. The second step considers the substrate with a hole with the radius of b.
Both steps are under the same temperature load ΔT which is defined as the temperature difference between operating temperature and annealing temperature. In the original problem, the interface between insulation layer and substrate are not free to expand as the temperature changes. To correctly reflect the displacement along the interface between insulation layer and substrate, a stress P is exerted on the edge of the insulation layer in step 1 and the edge of the hole in step 2, respectively.
In this paper, the material of TSV is Cu and the insulation layer is SiO2. Cu and SiO2 are still assumed to be isotropic material while Silicon is anisotropic material.
3.
Modeling of anisotropic TSV thermal-stress
3.1
Modeling without considering traction-free
The first step in Fig. 1 is considered. According to plane strain assumption, the displacement uri and normal stress of TSV and insulation layer along radial direction σrri can be given under isotropic assumption[13]:
$$begin{split}&{u_r^i = C_1^ir + displaystyle frac{{C_2^i}}{r},}& {sigma _{rr}^i = frac{{{E^i}}}{{(1 + {upsilon ^i})(1 - 2{upsilon ^i})}}left[ {C_1^i - frac{{C_2^i(1 - 2{upsilon ^i})}}{{{r^2}}} - left( {1 + {upsilon ^i}} ight){alpha ^i}Delta T} ight],}&{qquad :qquad :qquad :i in { { m{Cu,Si}}{{ m{O}}_2}} ,}end{split}$$  | (1) |
where υi is Poisson’s ratio, αi is CTE, Ei is the Young module and r is the radial distance from the center of TSV.
C1i and C2i are constants that can be solved by the specific boundary conditions as follows:
$$left{ begin{array}{*{35}{l}} r=0, ,,,,, u_{r}^{ m cu}=0, r=a,,,,, u_{r}^{ m cu}=u_{r}^{{ m sio}_2},,,,, sigma _{rr}^{ m cu}=sigma _{rr}^{{ m sio}_2}, r=b,,,,,,u_{r}^{{ m sio}_2}=u_{r}^{ m si},,,,,displaystyle sigma _{rr}^{{ m sio}_2}=P. end{array} ight.$$  | (2) |
Under the above boundary conditions, the normal stress of TSV and insulation layer along the radial direction can be expressed as the function of radius r. However, there is an unknown constant P in the expression of uri and σrri, i∈{SiO2, Cu}.
Then, we consider the second step, which is a 2-D plane stress problem under the anisotropic condition. In such case, the stress of the substrate should be solved by introducing a complex variable approach. Lekhnitskii has given the solution of step 2 in the Cartesian coordinate system[14]. The displacement and stress of substrate can be expressed by the complex function Ф(z1), Ф(z2) and their derivatives as follows:
$$begin{split}& sigma _{xx}^{{ m {s_i}}} = 2{mathop{ m Re}nolimits} [mu _1^2Phi '({z_1}) + mu _2^2Phi '({z_2})],& sigma _{yy}^{ m {s_i}} = 2{mathop{ m Re}nolimits} [Phi '({z_1}) + Phi '({z_2})],& u_x^{ m {s_i}} = 2{mathop{ m Re}nolimits} [{p_1}Phi ({z_1}) + {p_2}Phi ({z_2})],& u_y^{ m {s_i}} = 2{mathop{ m Re}nolimits} [{q_1}Phi ({z_1}) + {q_2}Phi ({z_2})],& {q_1} = {a_{12}}mu _1^{} + frac{{{a_{22}}}}{{{mu _1}}},{q_2} = {a_{12}}mu _2^{} + frac{{{a_{22}}}}{{{mu _2}}},& {p_1} = {a_{11}}mu _1^2 + {a_{12}},{p_2} = {a_{11}}mu _2^2 + {a_{12}}.end{split}$$  | (3) |
z1 and z2 are two complex variables z1 = x + μ1y and z2 = x + μ2y, where μ1 and μ2 are the roots of the following characteristic equations, which are only with positive imaginary.
$${a_{11}}{mu ^4} + (2{a_{12}} + {a_{66}}){mu ^2} + {a_{22}} = 0,$$  | (4) |
where a11, a12, a66, and a22 are the elastic constants of silicon substrate.
The complex function Ф(z1) and Ф(z2) are given as follows:
$$begin{split}&{Phi _1}({z_1}) = frac{{overline {{beta _1}} - {mu _2}overline {{alpha _1}} }}{{left( {{mu _1} - {mu _2}} ight){xi _1}}},,,,,, {Phi _1}({z_1}) = - frac{{overline {{beta _1}} - {mu _1}overline {{alpha _1}} }}{{left( {{mu _1} - {mu _2}} ight){xi _2}}},&overline {{alpha _1}} = frac{{Pr}}{2},,,,,, overline {{beta _1}} = frac{{iPr}}{2},&{xi _1} = frac{{{z_1} + sqrt {z_1^2 - {r^2} - mu _1^2{r^2}} }}{{r(1 - i{mu _1})}},,,,,, {xi _2} = frac{{{z_2} + sqrt {z_2^2 - {r^2} - mu _2^2{r^2}} }}{{r(1 - i{mu _2})}},end{split}$$  | (5) |
where r is the radius of hole in the substrate. As depicted in Fig. 1, the radius of hole in the substrate should be b.
It can be seen from Eq. (5) that the complex functions Ф(z1) and Ф(z2) also contain the stress P. Substituting Eq. (5) into Eq. (3), the displacement of step 2 at r = b can be obtained as follows:
$$u_{r = b}^{{ m si}} = (1 + {upsilon _{ m si}}){alpha _{ m si}}Delta T cdot b + left{ {[1 + i({mu _1} + {mu _2})]{a_{11}} + {a_{12}}} ight}P·b,$$  | (6) |
where υsi and αsi are Poisson’s ratio and CTE of silicon, respectively. The first item in Eq. (6) represents the displacement due to thermal expansion and the second item represents the displacement due to stress P.
The displacement of insulation layer
m{si}}{{
m{o}}_{ 2}}} $
$$begin{split} P =&left[ frac{{{b}^{2}}{{e}_{1}}{{f}_{1}}-{{e}_{1}}{{g}_{2}}-{{b}^{2}} {{g}_{1}}{{e}_{2}}+{{g}_{1}}{{f}_{2}}}{({{f}_{2}}{{f}_{1}}-{{g}_{2}}{{e}_{2}}) b}-(1+{{upsilon }_{ m si}}){{alpha }_{ m si}}Delta Tcdot b ight] & times frac{1}{{[1+i({{mu }_{1}}+{{mu }_{2}})]{{alpha }_{11}}+{{alpha }_{12}}}b}, end{split}$$  | (7) |
where e1, f1, g1, e2, f2 and g2 are constants related to the material property, which can be solved from Eq. (1) under the boundary condition Eq. (2). Here, we directly give the above constants as follows:
$$begin{split}& {e_1} = frac{{{E_{ m Cu}}}}{{(1 - 2{upsilon _{ m Cu}})}}{alpha _{ m Cu}}Delta T - frac{{{E_{{ m sio}_{2}}}}}{{(1 - 2{upsilon _{{ m sio}_2}})}}{alpha _{{ m sio}_2}}Delta T,& {f_1} = frac{{ - {E_{{ m sio}_2}}}}{{(1 + {upsilon _{{ m sio}_2}})}}{alpha _{{ m sio}_2}}Delta T cdot frac{1}{{{b^2}}},& {g_1} = frac{{{E_{{ m sio}_2}}}}{{(1 - 2{upsilon _{{ m sio}_2}})}}{alpha _{{ m sio}_2}}Delta T,& {e_2} = frac{{{E_{ m Cu}}}}{{(1 - 2{upsilon _{ m Cu}})(1 + {upsilon _{ m Cu}})}}frac{1}{{{a^2}}} + frac{{{E_{{ m sio}_2}}}}{{(1 + {upsilon _{{ m sio}_2}})}}frac{1}{{{a^2}}},& {f_2} = frac{{{E_{ m Cu}}}}{{(1 - 2{upsilon _{ m Cu}})(1 + {upsilon _{ m Cu}})}} - frac{{{E_{{ m sio}_2}}}}{{(1 - 2{upsilon _{{ m sio}_2}})(1 + {upsilon _{{ m sio}_2}})}},& {g_2} = frac{{{E_{{ m sio}_2}}}}{{(1 - 2{upsilon _{{ m sio}_2}})(1 + {upsilon _{{ m sio}_2}})}}.end{split}$$  | (8) |
Substituting the stress P in Eq. (5) and in Eq. (3), the analytical expressions of normal stress σxxSi and σyySi will be obtained. To facilitate the calculation of the stress in the cylinder coordinate system, the normal stress σxxSi and σyySi can be transformed to the following expressions, as in Ref. [15]:
$$begin{split}& sigma _{rr}^{{ m {s_i}}} = - P cdot {mathop{ m Re}nolimits} [Q({mu _1},{mu _2})F(frac{r}{b},{mu _1}) + Q({mu _2},{mu _1})F(frac{r}{b},{mu _2})],& sigma _{theta theta }^{ m {s_i}} = - P cdot {mathop{ m Re}nolimits} [S({mu _1},{mu _2})F(frac{r}{b},{mu _1}) + S({mu _2},{mu _1})F(frac{r}{b},{mu _2})],end{split}$$  | (9) |
where the functions Q(m, n), F(m, n), R(m, n) and S(m, n) are defined as:
$$begin{split}& F(m,n) = frac{1}{{sqrt {{m^2} - 1 - {n^2}} (m + sqrt {{m^2} - 1 - {n^2}} )}},& S(m,n) = frac{{(i - n)(1 - im)}}{{m - n}},& Q(m,n) = {m^2}S(m,n),R(m,n) = - mS(m,n).end{split}$$  | (10) |
3.2
Modeling refinement considering traction-free
Based on the plane strain solution, the normal stress σzz along the z-axis is non-zero. In order to meet the plane strain assumption, there are two vertical stresses exerting on the TSV and sidewall insulation layer, respectively. The expressions of these two stresses are given as follows:
$$begin{split}sigma _{zz}^{{ m{Cu}}} =,, & 2{upsilon _{{ m{Cu}}}}frac{{{E_{{ m{Cu}}}}}}{{(1 + {upsilon _{{ m{Cu}}}})(1 - 2{upsilon _{{ m{Cu}}}})}}& times left[ {C_1^{{ m{si}}{{ m{o}}_{ m{2}}}} + frac{{C_2^{{ m{si}}{{ m{o}}_{ m{2}}}}}}{{{a^2}}} - (1 + {upsilon _{{ m{Cu}}}}){alpha _{{ m{Cu}}}}Delta T} ight],sigma _{zz}^{{ m{si}}{{ m{o}}_{ m{2}}}} =,, & 2{upsilon _{{ m{si}}{{ m{o}}_{ m{2}}}}}frac{{{E_{{ m{si}}{{ m{o}}_{ m{2}}}}}}}{{(1 + {upsilon _{{ m{si}}{{ m{o}}_{ m{2}}}}})(1 - 2{upsilon _{{ m{si}}{{ m{o}}_{ m{2}}}}})}}& times left[ {C_1^{{ m{si}}{{ m{o}}_{ m{2}}}} - (1 + {upsilon _{{ m{si}}{{ m{o}}_{ m{2}}}}}){alpha _{{ m{si}}{{ m{o}}_{ m{2}}}}}Delta T} ight],C_1^{{ m{si}}{{ m{o}}_{ m{2}}}} = & frac{{{e_1}{f_1} - (P + {g_1}){e_2}}}{{{f_2}{f_1} - {g_2}{e_2}}},{kern 1pt} {kern 1pt} {kern 1pt} {kern 1pt} C_2^{{ m{si}}{{ m{o}}_{ m{2}}}} = frac{{{e_1}{g_2} - (P + {g_1}){f_2}}}{{{e_2}{g_2} - {f_1}{f_2}}}.end{split}$$  | (11) |
However, in our original issue, there is no stress on the TSV and insulation layer along the z-axis, which makes a mismatch between the original vertical boundary conditions and the current vertical boundary conditions. To correct this boundary condition mismatch, the common strategy is exerting two compensating stresses on the top surface of the TSV and insulation layer. These two compensating stresses have the same values with the stresses obtained by Eq. (11), but opposite in direction, as shown in Fig. 2.
onerror="this.onerror=null;this.src='http://www.jos.ac.cn/fileBDTXB/journal/article/jos/2018/2/PIC/17050014-2.jpg'" class="figure_img" id="Figure2"/>
Download
Larger image
PowerPoint slide
Figure2.
(Color online) The schematic of traction-free.
The additional normal stresses will be generated on the surface of silicon substrate under the vertical stress
m{Cu}}} $
m{si}}{{
m{o}}_2}} $
m{si - h}}} $
m{si - h}}} $
m{Cu}}} $
$$begin{split}& sigma _{rr}^{ m {si - h}} = {I_1}frac{r}{{5.5{e^{ - 6}}}} left[ {sigma _{zz}^{ m Cu}frac{{{a^2}}}{2}frac{{1 - {upsilon _{ m si}}}}{{{r^2}}} + sigma _{zz}^{ m Cu}frac{{{b^2} - {a^2}}}{2}frac{{1 - {upsilon _{ m si}}}}{{{r^2}}}} ight],& sigma _{theta theta }^{ m {si - h}} = {I_1}left( {1 + frac{{2b - 2a}}{r}} ight)left[ {sigma _{zz}^{ m Cu}frac{{{a^2}}}{2}frac{{1 - {upsilon _{ m si}}}}{{{r^2}}} + sigma _{zz}^{ m Cu}frac{{{b^2} - {a^2}}}{2}frac{{1 - {upsilon _{ m si}}}}{{{r^2}}}} ight],end{split}$$  | (12) |
where I1 is the modified Bessel functions of the first kind as the empirical factor.
With superposition method, the final effective normal stress σrr and σ?? on the substrate surface can be obtained as follows:
$$sigma _{rr}^{ m {si - eff}} = sigma _{rr}^{ m si} + sigma _{rr}^{ m {si - h}},sigma _{theta theta }^{ m {si - eff}} = sigma _{theta theta }^{ m {si}} + sigma _{theta theta }^{ m {si - h}}.$$  | (13) |
4.
Model verification
In this paper, the Synopsys TCAD-SInterconnect[17] is used to verify the accuracy of the proposed model. To reduce the computation time, the 3-D simulation structure in TCAD is a quarter of the original structure, as shown in Fig. 3. The top mesh is denser than the bottom mesh because we are concerned with the stress on the substrate surface more than other positions. To ensure the accuracy of simulation, the proper boundary conditions are exerted on the simplified simulation structure. The Dirichlet boundary conditions are exerted on the left, right, front and back faces of the simulation structure. In addition, the displacement at the bottom surface of the simulation structure is restricted.
onerror="this.onerror=null;this.src='http://www.jos.ac.cn/fileBDTXB/journal/article/jos/2018/2/PIC/17050014-3.jpg'" class="figure_img" id="Figure3"/>
Download
Larger image
PowerPoint slide
Figure3.
(Color online) 3-D simulation structure in TCAD. (a) A quarter of the original TSV structure. (b) Simulation structure in mesh.
The material properties and parameters involved in simulation are given in Table 1. The operating temperature is 25 °C and the annealing temperature is 250 °C. Therefore, the temperature load in simulation is –225 °C. In addition, the wafer orientation is assumed to be (110).
| Name | Description | Value | ||
| υCu | Poisson’s ratio of Cu | 0.343 | ||
| αCu | CTE of Cu (ppm/°C) | 17.7 | ||
| ECu | Young’s module of Cu (GPa) | 115.5 | ||
| υ$_{ m{sio}{_2}}$ | Poisson’s ratio of SiO2 | 0.16 | ||
| α$_{ m{sio}{_2}}$ | CTE of SiO2 (ppm/°C) | 0.51 | ||
| E$_{ m{sio}{_2}}$ | Young’s module of SiO2 (GPa) | 71.7 | ||
| υsi | Poisson’s ratio of Si | 0.28 | ||
| αsi | CTE of Si (ppm/°C) | 3.05 | ||
| Cij | Stiffness of silicon with orthotropic material property (GPa) | C11 = 165, C12 = 63.9, C44 = 79.6 | ||
| a | Radius of TSV (μm) | 2.5 or 5 | ||
| b-a | Thickness of insulation layer | 0.5 | ||
| ΔT | Temperature load (°C) | ?225 |
Table1.
The parameters involved in the TCAD simulation.
Table options
-->
Download as CSV
| Name | Description | Value | ||
| υCu | Poisson’s ratio of Cu | 0.343 | ||
| αCu | CTE of Cu (ppm/°C) | 17.7 | ||
| ECu | Young’s module of Cu (GPa) | 115.5 | ||
| υ$_{ m{sio}{_2}}$ | Poisson’s ratio of SiO2 | 0.16 | ||
| α$_{ m{sio}{_2}}$ | CTE of SiO2 (ppm/°C) | 0.51 | ||
| E$_{ m{sio}{_2}}$ | Young’s module of SiO2 (GPa) | 71.7 | ||
| υsi | Poisson’s ratio of Si | 0.28 | ||
| αsi | CTE of Si (ppm/°C) | 3.05 | ||
| Cij | Stiffness of silicon with orthotropic material property (GPa) | C11 = 165, C12 = 63.9, C44 = 79.6 | ||
| a | Radius of TSV (μm) | 2.5 or 5 | ||
| b-a | Thickness of insulation layer | 0.5 | ||
| ΔT | Temperature load (°C) | ?225 |
SInterconnect does not support editing elastic constants. In such case, we use the stiffness coefficient in modeling and simulation. The relationships between elastic constants and stiffness coefficient are given as follows:
$$begin{split}& {a_{11}} = {a_{22}} = frac{{{C_{11}} + {C_{12}}}}{{C_{11}^2 + {C_{11}}{C_{12}} - 2C_{12}^2}},& {a_{12}} = - frac{{{C_{12}}}}{{C_{11}^2 + {C_{11}}{C_{12}} - 2C_{12}^2}},{a_{66}} = frac{1}{{{C_{44}}}}.end{split}$$  | (14) |
The comparison of normal stress
m{si - eff}}} $
m{si - eff}}} $
m{si - eff}}} $
m{si - eff}}} $
onerror="this.onerror=null;this.src='http://www.jos.ac.cn/fileBDTXB/journal/article/jos/2018/2/PIC/17050014-4.jpg'" class="figure_img" id="Figure4"/>
Download
Larger image
PowerPoint slide
Figure4.
(Color online) The comparison of normal stress
m{si - eff}}} $
m{si - eff}}} $
5.
Conclusion
A more accurate analytical thermal-stress model of TSV is proposed both considering the anisotropic property and being traction-free. It is a feasible way to solve the anisotropic thermal-stress issue for TSV by using the Lekhnitskii complex method. The comparison results show that the maximum error between the model and the TCAD 3-D simulation is below 5%. With the proposed model, the thermal-stress distribution of substrate caused by TSV will be predicted more accurately than the isotropic stress model. In future works, the proposed model can be applied in TSV placement driven by thermal-stress, leading to the faster and more accurate optimization for 3-D IC physical design.
