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 SiO₂. Cu and SiO₂ 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 u_rⁱ and normal stress of TSV and insulation layer along radial direction σ_rrⁱ 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 υⁱ is Poisson’s ratio, αⁱ is CTE, Eⁱ is the Young module and r is the radial distance from the center of TSV.



C₁ⁱ and C₂ⁱ 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 u_rⁱ and σ_rrⁱ, i∈{SiO₂, 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 Ф(z₁), Ф(z₂) 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)

z₁ and z₂ are two complex variables z₁ = x + μ₁y and z₂ = x + μ₂y, where μ₁ and μ₂ 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 a₁₁, a₁₂, a₆₆, and a₂₂ are the elastic constants of silicon substrate.



The complex function Ф(z₁) and Ф(z₂) 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 Ф(z₁) and Ф(z₂) 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 $u_r^{{
m{si}}{{
m{o}}_{ 2}}} $ is equal to the displacement of silicon substrate $ u_{r}^{text{si}{{text{o}}_{ 2}}} $ at r = b. Consequently, the stress P can be expressed from Eq. (1), Eq. (2) and Eq. (6), as follows:

$$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 e₁, f₁, g₁, e₂, f₂ and g₂ 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 σ_xx^Si and σ_yy^Si will be obtained. To facilitate the calculation of the stress in the cylinder coordinate system, the normal stress σ_xx^Si and σ_yy^Si 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 $ sigma _{zz}^{{
m{Cu}}} $ and$sigma _{zz}^{{
m{si}}{{
m{o}}_2}} $. These additional normal stresses can be solved based on the classic half-plane problem. The solution of half-space under vertical loads has been given in Ref. [16], but this solution is only suitable for the transversely isotropic material and homogeneous materials. To eliminate the error brought by the material property mismatch and homogeneous assumption, we introduced the empirical factors by curve-fitting from the simulation results. The normal stress $sigma _{rr}^{{
m{si - h}}} $ and $sigma _{{theta theta }}^{{
m{si - h}}} $ caused by vertical stress $ sigma _{zz}^{{
m{Cu}}} $ can be given as follows:

$$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 I₁ 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 $sigma _{theta theta }^{{
m{si - eff}}} $ and $sigma _{rr}^{{
m{si - eff}}} $ between the proposed model and the TCAD simulation is shown in Fig. 4. Two different radiuses (r = 2.5, 5 μm) are adopted to verify the scalability of the proposed model. The normal stresses $sigma _{theta theta }^{{
m{si - eff}}} $ and $sigma _{rr}^{{
m{si - eff}}} $ are both oriented along the [110] direction. Although the proposed model is 2-D, the analytical results match well with the 3-D simulation results. The maximum error between model and simulation is below 5%. The models proposed in Refs. [6, 8] are also compared. Apparently, there is an error between the above two models and simulation since their proposed models assume the substrate as isotropic material. As the radius of TSV increase from 2.5 to 5 μm, the error induced by isotropic assumption becomes greater.






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 $sigma _{theta theta }^{{
m{si - eff}}} $ and $sigma _{rr}^{{
m{si - eff}}} $ between proposed model and TCAD simulations. (a) r = 2.5 μm. (b) r = 5 μm.

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.