The thermal fatigue plays an important role during of degradation of interconnection compartments of power electronic devices. The temperature variations resulting from the power cycling has as consequences the stresses and plastic deformations that can affect the microstructure of the materials at the interconnection interfaces of upper metallic parts. Wires and metallization layers more solicited than silicon layers lead to the distortion of material interfaces when the temperature increases, leading to the deformation or degradation of the material surfaces. This will decrease the composite life and leads to an accelerated degradation. The arrangement of grains and grain boundaries is key to understanding the microstructure of metals and composites. When subjected to thermal and mechanical stresses, the variation in surface energies between adjacent grains, confined by the grain boundary, can cause the grains to separate. This phenomenon occurs due to the thermal and mechanical deformation of the grain boundary and the grain groove profile. Such occurrences are commonly observed in the bonding wires utilized in electronic devices.

The studies [1,2,3] have focused on examining the impact of microstructure and physicochemical properties on degradation processes. In literature [4,5,6], three effects were investigated. The first two effects examined the influence of bonding procedures and temperature on crack formation and the microstructure of the interconnection zone. Meanwhile, the third effect explored the relationship between material purity, grain size, and hardness during cycling. The metallization layer, typically around 5μm thick, deposited on the chips undergoes significant distortion compared to materials like silicon when exposed to high temperature. This distortion results in substantial tensile and compressive stresses, leading to notable inelastic strains [7]. It has been reported that thermomechanical cycling can cause two main types of degradation on the topside of power chips: metallization reconstruction and degradation of bonding contacts [7,8,9]. It is assumed that during cyclic aging, a progressive effect of condensation-evaporation occurs, leading to structural degradation and grooving of the film. However, the precise mechanism of this degradation is not yet fully understood, and further efforts are required to better comprehend the effects of stress parameters on the degradation of contacts between metallization and bond wire. This involves finding a mathematical solution to describe the formation of grain boundary grooving in polycrystalline thin films. Several solutions to this mathematical problem have been proposed in the literature [10,11,12,13,14,15,16,17,18,19,20]. In 1957, Mullins [10] conducted a study on the thermal effect on the profile of grain boundary grooving, laying the foundation for subsequent research on this phenomenon [13,14,15,16,17,18,19,20]. Various studies have focused on the development of this phenomenon, particularly exploring evaporation-condensation, surface diffusion, and formulating the mathematical problem that describes the profile of grain boundary grooving [10,11,12].Some authors [22] tried to adapt integrable nonlinear evolution equations related to the well-known linearizable diffusion equation to derive a new integrable nonlinear equation which models the surface evolution of anisotropic material accompanying the action of evaporation-condensation and surface diffusion [22].

A multiple integration technique allowing to solve high-order diffusion equations was proposed by Hristov [23] based on multiple integration procedures by applying the heat-balance integral method of Goodman and the double integration method of Volkov. Hristov [24] presented a solution for the linear diffusion models of Mullins’ thermal grooving [10,11,12].

Fourth-order diffusion equations are commonly encountered in various applications, including surface diffusion on solids [10,11,12,25,26,27,28] and thin film theory [27,28]. Unlike second-order diffusion equations, fourth-order equations generally do not satisfy any known maximum principle. Even with simple time-independent linear boundary conditions, evolving solutions tend to generate additional extrema from initially smooth conditions [29].

Broadbridge [30] studied the problem of a surface groove by evaporation-condensation governed by   $\frac{\partial y}{\partial t}=\frac{\frac{{\partial }^{2}y}{\partial x^2}}{1+{\left(\frac{\partial y}{\partial x}\right)}^2}$. The depth of a groove at a grain boundary was predicted without any approximation [30].

Chugunova and Taranets [31] studied the initial–boundary value problem associated with the fourth-order Mullins equation with initial data. They considered this problem by assuming that the specific free energy of the boundary is lower than the surface free energy. The Mullins equation, originally introduced by Mullins in 1957 [10], is a model used to analyze the evolution of surface grooves at the grain boundaries of heated polycrystals. Chugunova and Taranets [31] successfully demonstrated the global existence of weak solutions over time and established that the energy minimizing steady state serves as the global attractor.

