The increasing demand of electronic products requires highly miniaturized integrated circuits (ICs). These ICs cause the chip integration density to grow sharply and lead to a higher current density in the solder joints, which will induce the phenomenon of electromigration (EM) and an enhanced formation of intermetallic compound (IMCs) in the solder joints. For the widely used Sn-Ag-Cu (SAC) series of solder, Cu atoms in Cu pads and Sn atoms in SAC solder joints may combine together to form layer-shaped IMCs (mainly Cu₃Sn and Cu₆Sn₅) during EM and the thermal ageing process. As miniaturized solder joints contain only several large bulk grains, the performance of the solder grains and IMCs has a great effect on the solder reliability, which has attracted great attention from the research community (Hau-Riege and Thompson, 2001). At the working temperature and current density, the dominant failure mode under EM is the mechanical-loading-induced crack initiation and extension, and experimental results have shown that the most common IMC is Cu₃Sn under such working conditions (Zhang et al., 2015). Cu₃Sn has great prospect in power electronic device applications due to its high melting temperature (Liu et al., 2020), and the solder joints with the interfacial structure of Cu/Cu₃Sn/Cu is thought to be a substitute due to the ability of servicing under high temperature (Yao et al., 2017). While for this kind of structure, cracks can easily form between Cu₃Sn and Cu (Tseng et al., 2010), or within Cu₃Sn (Yao et al., 2019), which affects the reliability of the solder joints. Thus, it is important to investigate the mechanical reliability of the Cu₃Sn/Cu interface.

Due to IC manufacturing, the Cu pad always contains heterogeneous bulk single grains or polycrystalline Cu (Hu et al., 2007), and the mechanical properties vary depending on the grain structure. It is possible to control the grain orientation and distribution in crystalline materials (Todaro et al., 2020) for the purpose of improving the reliability of packaging. Works about mechanical properties of Cu and Cu3Sn with different combination of grain structure and orientation have been reported by experimental research separately (Yang et al., 2017; Juang et al., 2021; Shen et al., 2022), both of them show the great difference between them. while the manner in which the grain structure could affect the reliability of the interface between the pad and the solder joints is an important research topic in advanced electronic packaging (Zhang et al., 2018). The development of the molecular dynamics (MD) theory has enabled the microscale simulation of phenomena regarding the influence of the degree of crystallinity on the interface reliability (Zhang and Asle Zaeem, 2018), and uniaxial tension simulations can be used to evaluate the mechanical properties of materials (Cheng et al., 2012; Mackenchery et al., 2016). In this work, the mechanical reliability of the Cu₃Sn/Cu interface is investigated through MD tension simulations.

The significant differences in the mechanical properties between single crystal and polycrystalline Cu (Liu et al., 2005; Qiu et al., 2017) affect the reliability of the Cu₃Sn/Cu interface considerably. On the other hand, nanotwinned Cu has many unique mechanical characteristics, which can have a dramatic effect on the solder reliability (Fu et al., 2018). Furthermore, the (111) grain orientation has been proved to be beneficial in improving EM (Cheng et al., 2018). Hence, models based on single-grain Cu, twin-crystal Cu with (111) twin grain interface, and polycrystalline Cu were implemented and compared. Additionally, three different strain rates were applied to the interface models to evaluate the interface strength and strain rate response under uniaxial tension. Finally, the failure mechanism of the Cu₃Sn/Cu interface structure was analyzed from an atomic perspective.

The MD simulations were performed using a Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) (Plimpton, 1995). The modified embedded atom method (MEAM) interatomic potential was used to describe the interactions between the Cu and Sn atoms. The MEAM interatomic potential has obtained good results in many fields (Xu and Kim, 2019; Attarian and Xiao, 2022). The parameters of the MEAM potential and the lattice structure of Cu₃Sn and Cu can be found in Refs. (Baskes et al., 1994; Aguilar et al., 2000; Cheng et al., 2012). and are shown in Table 1.

The crystal structure of Cu₃Sn was reported by Watanabe et al. (Watanabe et al., 1983) that it is a hexagonal close-packed Cu₃Ti-type transition-metal IMC. Villars and Calvet (Villars and Calvet, 1985) gave the complete description of the atomic positions within the crystal structure of Cu₃Sn. its structure could be treated as a supercell consisted by eight atoms with orthorhombic symmetry (space group Pmmn), which was proved to have a better measurement from crystallographic studies (Pang et al., 2008). In order to validate the accuracy of MEAM interatomic potentials used for the Cu₃Sn structure, the orthorhombic elastic constants of Cu₃Sn was calculated by applying a strain field to the simulation cell, as shown in Table 2. It can be seen our MD results are in good agreement with the ab initio calculations.

For the problems discussed in this domain, to better describe the different grain structures of Cu, single-crystal Cu, twin-crystal Cu, and polycrystalline Cu as the pad layer were investigated. Two types of twin copper models with one twin grain boundary (has two grains) and four twin grain boundaries (has five grains) were built with the (111) surface as the twin plane to study the effect of grain size. Furthermore, the polycrystalline Cu model contained 10 grains with random orientation. Both of these four types of Cu combine with single-crystal Cu₃Sn to form the interface structure for the axial tension simulations; the interface consisted of the (001) surface of Cu₃Sn combined with the (112) surface of single-crystal Cu and twin-crystal Cu. The total size of the interface structure built for the MD simulations was 110.58 × 90.6 × 180.72 Å with the same heights of Cu and Cu₃Sn, as shown in Figure 1.

For the axial tension simulations performed in this study, asymmetric boundary conditions were used in each direction. The Cu₃Sn and Cu layers were placed next to each other at a distance of 3 Å, and all structures were then equilibrated to obtain the balanced models, structures were preheated from 100 to 300 K for 50 ps under the condition of NVE (constant Number of atoms, Volume and Energy) with a berendsen thermostat. After that 34.6 Å for both edges along z-axial of models were set as rigid bodies; then, uniform velocities at strain rates of 0.01, 0.1, and 0.5% ps⁻¹ were applied to the rigid region of Cu₃Sn along the direction perpendicular to the interface. Virial stress was used to calculate the tensile stress at the atomic level according to (Spearot et al., 2004; Dandekar and Shin, 2011), which could be of great help to get the von Mises stress as introduced in (Rahmati et al., 2020). The expression of virial stress is defined as (Spearot et al., 2004): Π α β i = 1 Ω i [ 1 2 ∑ j ≠ i ( 1 r i j ∂ U ∂ r ) r α i j r β i j − m i v α i v β i ] , (1) With σ α β = 1 N ∗ ∑ i N ∗ Π α β i , (2) Where the first term within the square brackets in Eq. 1 accounts for the force between two atoms, as a function of the potential energy U, and the second term accounts for the kinetic contribution to the virial stress. In the above equations, m i , v α i and v β i are the mass and velocity of the ith atom in the α and β directions respectively; r α i j and r β i j are α and β components of total distance r i j between atoms i and j. the atomic dipole force tensor Π α β i is averaged over volume that the atom occupies in a unit cell. Then the global stress tensor σ α β is defined as the average over the number of atoms that in the region to be calculated. As the tension simulation was implemented in this paper, virial stress tensor of σ 33 was extracted for comparison which is along the loading direction.

After the simulations, virial stress versus displacement curves were extracted for the different cases, as shown in Figure 2. It can be directly seen that the polycrystalline Cu model show a great ductile character under all strain rates, and the maximum stress is considerably lower than that obtained for the other models. Furthermore, the maximum stress obtained for the single-crystal Cu models is slightly larger than that obtained for the twin-crystal Cu models at a strain rate of 0.5% ps⁻¹ (8,150 MPa for the single-crystal Cu models and 7,700 MPa for the twin-crystal Cu models). However, as the strain rate is reduced to 0.1 and 0.01% ps⁻¹, this difference decreases, and Figure 2 shows that the number of twin grain boundaries does not change significantly as the virial stress exhibits the same trend for all strain rates investigated in this work.

