，，

State Key Laboratory of Mechanics and Control of Mechanical Structures，College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，P.R.China

Abstract:Based on the bulk free energy density and the degenerate mobility constructed by the quartic double-well potential function，a phase field model is established to simulate the evolution of intragranular microvoids due to surface diffusion in a stress field.The corresponding phase field governing equations are derived.The evolution of elliptical microvoids with different stresses Λ，aspect ratios β and linewidths is calculated using the mesh adaptation finite element method and the reliability of the procedure is verified.The results show that there exist critical values of the stress Λc，the aspect ratio βc and the linewidth of intragranular microvoids under equivalent biaxial tensile stress.When ，the elliptical microvoids are instable with an extending crack tip.When Λ<Λc，β<βc or ，the elliptical microvoids gradually cylindricalize and remain a stable shape.The instability time decreases with increasing the stress or the aspect ratio，while increases with increasing the linewidth.In addition，for the interconnects containing two elliptical voids not far apart，the stress will promote the merging of the voids.

Key words：phase field method；stress migration；surface diffusion；finite element method；intragranular microvoid

0 Introduction

With the continuous progress of microelectron?ics technology，the size of typical integrated circuits has been drastically reduced.The miniaturization poses many challenges，especially when the line?width of the interconnect lines reaches the submi?cron or nanometer level.The thermal mismatch stress induced by the difference of thermal expan?sion coefficients between the interconnect lines and the passivation layer becomes more pronounced.As the thermal mismatch stress increases，the micro?void nucleates，grows，migrates and changes its shape.When the microvoid grows into a crack，it will cut off the entire interconnect line and cause a break.Therefore，it is of practical importance to un?derstand the evolution of stress induced microvoids.

The early research on the failure of intercon?nect lines is mainly by experimental observation.Hull and Rimmer［1］observed the growth of inter?granular microvoids under different stress states within metallic materials and concluded that the mi?crovoids may be generated due to stress migration.Blech and Herring［2］measured the gradient stresses generated by electromigration in aluminium thinfilm interconnect lines on silicon nitride substrates using an X-ray instrument.Wilson et al.［3］measured the stresses in copper interconnect lines with line?width of 50—500 nm using X-ray diffraction（XRD）and found that the stresses increased with decreasing linewidth and that the grain structure of the interconnect lines had a significant effect on the stresses and stress evolution.

Experimental observations form the basis for theoretical analysis and numerical simulation.The first theoretical investigation of the stress driven morphological instability in solids is given by Asaro and Tiller［4］.They discussed the linear stability of a planar surface separating a stressed，two-dimension?al semi-infinite solid from a fluid and found that the planar surface was unstable to small disturbances with wavelength greater than a characteristic length.The same result was later rediscovered by Srolo?vitz［5］.The starting point of these theoretical studies is the chemical potential，which consists of strain en?ergy density and surface energy density.Based on the above theory，Sun and Suo et al.［6-7］developed a theoretical model of microvoids evolution due to sur?face diffusion and derived the thermodynamic poten?tial with shape function to calculate the rate of micro?voids evolution using the variational principle.Wang and Suo［8］analyzed the shape instability threshold and the forming time of crack tip using conformal mapping method.Wang and Li［9-13］theoretically ana?lyzed the stability of two-and three-dimensional in?tra-and interg?ranular microvoids containing inter?nal pressure and providing a theoretical basis for the numerical simulation.In addition，the elliptical voids and inclusions morphological evolution under a gradient stress field［14-16］and the circular void mor?phological evolution under high current density was also analyzed theoretically［17-18］.

Due to the complexity of the material system and the limitation of the two methods，it is difficult for the experimental observation and theoretical analysis to completely solve the problem of micro?void evolution mechanisms.Therefore，numerical simulation becomes an effective method.The com?monly used numerical methods to analyze the evolu?tion of microvoids are the sharp interface method and the phase field method.The sharp interface method is a mature numerical simulation method.So far，the evolution of microvoids in a stress field is mostly based on the sharp interface model.For ex?ample，based on the weak statement proposed by Sun and Suo［19］，the morphological evolution pro?cess of microvoids in a stress field，electric field and gradient stress field［20-23］had been studied in the past few years.Moreover，the mesh adaptation finite ele?ment method was also used to study the evolution of intra-and inter?granular microcracks under electro?migration and stress migration［24-25］.

However，when it comes to complex interface problems（e.g.topological changes）or expansion to three dimensions，the sharp interface method for moving boundaries will become more difficult to handle.In recent years，the rapidly developing phase field method has provided a powerful tool for simulating the evolution of microvoids［26-28］.The phase field method overcomes the limitations of the sharp interface method.The morphological informa?tion of microvoids is implicitly included in the phase field equation，which avoids the explicit tracking of its interface location and thus has more advantages when dealing with more complex moving boundary problems.It is also easier to extend the phase field method from two dimensions to three dimensions.

The bulk free energy plays an important role in the phase field model and directly determines the construction and composition of the total free energy of the system［26］.At present，for the analysis of the evolution of microvoids morphology in the intercon?nect lines，there are two principal forms of the bulk free energy density，which are the quartic doublewell potential function and the nonsmooth doubleobstacle potential function.Using the quartic doublewell potential function，Yeon et al.［29］presented a phase field model for the surface corrugation of elas?tically stressed films where surface diffusion is a dominant mass transport mechanism.The method was also used to analyze such problems as void mi?gration driven by temperature gradients［30］，electro?migration-induced thermal grooving evolution on grain boundaries［31-32］and non-isothermal coales?cence in polycrystalline sintering［33］.Meanwhile，us?ing the nonsmooth double-obstacle potential func?tion，Barrett et al.［34］derived the phase field equa?tions for stress migration induced surface diffusion based on degenerate mobility and solved the govern?ing equations using the finite element method.Bhate et al.［35-36］proposed a phase field model and used the mesh adaptation method to capture the interfacial layer information to analyze the morphological evo?lution of microvoids under electro-and stress migra?tion induced surface diffusion and bulk diffusion，re?spectively.Li et al.［37］analyzed the effects of electri?cal conductivity and anisotropic surface energy on the migration and evolution of inclusions under the electric field.Subsequently，this method was ex?tended to research the morphological evolution of in?clusions in piezoelectric films under mechanical and electric loads［38］.

For the phase field method using the nons?mooth double-obstacle potential function，the order parameter only needs to be calculated within the in?terface layer，which improves computational effi?ciency.But it brings some disadvantages.The first is the need to introduce algorithms for capturing the interface layer，such as dynamic mesh method or narrowband interface method，which needs to cap?ture the node information of the interface layer and update the mesh after a certain time of calculation.Secondly，when coupling the stress field，different mesh densities need to be consistence with different stress conditions.However，the quartic double-well potential function is capable of continuously and syn?chronously coupling the external fields.And its com?putational effort and efficiency are rapidly improving with the development of adaptive grid technolo?gy［39］.The evolution of microvoids in a stress field based on a phase field model with a quartic doublewell potential function and degenerate mobility has not been reported in the literature.Therefore，in this paper，the evolution of elliptical microvoids in a stress field is investigated based on the modified phase field method using the quartic double-well po?tential function and the degenerate mobility.

The plan for the rest of the article is as fol?lows：in the first part，the phase field model and the governing equations are introduced；in the second and third parts，the results and discussions are pre?sented，and the effects of different stresses，aspect ratios and linewidths on the evolution of elliptical microvoids are investigated.

1 Phase Field Model and Govern?ing Equations

1.1 Phase field model

The phase field model is illustrated in Fig.1 and the interconnect line is idealized as an isotropic lin?ear elastic model of a two-dimensional single crys?tal，which is subjected to biaxial tensile stressesσxandσyon its external boundary.aandh0are the long semi-axis and short semi-axis of the elliptical micro?void，respectively.LandHare the length and the width of the interconnect line.As shown in Fig.1，the order parameters?=+1 and?=-1 represent the two equilibrium phases with the minimum bulk free energy.In this paper，they correspond to the metal?lic conductor and the microvoid respectively.For physically reasonable values of the material parame?ters，we may divide the regionRinto three por?tions.LetR+denote the region occupied by solid material（?=+1），R-denote the region occupied by void（?=-1），andRIdenote the interface re?gion.We denote the two extreme contours，associ?ated with?=-1 and?=+1 byΓ-andΓ+，re?spectively.The bisecting contour?=0 is denoted byΓ.srepresents the local coordinate direction along the hole tangent direction，andrthe local co?ordinate direction perpendicular to the hole tangent direction.

