# Phase field modelling combined with data-driven approach to unravel the orientation influenced growth of interfacial Cu6Sn5 intermetallics under electric current stressing

Shuibao Liang, Cheng Wei, Anil Kunwar, Upadesh Subedi, Han Jiang, Haoran

Ma, Changbo Ke

► To cite this version:

Shuibao Liang, Cheng Wei, Anil Kunwar, Upadesh Subedi, Han Jiang, et al.. Phase field modelling combined with data-driven approach to unravel the orientation influenced growth of interfacial Cu6Sn5 intermetallics under electric current stressing. Surfaces and Interfaces, 2023, 37, pp.102728. 10.1016/j.surfin.2023.102728. hal-04085593

HAL Id: hal-04085593

https://hal.science/hal-04085593

Submitted on 29 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Phase field modelling combined with data-driven approach to unravel the orientation influenced growth of interfacial Cu₆Sn₅ intermetallics under electric current stressing

Shuibao Liang<sup>a,b</sup>, Cheng Wei<sup>c</sup>, Anil Kunwar<sup>d,*</sup>, Upadesh Subedi<sup>d</sup>, Han Jiang<sup>a</sup>, Haoran Ma<sup>e</sup>, Changbo Ke<sup>b,**</sup>

<sup>a</sup>Wolfson School of Mechanical, Electrical and Manufacturing Engineering, Loughborough University,
Loughborough LE11 3TU, UK

<sup>b</sup>School of Materials Science and Engineering, South China University of Technology,
Guangzhou 510640, China

<sup>c</sup>School of Resource Engineering, Longyan University, Longyan 364000, China

<sup>d</sup>Faculty of Mechanical Engineering, Silesian University of Technology, Konarskiego 18A, 44-100
Gliwice, Poland

<sup>e</sup> School of Microelectronics, Dalian University of Technology, Dalian 116024, China

*Corresponding author

**Corresponding author

E-mail addresses: anil.kunwar@polsl.pl (A. Kunwar); mecbke@scut.edu.cn (C. Ke)

**Abstract** Microstructure of interfacial intermetallics plays an important role in determining the service performance and reliability of interconnects especially due to the anisotropic properties of intermetallic grains, preferential growth of intermetallics induced by the electric current is observed in experiments, but the exact mechanisms for this have not been understood completely. For this endeavor, a phase field model considering the free energy arisen from the applied electric current is developed to tackle the intermetallic grain growth behavior in Sn/Cu interconnects, the focus falls on the influence of anisotropic electric conductivity. Simulation results show that electric current stressing preferentially accelerates the intermetallic growth, and the orientation effects are more pronounced with higher value of electric current density. Due to the competitive growth of multiple Cu₆Sn₅ grains in the presence of electric current, most region of the intermetallic layer is occupied by the grains whose c-axis is along the direction of the electric current. The intermetallic grain with higher electric conductivity along the electron flow holds smaller
current field, while the phase boundaries own higher electric field. It is found that the higher local electric
field concentrating at phase boundary drives faster phase boundary migration. The data generated from
physics-based phase field simulations are utilized to train a neural network model that functionally maps
the area of a particular grain with the applied current density, simulation time, grain identification,
orientation of the grain and its neighboring grains. The prediction model reveals that that electric current
stressing preferentially accelerates the intermetallic growth, and the orientation effects are more
pronounced with higher value of electric current density. The results from phase field simulation and
physically informed machine learning model further deepen understandings of the microstructure
evolution and selective intermetallic growth in the context of electric current, and shed light on the in-
silico studies and design route of interconnects under other types of loadings.

**Keywords:** Crystal orientation, Preferential grain growth, Intermetallic, Electromigration, Phase field model, Machine learning

**1. Introduction**

The formation and existing of intermetallics are one of the most important evidences that the
interconnects have a reliable performance, whether through soldering [1-3], welding [4], and self-
propagating exothermic reaction [5, 6]. Intermetallics occupying a suitable volume at the bonding
interface could ensure a good service reliability of the interconnects [7]. Noticing that the physical
properties of intermetallic compounds are usually anisotropic [8-11], such as mechanical properties,
thermal and electric conductivities. For example, Young's modulus of Cu₆Sn₅ intermetallic grains with
[010] direction parallel to loading direction is lower than other grain orientations [11, 12], and strength of
Cu₆Sn₅ grains closing to the c-axis has 20% increase compared with that normal to the c-axis [13].
Therefore, the microstructure characteristics of the intermetallic compound phase, such as the distribution and size difference of the grains with different orientations, can determine the performance and lifetime of the interconnects [1, 7, 14-17], particularly for the micro- and nano-scale integrated structures.

Interconnects in service are frequently subjected to multiple mechanical and physical loadings [18-21], to name a few, thermal stress, electric current stressing, and temperature gradient. Meanwhile, field-assisted manufacturing is being an attractive means to optimize or accelerate the bonding process to achieve the interconnects with better properties [22-24]. However, the mechanical and theoretical mechanisms behind the physical fields acting on microstructure and bonding process of interconnects remains to be further clarified. Generally, the physical loadings can promote or hinder the growth of intermetallic, depending on the type and direction of loadings [25, 26]. Moreover, the grain distribution and orientation are affected by the applied loadings [27, 28]. The studies of transient liquid phase bonding of Sn/Cu interconnects under temperature gradient showed that, in addition to accelerating the growth of intermetallic phase [29], the c-axis of the Cu₆Sn₅ intermetallics after bonding tend to be parallel to the temperature gradient direction [22]. The orientation of Co₃W intermetallic grains at the W/Co bonding interface after spark plasma sintering was observed to be associated with the density of electric current, and intermetallic grains preferentially grew in a specific direction along current direction [30]. The favor growth of Cu-Sn intermetallic grains with preferred orientations was found in Cu-Sn intermetallic interconnects [31, 32], and grain orientation of intermetallics were greatly influenced due to the applied high-density high electric current. However, the behind mechanism of preferential growth of intermetallic grains induced by the applied temperature gradient or electric current still needs to be clarified, thus further understanding and studies are imperative.

A whole range of researches show the capability and advantages of phase field modelling in the study of formation and growth of intermetallics in Sn/Cu [33-37], Mg/Al [38], and Al/Ni [39, 40] systems. For instance, Huh et al [33, 34] developed a phase field model to study the interfacial reaction in Sn/Cu system
during the conditions of soldering and aging, the simulated results are in line with previous observation.

Park et al [35] has studied the formation characteristics of intermetallic compounds at the bonding layer Sn-Cu system during the early stages of soldering with considering the phase nucleation, providing a better understand of the Cu-Sn soldering reactions. Moelans [36] proposed a new type of interpolation functions for describing the multi-phase systems, and applied to study the intermetallic growth. Ke et al [37] investigated the microstructure evolution in the region between Sn and unidirectional Cu substrates by phase field modelling incorporating micro-elasticity, the simulated results revealed the mechanism of formation of different patterns on different Cu substructures.

Moreover, the evolution of intermetallic compounds in interconnects under external loadings can also be simulated by coupling the phase field model with extra energy density or driving force terms, such as the energy terms exerted by electric current [41], stress field [42], and temperature gradient [43, 44]. Noteworthily, the phase field model was previously attempted by many researchers to reveal the microstructure evolution behaviour under external loadings [45-48], which were consistent with experimental observations. The previous study by Attari et al [47] taken into account the electromigration and non-equilibrium formation of vacancies, the simulation result gave some insightful understandings of the intermetallic growth behaviour and vacancy distribution influence on void formation. However, anisotropy in electric properties of intermetallic compounds would lead to the preferential growth of intermetallics, rare attentions were paid to this issue. For sub-100 nm structures of nanoscale interconnects transmitting electric currents, the operating current density can reach as high as several MA/cm² thereby causing reliability issues in these devices [48, 49]. Thus the knowledge of electromigration physics is of utmost importance for safety of nano-electronic devices. Our recent study [50] has mechanistically explored the electric current stressing induced preferential growth Cu₆Sn₅ intermetallic grains at the solder/substrate interface, but the study only used two order parameters for the intermetallic phase and
considered intermetallic grains with two different orientations, it can be of limited understandings of real microelectronic interconnects consisting of randomly oriented multiple intermetallic grains.