Gurtin and Jabbour [32] developed a regularization theory that incorporates curvature effects, including surface diffusion and bulk-surface interactions. They investigated two specific cases: (i) the interface considered as a boundary between bulk phases or grains, and (ii) the interface between an elastic thin film bonded to a rigid substrate and a vapor phase depositing atoms on the surface [32].

Huang [33] conducted isothermal stress relaxation tests on electroplated Cu thin films, considering both passivated and unpassivated films. Based on a kinetic model, Huang [33] deduced grain-boundary and interface diffusivities and provided numerical and analytical solutions for the coupled diffusion problems. The study also analyzed the impact of surface and interface diffusivities on stress relaxation in polycrystalline thin films, comparing the results to experimental data.

Asai and Giga [34] considered the surface diffusion flow equation under specific boundary conditions. The problem of Mullins (1957) was proposed to model the formation of surface grooves on the grain boundaries, where the second boundary condition   $y^{\text{'}\text{'}\text{'}}\left(0\right)=0$ is replaced by zero slope condition on the curvature of the graph. Asai and Giga solved the initial-boundary problem with homogeneous initial data for construction of a self-similar solution and a solution was proposed by using a semi-divergence structure.

Escher et al. [35] demonstrated the existence and uniqueness of classical solutions for the motion of immersed hypersurfaces driven by surface diffusion. They focused the surface diffusion proposed by Mullins [10,11,12] to model surface dynamics for phase interfaces when the evolution is governed solely by mass diffusion within the interface. Other studies were devoted to the diffusion problems, grain boundary migration and grain dynamics evolution in materials [36,37,38,39,40,41,42].

Mullins et al. [43] have linearized the differential equation by assuming a very small slope at any point of the grain profile. In 1975, Brailsford & Gjostein [44] derived approximate solutions by studying the influence of surface energy anisotropy on morphological changes occurring by surface diffusion on simply shaped bodies. Wherever a grain boundary intersects the surface of a polycrystalline material, a groove develops. At the root of the groove, a balance between grain-boundary tension and surface tension produces an equilibrium angle [45]. The difference in chemical potential between the curved surface near the groove’s root and the smoother surface farther away results in material drift.

Tritscher [46] considered the boundary-value problem concerning the formation of a single groove due to surface diffusion at the junction of a bicrystal, assuming that the grain boundary remains planar.

Martin [47] extended the original Mullins theory of surface grooving due to a single interface to multiple interacting grooves formed by closely spaced flat interfaces. Martin considered two cases: the first involved simplifying Mullins’ analysis using Fourier cosine transforms instead of Laplace transforms, while the second dealt with an infinite periodic row of grooves. Martin [40] also solved the problem for two interacting grooves. Analytical solutions for the fourth partial differential equation governing the groove profile in metals have not been found in the literature.

In a previous study [48], we addressed the mathematical problem associated with the second non-linear partial differential equation in Mullin’s problem. We focused on the case of the evaporation-condensation and provided an exact solution for the geometric profile of grain boundary grooving when materials are subjected to thermal and mechanical stress, as well as fatigue effects.

This paper is devoted to model the grain groove profile governed by the fourth-order partial differential equation in the case of diffusion in thin polycrystalline films. An analytical and exact solution to the Mullins approximated problem,   $\frac{\partial y}{\partial t}+B\frac{{\partial }^{4}y}{\partial x^4}=0$, was given.

In this section, we were interested to the derivation of the differential equation that describes the evolution of a two-dimensional surface of small slope under capillary driving forces and surface diffusion transport. Surface properties are assumed to be independent of orientation. For a point on the surface at which the mean curvature is c, the chemical potential μ (c) per atom can be written as

where   ${\mu }_0$ is the chemical potential of reference for a flat surface (c = 0), γ is the surface tension of the metal/vapor interface and ω is the atomic volume of the film material. A gradient of surface curvature will therefore create a gradient of the chemical potential, which will produce a drift of atoms on the surface with an average velocity v given by the Nernst-Einstein relation.

where Ds is the surface diffusivity, k is the Boltzmann constant and T the absolute temperature.

The surface current of atoms J_S is defined by the product of the average velocity v by the atom number N_S per unit surface area S, it is given by the following equation:

The evolution of the surface may finally be described by the speed of movement v_n, of the surface element along its normal:

Notice that   $N_S$ is the number of diffusing atoms per unit area,   $J_S$ the surface current of atoms and B a rate constant given by the following equation:

Equation (7) can be written in the general case as:

If y is the coordinate of a point at the surface along the axis normal to the initial flat surface, the speed of motion of the point along this axis v_n is obtained by projection on the y-axis and one obtains:

Combining equations (9) and (10), one obtains:

Knowing that the curvature c is given by the following expression:

and

One obtains:

Using the same method for   $\frac{{\partial }^{2}c}{\partial s^2}$ , one obtains:

Therefore:

With   $y^{\text{'}}=\frac{\partial y}{\partial x}$ and   $y^{\text{'}\text{'}}=\frac{{\partial }^{2}y}{\partial x^2}$ , previous equation can be written as:

With the following boundary conditions:

Knowing that

One obtains:

By taking the following variable changes:

One obtains the different derivatives of   $y\left(x, t\right)$and   $u\left(x, t\right)$:

with

Equation (25) becomes:

or

By using the previous equations, one obtains:

With equation (20):

One writes:

$\frac{\partial y}{\partial t}=-B\frac{m}{{\left(Bt\right)}^{3/4 }}\frac{ \frac{{\partial }^{4}g}{\partial u^4}{\left(1+{m^2\left(\frac{\partial g}{\partial u}\right)}^2\right)}^2-m^2\left( {\left(\frac{{\partial }^{2}g}{\partial u^2}\right)}^3+10 \frac{\partial g}{\partial u} \frac{{\partial }^{2}g}{\partial u^2} \frac{{\partial }^{3}g}{\partial u^3}\right)\left(1+{m^2\left(\frac{\partial g}{\partial u}\right)}^2\right)+18m^4 {\left(\frac{\partial g}{\partial u}\right)}^2{\left(\frac{{\partial }^{2}g}{\partial u^2}\right)}^3}{{\left(1+{m^2\left(\frac{\partial g}{\partial u}\right)}^2\right)}^4 }$Let us put:

Using equation (32), one obtains:

If we suppose a second order approximation of the derivative,   $y^{\text{'}2} \ll 1$ , it is easy to deduce the following equation:

With the new boundary conditions:

In this section, we propose a new method to resolve the equation (38) by using the following equation:

and by considering the different solutions r in function of u.

Let us consider the following equation valid for all values of :

To resolve equation (38), we begin by transforming equation (39’) into difference between two perfect squares, therefore, the expression   ${(8\lambda r}^2+u r+{4\lambda }^2-1)$ will be transformed into perfect square, if it has a double solution and then his discriminant has to be cancelled.

Now, let us consider the equation:

The discriminant   $∆$ of this second-degree equation (40) function in r can be written as:

Putting   $∆=0$, one has:

Equation (42) can be written as:

With   $p=-\frac{1}{4}$ and   $q=-\frac{u^2}{128}$Putting   $\lambda =\alpha +\beta$ and taking   $\alpha \beta =-\frac{p}{3}=\frac{1}{12}$ or   ${\alpha }^3{\beta }^3=\frac{1}{{12}^3}$ one obtains   ${\alpha }^3+{\beta }^3=-q=\frac{u^2}{128}$ ; and   ${\alpha }^3$ et   ${\beta }^3$ will be the two solutions of the following second-degree equation:

or

The discriminant of equation (44′):

Can be calculated as a function of u:

Two cases have to be distinguished:

In this case, the solutions of equation (44′) will be given by:

This leads to the solution of equation (43):

This value of   ${\lambda }_2$ will cancel the discriminant of equation (40)

Therefore, the solution r is given by: and then:

Consequently, one obtains:

or

The four solutions of equation (37) can be then obtained from the solutions of the two following 2^nd degree equations:

The discriminants of equations (53) and (54) are given by the respective following expressions:

Two cases can be studied:

Solutions of   ${\mathit{r}}^2+\sqrt{2{\mathit{\lambda }}_2}\mathit{r}+{\mathit{\lambda }}_2+\frac{\mathit{u}}{8\sqrt{2{\mathit{\lambda }}_2}}=0$

Knowing that   ${∆}_1=-2{\lambda }_2-\frac{u}{2\sqrt{2{\lambda }_2}}$ is negative because of the condition   $u>\frac{2^{5/2}}{3^{3/4}}$, one obtains two conjugate complex solutions:

Solutions of   ${\mathit{r}}^2-\sqrt{2{\mathit{\lambda }}_2}\mathit{r}+{\mathit{\lambda }}_2-\frac{\mathit{u}}{8\sqrt{2\mathit{\lambda }}}=0$

Where   ${∆}_2=-2{\lambda }_2+\frac{u}{2\sqrt{2{\lambda }_2}}$

Let us prove that   ${∆}_2>0$

${∆}_2$ can be written as:   ${∆}_2={\lambda }_2\left(-2+\frac{u}{{\left(2{\lambda }_2\right)}^{3/2}}\right)$, To obtain the sign of   ${∆}_2$ , we have to study the sign of   $2\left(-1+\frac{u}{{2\left(2{\lambda }_2\right)}^{3/2}}\right)$ and then to compare between 1 and   $\frac{u}{{2\left(2\lambda \right)}^{3/2}}$ or between   $2{\lambda }_2$ and   $\frac{u^{2/3}}{2^{2/3}}$.

Let us put   $=\sqrt{1-\frac{2^{10}}{3^3{u}^4}}$ , one obtains:

Equation (61) shows that   $\frac{\partial Z}{\partial X}\le 0$ , this implies that Z decreases for all values of   $X\ge 0$ and   $Z<1$ for   $X>0$ and therefore   $\frac{2^{2/3}\left(2{\lambda }_2\right)}{u^{2/3}}<1$ or   $-2{\lambda }_2+\frac{u}{2\sqrt{2{\lambda }_2}}>0$ and   ${∆}_2>0$.

Therefore, the two other solutions are then given by equations (62) and (63):

Solution of equation (38) for   $\mathit{u}\ge \frac{2^{5/2}}{3^{3/4}}$

Now, the final solution, in the case of   $\ge \frac{2^{5/2}}{3^{3/4}}$ , is given by equation (64):

with

The solution in function of   $y\left(x, t\right)$ will be written as:

Using of the boundary conditions:

The first condition implies necessary: A₃ = A₄ = 0 and therefore the solution will be given by the following form:

This solution can be written as:

With

In this case, one obtains two conjugate complex solutions   ${\alpha }^3$ and   ${\beta }^3$:

Where   $A^2={\left|{\alpha }^3\right|}^2={\left|{\beta }^3\right|}^2=\frac{\frac{2^{10}}{3^3}}{2^{16}}=\frac{1}{2^6{*3}^3}$ and   $A=\frac{1}{2^3{*3}^{3/2}}$ and finally   $A^{1/3}=\frac{1}{2 \sqrt{3}}$With:

The real solution is given by: or

This leads to the solution of equation (43):

This value of λ will cancel the discriminant of equation (39)

The solution is given here by:

Remember that:

And therefore:

Their respective discriminants are given below:

Two cases can be studied for   $u<\frac{2^{5/2}}{3^{3/4}}$:

First case:   ${\mathit{r}}^2+\sqrt{2{\lambda }_1}\mathit{r}+{\lambda }_1+\frac{\mathit{u}}{8\sqrt{2{\lambda }_1}}=0$

Here one has   ${\lambda }_1=\frac{1}{\sqrt{3}}cos\left(\frac{\theta }{3}\right)$

Where   ${∆}_1=-2{\lambda }_1-\frac{u}{2\sqrt{2{\lambda }_1}}$ is negative. The two conjugate complex solutions of equation (53) are given below: