Engineering Fracture Mechanics 74 (2007) 416–430

[Link]/locate/engfracmech

A study of multiple singularities in multi-material wedges

and their use in analysis of microelectronic
interconnect structures
Yaping Luo, Ganesh Subbarayan *

School of Mechanical Engineering, Purdue University, 585 Purdue Mall, West Lafayette, IN 47906, United States

Received 16 February 2005; received in revised form 27 March 2006; accepted 21 April 2006
Available online 1 August 2006

Abstract

The nature of the singular stress field at bi-material and multi-material wedges with de-bonded or bonded interfaces is
investigated in this paper with the aim of relating singular stress field to failures in microelectronic interconnect structures.
The dependence of the strength and the number of the singularities on the relative material behavior as well as the wedge
angle is extensively analyzed using planar isotropic elasticity. The procedure is used to analyze modern copper interconnect
structures; the orders of singularities and angular functions are asymptotically determined, based on which the stress inten-
sity factors are calculated.
 2006 Elsevier Ltd. All rights reserved.

Keywords: Multiple singularities; Multi-material wedge; Interconnects; Microelectronics

1. Introduction

Copper interconnect technology based on the dual-damascene process is increasingly popular at the present
time. The common cause for the failure of devices constructed using this technology include electromigration
and stress migration [1]. For analyzing both of these failure causes, accurate determination of the stress state
in the interconnect is critical. The interconnect structure is complex and involves many interfaces between dis-
similar materials (Fig. 1) and therefore, stress concentration under temperature change is induced. This leads
to failure during fabrication or test. It is well known that the stress fields at the interface corners exhibit a sin-
gular behavior for a linear elastic material that is proportional to rk [2–5], where r represents the distance
from the interface corner, and k is the order of the singularity. Accounting for the singular behavior of stress
distribution during electromigration simulation would yield more accurate results. In addition, since the inter-
layer dielectric tends to be weak mechanically, they are more prone to thin film cracking or delamination [6].

*
Corresponding author. Tel.: +1 765 494 9770; fax: +1 765 494 0539.
E-mail address: ganeshs@[Link] (G. Subbarayan).

0013-7944/$ - see front matter  2006 Elsevier Ltd. All rights reserved.
doi:10.1016/[Link].2006.04.032
 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430 417

Fig. 1. An illustration of the structure and scaling in copper interconnects (International technology roadmap for semiconductors, 2001).

The failure due to stress concentration usually originates at the singular stress point, and the knowledge of
interfacial fracture mechanics is important to understand the failure mechanisms and to achieve better inter-
connect structure design. To date, detailed analysis of the singular stress field and the possible fracture modes
during processing do not appear to have been systematically studied. This is the goal of the present work.

2. A Survey of asymptotic analysis in multi-material wedges

Singular stress fields arise in electronic devices and packages consisting of dissimilar materials or geometric
discontinuities such as interfaces, corners, and free edges. The stress singularities arising in microelectronics
can be generally grouped into three categories (Fig. 2):

1. Angular corner of a homogeneous material (Fig. 2(a)).

2. Angular corner of a bi-material wedge (Fig. 2(b)).
3. A multi-material wedge with adhesion at all interfaces (Fig. 2(c) and 2(d)).

In general, the stresses near the edge of an interface can be described asymptotically as [7]:
X
N
rij ¼ K n rkn fijn ðhÞ þ rij0 ðhÞ ð1Þ
n¼1

where r and h are polar coordinates, N is the total number of singularities. kn is the order of singularities. The
angular function fijn(h) can be completely determined by the asymptotic analysis. The stress term rij0(h) is pro-
portional to the temperature change, and is also determined by asymptotic analysis; it is dependent on the
geometry as well as the mechanical and thermal properties of the composite structure.
The angular corner of a bi-material wedge (Fig. 2(b)) has been studied quite extensively in the literature.
Bogy [2,3] and Hein and Erdogan [4] were among the earliest researchers to determine the stresses in a bi-
material wedge-shaped region. Analytical results of the order of singularities on different geometries with dif-
fering material properties were explicitly derived in their work. Since the resulting equation for the solution of
the order of singularity is transcendental, exhaustive search using an iterative method such as the Newton–
Raphson method or the Bisection method is the commonly used technique to compute all the solutions. The-
ocaris [8] first analyzed the order of singularity in a perfectly bonded multi-material wedge (Fig. 2(d)). He
explicitly derived an analytical solution for the homogeneous material case. Although, the general form in
Eq. (1) contains multiple singularities, it is often assumed that only the largest singularity is important, and
therefore only the largest singularity was considered in earlier described research. Munz and Yang [9–11] ana-
lytically determined the non singular stress terms in Eq. (1) caused by a change of temperature or edge trac-
tions using Airy stress functions in manner similar to Williams [12]. In general, the angular variation of
stresses in Eq. (1) does not appear to have received a significant attention in the above literature.
418 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430