There has been an ongoing trend among the researchers to implement physically informed machine learning models at multiscale. The first approach is towards incorporating the laws of physics within the machine learning model itself, for an example solving the partial differential equations using neural network as described in the study by Raissi et al. [51]. The second approach is towards collecting the datasets with features that have variance in accordance with the governing equations of physics, and building the machine learning models upon these datasets. One way to generate physically intuitive datasets is by performing a large set of physics based multi-scale numerical simulations (molecular dynamics, density functional theory and finite element method) through variation of the parameters in the governing equations. The set of partial differential equations have been solved using finite element method in Kunwar et al. [44, 52] to generate data required by artificial neural network (ANN). The major advantage of this approach is the unlimited flexibility in the architecture of the machine learning model. The third approach is to develop multi-modal machine learning models capable enough of taking in the physics-based information in addition to the empirical datasets. This approach has been implemented recently in the development of physically informed machine learning potentials [53, 54].

In this work, the grain growth and electromigration behaviour of intermetallics in the Sn/Cu interconnects with emphasizing on effect of anisotropic electric conductivity are studied by developing a phase field model taking into account the electrostatic free energy. The cases of bi-crystal and multiple grains in the intermetallic phase are respectively investigated under electric current stressing, compared with the cases without electric current stressing. A set of phase-field simulations will be performed by varying the parameters such as current density, orientation angle etc., and data of these features along with time, grain identification and grain size etc. will be collected. The generated datasets will be utilized to train the ANN based on second approach mentioned above and then subsequently build a prediction model.
Thus, by integrating the physically intuitive computational method (phase field method) with machine learning technique, the present work will be able to unravel the hidden details of the mechanistic aspects associated with the electromigration-induced grain growth of interfacial intermetallics affected by orientation.

## 2. Numerical model

For the Cu/Sn bonding/soldering system, usually there are Cu₆Sn₅ and Cu₃Sn intermetallics at the bonding interface. Cu₆Sn₅ is merely considered in the present study since that it is the main phase in the reflowing process[35, 55], besides, the amount of Cu₃Sn intermetallic is nearly negligible compared with that of Cu₆Sn₅ intermetallic going through a considerable time of aging. The simulation domain is thus composed of three individual phases, i.e., the phases of substrate, Cu₆Sn₅ intermetallic and solder. The N number of orientations of grains could be represented by a set of order parameters {$\eta_i(x, y, z, t)$} ($i=1, 2, ..., N$). The variable $\eta_i$ and $\eta_N$ stand for substrate and solder, respectively. The other N-2 number of order parameters denote the intermetallic grains with different orientations. The grains labeled as same orientation are presented using one of order parameters with the value of 1, while other order parameters are 0. To further describe the different states for the intermetallic compared with the solder and substrate, another order parameter $\phi$ is introduced to represent the intermetallic phase, $\phi=1$ in the intermetallic and $\phi=0$ in the solder or substrate. The total free energy with the order parameters can be constructed as

$$F = \int_V \left[ \phi(\phi) + \xi(\eta_1, \eta_2, \dots, \eta_N, \phi) + \frac{k_\phi}{2} (\nabla\phi)^2 + \frac{k_\eta}{2} \sum_{i=1}^N (\nabla\eta_i)^2 \right] dV, \quad (1)$$

where $k_\phi$ and $k_\eta$ are gradient energy coefficients. The local free energy $\phi(\phi)$ is given by

$$\phi(\phi) = \left(\frac{a}{2}\phi^2 - \frac{b}{3}\phi^3 + \frac{c}{4}\phi^4\right)\Delta f, \quad (2)$$
where $\Delta f$ is determined by $\Delta f = (\eta_1^2\phi^2\Delta f_1 + \eta_N^2\phi^2\Delta f_2)/(\eta_1^2\phi^2 + \eta_N^2\phi^2)$, where $\Delta f_1$ is the chemical energy difference between the intermetallic and substrate phases, $\Delta f_2$ is the chemical energy difference between the intermetallic and solder phases. The constants are chosen as $a = 0.14, b = 12.42, c = 12.28$ [56]. The local free energy density $\xi(\eta_1, \eta_2, \dots, \eta_N, \phi)$ is given by

$$ \xi(\eta_1, \eta_2, \dots, \eta_N, \phi) = \xi_1(\eta_1, \eta_2, \dots, \eta_N, \phi) + \xi_2(\eta_1, \eta_2, \dots, \eta_N) + \xi_3(\eta_1, \eta_2, \dots, \eta_N, \phi), \quad (3) $$

where

$$ \xi_1(\eta_1, \eta_2, \dots, \eta_N, \phi) = m \left[ \phi^4/4 + (1-\phi)^4/4 - (\phi^2/2) \sum_{i=2}^{N-1} (2\eta_i^3 - \eta_i^4) + \sum_{i=1}^{N} \eta_i^4/4 - (1-\phi)^2 (2\eta_1^3 - \eta_1^4 + 2\eta_N^3 - \eta_N^4)/2 \right]; $$

$$ \xi_2(\eta_1, \eta_2, \dots, \eta_N) = m \left( 1.5 \sum_{i=2}^{N-1} \sum_{j=2, j \neq i}^{N-1} \eta_i^2 \eta_j^2 + 1.5 \eta_1^2 \eta_N^2 \right); $$

$$ \xi_3(\eta_1, \eta_2, \dots, \eta_N, \phi) = \left[ \phi^2(\eta_1^2 + \eta_N^2) + (1-\phi)^2 \sum_{i=2}^{N-1} \eta_i^2 \right] \Delta f / 2, $$

where $\xi_1(\eta_1, \eta_2, \dots, \eta_N, \phi)$ describes the coupling between the phase order parameter $\phi$ and the orientation order parameter $\eta_i$, $\xi_2(\eta_1, \eta_2, \dots, \eta_N)$ is used to account for the interaction between grains with different orientations, $\xi_3(\eta_1, \eta_2, \dots, \eta_N, \phi)$ denotes the coupling between the different phases. The parameter $m = 6\sigma_{gb} / l_{gb}$ [57, 58], in which $\sigma_{gb}$ and $l_{gb}$ stand for the grain boundary energy and width, respectively.

Based on the above free energy formulation, the “multi-well” can be formed between different grains in the intermetallic phase, and the thermodynamic stability between different phases is determined by $\Delta f$ (See Fig. 1).
<figure><img src="image_2.png" /><figcaption>Fig. 1. (a) Landscape of the free energy density $\xi(\eta_1, \eta_2, \dots, \eta_N, \phi)$, where red surface plot represents the free energy density changes from solder to intermetallic phase, and the blue surface represents the free energy density changes from substrate to intermetallic phase. Free energy density changes across the phase boundaries (b) and the intermetallic grain boundaries (c).</figcaption></figure>

When the electric current flowing through the system, the electrostatic free energy density difference across the grain boundaries will drive the atom to migrate. This effect can be accounted by coupling the electric potential gradient to Allen-Cahn equation [59], which is similar to the consideration of the temperature gradient influence on grain growth [60, 61]. The microstructure evolution of intermetallic grains (with specific orientation $\eta_i$) can be given by

$$ \frac{\partial \eta_i}{\partial t} = -L_{\eta} \left[ \frac{\partial \xi(\eta_1, \eta_2, \dots, \eta_N, \phi)}{\partial \eta_i} - k_{\eta} \nabla^2 \eta_i - \mu_{\eta} \nabla \eta_i \cdot \nabla \psi \right], \quad (4) $$

where $L_\eta$ is a constant relating to the grain boundary mobility, $\mu_\eta = N_a e Z_i \delta(\eta_i, \eta_j)/V_m |\nabla \eta_i|^2$ is the electric potential gradient coefficient, $N_a$ represents Avogadro's number, $e$ and $Z_i$ are electronic charge and effective charge number, respectively, $V_m$ denotes the molar volume, and $\psi$ is the electric potential, the
Dirac-δ function $\delta(\eta_i, \eta_j)$ is employed to apply the electromigration driving force at the grain boundaries in intermetallic phase, with the form of $\delta(\eta_i, \eta_j) = 1 - \left(1 - \tanh^2\left(\sum_{i=2}^{N-1} \sum_{j=2, j \neq i}^{N-1} \frac{\eta_i^2 \eta_j^2}{d}\right)\right)^2$ and d is a constant across the grain boundaries.

The governing equation of order parameter can be written as

$$ \frac{\partial \phi}{\partial t} = -L_{\phi} \left[ \frac{\partial \phi(\phi)}{\partial \phi} + \frac{\partial \xi(\eta_1, \eta_2, \dots, \eta_N, \phi)}{\partial \phi} - k_{\phi} \nabla^2 \phi - \mu_{\phi} \nabla \phi \cdot \nabla \psi \right], \quad (5) $$

where $L_φ$ is a kinetics mobility constant, the last term $\mu_φ = N_a e Z_φ w_{gb} M_{gb} / L_φ V_m$ is used to reflect the driving force induced by electric potential at phase boundaries, and $Z_φ = \left[ \sum_{i=2}^{N-1} \eta_i^2 Z_i - (\eta_1^2 Z_1 + \eta_N^2 Z_N) \right] / \sum_{i=1}^{N} \eta_i^2$, $w_{gb}$ is the actual physical width of the grain boundary, $M_{gb}$ is the grain boundary mobility.

In Eqs. (4) and (5), $\nabla\psi(\eta_1, \eta_2, \dots, \eta_N) = -\mathbf{J} \cdot \boldsymbol{\sigma}^{-1}(\eta_1, \eta_2, \dots, \eta_N)$, where $\mathbf{J}$ is the electric current density and $\boldsymbol{\sigma}(\eta_1, \eta_2, \dots, \eta_N)$ the electric conductivity tensor which is calculated through $\boldsymbol{\sigma}(\eta_1, \eta_2, \dots, \eta_N) = \sum_{i=1}^{N} \eta_i^2 \boldsymbol{\sigma}_i / \sum_{i=1}^{N} \eta_i^2$, the electric conductivity tensor of intermetallic grains with different orientations is given by $\boldsymbol{\sigma}_i^{\text{IMC}} = \mathbf{M} \boldsymbol{\sigma}_0^{\text{IMC}} \mathbf{M}^T$, with $\mathbf{M}$ the transformation matrix [62].

The model described above is applied to investigate the Cu₆Sn₅ intermetallic growth at the presence of electric current at 150°C. The multiphase system including Sn solder, Cu substrate and multiple Cu₆Sn₅ grains (N-2) are considered based on the experimental study [63] and simulation work launched by Ref. [42], and $\Delta x = \Delta z = 4.15 \times 10^{-9}$ m, as shown in Fig. 2. The schematic sketch of Fig. 2(a) is utilized to define the crystal orientation angle ($\theta$) of Cu₆Sn₅ intermetallic grains, and the definition of N-2 order parameters for intermetallic grains in Fig. 2(b) which allows the allocation of these different values of crystal orientation angles. The periodic conditions are adopted for $\eta_i$ and $\phi$ along the x axis, and Neumann boundaries are applied for the top and bottom edges. The periodic boundary condition allows the half
grains at the left and right boundaries to be defined by a single order parameter. Constant current density ($j_a$) values are ensured using these Neumann boundary conditions. The electric conductivities of substrate and solder are assumed to be isotropic, with the quantity of $\sigma^{Cu} = 5.85 \times 10^7$ $\Omega^{-1}m^{-1}$ and $\sigma^{Sn} = 9.09 \times 10^6$ $\Omega^{-1}m^{-1}$ [41, 64], respectively. The diagonal elements of electric conductivity tensor $\sigma_0^{IMC}$ of intermetallic grains are assigned as: $\sigma_{11}^{IMC} = \sigma_{22}^{IMC} = 2.38 \times 10^6$ $\Omega^{-1}m^{-1}$, $\sigma_{33}^{IMC} = 3.75 \times 10^6$ $\Omega^{-1}m^{-1}$ [65]. Other simulation parameters are: $L_{\eta} = 2.24 \times 10^{-7}$ m·s/kg, $L_{Sn/IMC} = 4.41 \times 10^{-12}$ m·s/kg, $L_{IMC/Cu} = 2.34 \times 10^{-13}$ m·s/kg, $Z_{Cu} = 5$ [66], $Z_{IMC} = 30$ [66], $Z_{Sn} = 18$ [67], $\sigma_{gb} = 0.5$ J/m$^2$ [42], $V_m = 1.629 \times 10^{-5}$ m$^3$/mol, $l_{gb} = 3 \times 10^{-8}$ m, $w_{gb} = 5 \times 10^{-10}$ m [68], $\Delta f_1 = 9.5 \times 10^8$ J/m$^3$ and $\Delta f_2 = 2.4 \times 10^8$ J/m$^3$ [69]. Then, the values of parameters $k_{\phi}$, $k_{\eta}$ and $m$ can be determined based on the assigned $\sigma_{gb}$ and $l_{gb}$ [58]. In the simulation procedure, the implicit backward difference method and the finite element scheme are used for the temporal and spatial discretization for solving the equations governing the grain morphology evolution and electric potential distribution.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Schematic configurations are illustrated in (a) to provide 3D perspective of angle $\theta$ between the z direction (i.e., electric current direction) and the c-axis of intermetallic grains, and (b) to show the view from y-axis direction of the simulation domain in x-z plane, where three phases of solder, intermetallic grains and substrate are considered. The sketch shown in (b) is used for the construction of the equivalent 2D geometry of the finite element model, the numbers 1, 2, ..., 8 tagged on the grains are the grain identifications.</figcaption></figure>
### 3. Phase field simulation of intermetallic growth under electric current stressing

We first consider the interconnect including the bicrystal intermetallic phase under the electric current stressing with different densities, and the growth rate and characteristics of intermetallic grains are compared with previous experimental studies to illustrate and verify the feasibility of the developed model. Subsequently, multiple grains with different orientations in the intermetallic phase are considered and to investigate the competitive and preferential growth of different grains though the grain size comparison. Finally, the typical simulation case in the last subsection is selected for analysing current density distribution, and discussion about the mechanism behind the preferential growth, mainly from the perspective of the electromigration induced driving force.

#### 3.1. Preferential grain growth of intermetallics

For simplicity, two grains are firstly considered in the intermetallic phase, and the orientation angle $\theta$ for grain 1 is set as 0°, the angle for grain 2 is 90°. Figure 3(a1-a3) show the simulated grain morphology at $t=1500$ s with different values of current density, 0, $2 \times 10^4$ and $2 \times 10^5$ A/cm$^2$. It is found intermetallic grains grows faster under the higher electric current density, the quantitative changes of the intermetallic phase layer can be clearly seen in Fig. 4(a). The increase velocity of intermetallic layer under current conditions of 0 and $2 \times 10^4$ is 0.024 and 0.034 nm/s, respectively, which is consistent with the experimental studies [70].

More importantly, from the simulation results, it is observed that the grain 1 grows faster than grain 2, and the difference is more pronounced under the electric current with higher density. As shown in Fig. 3(c), grain 1 not only occupies the spaces of solder and substrates, but also it can be seen that the width of grain 1 increases significantly especially adjacent to the solder region, while the width of the grain 2 decreases. The size comparison of grain 1 and 2 is shown in Fig. 4(b), and the grain size difference is more obvious when higher current density is applied. Noticeably, high current electric densities have been found
in the experimental studies on Cu/Sn and Co/W interconnects to induce the rapid growth or preferential
presence of intermetallic grains with certain orientations [30, 32]. The study by Deng et al [30] found that
the intermetallic grains with specific orientations are easier to grow and prevail after applying electric
current, The present simulation further confirms this behaviour, and it is in good agreement with the
findings by experiments.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Comparison of grain morphologies at *t*=1500 s for conditions of 0 A/cm² (a1), 2×10⁴ A/cm² (a2) and 2×10⁵ A/cm² (a3),
the dashed lines denote the initial position of grain boundaries.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** (a) Variation in intermetallic thickness with time under different current densities. (b) Comparison of the diameters of grain 1 and grain 2 with time under different current densities.</figcaption></figure>
## 3.2. Competitive growth of multiple intermetallic grains under electric current stressing

As intermetallic grains in the Sn/Cu interconnects usually have random crystal orientations [63], it is very important from the perspective of materials design engineering to consider simulation cases that account into the effect of multiple grain orientations. Two cases of systems with multiple orientation angles (8 number of $\theta$ in the range 0°-90°) for intermetallic grains of Fig. 2(b) are performed for studying the competitive growth under multiple intermetallic grains with different orientations. Case I is designed with an ascending order in the magnitude of orientation angles of the intermetallic grains from left to right of the intermetallic phase layer (0°, 20°, 30°, 40°, 50°, 60°, 70°, 90°), while case II assigns a more randomized values for the corresponding grains 50°, 60°, 20°, 90°, 0°, 70°, 30°, 40° for intermetallic phase from left to right. Fig. 5(a) shows the simulated morphology at $t$=500, 1000 and 1500 s without applied electric current stressing. Fig. 5(b) and Fig. 5(c) respectively present the simulated grain morphology at times $t$=500, 1000 and 1500 s of case I and case II under current density condition of $2 \times 10^4$ A/cm$^2$. It is apparent that the presence of electric current stressing induces the intermetallic phase grow faster, which was also found in some previous experimental studies [31, 32]. Moreover, the comparison between case I and case II shows that there is almost no morphological difference in terms of height and width of the intermetallic grains under current density condition of $2 \times 10^4$ A/cm$^2$, even though the grain orientations are evidently different.

When the applied electric current density increases to $2 \times 10^5$ A/cm$^2$, the simulated grain morphology at times $t$=500, 1000 and 1500 s of case I and case II is shown in Fig. 5(d-e). The intermetallic grains grow much faster, especially along the direction towards the solder matrix. At $t$= 500 s, the sizes of intermetallic grains with different orientations are no longer the same. Transient grain evolution is influenced by not only the magnitude of current density values but also the net flow of electrons along a grain of certain orientation. As time proceeds to 1000 s in Fig. 5(d-e), the size of the grain with orientation angle $\theta$=90° is very small, while the size of grain with $\theta$=0° increases sharply. At $t$= 1500 s, the grain with $\theta$=90° in case
I have been swallowed by the surrounding intermetallic grains, similar event also occurs in case II, and the grains with orientation angles (θ) of 90° and 70° have also disappeared. The above-mentioned evolution behavior results in the intermetallic grains with θ=0° occupying a relatively large area volume fraction of the Cu₆Sn₅ intermetallic phase region.

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** The simulated microstructure of solder interconnects under current density conditions of 0 A/cm² (a), 2×10⁴ A/cm² (b-c), 2×10⁵ A/cm² (d-e), at different times for case I (a, b, d) and case II (a, c, e). The dashed lines indicate the initial position of intermetallic grain boundaries.</figcaption></figure>
The analysis of intermetallic phase thickness over time is shown in Fig. 6. Clearly, the intermetallic phase thickness increases faster when a higher electric current stressing is applied. Moreover, the thickness of intermetallic phase in case I is almost the same with that in case II under the same electric current density.

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Average thickness of the intermetallic phase as a function of the intensity of applied current density for two different designed cases of grain orientations. The layer thickness of intermetallic phase represents the value obtained by dividing the total area of all the intermetallic grains by the base width of the interface.</figcaption></figure>

To quantitatively understand the growth and competitive behaviour of intermetallic grains with different orientations, the grain sizes with different orientations after electric current stressing at different times are plotted, as shown in Fig. 7. It is clearly seen that all the grains share the same trend, that the size of all grains increase with time regardless of grain orientations, this trend is particularly prominent in grains with $\theta=0°$. In case I, except that the sizes of grains with orientation angles ($\theta$) of 90° and 0° are very different from the grains with other orientations, and the size of grains with other orientations do not differ greatly. While in case II, it can be found that the grains with larger orientation angle ($\theta$) have smaller grain size. This result indicates that the intermetallic phase tends to favor more grains with orientation angle equal to or very close to 0°. The quantitative simulation results have revealed that the changes in the morphology of the intermetallic grains caused by the competitive grain growth are consistent with the experimental phenomenon that the Cu₆Sn₅ phase often have the grains with c-axis along the direction of
the electric current if a high-density electric current is applied [31, 32]. Moreover, the competitive growth becomes more obvious with the increase of electric current density, and the grains with c-axis along the electric current dominate the growth and exist in the intermetallic phase.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Size distribution pattern of intermetallic grains with different orientations for case I (a) and case II (b) at different times. The orientation angles in the x-axes of the histograms also serve as the grain identity in the intermetallic phase.</figcaption></figure>

**3.3. Current density distribution and driving force analysis**

Noting that the intermetallic grains with different orientations have different electric conductivity tensors consider the anisotropy of intermetallics, characterization of the correlation between the grain morphology and electric current density is critical for revealing the mechanism of the preferential growth of intermetallic grains. In general, considering that the difference of the electric conductivity tensor in different intermetallic grains, when the effective component of the electric conductivity of the grains along the electric current direction is larger, the electron flow tends to pass through these grains. The effective electric conductivity along the direction of the applied electric current could be given by

$$ \sigma_e^{\text{IMC}} = \sigma_{33}^{\text{IMC}} \cos^2 \theta + \sigma_{11}^{\text{IMC}} \sin^2 \theta. $$

Due to the conservation of electric current, the current density distribution closely depends on the intermetallic morphology, thus assessment of electric current density distribution in evolving interfacial intermetallic grains can be largely insightful in understanding the dynamics of a multi-phase system.
The quantitative distribution of current density is shown in Fig. 8, in which the grain morphology shown in Fig. 5(e) is considered. It is observed that most of the electrons flow through the grains with orientation angles of 0°, 20° and 30°, and other grains are observed with the smaller values of current density. Comparing Fig. 8(a) with 8(b) and 8(c), we see that the intermetallic grains with larger electric current density occupy the original spaces of the intermetallic grains with smaller one, which also means that those grains having higher electric conductivity preferentially grow and prevail in the intermetallic phase region.

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** The current density distributions in accordance with the grain morphology shown in Fig. 5(e) (case II) at $t=500$ (a), 1000 (b) and 1500 s (c). The red arrows indicate the directions of the electron flow.</figcaption></figure>

The driving force acting on the grain boundary due to the applied electric current can be expressed as $F_{em} = eZ^*Ew_{gb}$, where $E$ is the electric field [71, 72], $Z^*$ stands for the effective charge number, it can be seen that the driving force is proportional to the magnitude of electric field. Fig. 9(a) shows the electric field distribution of the simulation domain with morphology corresponding to $t=1000$ s of case II shown in Fig. 5(e). It is observed that the electric field in the grains with small orientation angles (e.g., grain with 0°) is lower than that in grains with large orientation angles (e.g., grains with 70° and 90°) due to the higher electric resistivity of large-orientated grains along the electron flow direction. However, the electric field in the micro-domains of solder near the grains with small orientation angles is higher than that near grains with large orientation angles.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** (a) The electric field distribution corresponding to *t*=1000 s of Fig. 5(e) (case II) is presented. (b) Plot of electric field and its vector direction at the phase boundaries and grain boundaries, color of the vector flows at the phase boundaries of solder/intermetallic and intermetallic/Cu depends on the magnitude of the electric field.</figcaption></figure>

Additionally, from the plot of electric field and its vector direction at the phase boundaries and grain boundaries shown in Fig. 9(b), it is found that electric field acting on the phase boundaries of grain A/solder and grain A/substrate is higher than that of grain B(C)/solder and grain B(C)/substrate, respectively. This could be the main reason for having different growth velocities of grain A and B(C). When the driving force for the boundary changes, the boundary between A and B(C) is no more the same with the electric current direction. Electric field force also exists at the boundary region between neighboring grains in the intermetallic phase, as shown in the gray region in Fig. 9(b), which was reported in previous study [59] and observed in the experiment [27]. This extra driving force at the grain boundaries in intermetallic layer will induce the grain boundaries to migrate along the direction parallel to the electric field, further leading to a more rapid growth of grain A by shrinkage of grain B(C). Base on the experimental studies [32, 65], it was indicated that the preferred orientation of intermetallic grains [0001] direction paralleling to the electric current direction is closely related to the electric resistance of the intermetallic phase, which is also shown in the present simulation. Moreover, from the perspective of the relationship between electric conductivity and carrier mobility, σ = neμ, where n is carrier concentration, μ is the carrier mobility and σ is the electric conductivity, it is seen that lower electric resistance induces higher carrier mobility since n is relatively constant [65]. When the c-axis of the grain is along the direction
of electric current, i.e., grain A with $\theta=0°$, the scattering obeyed by the electrons is the lowest due to the smallest atomic density through the Cu₆Sn₅ (0001) plane [32], then may lead to the higher electromigration and then faster growth of grain A.

To further explain our findings in the above discussion, we can regard the computational domain as the one consisting of two rectangular domains with the same size, domain 1 has three layers including solder, substrate and whole grain A, and the three-layered domain 2 owns the whole grain B. At the initial time, because the volume of the intermetallic grains is smaller compared with that of the solder and the substrate, and the resistivity of the solder/substrate is smaller than the intermetallic phase, the voltage across the two rectangular domains is almost the same. However, the voltage across grain B occupies a larger ratio of the total voltage than grain A due to the higher resistance of grain B, which leaves the voltage acting on the phase boundaries of grain B smaller than that on grain A. Hence, the electric field at phase boundaries of grain A will be greater, which in turn leads to the phase boundary adjacent to grain A growing faster than that of grain B.

4. Physically informed neural network for learning the differential grain growth

As discussed in previous sections, multiple number of phase-field simulations have been performed by varying the magnitude of applied current density and the orientation angles of the grain system. Thus, these factors are very significant in relation to the intermetallic grain growth kinetics. The grains will have different areas at different simulation time, and thus the time is also an important feature during grain growth. If a grain of grain identification $p$ ($p$ = 1, 2, ..., 8) is chosen for the intermetallic grains from left to right of the intermetallic phase layer shown in Fig. 2(b), it is already found that its orientation angle ($\theta^p$) will affect the evolution rate of its area ($A^p$). The adjacent grain also will affect the relative growth of the grain; and in this respect the orientation angle of the grain at vicinity also needs to be considered as another significant feature affecting the growth kinetics of the grain $p$. In the present study, as shown in
Fig. 2(b), the grains are linearly arranged at the interface. Thus, for a grain *p*, the grain at its left side will be referred to as grain (*p*-1) and the grain at its right side will be known as grain (*p*+1). The orientation angles $\theta^{p-1}$ and $\theta^{p+1}$ of these two neighboring grains subsequently matter to the evolution rate of grain *p*. In summary, seven features, namely, grain identification (*p*), $\theta^p$, $\theta^{p-1}$, $\theta^{p+1}$, applied current density (*j*), *t* and $A^p$ are chosen as the variables for the data-driven model of the present study. From the phase field simulations, 144 observations were made regarding the values of these 7 features. Since the data were generated from the solutions of partial differential equations, the neural network model will be adequately informed with the underlying physics.

It is always important to assess these features based upon the principle of maximum relevance and minimum redundancy. Pearson's correlation coefficient (PCC) heat map representation of the dataset is a good way to get insight about the linear relationship between any two features in the list of these 7 features - *p*, $\theta^p$, $\theta^{p-1}$, $\theta^{p+1}$, *j*, *t* and $A^p$. As it is aimed to define $A^p$ as the output feature and the features *p*, *t*, $\theta^p$, $\theta^{p-1}$, $\theta^{p+1}$, and *j* as the input features in the ANN in Section 4.1, the minimum redundancy and maximum relevance criteria will be established based upon this input-output structure. In the context of this study, an input feature or attribute is considered to be relevant and not redundant if it has a correlation coefficient of absolute magnitude less than 0.5 with other input feature whereas it has the coefficient of magnitude more than 0.05 with the output feature. A PCC heatmap for the dataset is shown in Fig. 10(a). It should be noted that the coefficient with the absolute magnitude less than $1\times10^{-9}$ is established as 0.0 in the heatmap. It should be noted that the feature related to *time* (*t*) has correlation coefficient of 0.0 with every other remaining 5 input features and has a correlation coefficient of +0.25 with $A^p$ (output feature). Thus *t* is a very relevant feature and does not show redundancy. Compared to *t*, *j* is even more relevant as its correlation coefficient (+0.52) with $A^p$ is the highest, and also it has no observable linear relationship with the other five features. Thus, the current density is considered as the most relevant input feature of this
study. After *j*, the orientation angle of *p*th grain (θP) is the second most relevant feature as it has the
correlationship coefficient of -0.26 with the output feature. None of the features are seriously redundant
and irrelevant as all of the input features pairs considered will have a correlationship coefficient well below
an absolute magnitude of 0.5. The input feature pair *p* and θP showing a correlationship coefficient of
+0.42 show a little redundancy when considered together, but since the absolute magnitude 0.42 is smaller
than the threshold of 0.5, both of these features will be considered relevant for this study. Hence, all of
the 6 input features will be considered relevant and non-redundant in the present study.

Data visualization can provide a good summary regarding the distribution of the variables. As shown
in Fig. 10(b), a pairplot diagram was constructed upon the dataset for the 5 features - grain identification
(*p*), θP, *j*, *t* and *A*P. The diagram clearly reveals that the grain growth represented by *A*P is affected
importantly by the choice of θP, time or current density, and thus provides a green signal to proceed with
the development of artificial neural network.
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Pearson's correlation coefficient (PCC) heatmap is constructed in (a) for the correlation analysis among the 7 features - *p*, *θ*^p, *θ*^p-1, *θ*^p+1, *j*, *t* and *A*^p. It is evident from the image that the grain area (*A*^p) has the strongest correlation with the current density *j* (correlation coefficient = +0.52) whereas the mildest correlation with *p* (correlation coefficient = -0.1). The correlation between the four features: grain identification (*p*), orientation angle for the *p*^th grain, applied current density (*j*) and area of the *p*^th grain (*A*^p) have been visualized using pairplot diagram in (b). It is important to note in the figure that some grains with non-zero initial area have a zero area at *t* = 1500 s, these grains correspond to larger orientation angle and larger current density. For larger *j*, the grains with orientation angle of 0 have exceptionally large areas at *t* =1500 s.</figcaption></figure>
## 4.1 Architecture of the artificial neural network

Among the seven features, the grain identification *p* of intermetallic grain, θ<sup>p</sup>, θ<sup>p-1</sup>, θ<sup>p+1</sup>, *j* and *t* were taken as the input features whereas the area (A<sup>P</sup>) of the *p*<sup>th</sup> grain was considered as the output feature for the ANN shown in Fig. 11(a). Corresponding to the 6 input features at the input layer, the size or dimension of the input was 6. Three hidden layers, namely *hl1*, *hl2* and *hl3* were defined in between the input and output layer. The numbers of neurons in *hl1*, *hl2* and *hl3* were respectively assigned as 20, 10 and 5. While ReLU activation function was defined for the first hidden layer, LeakyReLU activation function was chosen for the second and third hidden layers. The output layer was assigned a single neuron, and to account for the non-linearity associated with the data consisting of both growing and shrinking intermetallic grains, sigmoid activation function was defined for the output layer.

The original data supplied by the phase field simulations is marked by the differences in ranges of the datasets of the features. The grain identification varies in the range 1-8, whereas the orientation angle features (θ<sup>P</sup>, θ<sup>P-1</sup>, θ<sup>P+1</sup>) are in the values lying from 0 to 90°. The current density features range from 0 to 2×10<sup>5</sup> A/cm<sup>2</sup> and time is in the range 500-1500 s. The use of datasets with features of different ranges of variation in the machine learning model will cause the result to be dependent mainly on the input feature having the largest range - current density in the present context. To avoid the erroneous prediction caused by the hegemony or dominance of a single input feature, a technique called normalization [52, 73] was implemented on the dataset in accordance to which all of the input and output features were scaled within the range 0-1. With this procedure, all the input features can have their proportional claim on the prediction result.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 11.** A neural network model with an input layer (i), three hidden layers - hl1, hl2, hl3 and an output layer (o) is shown schematically in (a). The grain identification (p) of an intermetallic grain, orientation angles of the grain as well as its adjacent neighbors, current density and time are chosen as the input features whereas the area of the grain (A^p) is selected as the output feature for the ANN. The image in (b) shows the mean square error (MSE) values at 199th epoch for 10 different ANN models compiled with different values of constant learning rates. In image (c), the variation of MSE values for training and validation data is plotted against the number of Epochs for ANN model corresponding to learning rate = 2.5E-04.</figcaption></figure>

To enable the ANN to predict well beyond the training data, 25% of the total observations were designated as validation data. Mean square error (MSE) is chosen as the metrics for performance assessment of the ANN, and is defined mathematically as following:

$$MSE = \frac{1}{2n} \sum_{n=1}^{N} (y_e - A^{p,c})^2, \qquad (6)$$

where A^p,c is the area of the p^th grain computed from phase field simulations and y_e is the corresponding value estimated by the neural network. With Adam optimizer being selected as the optimizer for model, the neural network model was run and compiled in TensorFlow software [74]. Learning rate (lr) is the
hyperparameter that has a significant role on the performance of the neural network. 10 ANN models were built in this study, and each model was assigned a unique constant value of learning rate. The learning rates used in this work range from 1.0E-05 to 1.0E-02. The training and validation MSE values at the 199th epoch for the 10 models are plotted graphically in Fig. 11(b). As revealed by the black colored dots in the figure, the models with $lr = 2.5E-04$, $lr = 5.0E-04$ and $lr = 1.0E-02$ have the low values of MSE for the training data. Among the three models, the model with $lr=5.0E-04$ has the lowest MSE value for training data (MSE train =4.56E-05), followed by the model with $lr = 2.5E-04$ (MSE train =4.97E-05). The model corresponding to $lr$ of 1.0E-02 has MSE train =5.42E-05. In terms of MSE for validation data shown by red colored dots in the image, the model with $lr = 2.5E-04$ has the lowest error (MSE validation = 3.3E-03). The validation mean square error for models with $lr=5.0E-04$ and $lr=1.0E-02$ are 3.8E-03 and 8.2E-03 respectively. Owing to the reason that it has the lowest validation MSE, the model corresponding to $lr = 2.5E-04$ is chosen subsequently as the prediction model. The training and validation MSE curves for this model ($lr=2.5E-04$) for the entire 199 epochs are presented in Fig. 11(c).

## 4.2 Prediction of grain area as a function of its orientation

As the ANN was trained with datasets generated from phase field simulations, the machine learning model is physically informed and so it is able not only to interpolate the data but also to extrapolate it without making mistakes. Physically informed neural network models are a suitable choice to make predictions as they are computationally less expensive than the phase field simulations. That the ANN prediction model as applied to make estimations beyond the training data, is essentially acting as the sole model to make predictions for all the new values of input feature sets.

The main purpose of this study is to assess the influence of area of Cu₆Sn₅ intermetallic grain *p* by a single variable, namely its crystal orientation angle ($\theta^P$). The prediction ANN model can be utilized for this purpose by predicting $A^P$ as a function of $\theta^P$ solely and keeping other 5 input features as constant. The
curve of the grain area plotted as a function of θP by defining constant values of other input variables is known as the Ceteris Paribus plot, and several such plots are presented in Fig. 12. The time feature is selected as t=1250 s. The fourth intermetallic grain from left is chosen, i.e. grain identification p= 4. Thus, the neighboring intermetallic grains (p-1) and (p+1) are 3 and 5 respectively. The red colored curves correspond to a constant current density of 5×10⁵ A/cm² whereas the black colored lines represent j of magnitude of 5×10⁴ A/cm². The image in Fig. 12(a) correspond to the orientation angle of both the neighboring grains kept at constant value of 0° whereas Fig. 12(a) reveals that these two grains 3 and 5 have orientation angle of 90°. The reference area Aᵢ =1.70×10⁻¹⁴ m² shown in the images is the original area of grain 4 (as defined in the initial condition of the phase field simulation).

<figure><img src="image_16.png" /><figcaption>**Fig. 12.** The estimation results from the prediction model are presented in the images. Inspired from the concept of Ceteris Paribus plots, the prediction model tries to study the variation in area of 4th Cu₆Sn₅ intermetallic grain (grain identification p with p=4) with the change in its orientation angle (θP; p=4) when other input features are held constant. The orientation angles of neighboring grain p-1 on the left grain p is considered equal to that of another neighboring grain p+1 on the right of grain p. That is, θp-1 = θp+1 = 0° in (a) whereas θp-1 = θp+1 =90° in (b). The value of time (t) in both images is 1250 s. The red colored predicted curves in the two images correspond to the constant current density of magnitude 5×10⁵ A/cm² whereas the black colored curve is predicted for j value kept constant at 5×10⁴ A/cm². The dotted green line corresponds to the numerical value of area (Aᵢ) of the grain p at t =0 (initial condition of the phase field simulation).</figcaption></figure>

As shown in Fig. 12(a) and 12(b), at t= 1250 s, the area of grain 4 is predicted to become smaller and smaller as its orientation angle (θP) increases from 0 to 90°. It is important to note the effect of the orientation angles of the neighboring grains as well. As revealed in the Fig. 12(a), when the neighboring grains have orientation angle of 0°, these grains will not only grow faster than their adjacent intermetallic
grains with larger orientation angle but also grow at the expenses of the later. So, when the grain 3 and grain 5 have orientation angle of 0°, the grain 4 for $\theta^p = 90°$ is predicted to have an area of $1.7055 \times 10^{-14}$ m² at $t=1250$ s and $j=5 \times 10^4$ A/cm². It could be understood that though the grain 4 has first grown due to combined electromigration and diffusion, but has again reduced to an area (nearly equal to its initial area) due to the presence of larger neighboring grains. The reduction in area is more pronounced at larger current density of $5 \times 10^5$ A/cm². In Fig. 12(a), the grain 4 with $\theta^p = 90°$ is estimated to be reduced to an area of $2.43 \times 10^{-15}$ m² at $t=1250$ s, and this is very small compared to the initial area $A_i$.

The prediction results in Fig. 12(b) illustrate another scenario regarding the roles of the orientation angles of the surrounding intermetallic grains. When 3rd and 5th grains have orientation angle of 90°, the grain 4 is predicted to have an area always larger than $A_i$ at all values of orientation angles in the range 0-90° at $t=1250$ s for both values of current densities. The grain 4 grows at the expense of its neighboring grains, and this phenomenon is more pronounced at larger current density and smaller $\theta^p$. With the increase in $\theta^p$, the grain's area gradually decreases, and this is more severe at higher magnitude of current density. At $\theta^p=0°$, the areas of the grain 4 corresponding to $t=1250$ s for $j=5 \times 10^4$ A/cm² and $j=5 \times 10^5$ A/cm² are respectively $2.016 \times 10^{-14}$ m² and $3.899 \times 10^{-14}$ m². When the orientation angle of this grain changes to 90°, the area values predicted at $j=5 \times 10^4$ A/cm² and $j=5 \times 10^5$ A/cm² are $1.90 \times 10^{-14}$ m² and $2.70 \times 10^{-14}$ m² respectively.

The physically informed neural network has clearly revealed that the orientation of a grain, and its neighboring grains has a significant effect on the growth kinetics of the Cu₆Sn₅ intermetallic grains at the interface of Cu-Sn interconnects subjected to electric current stressing. Moreover, the computational efficiency of these machine learning models are expected to enable the accelerated design of material interfaces characterized with multi-variate phenomena and mechanisms.
**5. Limitations of the simulation**

The above findings reveals that the morphology of intermetallic layer is significantly related to the anisotropy nature of electric conductivity, since the driving forces of electromigration are different for grains with different orientations. This consolidates the conclusions proposed based on experimental study that the low-resistance route for electron flow may be formed by the preferential epitaxial growth of intermetallics [65]. The developed phase field model is beneficial for digging the mechanism of materials’ selective growth behaviour in the environment of electric current, as well as helpful for providing some inspiration to understand the microstructure evolution and preferential grain growth under various types of loadings, such as temperature gradient [22] and magnetic field [75]. However, it should be mentioned that the present model do not take into account the concentration field since that the concentration gradients in the same phase region are small based on the experimental studies [76, 77], the ruling out of concentration field would better inspect the effect arising from current density at the cost of negligible accuracy. Despite it was reported that the electric current drives the dominant growth of intermetallics and it is mainly related to the electric properties of intermetallic grains [30, 31, 78], which is confirmed and further clarified in this study, the anisotropies of the interfacial energy [37] and diffusion coefficient [66, 79] may affect the growth behavior of the grains. It is anticipated that the above-mentioned limitations of the present model will be tackled in the forthcoming study together with considering effects of other factors.

**6. Conclusions**

In this study, a phase field model taking into account the electric field is developed to study the growth behaviour of multiple Cu<sub>6</sub>Sn<sub>5</sub> grains in the intermetallic layer under electric current stressing, and the generated datasets is utilized to train the ANN and subsequently build a prediction model, with focuses
on clarifying how the intermetallic grain grows when considering the anisotropy of electric conductivity.

The following conclusions can be drawn from this work:

(1) Simulation results show that the anisotropy of electric conductivity induces different growth rates for grains with different orientations. The electromigration investigation cases of multiple intermetallic grains further confirm that grain grows faster and more preferentially prevails in the intermetallic phase when the direction of the applied electric current is along the *c*-axis of the intermetallic grain, which is qualitatively consistent with the previous experimental study.

(2) The initial distribution characteristics of intermetallic grain orientations are found to seriously affect the competitive and preferential growth of grains. Compared with the case having an ascending order in the magnitude of orientation angles of the grains, the grain growth is more obvious for the case with relatively random orientated intermetallic grains.

(3) The grains with higher electric conductivity along the current direction grow faster, and those grains dominate the thickening and coarsening process in intermetallic phase. It is revealed by the simulation that the faster growth of grains mainly depends on the higher local electric field near the phase interfaces.

(4) The physically informed neural network model was utilized to predict the area of a grain as a function of its orientation angle, while maintaining the other features as constant. This grain area function was then determined for two different constant values of other three features orientation angle of a neighboring grain, orientation angle of another neighboring grain, and applied current density. It is revealed that the increase of orientation angle is responsible for the decrease in its area during current stressing. When both neighboring grains have orientation angle of 0°, these neighbors grow at the expense of the area of the intermetallic grain. On the other hand, when both the neighboring intermetallic grains are oriented at 90°, the grain expands in areas at their
expense. The area expansion and reduction of the grain is more pronounced at larger current
density.

**Declaration of Competing Interests**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Acknowledgments**

This work is supported by National Natural Science Foundation of China (No. 51205135 and No. 52101035) and the National Science Centre, Poland (SONATA BIS grant, Grant Number: 2021/42/E/ST5/00339).

**Data Availability**

Data will be made available on request.

**References**

[1] Annuar S, Mahmoodian R, Hamdi M, Tu K-N (2017) Intermetallic compounds in 3D integrated circuits technology: a brief review. Sci Technol Adv Mater 18: 693-703

[2] Sun L, Chen M-h, Zhang L, He P, Xie L-s (2020) Recent progress in SLID bonding in novel 3D-IC technologies. J Alloys Compd 818: 152825

[3] Kannojia HK, Dixit P (2021) A review of intermetallic compound growth and void formation in electrodeposited Cu–Sn Layers for microsystems packaging. J Mater Sci: Mater Electron 1-36

[4] Liu B, Tian Y, Wang C, An R, Liu Y (2016) Extremely fast formation of CuSn intermetallic compounds in Cu/Sn/Cu system via a micro-resistance spot welding process. J Alloys Compd 687: 667-673
[5] Liang S, Zhong Y, Robertson S, Liu A, Zhou Z, Liu C (2021) Thermo-mechanical Characteristics and Reliability of Die-attach through Self-propagating Exothermic Reaction Bonding. IEEE Trans Compon Packaging Manuf Technol

[6] Dombroski DM, Wang A, Wen JZ, Alfano M (2022) Joining and welding with a nanothermite and exothermic bonding using reactive multi-nanolayers–A review. Journal of Manufacturing Processes 75: 280-300

[7] Tu K-N, Hsiao H-Y, Chen C (2013) Transition from flip chip solder joint to 3D IC microbump: Its effect on microstructure anisotropy. Microelectron Reliab 53: 2-6

[8] Choudhury SF, Ladani L (2015) Single Crystal Plasticity Finite Element Analysis of Cu 6 Sn 5 Intermetallic. Metallurgical and Materials Transactions A 46: 1108-1118

[9] Zhang Z, Li M, Liu Z, Yang S (2016) Growth characteristics and formation mechanisms of Cu6Sn5 phase at the liquid-Sn0.7Cu/(111) Cu and liquid-Sn0.7Cu/(001) Cu joint interfaces. Acta Mater 104: 1-8

[10] Daeumer M, Sandoval ED, Azizi A et al (2022) Orientation-Dependent Transport Properties of Cu3Sn. Acta Mater 117671

[11] Choudhury SF, Ladani L (2014) Grain growth orientation and anisotropy in Cu 6 Sn 5 intermetallic: nanoindentation and electron backscatter diffraction analysis. J Electron Mater 43: 996-1004

[12] Mu D, Huang H, Nogita K (2012) Anisotropic mechanical properties of Cu6Sn5 and (Cu, Ni) 6Sn5. Mater Lett 86: 46-49

[13] Jiang L, Jiang H, Chawla N (2012) The effect of crystallographic orientation on the mechanical behavior of Cu 6 Sn 5 by micropillar compression testing. J Electron Mater 41: 2083-2088

[14] Tu K-N, Liu Y (2019) Recent advances on kinetic analysis of solder joint reactions in 3D IC packaging technology. Mater Sci Eng R Rep 136: 1-12
[15] Kelly MB, Niverty S, Chawla N (2020) Four dimensional (4D) microstructural evolution of Cu6Sn5 intermetallic and voids under electromigration in bi-crystal pure Sn solder joints. Acta Mater 189: 118-128

[16] Ahmed MT, Motalab M, Suhling JC (2021) Impact of Mechanical Property Degradation and Intermetallic Compound Formation on Electromigration-Oriented Failure of a Flip-Chip Solder Joint. J Electron Mater 50: 233-248

[17] Lee C-C, Lin Y-M, Liu H-C, Syu J-Y, Huang Y-C, Chang T-C (2021) Reliability evaluation of ultra thin 3D-IC package under the coupling load effects of the manufacturing process and temperature cycling test. Microelectron Eng 244: 111572

[18] Qi H, Osterman M, Pecht M (2009) A rapid life-prediction approach for PBGA solder joints under combined thermal cycling and vibration loading conditions. IEEE Trans Compon Packag Technol 32: 283-292

[19] Wong EH, Wong CK (2009) Approximate solutions for the stresses in the solder joints of a printed circuit board subjected to mechanical bending. Int J Mech Sci 51: 152-158

[20] Chen H, Wang L, Han J, Li M, Liu H (2012) Microstructure, orientation and damage evolution in SnPb, SnAgCu, and mixed solder interconnects under thermomechanical stress. Microelectron Eng 96: 82-91

[21] Liang S, Ke C, Huang J, Zhou M, Zhang X (2019) Phase field simulation of microstructural evolution and thermomigration-induced phase segregation in Cu/Sn58Bi/Cu interconnects under isothermal aging and temperature gradient. Microelectron Reliab 92: 1-11

[22] Yang TL, Aoki T, Matsumoto K et al (2016) Full intermetallic joints for chip stacking by using thermal gradient bonding. Acta Mater 113: 90-97
[23]Liu B, Tian Y, Feng J, Wang C (2017) Enhanced shear strength of Cu–Sn intermetallic interconnects with interlocking dendrites under fluxless electric current-assisted bonding process. J Mater Sci 52: 1943-1954

[24]Jung D, Sharma A, Mayer M, Jung J (2018) A review on recent advances in transient liquid phase (TLP) bonding for thermoelectric power module. Rev Adv Mater 53: 147-160

[25]Chen C, Tong H, Tu K-N (2010) Electromigration and thermomigration in Pb-free flip-chip solder joints. Annu Rev Mater Res 40: 531-555

[26]Chen C, Hsiao H-Y, Chang Y-W, Ouyang F, Tu K-N (2012) Thermomigration in solder joints. Mater Sci Eng R Rep 73: 85-100

[27]Wu AT, Tu K-N, Lloyd J, Tamura N, Valek B, Kao C (2004) Electromigration-induced microstructure evolution in tin studied by synchrotron x-ray microdiffraction. Appl Phys Lett 85: 2490-2492

[28]Wu AT, Hsieh Y (2008) Direct observation and kinetic analysis of grain rotation in anisotropic tin under electromigration. Appl Phys Lett 92: 121921

[29]Zhao N, Zhong Y, Huang M, Ma H, Dong W (2016) Dissolution and precipitation kinetics of Cu6Sn5 intermetallics in Cu/Sn/Cu micro interconnects under temperature gradient. Intermetallics 79: 28-34

[30]Deng S, Li R, Yuan T, Yang P, Xie S, Li J (2020) The preferential growth behaviors of the intermetallics at the W/Co interface during spark plasma sintering. Appl Phys Lett 117: 194102

[31]Liu B, Tian Y, Wang C, An R, Wang C (2017) Ultrafast formation of unidirectional and reliable Cu3Sn-based intermetallic joints assisted by electric current. Intermetallics 80: 26-32

[32]Feng J, Hang C, Tian Y, Wang C, Liu B (2018) Effect of electric current on grain orientation and mechanical properties of Cu-Sn intermetallic compounds joints. J Alloys Compd 753: 203-211

[33]Huh JY, Hong K, Kim Y, Kim K (2004) Phase field simulations of intermetallic compound growth during soldering reactions. J Electron Mater 33: 1161-1170
[34] Hong K, Huh JY (2006) Phase field simulations of morphological evolution and growth kinetics of solder reaction products. J Electron Mater 35: 56-64

[35] Park M, Arróyave R (2010) Early stages of intermetallic compound formation and growth during lead-free soldering. Acta Mater 58: 4900-4910

[36] Moelans N (2011) A quantitative and thermodynamically consistent phase-field interpolation function for multi-phase systems. Acta Mater 59: 1077-1086

