A New Solution to the Grain Boundary Grooving Problem in Polycrystalline Thin Films When Evaporation and Diffusion Meet in Power Electronic Devices

Power semiconductor devices and modules are critical components in various applications, including electric vehicles, renewable energy systems, and industrial automation. These devices often employ thin film technologies to achieve high performance and compact designs [1,2,3]. However, their reliability can be compromised over time due to fatigue phenomena, particularly in environments characterized by high temperatures and cyclic loading. Understanding fatigue mechanisms in thin films within the context of power semiconductor devices is essential for ensuring their long-term stability and efficiency [4,5,6].

In power semiconductor devices and modules, such as insulated gate bipolar transistors (IGBTs) and metal-oxide-semiconductor field-effect transistors (MOSFETs), thermal stress arises from the operation and dissipation of heat during switching events and continuous operation [7,8,9,10,11,12]. The cyclic nature of these temperature fluctuations subjects thin films within the devices to mechanical strain, leading to fatigue-related phenomena [13]. Groove formation, a characteristic manifestation of thermal stress, can occur on the surfaces of critical components such as substrates, interconnects, and passivation layers [14,15,16]. These grooves typically emerge along regions subjected to maximum thermal stress, such as the interfaces between different material layers or near localized heat sources. Over time, the repetitive stress cycles deepen and widen these grooves, potentially compromising the electrical and thermal performance of the semiconductor device. Additionally, the presence of grooves can increase the risk of electrical breakdown and thermal hotspots, further exacerbating reliability issues [17,18].

IGBT power modules are intricately structured with multiple material layers, subject to diverse sources of stress that can induce gradual degradation through intrinsic (chip-related) and extrinsic (package-related) failure mechanisms. Figure 1 provides illustrations of several of these mechanisms. However, the predominant causes of degradation primarily stem from thermomechanical factors, owing to temperature fluctuations experienced during the module’s operational lifespan. These fluctuations lead to progressive deterioration between successive layers, particularly pronounced in layers exhibiting significant mismatches in the coefficient of thermal expansion (CTE). Additionally, the degradation of metallization and the interface between metallization and bond wires represent prominent contributors to power module failures, further underscoring the multifaceted nature of degradation processes in these modules.

In the context of power semiconductor devices, groove formation due to thermal stress is intricately linked to evaporation and diffusion phenomena. Elevated temperatures within the device can accelerate the evaporation of volatile species from thin film materials or encapsulants, altering their composition and mechanical properties [19,20]. This loss of material may exacerbate stress concentrations and promote the initiation of grooves along vulnerable regions of the device [21,22,23,24]. Furthermore, diffusion processes play a crucial role in redistributing material within the thin films of power semiconductor devices. Diffusion-induced phenomena, such as interdiffusion at material interfaces or dopant migration within semiconductor layers, can influence the mechanical and electrical properties of the device [25,26]. The interplay between diffusion and thermal stress contributes to the evolution of grooves and the overall degradation of device performance over time [27,28,29].

Fatigue in thin films within power semiconductor devices and modules poses significant challenges to their reliability and performance under demanding operating conditions [30,31]. The formation of grooves due to thermal stress, coupled with evaporation and diffusion phenomena, represents complex interdependencies that must be addressed to enhance device durability [31,32,33,34,35].

The problem of grain boundary grooving in polycrystalline thin films was studied by many researchers [14,15,16,36,37,38,39,40,41,42,43,44,45,46,47,48,49]. The works of Mullins [14,15,16] were devoted to solving the problem governing the profile of grain boundary grooving. The mathematical formulation of this problem was developed by several scientists by focusing on the evaporation–condensation and surface diffusion by adopting several non-linear methods [37,38,39,40,41,42,43,44,45,46,47,48] due to the non-linear partial differential equation previously formulated by Mullins [14,15,16,50].

Li et al. [51] studied the interaction between grain boundary and cascade in some composites with reference to single-element Ni metal by applying molecular dynamics. They showed that the capture of point defects by grain boundary was effective, and the number of captured defects was related to the distance between primary knock-on atoms and grain boundary.