Fig. 2. Classification of multi-material wedges. (a) Angular corner of a homogeneous material. (b) Angular corner of a bi-material wedge.
(c) Bi-material wedge with adhesion at both interfaces. (d) Multi-material wedge with adhesion at all interfaces.

Pageau et al. [13] proposed a definition of generalized stress intensity factors based on a single complex root
for n-material in-plane problem, consistent with that for bi-material interface crack. They demonstrated the
occurrence of logarithmic singularity at bifurcation points where two complex conjugate roots change to
two real roots, corresponding to a variation in wedge geometry or material properties. Pageau et al. also stud-
ied the angular distribution of stresses associated with the complex roots corresponding to specific wedge
geometries, which included tri-material junction with a de-bonded interface and a bi-material junction.
The angular corner of a homogeneous material (Fig. 2(a)) was extensively analyzed by Seweryn and Molski
[14] for combinations of free, clamped, roller support boundary conditions applied to the wedge faces, and
under opening, sliding, and anti-plane loading modes. Instead of Airy stress function, Seweryn and Molski
substituted displacement fields of the following form into the Lamê equations:
uri ðr; hÞ ¼ rk f ðhÞ; uhi ðr; hÞ ¼ rk gðhÞ ð2Þ
They characterized the complex order of singularity as a function of the vertex angle. In addition, they dis-
cussed the possibility of applying the generalized stress intensity factors to practical problems. Analysis follow-
ing the solution approach of Seweryn and Molski for multi-material wedges appears to be missing in the
literature at the present time.
There appear to be inconsistencies in the literature regarding the nature and order of the singularities. The
magnitude of the order of singularity was assumed to be in the range [0, 1] in most of the papers in the liter-
ature, while in Munz and Yang’s work, the order of singularity is assumed to be valid only when it has a value
below 0.5. There is no agreement on the total number of the singularities also. In the homogeneous material
case, two singularities are expected since symmetric and anti-symmetric boundary conditions will clearly lead
to two different singularities in this case. Also, in analyzing the state of stress at the apex of a composite wedge
consisting of more than one material, most studies have considered only the effect of the largest singularity.
The complexity of the analysis is further aggravated by the oscillatory nature of stress near the point of
 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430 419

singularity; generally, when solving the governing equation for order of singularity, the result is a complex
value. The real part characterizes the strength of the singularity, while the imaginary part manifests as an oscil-
latory stress.
Van Vroonhoven [15] appears to be among the first to study stress singularities in bi-material wedges occur-
ring in Microelectronic structures. He carried out a parametric study of the bi-material wedge with two
bonded interfaces (Fig. 2(c)), and compared the solution to the case when the interface was delaminated
and concluded that cracking occurs only in the presence of delamination. Only the largest singularity was com-
puted in Van Vroonhoven’s work; thus, studies that have applied the results of Ref. [15] have not considered
the impact of multiple singularities [16,17].
Liu et al. [18] were among the first to demonstrate the importance of multiple singularities in microelec-
tronic structures. Liu et al. obtained a better fit of stress field to the finite element results by taking into
account two singularities, thus showing that it is critical to include multiple singularities to better characterize
the singular stress field. But, Liu et al. did not study failure and their analysis was limited to the case of an
angular corner of a bi-material wedge.
In the current paper, we parametrically analyze the number and the order of the singularities in multi-mate-
rial wedges that are of relevance to microelectronic copper interconnect structures. We extend the derivation
of Seweryn and Molski made for homogeneous angular corners to multi-material wedges. We exhaustively
study the effects of varying both the wedge angle as well as the properties of the materials that constitute
the wedge. We modify the numerical procedure applied by Liu et al. to develop a consistent procedure for
calculating the stress intensity factor in the presence of multiple singularities. The analysis is then applied
to copper interconnect structures of modern microelectronic components to identify critical corners.

3. Asymptotic analysis of stress singularities in multi-material wedges

Consider a multi-material wedge in a polar coordinate system (r,h). There are two independent displace-
ment components ur and uh. The subscript i hereafter refers to the different materials present in the wedge.
The constant stress terms due to thermal expansion are included in the Appendix for completeness. The solu-
tion procedure is extended from Seweryn and Molski’s [14] formulation for the homogenous wedge problem.
The three strain components in polar coordinates are as follows:
ouri 1 ouhi uri 1 ouri ouhi uhi
er i ¼ ; eh i ¼ þ ; erhi ¼ þ  ð3Þ
ori ri ohi ri ri ohi ori ri
Assuming plane strain conditions, the stress components are then given by
rri ¼ Ki ðeri þ ehi Þ þ 2li eri ; srhi ¼ li erhi ; rhi ¼ Ki ðeri þ ehi Þ þ 2li ehi ð4Þ
The equations of equilibrium in polar coordinates are
orri 1 osrhi rri  rhi osrhi 1 orhi srh
þ þ ¼ 0; þ þ2 i ¼0 ð5Þ
ori r ohi ri ori r ohi ri
Combining above three sets of equations, the equation of equilibrium in terms of displacements can be written
as (Lame equations)
   
o ouri 1 ouhi uri 1 o ouhi 1 ouri uhi
ðKi þ 2li Þ þ þ  li  þ ¼0
ori ori ri ohi ri ri ohi ori ri ohi ri
    ð6Þ
1 o ouri 1 ouhi uri o ouhi 1 ouri uhi
ðKi þ 2li Þ þ þ þ li  þ ¼0
r ohi ori ri ohi ri ori ori ri ohi ri

The asymptotic solution is determined by assuming that the displacement has the form shown in Eq. (2). In the
present study, the oscillatory singularity is not considered because the stress intensity factor is defined based
only on the real part of the singularity. Further, as the comparisons with finite element analysis (presented
later in the paper) reveals, the complex part is negligibly small in the application considered in the present
paper. Substituting Eq. (6) into Eq. (5), we get
420 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430

ðji  1Þf 00 ðhi Þ þ ðji þ 1Þðk2  1Þf ðhi Þ þ 2ðk  ji Þg0 ðhi Þ ¼ 0;
ð7Þ
ðji þ 1Þg00 ðhi Þ þ ðji  1Þðk2  1Þgðhi Þ þ 2ðk  ji Þf 0 ðhi Þ ¼ 0;
where
8
< 3  mi for plane stress
ji ¼ 1 þ mi
:
3  4mi for plane strain
The general form of the two unknown functions f(h) and g(h) is given by [19]
fi ¼ Ai cosð1 þ kÞhi þ Bi sinð1 þ kÞhi þ C i cosð1  kÞhi þ Di sinð1  kÞhi ;
ji þ k ji þ k ð8Þ
gi ¼ Bi cosð1 þ kÞhi  Ai sinð1 þ kÞhi þ Di cosð1  kÞhi  C i sinð1  kÞhi
ji  k ji  k
The general solution for displacements and stresses can then be expressed as
uri ¼ rki ½Ai cosð1 þ kÞhi þ Bi sinð1 þ kÞhi þ C i cosð1  kÞhi þ Di sinð1  kÞhi ;
uhi ¼ rki ½Bi ðji  kÞ cosð1 þ kÞhi  Ai ðji  kÞ sinð1 þ kÞhi þ ðji þ kÞDi cosð1  kÞhi
 C i ðji þ kÞ sinð1  kÞhi =ðji  kÞ;
rri ¼ 2krk1
i li ½ðji  kÞAi cosð1 þ kÞhi þ ðji  kÞBi sinð1 þ kÞhi
þ ð3  kÞC i cosð1  kÞhi þ ð3  kÞDi sinð1  kÞhi =ðji  kÞ; ð9Þ
rhi ¼ 2krk1
i li ½ðji  kÞAi cosð1 þ kÞhi  ðji  kÞBi sinð1 þ kÞhi
þ ð1 þ kÞC i cosð1  kÞhi þ ð1 þ kÞDi sinð1  kÞhi =ðji  kÞ;
srhi ¼ 2krk1
i li ½ðji  kÞAi sinð1 þ kÞhi þ ðji  kÞBi cosð1 þ kÞhi
þ ð1  kÞC i sinð1  kÞhi  ð1  kÞDi cosð1  kÞhi =ðji  kÞ
For the case of a bi-material wedge bonded at both interfaces (Fig. 2(c)), the wedge is symmetric with
respect to x-axis; the symmetric (opening) loading will lead to the following boundary conditions:
uh1 ¼ srh1 ¼ 0 for h ¼ 0;
ð10Þ
uh2 ¼ srh2 ¼ 0 for h ¼ p
Continuity condition is enforced at the material interface
ur 1 ¼ ur 2 ; uh1 ¼ uh2 ; rh1 ¼ rh2 ; srh1 ¼ srh2 at h ¼ a ð11Þ
When the load is anti-symmetric (shear), the conditions become
ur1 ¼ rh1 ¼ 0 for h ¼ 0;
ð12Þ
ur2 ¼ rh2 ¼ 0 for h ¼ p
We apply the same continuity boundary condition at the interface
ur 1 ¼ ur 2 ; uh1 ¼ uh2 ; rh1 ¼ rh2 ; srh1 ¼ srh2 at h ¼ a ð13Þ
Another way to describe the above boundary conditions is to set the displacement field and stress field con-
tinuous at both interfaces