Fig.1 Schematic diagram of phase field model of a microvoid in a stress field

This paper adopts some basic assumptions as follows：

（1）The isotropic surface diffusion is the only material transport mechanism in the microvoids evo?lution.

（2）The model studied in this paper is a plane strain problem.

（3）The stress is uniformly distributed on the stress boundary and the deformation of the line satis?fies the line elastic deformation.

1.2 Governing equations

For the surface diffusion mechanism，the order parameter of the conserved field needs to satisfy the law of mass conservation and energy dissipation.Ac?cording to the Gurtin principle of microcode equilib?rium［40］，the expression for the total free energy func?tion of the systemFin a stress field is

whereγsis the surface energy per unit area，?the parameter controlling the thickness of the interface，the dimensionless parameter［41］andW(ε，?) the strain energy density of the system，which can be expressed as

whereC(?) is the second-order elasticity tensor with order parameter，σthe stress tensor，εthe strain tensor，and the bulk free energy density func?tionfb(?) is

The chemical potentialμis defined by the varia?tional derivative ofFwith respect to?

where the factor 2 is the change of the order parame?ter from +1 to -1，i.e.，the transition from solid material to voids across the interfacial layer，F(xiàn)the free energy functional density of the system andΩthe volume of the atom.

The driving force of atomic migrationfis relat?ed to the gradient ofμ

where the negative sign indicates that the direction of driving force is opposite to the chemical potential gradient，so that atoms migrate from the high chemi?cal potential region to the low chemical potential re?gion.

The gradients of strain energy and surface cur?vature cause material to be deposited or removed lo?cally，and the surface motion will occur where un?balanced flux exists.The number of atoms per time crossing a unit length on the surface is determined by

whereB=is the dimensionless parame?ter［41］andDsthe atomic mobility for surface diffu?sion.In addition，M=M(?) denotes the degener?ate coefficient［42］associated with?，shown as

Based on the law of conservation of mass，the derivative of?with respect to timetis related to the divergence of the diffusion fluxJ，shown as

Coupling Eq.（4）to Eq.（8）yields the phase field equations combining the thermodynamic and kinetic laws，which are the modified fourth-order nonlinear Cahn-Hilliard equations

Based on the phase-field governing equation，the morphological evolution and migration of void depend on the gradient of strain energy along the surface.At each time，we solve the distribution of strain energy for a given void morphology.

The equations to solve the stress field are

whereG(?) andλ(?) are the shear modulus func?tion and the Lame function of the order parameter，respectively.trεis the trace of the second?order strain tensor andIthe unit tensor.In the microvoid region，the shear modulus and the Lame constant are zero.In the solid region，the shear modulus isG0，and the Lame constant isλ0.

Based on the governing equations of the phase field model，the corresponding program is devel?oped through the open-source finite element frame?work Moose［43］.For the convenience of description，the dimensionless process isThen，the dimension?less phase field method governing equation can be written as

For each time step Δτ，the calculations proceed as follows：

（1）Solve the elasticity problem on the current configuration，obtaining the strain energy density at all of the element nodes.

（2）Compute the new configuration by itera?tive calculation in space to obtain the values of the order parameters that satisfy the residual require?ments.

（3）Update the adaptive time step.

2 Numerical Simulation and Dis?cussion

It is known from diffusion theory that atoms mi?grate from a region of high chemical potential to a re?gion of low chemical potential due to the presence of chemical potential differences.The chemical poten?tial of each point on the surface of the intragranular microvoid is related to the external stress，aspect ra?tio and linewidth.Therefore，in this section，the phase field method proposed in this paper is used to analyze the evolution of intragranular microvoid un?der equivalent biaxial tensile stress with different stresses，different aspect ratios and different line?widths.