Chepak-Gizbrekht and Knyazeva [52] investigated a two-dimensional model of grain-boundary diffusion with consideration for triple junctions. They studied the effect of grain size on diffusion in a polycrystal under varying temperature conditions; found that the diffusant behaved non-monotonously, accumulating and consuming at the triple junction during diffusion; and showed a dependency of the diffusion coefficient on the grain size.

Moshtaghi et al. [53] studied the effect of ferrite grain boundary in some materials on hydrogen trapping, desorption, diffusion, and grain size. Their modeling results showed that in grain sizes smaller than a certain critical grain size, the hydrogen diffusion coefficient decreases with decreasing grain size, indicating that H trapping dominates short-circuit diffusion mechanisms along high-angle grain boundaries (HAGBs). They suggested that the fraction of acicular ferrite grains can define the HAGB content in the alloy and can be a determinant factor in the behavior of weldments in H-containing media.

Schweizer et al. [54] studied the volume and grain boundary diffusion processes of tungsten atoms in a metal matrix on the atomic scale. Using in situ high-resolution scanning transmission electron microscopy, they followed the random movement of single atoms within a lattice at elevated temperatures. They used the STEM to confirm random walk processes, quantify diffusion kinetics, and distinctly separate diffusion in the volume from diffusion along defects. They also confirmed that the mechanism of grain boundary diffusion can rely on both lattice sites and excess volume depending on the type of boundary.

Wang et al. [55] modified the microstructure of the grain boundary phase through the diffusion method, improved the coercivity of magnetic powders, and then their thermal stability. They studied the effect of GBD on the microstructures and magnetic properties of some magnetic powders and the modification of different content of the alloy to grain boundary and microstructure.

Zhao et al. [56] used a grain growth model based on a Monte Carlo algorithm integrated with the diffusion model by considering grain growth and diffusion in nanocrystalline materials. They studied the effects of grain growth and impurity on diffusion and provided a numerical strategy to study microstructure evolution in nanocrystalline materials.

Chatzimichail et al. [57] determined the temperature dependence of the surface energy and the grain boundary energy of polycrystalline MgO, as well as of the interfacial energy of MgO in contact with liquid Ag and Cu. They studied the mechanism and the kinetics of mass transport in the grooving process at the grain boundaries/free surface intersection of MgO. They found that surface diffusion as a mass transport mechanism is responsible for grain boundary grooving in polycrystalline MgO.

The kinetics of atom adsorption, desorption, and diffusion in polycrystalline materials were analyzed [58] to understand the influence of grain boundaries and grain size. The effect of the ratio of grain boundary and grain diffusion coefficients in concentration profiles was studied, as well as the influence of grain size on diffusion processes and the impact of kinetic processes taking place on the surface.

Gheno et al. [59] studied the grain boundary diffusion of chromium in polycrystalline nickel by SIMS and examined the shape of individual grain boundary profiles, and the applicability of the existing methods to process grain boundary diffusion profiles.

By advancing our understanding of these mechanisms and developing robust materials, designs, and manufacturing processes, it is possible to mitigate the effects of fatigue and improve the reliability of power semiconductor devices. Through interdisciplinary research and innovation, we can ensure the continued advancement and widespread adoption of power electronics technology in critical applications. Peyghami et al. [60] proposed a long-term performance indicator for power electronic converters based on their reliability. The converter reliability was represented by constant lifetime curves, which have been developed using an Artificial Neural Network (ANN) under different operating conditions. These proposed lifetime curves presented the long-term performance of converters facilitating optimal system-level design for reliability, reliable operation, and maintenance planning in power electronic systems.

We proposed, in previous works [61,62,63], analytical solutions to the mathematical problem in the case of the evaporation–condensation in polycrystalline thin films by resolving the corresponding second non-linear partial differential equation [61,62] without any approximation when materials are submitted to thermal and mechanical stress and fatigue effects. One proved the non-validity of Mullins approximation that neglected the first derivative in the mathematical equation associated with the evaporation case. In a recent study, we studied the problem of diffusion in thin polycrystalline films, resolved the fourth-order partial differential equation of diffusion, and provided an analytical solution to the equation $\frac{\partial y}{\partial t}+B\frac{{\partial }^{4}y}{\partial x^4}=0$, obtained by admitting the approximation ${y\prime }^2\ll 1$. However, we did not find any analytical study of the association of evaporation and diffusion phenomena in polycrystalline films.

In this paper, we studied the effect of the simultaneous combination of evaporation–condensation and diffusion on the grain groove profile by resolving the fourth-order partial differential equation associated with these two combined effects. An analytical solution was proposed and detailed in this present work and compared to the solution obtained separately in the diffusion and evaporation–condensation cases.

The differential equation describing the grain boundary grooving when evaporation and diffusion simultaneously act in polycrystalline thin films is given by Equation (1) or (2): or where $y^{\prime }=\frac{\partial y}{\partial x}$ and $y^{\prime \prime }=\frac{{\partial }^{2}y}{\partial x^2}$ and the boundary conditions are given by Equation (3): where $x$ and $y(x,t)$ represent the coordinates of a point at the surface along the axis normal to the initial flat surface at a time $t$. $C$ and $B$ are two constants relative, respectively, to evaporation and diffusion phenomena. $C$ and $B$ are given by Equations (4) and (5): where $\gamma$ is the isotropic surface energy or tension of the metal/vapor interface, $P_0\left(T\right)$is the vapor pressure at temperature $T$ in equilibrium with the plane surface of the metal, $\omega$ is the atomic volume, $M$ is the molar mass of the metal, $R$ is the perfect gas constant, $D_s$ is the surface diffusivity, and $N_S$ is the number of diffusing atoms per unit area.

Equation (2) can be written as follows:

By taking the following variable changes Equation (7):

The different derivatives of $y\left(x, t\right)$ and $u\left(x, t\right)$ are given by Equation (8):

One deduced the following general equation of the grain groove profile:

The above partial differential equation cannot analytically be resolved without approximation. In the following sections, we resumed the essential results obtained in previous studies [61,62,63].

The evaporation–condensation problem is governed by the following mathematical equation:

Equation (9) can be transformed with the same notations to the following equation: or

To resolve the non-linear differential Equation (11), Mullins did the approximation of a small slope by supposing $\left|y\prime \right|\ll 1$. He wrote the following:

The integration of Equation (12) gave the following:

With A being a constant of the problem determined by the condition boundary $y\prime \left(0, t\right)=m$

The solution of the differential Equation (12) obtained by the approximated Mullins problem is given by Equation (13):

In previous studies [61,62], we corrected the solution given by Mullins by considering general Equation (10) without any approximation and obtained the following equation:

And the final solution is given by Equation (15): where $\theta$ is the groove contact angle.

The diffusion case previously studied [63] was relative to the mathematical solution of the formation of grain boundary grooving in polycrystalline thin films by taking the approximation ${y\prime }^2\ll 1$. The following fourth differential equation (Equation (16)) was obtained:

Equation (16) satisfied the previous boundary conditions.

The analytical solution of Equation (16) was obtained and given by Equation (17): where the explicit expressions $g_1\left(u\right)$ and $g_2\left(u\right)$ of the solution were given as follows:

It was proved that the solution given by Hamieh et al. [63] revealed a damped sinusoidal groove profile in the case of electronic power devices. The expressions of zeros, minima, and maxima of the profile as a function of the order number, as well as detailed information about the groove profile y(x) and its derivatives, were given. The comparison of this new solution with Mullins’ results showed that the approximated solution obtained by Mullins [14,15,16] overestimated the values of the groove geometric characteristics, exceeding the actual values by more than 2.5 times. Additionally, valuable insights into the diffusion behavior of various metals were gained through this study. The new expressions relative to the diffusion were used to study several metals such as Cu, Al, Sr, Li, Cs, Ti, Co, Ga, and Tl and give the geometric parameters such as the depth $h_{Max}$ and the width $w_{Max}$ of the grain groove.

The evaporation/condensation and diffusion phenomena in polycrystalline thin films were studied separately in the literature. We did not find a complete and rigorous development of the simultaneous combination of evaporation/condensation and diffusion. In this section, we developed a new mathematical solution to this case by supposing that ${y\prime }^2\ll 1$.