ur 1 ¼ ur 2 ; uh1 ¼ uh2 ; rh1 ¼ rh2 ; srh1 ¼ srh2 at h ¼ a ð14Þ

ur1 jh¼a ¼ ur2 jh¼2pa ; uh1 jh¼a ¼ uh2 jh¼2pa ; rh1 jh¼a ¼ rh2 jh¼2pa ; srh1 jh¼a ¼ srh2 jh¼2pa ð15Þ
The first set of boundary conditions is preferred because each assumption (symmetric or anti-symmetric
loading) will result in one definite value of the order of singularity. This will make the solution procedure easier
since the initial value could be randomly set in this case, and the resulting transcendental equations converge
much easier. Under both loading conditions, there are eight equations to be solved using eight boundary con-
ditions. For a nontrivial solution of constants: A1, B1, C1, D1, A2, B2, C2, D2, the determinant of the eigenvalue
 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430 421

problem must vanish, which results in a transcendental equation for k. Solving the transcendental equation will
yield all possible values of k. In the current study, an exhaustive search was used to identify all the solutions of k.
The results of k varying with l2/l1 are shown in Figs. 3 and 4. In both cases, the wedge angle a was set to be p/2,
the chosen j values are as shown in the figures. The largest possible value of k in this case is found to be 0.5. The
symmetric mode is dominant when l2/l1 > 1, and the anti-symmetric mode is dominant when l2/l1 < 1.
In the case of a wedge with more than two materials, with every interface is bonded (Fig. 2(d)), and i rang-
ing from 1 to n in Eq. (9), the boundary conditions that ensure the continuity of the displacement field and
stress field at all interfaces are
uri ¼ uriþ1 ; uhi ¼ uhiþ1 ; rhi ¼ rhiþ1 ; srhi ¼ srhiþ1 at h ¼ ai ð16Þ

0.4 Anti-symmetric
Symmetric

0.3
λ

0.2

0.1

0
10-3 10-2 10-1 100 101 102 103
μ2 / μ1

Fig. 3. k versus l2/l1 in a bi-material wedge with two bonded interfaces, j1 = 2, j2 = 1.5.

0.4

0.3
λ

0.2

0.1 Anti-symmetric
Symmetric

0
10-3 10-2 10-1 100 101 102 103
μ2 / μ1

Fig. 4. k versus l2/l1 in a bi-material wedge with two bonded interfaces, j1 = 1.8, j2 = 2.2.
422 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430

Thus we have 4n equations, and again for a nontrivial solution of constants: Ai, Bi, Ci, Di, i = 1, 2, . . . , n the
determinant of the eigenvalue problem must vanish. This gives rise to a transcendental equation for k. In gen-
eral there are up to two singularities in this case, and each singularity being associated with one mode. But,
now, the modes can no longer be classified as being symmetric or anti-symmetric.
In case of bi-material wedge with only one bonded interface (Fig. 2(b)), the boundary conditions that need
to be applied are different. The first set of boundary conditions are the continuity condition at the material
interface
ur 1 ¼ ur 2 ; uh1 ¼ uh2 ; rh1 ¼ rh2 ; srh1 ¼ srh2 at h ¼ 0 ð17Þ

The second set of boundary conditions are the traction free condition at the two free surfaces

rh1 ¼ 0; srh1 ¼ 0 at h ¼ a;
ð18Þ
rh2 ¼ 0; srh2 ¼ 0 at h ¼ b