2.1 Effect of stress field

Fig.2 gives a graphical representation of the evo?lution of a circular void subjected toy-directional tensile stress，the stress nephogram of the model at the initial moment and the schematic diagram of the adaptive grid at the initial time for a quarter model.It can be seen from Fig.2 that the curvature of each point on the microvoid surfaceκis equal at the ini?tial time andκC=κB=κA.However，the strain en?ergy density of each point is different because of the stress concentration.We can find thatWC>WA>WB，which indicate that the atoms will migrate from pointCto pointAand pointB.When the stress is relatively small（Fig.2（a）），the circular microvoid will become elliptical microvoid.At the same time，the curvature of the microvoid becomes differentκC>κB>κA，which impedes further ellipticaliza?tion of the shape.And the circular microvoid will re?main a stable elliptical void.When the stress is rela?tively high（Fig.2（b）），the strain energy dominates，and the instability is formed by the gradual forma?tion of crack tips at the end of the void due to the at?om migration.As can be seen from Fig.2，there ex?ists a critical stressΛc.When the stress is larger than the critical stress，the microvoid will collapse into a crack.In addition，F(xiàn)ig.2（d）gives the sche?matic diagram of the adaptive grid at the initial time for a quarter model.The adaptive grid is updated synchronously with the evolution of the microvoid.The mesh refinement level of the transition between the surface layer is 4—6.

Fig.2 Evolution of intragranular circular microvoid，the initial stress nephogram for Λ=0.4,σx=0 and the schematic diagram of the adaptive grid for a quarter model

Fig.3 shows the variation of the critical stress as a function of the stress ratioσx/σy.Below the curve，the voids remain stable；above the curve，the voids are unstable and the crack tips will appear.It can be seen that the solution of the phase field method due to surface diffusion in the stress field is basically consistent with the curve of the theoretical solution.Therefore，the algorithm in this paper is reliable.

Fig.3 Λc as a function of the stress ratio σx/σy

Fig.4 shows the evolution of intragranular mi?crovoid forβ=2，=20 as the stress increases（Λ=0.1，0.3）.As shown in Fig.4，the evolution of elliptical microvoid under different stresses is basi?cally the same as that of the circular microvoid.The only difference between the circular microvoid and the elliptical microvoid is that the curvature of the microvoid surface is different at the initial time.For the elliptical microvoid，atoms migrate from the rel?atively flat surface to the tip（Fig.4（a））.That is，the atoms migrate from pointAto pointBand then from pointBto pointC.When the stress is relative?ly small，the elliptical microvoid will gradually cylin?dricalize and maintain a stable form.But，the ellipti?cal microvoid will form crack tips at both ends and gradually destabilize when the stress is large.There?fore，there exists a critical stressΛcas shown in Fig.5.In addition，it can be seen from Fig.5 that the critical stress decreases with increasing the aspect ra?tio.The elliptical microvoid will destabilize and form a cusp whenΛ≥Λc.The elliptical microvoid will cylindricalize and remain stable whenΛ<Λc.Fig.6 shows the variation of instability time with stress.From Fig.6，it can be seen that the instabili?ty time decreases with the increase of the aspect ra?tio and the stress.WhenΛ>0.4 andβ>3，the mi?crovoid will destabilize rapidly.Fig.7 shows the evo?lution process of the microvoid aspect ratio with time under the two kinds of external loads in Fig.4.When the microvoid is unstable，the evolution of the aspect ratio with time can be seen as the evolu?tion of cracked tip growth with time.From Fig.7，it can be seen that whenΛ=0.3 andτ>1.5，the cracked tip expands rapidly and eventually leads to the void destabilization.

Fig.4 Evolution of the intragranular microvoid for β=2,=20

Fig.5 Λc as a function of β(=20)

Fig.6 τs as a function of Λ

Fig.7 β as a function of τ

2.2 Effect of the initial aspect ratio

Fig.8 shows the evolution of the intragranular microvoid forΛ=0.2，=20 as the initial aspect ratio increases（β=2，3）.As shown in Fig.8，there exists a critical aspect ratioβc.Whenβ≥βc，the el?liptical microvoid will form a cracked tip and cause instability.Whenβ<βc，the elliptical microvoid will gradually cylindricalize and maintain a relatively stable shape.The critical aspect ratio decreases rap?idly with increasing the stress（Fig.9）.Fig.10 shows the variation of the instability time with increasing the aspect ratio.It can be seen that the instability time decreases with increasing the aspect ratio or the stress.In addition，the instability time is approxi?mately 0 whenβ≥3.5，which means that the mi?crovoid will be unstable instantly.

Fig.8 Evolution of intragranular microvoid for Λ=0.2,=20

Fig.9 βc as a function of Λ(=20)

Fig.10 τs as a function of β

2.3 Effect of linewidth

Fig.11 shows the evolution of the intragranular microvoid forΛ=0.2，β=2 as the linewidth de?creases（=20，8）.As displayed in Fig.11，there exists a critical linewidth，the ellip?tical microvoid will form a cracked tip and be insta?bility.When，the elliptical microvoid will gradually cylindricalize and maintain a relatively sta?ble shape.As shown in Fig.12，the critical line?width increases as the aspect ratio increases.And the effect of linewidth can be eliminated when the as?pect ratio is greater than 3.6，which means that re?gardless of the linewidth，the microvoid will be uns?tabe.Fig.13 shows the variation of the instability time with increasing the linewidth.It is obvious that the instability time decreases as the aspect ratio in?creases and increases as the linewidth increases.When the linewidth is relatively small，the micro?void will quickly become unstable.

Fig.11 Evolution of the intragranular microvoid for Λ=0.2,β=2

Fig.12 as a function of β(Λ=0.2)

Fig.13 τs as a function of

In addition，the interconnect line does not al?ways contain just one void，and sometimes contains two or more voids.Fig.14 shows the morphological evolution of two intragranular microvoids forΛ=0.1，β=2.The vertical distance between the cen?ters of the two elliptical voids is 3.The atoms mi?grate from the region with high strain energy to the region with low strain energy（Fig.14（b）），and even?tually the two elliptical voids evolve into symmetri?cal“heart-shaped”voids on top and bottom.That is，the tensile stress can promote the merging of two voids that are not far apart.Therefore，it is un?derstandable that if the microvoids in the intercon?nect line evolve and combine to form larger voids，it is more likely to lead to the formation of micro?cracks and the open circuit failure of the interconnect line.

Fig.14 Evolution of two intragranular microvoids and the final stress nephogram

Because of the interaction between multiple voids，the interconnect line with multiple micro?voids is a complex system，which will be analyzed in detail in the next work.The phase field simula?tion of the double-void problem here is only to dem?onstrate the applicability of the method to the multivoid problem.Due to the microscale of the intercon?nect line，especially the harsh multi-physical field coupling of its service conditions，the variability of its service behavior is caused，which greatly increas?es the complexity of the characterization and evalua?tion of its service behavior.In this paper，we focus on the phase field simulation of the void evolution under stress induced surface diffusion in the intercon?nect line.The research on the failure behavior of multiple voids under multi-physical field will be car?ried out in the following work.

3 Conclusions

The governing equations of the phase filed method are derived to simulate the evolution of mi?crovoids due to surface diffusion in a stress field.The reliability of the procedure is verified.And the evolution of the elliptical microvoids under different stresses，different aspect ratios and different line?widths is studied in detail.The main conclusions are as follows：

（1）Due to surface diffusion induced by stress migration，the elliptical microvoids have two trends：Firstly，the elliptical microvoids gradually cylindricalize and maintain a stable form；secondly，the elliptical microvoids form a crack tip and will col?lapse into a crack.

（2）There exist critical values of the stress fieldΛc，the initial aspect ratioβcand the linewidth.WhenΛ≥Λc，β≥βcor≤，the intragranular microvoid can form crack tips and have a tendency to cause a failure of interconnects；conversely，the in?tragranular microvoid will gradually cylindricalize.

（3）The higher of the stressΛ，the larger of the aspect ratioβor the smaller of the linewidth，the more easily to cause instability of the intragranu?lar microvoid in interconnects，and the smaller time is required for destabilization.