In this case, Equation (6) can be written as

By using the reduced function $g\left(u\right)$, Equation (19) became as

The ratio $\frac{C}{B}$ is given by where M is the molar mass of the metal and R is the perfect gas constant.

Equation (21) can be resolved by using the following characteristic equation:

The thermo-diffusion coefficient $a$ is given by Equation (23):

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

Our method consisted of transforming the expression E into a perfect square by finding the double root of this second-degree equation. E is given by Equation (25):

This led to the study of the discriminant $∆$ of Equation (24) by distinguishing two cases relative to the positive and negative signs of $\left(u-\sqrt{a_1+a_2}\right)$ or $\left(u-u_0\right)$ where, $u_0$, $a_1$, and $a_2$ are given by Equation (26):

The analytical solution $g\left(u\right)$ was obtained by applying the boundary conditions to the studied cases, thus allowing the study of the profile variation of the grain groove (all mathematical details are given in Appendix A).

The final analytical solution of the grooving profile in the case of the combined cases of evaporation and diffusion is given as follows: where the functions $g_1\left(u\right)$ and $g_2\left(u\right)$ are given, respectively, by Equations (28) and (29):

The various expressions of the variables $\mu \left(u\right)$, $\theta \left(u\right)$, $\rho \left(u\right)$ and of the problem constants ($A_{11}, A_{21}, A_{12}, A_{22}, A_{32}, A_{42})$are provided in Appendix A.

The mathematical solution given by Equations (27)–(29) led to the drawing of the variations of the profile $y\left(x, t\right)$ in Figure 2 as a function of the distance x by taking the symmetric axis of the groove as the y-axis.

The evolution of $y\left(x, t\right)$ in Figure 2 showed similar variations with the curve obtained with the diffusion case alone. A damped sinusoidal profile of the groove with an infinity of maxima, minima, and zeros of the solutions was revealed (Figure 3) with smaller oscillations. Figure 3 clearly shows the first minimum and the second maximum of the profile with a rapidly decreasing amplitude with the distance $x$.

The first values of the different parameters of the grain groove shape, such as the coordinates of maxima ($x_{Max};$ $y_{Max}$) and minima ($x_{min};$ $y_{min}$) of the function $y(x,t)$, the differences between two consecutive maxima ($∆x_{Max}$) and two minima ($∆x_{min}$), and the zeros ($x_0$) of y, obtained from the solution $y\left(x, t\right)$ were presented in Table 1 for the first eight numbers by comparing the combined effect of evaporation and diffusion and the lone diffusion case.

Table 1 shows that the maximum $y_{Max}$ and the minimum $y_{min}$ of the solution $y(x,t)$ decreased towards zero when x increased to infinity. However, the difference between the two maxima $∆x_{Max }$ or two minima $∆x_{min}$ of the grain groove profile decreased and tended to zero at infinity. The different parameters of the grain groove profile presented some differences in their values when considering the combination of evaporation and diffusion phenomena compared to the lone diffusion case (Table 1). The curves relative to the comparison between the case of the diffusion alone and the combined evaporation and diffusion as a function of the distance were plotted in Figure 4. The maximum $y_{Max}$ of the grain profile in the combination of evaporation and diffusion effects showed a decrease in its value with respect to the case of diffusion. The same decrease was observed for all maxima $y_{Max}$ and minima $y_{min}$ of the function $y\left(x, t\right)$. However, the separation distance between two consecutive maxima $d_{Max}$ or minima $d_{min}$ increased relatively to the diffusion case.

To clarify the difference in the diffusion process and the combined effects, one gave, in Table 2, the ratios of the different parameters in the two examined cases. The case of the diffusion alone was symbolized by the letter $D$, whereas the evaporation/condensation combined with the diffusion was represented by $ECD$. The results in Table 2 showed that the values of $x_{Max }$, $x_{min}$, $∆x_{Max }$, and $x_{min}$ are practically the same in $ECD$ and $D$. However, an important decrease in the values of $y_{Max}$ and $y_{min}$ in a $ECD$ case was observed, and this decrease tended to zero when the distance $x$ increased.