The above boundary conditions result in an eigenvalue problem with eight equations, similar to previous dis-
cussion, from which k can be computed. The results indicate that the total number of singularities can be one,
two or three, suggesting a more complicated coupling of failure modes. Fig. 5 shows how the singularity
changes with respect to b, when a is fixed at p/2. The stresses are non-singular when b is under 0.35p, and
only one singularity exists when b is under p. The second singularity arises when b increases beyond p. The
singularity continues to exist as the sum of angles a and b approaches p. The case corresponding to eight
b = 3/2p with the second interface delaminated is described further in Fig. 6. In the figure, all singularities
are plotted as functions of l2/l1; it is seen that the first two singularities decrease as l2/l1 increases, while the
third singularity vanishes at l2/l1 = 1. The larger the difference between l2 and l1, the larger is the value of the
third singularity. The third singularity is probably due to out-of-plane shear generated for the particular com-
bination of l2 and l1. It should be noted that it is possible for the value of the largest singularity to be greater
than 0.5 depending on the chosen l2/l1 ratio. In Fig. 7 another case at a = p/4 is considered; the trends are
similar to those corresponding to a = p/2.
Based on the calculated order of singularities, the angular functions can be determined as a function of an
unknown scale parameter, which is incorporated into stress intensity factor. The angular functions depend not
only on material properties and geometries, but also on the loading conditions.

0.4

0.3
λ

0.2
Singularity 1
Singularity 2
Singularity 3
0.1 Singularity 3

0.4π 0.6π 0.8π 1.0π 1.2π 1.4π

β
Fig. 5. k versus b at the angular edge of a bi-material wedge with a = p/2, j1 = 2, j2 = 1.5, l1 = 1, l2 = 10.
 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430 423

0.7

Singularity 1
0.6 Singularity 2
Singularity 3

0.5

0.4
λ

0.3

0.2

0.1

0 -4
10 10-3 10-2 10-1 100 101 102 103 104
μ2 / μ1

Fig. 6. k versus l2/l1 in a bi-material wedge with a delaminated interface, a = p/2, b = 3p/2, j1 = 2, j2 = 1.5.

0.6

0.4
λ

Singularity 1
Singularity 2
0.2 Singularity 3

Singularity 3

0
0.8π 1.0π 1.2π 1.4π 1.6π 1.8π
β

Fig. 7. k versus b at angular edge of bi-material wedge with a = p/4, j1 = 2, j2 = 1.5, l1 = 10, l2 = 1.

4. Characterization of stress intensity factors

In the usual definition of the stress intensity factor, the angular function is set to one at the crack or mate-
rial interface. While in cases where multiple singularities and/or material interfaces are involved, a consistent
definition of the stress intensity factor is critical. Liu et al. and other studies have calculated stress intensity
factors by studying the variation of stress along prescribed interfaces, while arbitrarily setting the value of
the angular function to unity. That is, the simultaneous fit for the multiple stress intensity factors is through
an equation of the form
424 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430

K1 K2
s¼ þ ð19Þ
rk1 1=rk2

In this form it is difficult to make a relative comparison of the intensity factors since the two angular
functions have been assumed to simultaneously possess a value of unity. A further disadvantage of this pro-
cedure is that the value of the stress intensity factors depends on the interface which was chosen for the fit
(on which the angular functions are assumed to have unit values). To correct the ambiguity in the definition
of the stress intensity factor, in the current paper, we determine the angular functions as follows. We first
solve for the angular functions (to within an arbitrary constant) as the eigenvectors corresponding to the
eigenvalue problem defined by Eqs. (9)–(18) through which the multiple singularities are identified. The

Fig. 8. Cross-sectional TEM image of a six-layered dual damascene copper interconnect structure on a Si surface. The diameter of the
narrowest Cu via is 0.25 lm [6].

Fig. 9. A schematic illustration of the interconnect structure showing susceptible corners (after [21]).
 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430 425

angular functions are then normalized such that the maximum magnitude of the shear stress remains one.
Using the above procedure, the angular functions can be determined for any combination of wedge angle,
material properties and load conditions. More importantly, the values of the stress intensity factors will
remain constant for any choice of the numbers and locations of material interfaces or the origin used to
define the angular function.

Ta SiO2

SiO2 Ta
2. SiO2(90°)/Ta(270°)
1. Ta(90°)/SiO2(270°)

Cu Ta

Ta Cu
3. Cu(90°)/Ta(270°) 4. Ta(90°)/Cu(270°)

SiC SiC

SiO2 Ta Cu Ta