On the other hand, it can be seen from Figure 2 that the displacement at which the stress reaches the highest peak for the single-crystal Cu model is always larger than that of the other models for all strain rates. Since the difference in the maximum stress obtained for the single-crystal Cu and twin-crystal Cu models is not large at a low strain rate, we investigated the elastic stage of these two cases. It can be seen that the single-crystal Cu and polycrystalline Cu models are characterized by stress curves that always have the same slope. The same is true for the twin-grain-boundary models. Furthermore, the slope of the stress curve obtained for the twin-grain-boundary models is higher than that obtained for the other two models. At a strain rate of 0.5% ps⁻¹, the stress curve of the polycrystalline Cu model follows that of the twin-crystal Cu model at first, while, after a small peak, it follows that of the single-crystal Cu model. To investigate the mechanism of this phenomenon, the work of adhesion and dislocation distribution were studied for both cases.

The work of adhesion is typically used to describe the adhesion strength of an interface. In this case, the work of adhesion (W _Cu/Cu3Sn ) was obtained as the energy change when splitting the Cu/Cu₃Sn interface into the two Cu and Cu₃Sn surfaces. The value of W _Cu/Cu3Sn was calculated as follows (AlMotasem et al., 2016): W C u / C u 3 S n = ( E C u + E C u 3 S n − E t o t ) / S inter , (3) where E _Cu and E _Cu3Sn correspond to the potential energies of Cu and Cu₃Sn in a system with a free surface, respectively. E _tot represents the total potential energy of the whole model, and S _inter is the surface-interactive area of the model. The work of adhesion calculated for both models is shown in Table 3.

It can be seen that the work of adhesion between Cu and Cu₃Sn increases due to the existence of a twin crystal boundary, and more twin crystal boundaries in the same Cu region result in a higher work of adhesion. Considering the maximum stress illustrated in Figure 2, the single-crystal Cu model shows higher stresses than the other models, which appears in contradiction with our previous conclusions. However, it is noteworthy that the work of adhesion obtained for the polycrystalline Cu model is 2.89 J/m², which is twice larger than that of the other models, while the maximum stresses at different strain rates are substantially lower than those of the other models. This indicates that the existence of twin crystal boundaries or even ordinary grain boundaries has a great influence on the reliability of the Cu₃Sn/Cu interface. Thus, the work of adhesion alone cannot explain the reliability of interfaces that contain grain boundaries. Indeed, the mechanical strength of each layer changes considerably in the presence of grain boundaries as these cause a change in the defects in each layer during loading, which directly affects the stress evolution.

To investigate the mechanism of the phenomenon introduced above, the atomic position snapshot for all interface models was extracted at a displacement of 60 Å, where the structures are close to undergoing total failure. Additionally, to derive the defect evolution during the tensile process, the centro-symmetry parameter (CSP) introduced by Kelchner et al. (2011) was used to visualize the defects (dislocations, stacking faults, and nanotwins) during deformation. The CSP for a given atom is defined as follows: CSP = ∑ i = 1 N / 2 | R i − R i + N / 2 | 2 , (4) where N is the number of nearest neighbors (as Cu has the face-centered cubic (FCC) structure, N = 12) for the underlying lattice of atoms, and R _i and R _i+N/2 are the vectors corresponding to a specific pair of nearest neighbors opposite to the central atom. For a perfect lattice, the value of CSP is zero, whereas, for the regions characterized by defects or atoms close to the free surfaces, the CSP value is positive. A large CSP value indicates that the symmetry is broken. The OVITO software (Stukowski, 2009) was used for the CSP analysis and the dislocation extraction algorithm (DXA) (i.e., the dislocation analysis module in the software) was used to visualize the extracted dislocations and the corresponding Burgers vectors, which are indicated by lines in Figure 3.

It is noticeable from Figure 3A that the degeneration position is located at the Cu layer that is close to the interface for the single-crystal Cu model, while the existence of twin grain boundaries causes the degeneration position to move toward the Cu₃Sn layer: for the model that contains only one twin boundary, the degeneration position is located at the interface, while when four twin boundaries are considered, the degeneration position is located inside the Cu₃Sn layer, as shown in Figures 3B,C. Regarding the polycrystalline Cu model, although this model exhibits the maximum work of adhesion, its mechanical strength is lower than that of the other models; thus, the failure of this model is predominantly caused by the collapse of polycrystalline Cu, as shown in Figure 3D. Furthermore, the snapshot shows the great ductile character of this model, which explains why its stress evolution curve has a lower maximum stress than the other models at all strain rates.