The resulting variations of the different geometric characteristics of the grain profile in the diffusion alone and in the simultaneous case of evaporation/condensation and diffusion are certainly due to the thermal effect of the evaporation/condensation phenomenon. A decrease in the surface area of the groove profile was observed in the combined cases, proving the important competition between evaporation and diffusion, as it was shown in Section 8.

The effect of the slope $m$ at the origin or the contact angle on the shape of the groove profile was studied and is shown in Figure 5. It was shown that the groove depth increased when the contact angle of the groove $\theta$ increased while the other characteristics, such as $d_{Max}$ and the width $w_{Max}$ of the groove, remained the same.

The analytical rigorous solution was obtained for ${y\prime }^2\ll 1$. The variations of ${y^{\prime }}^2\left(x\right)$ versus the distance x were plotted in Figure 6. Two conclusions can be deduced from the curves of Figure 6. The first one was for the validity of the above approximation; it was concluded that the solution remains still valid even for all values of ${y\prime }^2$ with the condition of contact angle of the groove less than 27.6°. Knowing that the contact angle was defined as the angle formed between the y-axis and the tangent at the origin of the groove, it can be deduced that the total contact angle of the groove can reach 54.2° without any loss on the analytical solution. The second conclusion was that the analytical solution remains valid for distances $x\ge 6 \mathsf{\mu }\mathrm{m}$ even if ${y\prime }^2$ is not neglected behind 1.

Equation (20) can be also written as follows:

By considering the development of $g\left(u\right)$ in whole series as a function $u$, one writes the following:

The different derivatives of the above whole series are given by Equation (32):

Equation (30) can then be transformed into the following equation for all values of $u$ by replacing all derivatives with their developments: and one obtains a recurrent relation between the different terms $a_i$ of the series representing the function $g\left(u\right)$:

The complete determination of coefficients $a_i$ can be obtained by knowing the first four terms of the series $a_0$, $a_1$, $a_2$, and $a_3$. The other terms $a_i, \forall i\ge 4$can be deduced by recurrent relations. The boundary conditions of the grain grooving profile are used to determine the different terms of the series. Two types of boundary conditions can be applied: the first one was proposed by Mullins [14,15], and the second one by Amram et al. [64].

The use of the boundary conditions obtained by Mullins led to the following relations:

Using relation $a_2=\frac{a_0}{8}$ and Equation (33), one obtained the expression of $a_4$ as a function of the evaporation/diffusion coefficient $a$:

The calculations proved that the odd terms of the series are zero, and one writes:

The even terms are given by the following recurrent relation:

Therefore, the analytical solution $g\left(u\right)={\sum }_{i=0}^{\infty }{a}_i{u}^i$ was fully determined and obviously convergent towards zero.

By using the boundary conditions proposed by Amram et al. [54], one writes Equation (36):

This allowed the provision of the value of $a_4=-\frac{a_0}{96}$ and the recurrent formula giving all values $a_n$:

It was found that the boundary conditions proposed by Mullins [12,13] and Amram et al. [64] gave a very large difference between the two solutions. Figure 7 showed that the maximum of the groove with the boundary conditions of Amram et al. [64] was four times smaller than that obtained by using the boundary conditions of Mullins, whereas the grain groove deep was twice deeper than that of Mullins’s conditions. All geometric characteristics of the grain groove profile are affected by the change in the boundary conditions. It seems that the classical Mullins-type groove growing by surface diffusion alone cannot be satisfied by the presence of an evaporation/condensation process. The correction made by this present work represents a key microstructural element of thin metal polycrystalline films and, in many cases, participates in the determination of their physical, mechanical, and functional properties.

The general partial differential equation governing the grain boundary grooving submitted to evaporation/condensation and diffusion was given by the following:

The profile area of the groove ≙ was given by the following integral:

The integration of Equation (39) led to the following solution:

Using the boundary conditions (41): the integration of Equation (40), with respect to time $t$, led to Equation (42) giving the profile area of the grain boundary grooving as a function of time:

Equation (42) obtained from the general partial differential equation without any particular approximation is composed of two terms of the profile area relative to the evaporation/condensation and diffusion, respectively, given by the following relations:

Considering that $m=tan \theta$, where $\theta$ is the contact angle of the groove at origin, relations (42) can be written as follows:

To study the variations of the different profile areas of several metals heavily used in power electronic devices, one took the experimental data from previous works [62,63] and gave the values of $B$ and $C$ in Table 3.

Experimental results allowed us to draw, in Figure 8, the variations of the profile area of metals as a function of the parameter $m=tan \theta$, for a fixed time $t=12 \mathrm{h}$. The curves of Figure 1 show an increasing profile area of metals when the contact angle $\theta$ increases. The comparison between the results of the different metals led to the classification of these metals in order of increasing profile area as follows:

The above classification (45) is in good agreement with the evaporation and diffusion constants. The results showed that when these constants increase, the profile area increases, too.

The effect of the time on the profile of metals was shown in Figure 9 for a contact angle of the grain groove equal to $\theta =30°$ or for $m=0.58$. The results obviously confirmed the increase in the profile area when the time increased. This quantified the effect of time on the mechanical and thermal degradation of polycrystalline metals. The time effect gave the same classification order (given by relation (45)) of the profile area of the different metals.

The results relative to the two contributions of evaporation and diffusion in the values of the profile area of metals are represented in Figure 10 and Figure 11.

The results in Figure 10 showed similar variations of the profile area due to the evaporation process compared to those of Figure 8, and the same classification order of metals was obtained. However, the curves in Figure 11 relative to the diffusion part of the profile area presented non-monotonous variations. Indeed, the profile area due to the diffusion increased until a maximum of the parameter $m=m_{Max}=1$ corresponding to a contact angle $\theta =45°$ for all metals followed by a decrease for $\theta >45°$. The metals were classified in increasing order of the profile area as follows:

Relation (46) gave an order of metals completely different from that of relation (45) except for metals Sr and Cs, which always exhibited the highest profile area.

Now, in order to compare the two effects of evaporation and diffusion, we drew in Figure 12 the evolution of the profile area ratio of evaporation/diffusion for $t=12 \mathrm{h}$.

Figure 12a,b show that the ratio of evaporation on diffusion for all metals decreased until a minimum of the contact angle equal to $\theta =38°$. The results highlighted a higher effect of the diffusion for $\theta <38°$, while the evaporation effect predominated for $\theta >38°$. However, in all studied cases, it was shown that the ratio evaporation/diffusion is greater than 250 for all tested metals, proving that even if the diffusion exerted an important effect for $\theta <38°$, the profile area of the evaporation was at least 250 times higher than that corresponding to the diffusion.

A new analytical solution to the partial differential equation $\frac{\partial y}{\partial t}=C{y}^{\prime \prime }-B{y}^{\prime \prime \prime \prime }$ was proposed in the combined cases of evaporation/condensation and diffusion in polycrystalline thin films. A detailed mathematical method was given. The determination of the groove profile $y\left(x, t\right)$ and its derivatives allowed us to conclude that the validity of the new mathematical solution can be extended to the domain located out of that relative to the approximation ${y\prime }^2\ll 1$. Indeed, the obtained solution is valid even if the contact angle of the groove is greater than 30° for distance x greater than 6 μm. It was shown that the values of the parameters $x_{Max }$, $x_{min}$, $∆x_{Max }$, and $x_{min}$ of the groove profile were practically the same in the case of diffusion alone and combined evaporation/diffusion case. However, an important decrease in the values of $y_{Max}$ and $y_{min}$ in the combined case was observed, proving the capital role of the evaporation on the characteristics of the grain groove profile. The accurate determination of the profile area relative to the evaporation, diffusion, and combined case led to the interesting conclusion of the highest effect of the evaporation relative to that of the diffusion. This conclusion was confirmed by the experimental data of several metals, thus showing the important effect of evaporation. The effects of the time and contact angle on the profile area were highlighted in different metals. A maximum of the profile area relative to the diffusion was obtained for a contact angle $\theta =45°$, while a minimum in the ratio of the profile areas of evaporation on diffusion was observed for $\theta =38°$. Furthermore, another analytical solution of the partial differential equation using the series development method was proposed.