5. SiO2(90°)/Ta(90°)/SiC(180°)
(90°)/Ta(90°)/SiC(180° 6. Cu(90°)/Ta(90° 180°)
Cu(90°)/Ta(90°)/SiC(180°

SiC SiO2
Ta
Ta
Cu SiC

SiC(90°)/Ta(90°)/Cu(180°
7. SiC(90°)/Ta(90°)/Cu(180°) 8. SiO2(90°)/SiC(90°
(90°)/SiC(90°)/Ta(180°)
0°

Fig. 10. The material sets constituting the eight different corners identified in the interconnect structure.

Table 1
Properties of materials used in the copper interconnect structure
E(Pa) Possion’s ratio a (lm/m C)
SiC 4.0E + 11 0.16 4.50
Cu 1.17E + 11 0.34 16.56
Ta 1.86E + 11 0.35 6.48
SiO2 7.0E + 10 0.17 0.55

Table 2
The calculated strengths of singularities obtained through asymptotic analysis at the six corners of the copper interconnect structure
Singularity (I) Singularity (II)
Corner 1 0.143 0.0376
Corner 2 0.116 0.0217
Corner 3 0.0715 0.0980
Corner 4 0.0606 0.00969
Corner 5 0.200 0.00957
Corner 6 0.0872 0.00312
Corner 7 0 0
Corner 8 0.247 0.0869
426 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430

2 1
frθ2
fθ2
1.5

0.5
Angular functions

Angular functions
1

0.5 0

-0.5

-0.5 frθ1
fθ1
-1 -1
0 0.5π 1π 1.5π 2π 0 0.5π 1π 1.5π 2π
θ ( Corner 1 ) θ ( Corner 1 )

1 1

0.5
0.5

0
Angular functions

Angular functions

0
-0.5

-1
-0.5

-1.5 frθ2
frθ1 -1 fθ2
-2 fθ1

-2.5 -1.5
0 0.5π 1π 1.5π 2π 0 0.5π 1π 1.5π 2π
θ ( Corner 5 ) θ ( Corner 5 )

0.5 1 frθ2
fθ2
0
Angular functions
Angular functions

0.5
-0.5

-1 0

-1.5

frθ1 -0.5
-2 fθ1

-2.5 -1
0 0.5π 1π 1.5π 2π 0 0.5π 1π 1.5π 2π
θ ( Corner 8 ) θ ( Corner 8 )

Fig. 11. The angular functions at Corners 1, 5 and 8.

 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430 427

Thus, the stress intensity factors are obtained through a fit of the following form:
K1 K2
s¼ frh1 ðhÞ þ frh2 ðhÞ ð20Þ
r k1 1=rk2
where, the fit is carried out over all combinations of h and r values. To compute the generalized stress intensity
factors, asymptotic methods are most commonly adopted, in which the structure is analyzed using numerical
schemes such as the finite element method and the stress intensity factors are determined by comparing the
asymptotic solution with the numerical solution.
At least two methods exist for finite element analysis of singular stresses. One is to use special finite ele-
ments at the interface corner, where the order of the singularity is incorporated into the shape function of
these elements. Alternatively standard finite elements can be used over the entire domain. A very fine mesh
is generally needed near the singular point in order to capture the asymptotic behavior. In the current study,
a standard finite elements are used along with the asymptotic analysis and a convergence study is performed to
ensure accurate characterization of the stress intensity factors [20].

5. Application to interconnect stress analysis

In general, the interconnect structure is a three-dimensional network of conductor lines (copper) embedded
in a dielectric matrix, fabricated on the silicon wafer surface that contains the active devices (Fig. 8). The con-
ductor lines are on one of several layers that are linked by vias. A silicon nitride or a silicon carbide film that
serves as the ‘‘etch-stop’’ layer caps the interconnect structure. Tantalum or tantalum nitride is used as a bar-
rier to prevent diffusion of copper into the dielectric layer. The structure shown in Fig. 9 is analyzed in the
present study.
There are in total eight different geometric corners in the interconnect structure shown in Fig. 9 [21]. In gen-
eral different materials combine to form the corners. But at some corners, while the materials constituting the
corner are the same, the angles at which they combine vary. These combinations are illustrated in Fig. 10. The
material properties used in the present analysis are listed in Table 1, the asymptotic solution to all the singu-
larities at the corners identified in Fig. 9 are listed in Table 2. Clearly, Corners 1, 5 and 8 are the weakest as
identified by the asymptotic analysis. Plots of the angular functions at these three corners are shown in Fig. 11.
The stress intensity factor depends on the elastic properties of the materials as well as the loading conditions.
The applied loading condition was as follows. The temperature at the top surface was fixed at T = 20 C and
that at the bottom surface was set at T = 100 C. The remaining surfaces were enforced to be adiabatic. The
thermal-stress analysis was performed using ANSYS. The steady state temperature field solution and the mesh
for singular stress modeling are shown in Fig. 12. The thickness of the barrier and the etch stop was assumed to
be t = 0.1 lm everywhere. The overall length of the model was 5 lm, and the height was 4 lm.

Fig. 12. Temperature distribution in interconnect structure determined using finite element analysis. The bottom surface was held at
100 C and the top surface at 20 C.
428 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430

Corner 1 was selected for further detailed analysis. This corner is a case of bonded bi-material wedge
(see Fig. 2). The stress intensity factors at Corner 1 were calculated using the procedure described in
the previous section; the obtained values for the intensity factors were K1 = 1.28 · 106 and K2 = 6.39 ·
106. In Fig. 13, we compare the finite element solution against the fit obtained by using either one or
two singularities. We do this to assess the impact of including multiple singularities in describing the stress
field. As can be seen from the figure, the angular variation of the stress based on the dominant singularity
is a poor fit to the finite element solution. The fit is accurate when the angular functions associated with
both the singularities are included when describing the stress field. This suggests that it would be difficult, if
not impossible, to capture the angular variation of stress at the singular point when only the dominant
singularity is included. However, the variation in the radial direction (Fig. 13) is not significantly affected
by the inclusion of both singularities. This is because at h = 0 (the angle which was chosen for the
comparison) the angular functions have identical value. At other h values, the inaccuracy in the angular

-6.5

-7
Shear stress (x10 7 )

-7.5

-8 Fit using one singularity

Fit using two singularities
FE result

0 0.5π 1π 1.5π 2π
(a) θ (r=0.01, Corner 1)

-0.7

-0.8

-0.9
Shear stress (x10 )
8

-1

-1.1

-1.2
Fit using one singularity
Fit using two singularities
-1.3 FE result

-1.4
0.002 0.004 0.006 0.008 0.01
(b) r (distance from corner 1 at θ=0 o)

Fig. 13. Comparison of fit with FE result at Corner 1. (a) Stress as function of h at r = 0.01. (b) Stress as function of r at h = 0.
 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430 429

-0.7

-0.8

-0.9

Shear stress (x10 )

8 -1

-1.1

-1.2
Fit using one singularity
Fit using two singularities
-1.3 FE result

-1.4
0.002 0.004 0.006 0.008 0.01
r (distance from corner 1 at θ=45 )
o

Fig. 14. Comparison of fit with FE result at Corner 1 as function of r at h = 45.

variations plays a critical role in determining the accuracy of the stress in the radial direction. This is
illustrated in Fig. 14.

6. Summary

A systematic study of the nature of the singular stress field at bi-material and multi-material wedges was
presented in this paper. The analytical procedure for solving the order of singularities was extended from
Seweryn and Molski’s [14] formulation for homogenous wedge problems. The total number of singularities
as well as the order of the singularities was extensively investigated for a multi-material wedge with bonded
or delaminated interfaces as a function of wedge angle and material properties. It was found that in general
that the singularities range in number between two and three in multi-material wedges, and that their order
varies between 0 and 0.7 over three orders of magnitude change in shear modulus.
To relate the singular stress field to failures in microelectronic interconnect structures, the angular stress
variation and the stress intensity factors were calculated. The analysis showed that when only the leading
singularity is considered, the resulting stress field is considerably inaccurate in describing the variation of
stress in the angular direction. Therefore, accurate descriptions of stress state at singular corners require
the consideration of multiple singularities and an accurate derivation of the angular functions corresponding
to these singularities. Furthermore, a consistent definition of the stress intensity factors that is independent
of the angle at which the numerical fit to stress is carried out will enable the development of meaningful
failure criteria.

Appendix A. Constant stress term due to temperature change

To get the constant stress under thermal loading, k is chosen to be 0.

uri ¼ ri ½Ai cos 2hi þ Bi sin 2hi þ C i 

uhi ¼ ri ½Bi cos 2hi  Ai sin 2hi þ m2i Di 
rri ¼ li ½2Ai cos 2hi þ 2Bi sin 2hi  2ð1  m2i ÞC i  ðA1Þ
rhi ¼ li ½2Ai cos 2hi  2Bi sin 2hi  2ð1  m2i ÞC i 
srhi ¼ li ½2Ai sin 2hi þ 2Bi cos 2hi 
430 Y. Luo, G. Subbarayan / Engineering Fracture Mechanics 74 (2007) 416–430

The thermoelastic relation is

2li ðer  
ai DT Þ ¼ rr  ð3  ji Þðrr þ rh Þ=4
2li ðeh  
ai DT Þ ¼ rh  ð3  ji Þðrr þ rh Þ=4 ðA2Þ
li erh ¼ srh
where
3  mi
ji ¼ ; for plane stress
1 þ mi
ji ¼ 3  4mi ; for plane strain
The strains are expresses in displacements as
ou u 1 ov 1 ou ov v
er ¼ ; eh ¼ þ ; erh ¼ þ  ðA3Þ
or r r oh r oh or r
Substituting the stresses and strains into the thermoelastic relation, the displacements are found by integration
u 1
¼aDT þ Ai cos 2hi þ Bi sin 2hi þ ð1  m2i ÞC i  ð1  m2i ÞC i ji
r 2 ðA4Þ
v
¼ Ai sin 2hi  Bi cos 2hi
r
Using the same boundary conditions as used in calculating the orders of singularity, parameters Ai, Bi, Ci
could be found.

References

[1] Wilson, R. Failures plague 130-nanometer IC processes. Available from: [Link]

[2] Bogy DB. Edge-bonded dissimilar orthogonal elastic wedges under normal and shear loading. J Appl Mechanics 1968;35:460–6.
[3] Bogy DB. Two edge-bonded elastic wedges of different materials and wedge angles under surface tractions. J Appl Mechanics
1971;38:377–86.
[4] Hein VL, Erdogan F. Stress singularities in a two material wedge. Int J Fract Mech 1971;7:317–30.
[5] Rice JR. Elastic fracture mechanics concepts for interfacial cracks. J Appl Mech 1988;55:98–103.
[6] Tu KN. Recent advances on electromigration in very-large-scale-integration of interconnects. J Appl Phys 2003;19:5451–73.
[7] Munz D, Yang YY. Stresses near the edge of bonded dissimilar materials described by two stress intensity factors. Int J Fract
1993;60:169–77.
[8] Theocaris PS. The order of singularity at a multi-wedge corner in a composite strip. Int J Engng Sci 1974;12:102–20.
[9] Munz D, Yang YY. Stress intensities at the interface in bonded dissimilar materials under mechanical and thermal loading. J Appl
Mech 1992;59:857–61.
[10] Munz D, Fett T, Yang YY. The regular stress term in bonded dissimilar materials after a change in temperature. Engng Fract Mech
1993;44:185–94.
[11] Munz D, Yang YY. Stress singularities in a dissimilar materials joint with edge tractions under mechanical and thermal loadings. Int
J Solids Struct 1997;34:1199–216.
[12] Williams ML. Stress singularities resulting from various boundary conditions in angular corners of plates in extension. ASME J Appl
Mech 1952;19:526–8.
[13] Pageau SS, Gadi KS, Biggers SB, Joseph PF. Standardized complex and logarithmic eigensolutions for n-material wedges and
junctions. Int J Fract 1996;77:51–76.
[14] Seweryn A, Molski K. Elastic stress singularities and corresponding generalized stress intensity factors for angular corners under
various boundary conditions. Engng Fract Mech 1996;55:529–56.
[15] van Vroonhoven JCW. Stress singularity in bimaterial wedges with adhesion and delamination. Fatigue Fract Engng Mater Struct
1992;15:159–71.
[16] Fan XJ, Wang HB, Lim TB. Investigation of the underfill delamination and cracking in flip-chip modules under temperature cyclic
loading. IEEE Trans Comp Pack Technol 2001;24(1):84–91.
[17] Chue CH, Chen TH, Lee HT. Effect of stiffness and thickness ratios on popcorn cracking in IC packages. IEEE Trans Comp Pack
Technol 2003;26(2):340–8.
[18] Liu XH, Suo Z, Ma Q. Split singularities: stress field near the edge of a silicon die on a polymer substrate. Acta Mater 1999;47:67–76.
[19] Parton VA, Perlin PI. Mathematical methods of the theory of elasticity. Moscow: Mir Publishers; 1984.
[20] Akisanya AR. On the singular stress field near the edge of bonded joints. J Strain Anal 1997;32:301–11.
[21] Li BZ, Sullivan TD, Lee TC, Badami D. Reliab challenges copper interconnects 2004;44:365–80.