The model with one twin grain boundary has a lower work of adhesion than the model with four twin grain boundaries, which indicates that the adhesion between Cu and Cu₃Sn is weaker; the lower work of adhesion leads to the structure failure at the interface first. By contrast, the maximum stress of the single-crystal Cu model is slightly higher than that of the other models, as shown in Figure 2A; however, this model has the lowest work of adhesion, and its degeneration position is located in the Cu region, not at the interface. Since this phenomenon does not occur at strain rates of 0.1 and 0.01% ps⁻¹, it may be attributed to the fact that the rate dependence of mechanical properties that affects the atom distribution during deformation. Indeed, Figure 2A shows that the stress evolution curve of this model exhibits a peak before reaching the maximum stress, which indicates that several Cu atoms changed their position at this stage, thereby affecting the interface structure.

From the DXA results shown in Figure 3, it can be observed that a large number of Shockley dislocations with a Burgers vector of 1/6 [112] appeared in the single-crystal Cu model, and the existence of twin grain boundaries resulted in the main type of dislocation being the stair-rod dislocation with a Burgers vector of 1/6 [110]. Both dislocations types formed around the twin grain boundaries and near the interface; when considering four twin grain boundaries in the Cu layer, only a small number of dislocations could be found along the interface. As twin grain boundaries enhance the mechanical strength (Gao et al., 2012), once the mechanical strength exceeds the interface strength, a degeneration will occur at the interface. This is why a lower number of dislocations appeared in the twin-crystal Cu model, and the degeneration positions of each model are different.

On the other hand, the mechanical behavior of the Cu/Cu₃Sn interface depends considerably on the strain rate. Not only the stress evolution (shown in Figure 2) but also the degeneration position was significantly influenced by the strain rate, especially for the single-crystal Cu model. By comparing the degeneration positions of the single-crystal Cu model at strain rates of 0.01 and 0.1% ps⁻¹, as shown in Figure 4, it is clear that a lower strain rate leads the structure failure at the Cu layer near the interface, while a higher strain rate results in the structure failure inside the Cu₃Sn region. This indicates that Cu₃Sn is highly affected by the strain rate, and a higher strain rate causes the brittle failure of Cu₃Sn, while a lower strain rate results in the degeneration position being located at the interface for the twin-crystal Cu models and inside the Cu layer for the single-crystal Cu model.

The effect of the grain structure on the failure mechanism of the Cu₃Sn/Cu interface was investigated via MD simulations. Single-crystal Cu, twin-crystal Cu, and polycrystalline Cu models were built and compared. Axial tension simulations were performed at strain rates of 0.01, 0.1, and 0.5% ps⁻¹. The stress evolution and work of adhesion were then calculated and compared for each case to evaluate the mechanical reliability of the interface. Furthermore, both CSP and DXA were analyzed to investigate the interface failure mechanism from the perspective of defect evolution.

The stress versus displacement curves indicates that the ductile failure of polycrystalline Cu lead polycrystalline Cu model has a lower stress level and a great ductile character, which is induced by the sliding of grain boundaries. Additionally, the number of twin grain boundary has less influence on the interfacial character under the strain rates of 0.01 and 0.1% ps⁻¹. On the other hand, more twin crystal boundaries in the same Cu region would result in a higher work of adhesion, while with the existence of grain boundaries, only by comparing the work of adhesion cannot be used to evaluate the mechanical reliability of the interface. Investigation of the failure position and calculation of CSP and DXA could be of great help to gain a deeper understanding of this problem. By combining with the atom position, CSP and DXA results, it can be concluded that twin grain boundaries have the ability to enhance the mechanical strength of Cu layer, and when this mechanical strength exceeds the interfacial strength, degeneration will occur at the interface, this is why a lower number of dislocations appeared in the twin-crystal model. Results also show that the strain rate have great influence on mechanical character of Cu₃Sn. A high strain rate induces the brittle failure of Cu₃Sn, while a low strain rate causes the degeneration position to be located at the interface for the twin-crystal Cu models and inside the Cu layer for the single-crystal Cu model. In summary, our results provide information about the failure mechanism of the Cu₃Sn/Cu interface from MD simulation, which will be helpful for applications in the electronic packaging industry.