[37] Ke J, Gao Y, Kao C, Wang Y (2016) Pattern formation during interfacial reaction in-between liquid Sn and Cu substrates–A simulation study. Acta Mater 113: 245-258

[38] Raza SH, Klusemann B (2020) Multiphase-field modeling of temperature-driven intermetallic compound evolution in an Al–Mg system for application to solid-state joining processes. Model Simul Mat Sci Eng 28: 085003

[39] Tang S, Wang J, Yang G, Zhou Y (2011) Phase field modeling the growth of Ni3Al layer in the β/γ diffusion couple of Ni–Al binary system. Intermetallics 19: 229-233

[40] Kunwar A, Yousefi E, Zuo X, Sun Y, Seveno D, Guo M, Moelans N (2022) Multi-phase field simulation of Al3Ni2 intermetallic growth at liquid Al/solid Ni interface using MD computed interfacial energies. Int J Mech Sci 215: 106930

[41] Park M, Gibbons S, Arróyave R (2013) Phase-field simulations of intermetallic compound evolution in Cu/Sn solder joints under electromigration. Acta Mater 61: 7142-7154

[42] Hektor J, Ristinmaa M, Hallberg H, Hall SA, Iyengar S (2016) Coupled diffusion-deformation multiphase field model for elastoplastic materials applied to the growth of Cu6Sn5. Acta Mater 108: 98-109

[43] Attari V, Arróyave R (2020) Phase-Field Study of Thermomigration in 3-D IC Micro Interconnects. IEEE Trans Compon Packaging Manuf Technol 10: 1466-1473
[44] Kunwar A, Hektor J, Nomoto S, Coutinho YA, Moelans N (2020) Combining multi-phase field simulation with neural network analysis to unravel thermomigration accelerated growth behavior of Cu6Sn5 IMC at cold side Cu–Sn interface. Int J Mech Sci 184: 105843

[45] Zhou P, Johnson WC (2011) A diffuse interface model of intermediate-phase growth under the influence of electromigration. J Electron Mater 40: 1867-1875

[46] Xiong H, Huang Z, Conway P (2014) Effects of stress and electromigration on microstructural evolution in microbumps of three-dimensional integrated circuits. IEEE Trans Device Mater Reliab 14: 995-1004

[47] Attari V, Ghosh S, Duong T, Arroyave R (2018) On the interfacial phase growth and vacancy evolution during accelerated electromigration in Cu/Sn/Cu microjoints. Acta Mater 160: 185-198

[48] Kunwar A, Coutinho YA, Hektor J, Ma H, Moelans N (2020) Integration of machine learning with phase field method to model the electromigration induced Cu6Sn5 IMC growth at anode side Cu/Sn interface. J Mater Sci Technol 59: 203-219

[49] Beyne S, Pedreira OV, Wolf ID, Tökei Z, Croes K (2019) A novel electromigration characterization method based on low-frequency noise measurements. Semicond Sci Technol 34: 075002

[50] Liang S, Kunwar A, Wei C, Ke C (2021) Insight into the preferential grain growth of intermetallics under electric current stressing–A phase field modeling. Scr Mater 203: 114071

[51] Raissi M, Perdikaris P, Karniadakis GE (2019) Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J Comput Phys 378: 686-707

[52] Kunwar A, An L, Liu J, Shang S, Råback P, Ma H, Song X (2020) A data-driven framework to predict the morphology of interfacial Cu6Sn5 IMC in SAC/Cu system during laser soldering. J Mater Sci Technol 50: 115-127
[53] Pun G, Batra R, Ramprasad R, Mishin Y (2019) Physically informed artificial neural networks for atomistic modeling of materials. Nat Commun 10: 1-10

[54] Mishin Y (2021) Machine-learning interatomic potentials for materials science. Acta Mater 214: 116980

[55] Yao P, Li X, Liang X, Yu B, Jin F, Li Y (2017) A study on interfacial phase evolution during Cu/Sn/Cu soldering with a micro interconnected height. Mater Charact 131: 49-63

[56] Jin Y, Artemev A, Khachaturyan A (2001) Three-dimensional phase field model of low-symmetry martensitic transformation in polycrystal: simulation of ζ' 2 martensite in AuCd alloys. Acta Mater 49: 2309-2320

[57] Chen L-Q, Yang W (1994) Computer simulation of the domain dynamics of a quenched system with a large number of nonconserved order parameters: The grain-growth kinetics. Phys Rev B 50: 15752

[58] Moelans N, Blanpain B, Wollants P (2008) Quantitative analysis of grain boundary properties in a generalized phase field model for grain growth in anisotropic systems. Phys Rev B 78: 024113

[59] Liang S, Ke C, Wei C, Zhou M, Zhang X (2018) Phase field modeling of grain boundary migration and preferential grain growth driven by electric current stressing. J Appl Phys 124: 175109

[60] Tonks MR, Zhang Y, Bai X, Millett PC (2014) Demonstrating the temperature gradient impact on grain growth in UO2 using the phase field method. Mater Res Lett 2: 23-28

[61] Flint TF, Panwisawas C, Sovani Y, Smith MC, Basoalto H (2018) Prediction of grain structure evolution during rapid solidification of high energy density beam induced re-melting. Mater Des 147: 200-210

[62] Newnham RE (2005) Properties of materials: anisotropy, symmetry, structure. Oxford University Press, Oxford

[63] Wang S, Hsu L, Wang N, Ho C (2014) EBSD investigation of Cu-Sn IMC microstructural evolution in Cu/Sn-Ag/Cu microbumps during isothermal annealing. J Electron Mater 43: 219-228
[64] McCurdy A, Maris H, Elbaum C (1970) Anisotropic heat conduction in cubic crystals in the boundary scattering regime. Phys Rev B 2: 4077

[65] Liu C, Hu Y, Liu Y et al (2013) Epitaxial Cu–Sn bulk crystals grown by electric current. Acta Mater 61: 5713-5719

[66] Chao B, Chae S-H, Zhang X, Lu K-H, Im J, Ho PS (2007) Investigation of diffusion and electromigration parameters for Cu–Sn intermetallic compounds in Pb-free solders using simulated annealing. Acta Mater 55: 2805-2814

[67] Khosla A, Huntington H (1975) Electromigration in tin single crystals. J Phys Chem Solids 36: 395-399

[68] Gottstein G, Shvindlerman LS (2009) Grain boundary migration in metals: thermodynamics, kinetics, applications. CRC Press, Boca Raton

[69] Durga A, Wollants P, Moelans N (2020) Phase-field study of IMC growth in Sn–Cu/Cu solder joints including elastoplastic effects. Acta Mater 188: 241-258

[70] Gan H, Tu K-N (2005) Polarity effect of electromigration on kinetics of intermetallic compound formation in Pb-free solder V-groove samples. J Appl Phys 97: 063514

[71] Haessner F, Hofmann S, Seekel H (1974) Influence of high electric direct current densities on grain boundary migration in rolled gold single crystals. Scr Metall 8: 299-305

[72] Thouless M (1996) Stress evolution during electromigration in a bamboo structure. Scr Mater 34: 1825-1831

[73] Subedi U, Kunwar A, Coutinho YA, Gyanwali K (2021) pyMPEALab Toolkit for Accelerating Phase Design in Multi-principal Element Alloys. Metals and Materials International 1-13

[74] Abadi M, Barham P, Chen J et al (2016) TensorFlow: A System for Large-Scale Machine Learning. 12th USENIX symposium on operating systems design and implementation (OSDI 16), pp 265-283
[75] Chen J-Q, Guo J-D, Ma H-C, Liu K-L, Zhu Q-S, Shang JK (2015) Magnetic-field induced anisotropy in electromigration behavior of Sn-Ag-Cu solder interconnects. J Mater Res 30: 1065-1071

[76] Kumar S, Handwerker CA, Dayananda MA (2011) Intrinsic and interdiffusion in Cu-Sn system. J Phase Equilib Diffus 32: 309-319

[77] Yuan Y, Guan Y, Li D, Moelans N (2016) Investigation of diffusion behavior in Cu-Sn solid state diffusion couples. J Alloys Compd 661: 282-293

[78] Deng S, Li R, Yuan T, Cao P (2018) Effect of electric current on crystal orientation and its contribution to densification during spark plasma sintering. Mater Lett 229: 126-129

[79] Xian J, Belyakov S, Ollivier M, Nogita K, Yasuda H, Gourlay C (2017) Cu6Sn5 crystal growth mechanisms during solidification of electronic interconnections. Acta Mater 126: 540-551