Atomistic simulations of Ag–Cu–Sn alloys based on a new modified embedded-atom method interatomic potential

An interatomic potential for the ternary Ag–Cu–Sn system, an important material system related to the applications of lead-free solders, is developed on the basis of the second nearest-neighbor modified embedded-atom-method formalism. Potential parameters for the ternary and related binary systems are determined based on the recently improved unary description of pure Sn and the present improvements to the unary descriptions of pure Ag and Cu. To ensure the sufficient performance of atomistic simulations in various applications, the optimization of potential parameters is conducted based on the force-matching method that utilizes density functional theory predictions of energies and forces on various atomic configurations. We validate that the developed interatomic potential exhibits sufficient accuracy and transferability to various physical properties of pure metals, intermetallic compounds, solid solutions, and liquid solutions. The proposed interatomic potential can be straightforwardly used in future studies to investigate atomic-scale phenomena in soldering applications.


Introduction
Soldering is an important technique that critically determines the reliability of electronic devices. Among various types of solder materials, the most popular and traditional materials are the Sn-Pb-based alloys. However, due to environmental and toxicity concerns, the electronics industry has pursed the use of lead-free solders. Among alloy systems applicable for soldering applications, Sn-based alloys, especially binary Sn-Ag (SA) [1] Atomistic simulations of Ag-Cu-Sn alloys based on a new modified embedded-atom method interatomic potential formation of intermetallic compounds and Kirkendall voids [7]. In practical electronic devices, the SA and SAC alloys are usually bonded to Cu. During the soldering process, interdiffusion between the solder alloy and Cu occurs, which enables a strong connection between them. However, such interdiffusion phenomena often lead to the formation of intermetallic compounds at the interface during possible reflows and/or aging processes. In the SA and SAC alloys, two intermetallic compounds typically formed at the interface are Cu 6 Sn 5 on the solder side and Cu 3 Sn on the Cu side [8,9]. If the solder is the SA alloys, the Ag 3 Sn compound may also be observed depending on the processing conditions [7]. The formation of these intermetallic compounds often accompanies the Kirkendall voids that appear around or within those compounds [10]. This phenomenon critically deteriorates the mechanical reliability of electronic devices, but the underlying mechanisms depending on specific alloying and processing conditions are still questionable.
Moreover, there has been increasing attention on the properties of lead-free solders at small scales since the sizes of the electrical devices steadily decreased. For example, the possible melting point depression of nano-sized solder alloys is of interest [11] because the melting point of the usual bulk SAC alloys is ~ 34 K higher than that of the bulk Sn-Pb-based solder alloys [4]. The wettability of small-sized solders on copper substrates is also of particular interest because the synthesis of small-sized solder particles with a binary or ternary alloy of uniform composition is much more difficult than that of monolithic particles and larger solder balls [12,13].
To obtain a detailed understanding of related phenomena in SA and SAC alloys for soldering applications, simulations at small scales are an effective complement to experimental characterizations. A widely used simulation method at the small scale is the first-principles density functional theory (DFT) calculation. This method has a significant advantage of being able to easily access a variety of multi-component systems without the need for developing individual databases. Several DFT studies [14,15] were already conducted on the structural and elastic properties of related intermetallic compounds in the SA and SAC alloy systems. However, there are several target phenomena that cannot be solely covered by the DFT calculation due to the size limitation. Alternatively, atomistic simulations such as molecular dynamics (MD) and Monte Carlo (MC) simulations can be properly utilized when relatively large-scale phenomena are involved. However, such atomistic simulations require the development of individual interatomic potentials of target systems, and the accuracy of predictions are critically determined by the performance of interatomic potentials.
Although atomistic simulations on the SA and SAC alloys on the soldering applications are highly desirable, relevant interatomic potentials are currently not available. There were two reported potentials for the Ag-Cu-Sn system based on the modified embedded-atom method (MEAM) model [16] presented in the previous works by Sellers et al. [17] and Motalab et al. [18]. For the constituent binary systems, there were reported MEAM potentials for the Cu-Sn binary system suggested by Aguilar et al. [19], Cheng et al. [20], and Liang et al. [21]. A MEAM potential for the Ag-Sn binary system was suggested by Dong et al. [22]. In previous works, the binary and ternary MEAM potentials were developed based on the MEAM potential for the unary Sn system proposed by Ravelo and Baskes [23]. The reason that all the previous studies have attempted to develop the interatomic potential based on the MEAM model seems to be due to the difficulty of describing the properties of pure Sn. This unary system exhibits an allotropic phase transformation between two solid phases (diamond (α) ↔ body-centered tetragonal (β)), which is difficult to deal with using simpler and more computationally efficient potential models such as the embedded-atom method (EAM). Because the MEAM model was originally devised to describe various types of atomic bonds such as the metallic and covalent bonds in a single formalism, it seems reasonable to apply this model to Sn-based alloy systems.
However, the MEAM potential for pure Sn developed by Ravelo and Baskes [23], which underlies the previous development of alloy potentials [17][18][19][20][21][22], exhibits deficiencies in reproducing properties of the β-Sn phase that is stable at ambient conditions [24]. Moreover, the previous MEAM potentials for alloy systems [17][18][19][20][21][22] were not developed considering various physical properties of binary and ternary alloy systems, and the optimization of potential parameters was not systematically performed. For example, the binary potential parameters were determined by only focusing on selected properties (e.g., lattice constant and bulk modulus) of hypothetical L1 2 compounds, which are directly associated with the adjustable binary potential parameters of the MEAM model. In previous studies [17][18][19][20][21][22], the properties of various stable intermetallic compounds, solid solutions, and liquid solutions were not considered at all during either fitting or evaluation of potential parameters.
In the present study, we present a new interatomic potential for the Ag-Cu-Sn ternary system that enables atomistic simulations of SA and SAC alloys for soldering applications. The potential was developed based on the second nearestneighbor modified embedded-atom method (2NN MEAM) model [25][26][27][28][29]. This is an improved model that overcomes several shortcomings of the original MEAM approach [16] and has been reported to exhibit good accuracy and transferability for various alloy systems. The potential was developed based on a relevant optimization process of potential parameters to ensure sufficient reliability. Specifically, we used the forcematching method proposed by Ercolessi and Adams [30] for the optimization of potential parameters. According to this method, forces and energies associated with various atomic Invited Feature Paper configurations, including configurations at finite temperatures, should first be obtained through DFT calculations. The optimization of potential parameters was then performed by minimizing the difference between expected forces and energies from the potential parameters and the DFT calculations. It has been reported that this method can significantly improve the overall performance of developed interatomic potentials [31][32][33]. In this study, we particularly extend our recent work [24] on the development of an interatomic potential for pure Sn-based on the force-matching method. We will present a detailed process for the development of interatomic potential that realizes the atomistic simulations of the Ag-Cu-Sn ternary system. The performance of the developed potential will be presented by comparing the prediction by the optimized potential parameters with the experimental and DFT data.

Required parameters in the 2NN MEAM potential model
In the 2NN MEAM potential formalism, 14 independent parameters are required to describe a pure element. There are four parameters related to the universal equation of state, e.g., the cohesive energy ( E c ), the equilibrium nearest-neighbor distance ( r e ), the bulk modulus (B) of the reference structure, and the adjustable parameter d. Other 7 parameters are related to the electron density, e.g., the decay lengths ( β (0) , β (1) , β (2) , β (3) ) and the weighting factors ( t (1) , t (2) , t (3) ). Of the remaining three parameters, the parameter A belongs to the embedding function and the parameters C min and C max are responsible for the manybody screening. For a binary alloy system, the 2NN MEAM potential formalism requires 13 independent parameters in addition to the parameters of constituent unary systems. There are four parameters related to the universal equation of state, e.g., E c , r e , B, and d. Another parameter, the atomic electron density scaling factor ρ 0 , corresponds to the electron density. The remaining eight parameters are responsible for the manybody screening between two different types of atoms (four C min and four C max ). To describe a ternary alloy, additional 6 parameters are required in addition to the parameters for constituent unary and binary systems. These parameters are responsible for the many-body screening between three different types of atoms (three C min and three C max ). Detailed descriptions of the 2NN MEAM formalism are provided in the previous studies [25][26][27][28][29].

Optimization of potential parameters
We extended our previous optimization of the potential parameters for pure Sn [24] to obtain the potential parameters for the ternary Ag-Cu-Sn. We determined remaining unary (Ag and Cu), binary (Ag-Cu, Ag-Sn, and Cu-Sn), and ternary (Ag-Cu-Sn) parameters by fitting to the established DFT database of energies and forces. The optimization was conducted by comparing energies of target configurations and forces on each atom as expected by the DFT calculation and those expected by a candidate set of potential parameters. Fitting error, defined as the difference between expectations of the DFT and MEAM, was minimized using an in-house optimization code based on a genetic algorithm for efficient searches in high-dimensional spaces of potential parameters. An initial population of candidate solutions were randomly generated and evolved toward better solutions for further iterations. The population size (36) was greater than the required 2NN MEAM potential parameters for pure (14), binary (13), and ternary (6) systems. The optimization started by specifying a radial cutoff distance, reference structures for target unary and binary systems, and fitting weights of each atomic configuration. The radial cutoff distance of 5.0 Å, which was previously reported to be large enough to reproduce various physical properties of pure Sn [24], was also used in this study. The reference structures for unary (Ag and Cu) and binary systems (Ag-Cu, Ag-Sn, Cu-Sn) were determined to be fcc and CsCl-type B2 structures, respectively. If the obtained temporary set of potential parameters could not satisfactorily reproduce the overall properties of the target alloy system, another optimization process was considered using different fitting weights of atomic configurations.
During the optimization process, we realized that a hierarchical force-matching process, i.e., the sequential and separate fitting of parameters (unary → binary → ternary), is not very successful for obtaining reliable interatomic potentials for the target multi-component system. For example, if we conduct an optimization of binary parameters after the decision of constituent unary parameters, the performance of the potential to reproduce the properties of binary systems was much lower than expected compared to the performance of the potential to reproduce the properties of unary systems. As reported in the previous study [33], the alloy potential parameters presented in the MEAM formalism do not provide full flexibility for the force-matching process because the description of pair interactions is determined from the universal equation of state without involving a specific functional form.
We overcame this difficulty by a strategy of fitting unary and binary parameters simultaneously. During the optimization, we found that the potential model exhibits significantly different performances of fittings for each unary system. For example, we can obtain a number of relevant candidate parameter sets for pure Ag and Cu, but optimization of the potential for pure Sn was a relatively difficult task. Therefore, we used a previously determined potential parameter set for pure Sn [24], which was obtained from a separate force-matching process targeting a specific unary system. For pure Ag and Cu, we did not perform separate fitting for unary parameters. Instead, a simultaneous optimization of the binary and unary parameters was performed for the Ag-Sn and Cu-Sn systems. For example, the binary potential parameters for the Ag-Sn system and constituent unary parameters for pure Ag were obtained from one optimization process. This process makes the reproducibility of properties of constituent unary systems slightly worse compared to the result of a fitting focusing only on the target unary system. However, since this process provides more flexibility of fitting, the reproducibility of properties of alloy systems can be further improved. Based on this consideration, we can obtain reliable interatomic potentials for relatively demanding binary systems (Ag-Sn and Cu-Sn), which exhibit various stable intermetallic compounds in the phase diagram [34]. For the remaining binary (Ag-Cu) system without any stable intermetallic compounds in the phase diagram [34] and the ternary (Ag-Cu-Sn) system, subsequent force-matching processes were performed separately to obtain optimum potential parameters.
Optimization of the potential parameters for unary Ag and Cu systems was performed targeting various atomic configurations of the fcc, bcc, hcp, and liquid phases. During the fitting, a rescaling of the structural energies was performed to match the experimental cohesive energy of the fcc Ag (2.85 eV [35]) and fcc Cu (3.54 eV [35]). We confirmed that the rescaling only has a marginal effect on the overall performance of fitting because it involves an arbitrary selection of the reference state. After deriving the optimum set of potential parameters, the equilibrium nearest-neighbor distance of the reference structure (fcc) was also adjusted to match the experimental lattice constant.
The binary Ag-Sn phase diagram [34] exhibits a stable intermetallic compound for the Ag 3 Sn composition. The binary Ag-Cu phase diagram [34] exhibits only a eutectic reaction without any stable intermetallic compound. The binary Cu-Sn system exhibits a relatively complex characteristic of the phase diagram through the presence of various intermetallic compounds [34]. We particularly focused on the stable intermetallic compounds (e.g., Cu 20 Sn 6 , Cu 3 Sn, and Cu 6 Sn 5 ) for fitting targets that can be expressed in a stoichiometric composition without significant disordering at the sublattice sites. Cu 3 Sn and Cu 6 Sn 5 compounds are worth noting because these structures were often observed at the interface between the solder alloy and Cu [8,9].
Tables S1, S2, and S3 (Supplementary material) present the finally determined potential parameter sets for the unary (Ag and Cu), binary (Ag-Cu, Ag-Sn, and Cu-Sn), and ternary (Ag-Cu-Sn) systems. After optimization of the potential parameters, we again confirmed that the radial cutoff distance of 5.0 Å was large enough to reproduce various properties of target alloy systems. The simulation results in the following sections were performed using this value.

Fitting result: Statistical quality of energies and forces
In this section, the results of parameter optimization are demonstrated by the statistical correlation between the energies and forces from the DFT database and the present potential. Figures S1a and b-d (Supplementary material) indicate scatter plots for the statistical correlations of the energies and forces, respectively. Table S4 also lists the corresponding root mean square (RMS) errors of the energies and forces. These correlations were obtained using configurations obtained directly from the DFT calculation without additional atomic relaxations of structural parameters and atomic positions. Because a well converged atomic configuration at 0 K usually involves a force of almost zero on each atom, the effective comparison of forces can be obtained by utilizing the configurations at finite temperatures. These configurations were already implemented in the DFT database from the ab initio MD simulations at finite temperatures.
As presented in Fig. S1 and Table S4, the fitting quality of binary and ternary alloy systems is comparable to that of unary systems. The result demonstrates that the current strategy of fitting unary and binary parameters simultaneously is indeed effective for obtaining relevant potentials through the forcematching process. Since we used unary parameters for adjustable variables during the fitting of binary Ag-Sn and Cu-Sn potentials, the reproducibility in properties of unary systems could be worse than the fitting results focusing only on the specific unary system. However, as expected from the statistical correlations, the accuracy of fitting for unary Ag and Cu systems appears to be sufficient to reproduce various physical properties.
In this section, all the configurations examining the statistical correlations were used explicitly in the optimization process, and the comparisons were made without considering possible relaxations of structural parameters and atomic positions. We thus must further extend the evaluation of developed potentials even considering MD simulations at finite temperatures, as will be shown in the subsequent sections.

Performance of the developed unary potentials
In this section, the performance of developed unary potentials is examined by performing molecular statics simulations at 0 K and MD simulations at finite temperatures. A wide variety of physical properties were selected for the examination to clarify the strengths and weaknesses of the developed potentials. All validation simulations were performed using the LAMMPS code [36]. By default, the results in this section were obtained from molecular statics simulations, and cell dimensions and atomic positions were fully relaxed under periodic boundary conditions for all directions. Results based on other conditions, e.g., results obtained from MD simulations, and results obtained using different boundary conditions, will be specifically designated. The number of atoms in a simulation cell was at least 4000, ensuring a marginal size effect on the calculated physical properties. MD simulations were performed using a timestep of 2 fs, and temperature and pressure were controlled by the Nosé-Hoover thermostat and barostat [37,38], respectively.
The physical properties calculated in this section can be divided into two groups. The first group comprises properties calculated at 0 K ("Structural, elastic and defect properties" section) that are closely related to the atomic configurations used in the parameter optimization. The comparison of these properties indicates the accuracy of the fitting. The other group comprises properties at finite temperatures ("Physical properties at finite temperatures" section) that were not used directly in the parameter optimization process. The comparison of these properties indicates the transferability of the developed potential.

Structural, elastic, and defect properties
The performance of the developed unary potentials was first confirmed by comparing it with the physical properties reported in the previous experiments and/or those obtained through the preset DFT calculations. Table 1 lists the validation results on structural, elastic, and defect properties at 0 K. The properties predicted by the optimized potential parameters are compared with experimental and DFT results. The predictions by the new potentials are also compared with those by the previous potentials [27] developed based on same potential model (2NN MEAM) without the force-matching process. Experimental values of the cohesive energy and the lattice constant are exactly reproduced for both Ag and Cu because we employed the rescaling procedure for the cohesive energy and the equilibrium nearest-neighbor distance ("Optimization of potential parameters" section). The reproducibility of elastic properties appears to be somewhat inferior to the previous potential [27] which exhibits an exact match with the reported experimental values. However, the fitting of potential parameters focusing too much on the elastic properties usually leads to worse fits of other important properties. In fact, the reproducibility of elastic properties by the present potentials is satisfactory considering the balance of reproducibility of other properties and possible experimental deviations. The present potentials also well reproduce the order in structural energies of fcc, bcc, and hcp structures ( E fcc < E hcp < E bcc ). For defect properties related to the vacancy configurations, there is a significant discrepancy between the experimental and DFT values, especially for pure Ag. Predictions by the present unary potentials are in better agreement with the experimental values of vacancy formation and activation energies. A deviation between the experiments and DFT calculations is also observed for the surface energy, while only average values for various grains of a polycrystalline solid are available for the experiment. Results from both the present and previous [27] potentials better match the DFT predictions for individual fcc crystalline surfaces. The following quantities are listed: the cohesive energy E c (eV/ atom), the lattice constant a (Å), the bulk modulus B , the elastic constants C 11 , C 12 , and C 44 (10 12 dyne/cm 2 ), the structural energy differences E (eV/atom), the vacancy formation energy E vac f (eV), the vacancy migration energy E vac m (eV), the activation energy of vacancy diffusion Q vac (eV), and the surface energies E surf (erg/cm 2 ) for the orientations indicated by the superscript. The DFT and 2NN MEAM calculations were performed at 0 K, while the experimental data were obtained at finite temperatures. a Ref. [ Overall, the performance of the present unary potentials appears to be satisfactory in reproducing bulk, elastic, and defect properties. However, compared to the performance of previous unary potentials [27] based on the same potential model, the improvement by the force-matching process appears to be insignificant. Refitting of unary potentials is thus not justified at this stage. The feasibility of new optimizations will be presented in the next section, which focuses on the evaluation of properties at finite temperatures.

Physical properties at finite temperatures
Now, we address a group of physical properties that are important for applying the developed potentials for phenomena at finite temperatures. We first examine the performance of developed unary potentials in reproducing phonon properties. Figures 1a and b show phonon dispersion curves of fcc Ag and fcc Cu predicted by the present potentials as compared to the experimental data [39,40], present DFT calculations, and previous potential [27]. The predictions by the present potentials are in contrast with those by the previous potential [27], which significantly underestimate the frequencies at the high symmetry points. For fcc Ag, the present potential reproduces the overall location of phonon branches by slightly overestimating the experimental/ DFT frequencies. For fcc Cu, the overall location of phonon branches is also better reproduced by the developed potential as compared to the previous potential [27].
We further evaluated the performance of the present unary potentials in reproducing physical properties at finite temperatures associated with the phonon calculation. Validation of the thermal expansion coefficient and the heat capacity of fcc Ag and fcc Cu were performed based on the quasiharmonic approximation (QHA). The calculation based on the QHA can also be considered as an extension of the verification of phonon properties. The obtained thermal expansion coefficient and heat capacity values are compared with the reported experimental data [41][42][43][44] and present DFT results in Fig. 1c-f. For the thermal expansion coefficient of fcc Ag, DFT results predicted by the QHA somewhat deviates from the experimental results. The results of the present potential show better agreement with the experimental trend. For the thermal expansion coefficient of fcc Cu, the present potential well reproduces both experimental and DFT trends, while a noticeable overestimation is observed for the prediction of the previous potential [27]. For the heat capacity of fcc Ag and fcc Cu, the present and previous [20] potentials do not exhibit a clear difference in the reproducibility of the experimental and DFT trends.
The performance of developed potentials in reproducing additional thermal properties associated with the phase transformation between the solid and liquid phases is also presented in Table 2. All these results were obtained by MD simulations using the NPT ensemble at zero external pressure. The melting point was calculated based on the interface method which utilizes a simulation cell consisting of liquid and solid phases in contact with each other. This method is suitable to calculate the exact melting point, thus overcoming the overheating of the solid phase and undercooling of the liquid phase that typically appears during MD simulations. The volume change upon melting and the enthalpy of melting were obtained by comparing the atomic volume and enthalpy values of solid and liquid phases at the calculated melting point. As listed in Table 2, the performances of the present and previous [27] potentials differ considerably in reproducing these properties. For example, the present potential predicts only about a 3% overestimation for pure Ag and a -7% underestimation for pure Cu compared to the experimentally reported melting temperatures. For the previous potential [27], the discrepancies are even more pronounced (e.g., 9% and 18% overestimations for pure Ag and Cu, respectively).
Overall, it has been shown that the developed unary potentials successfully reproduce a wide range of fundamental properties of pure Ag and Cu. In particular, the results are satisfactory considering that we performed simultaneous optimization of unary (Ag or Cu) and binary parameters (Ag-Sn or Cu-Sn) to better describe binary systems. The good performance of the developed potentials in reproducing physical properties at finite temperatures makes them applicable to practical examples through relevant MD simulations.

Performance of the developed potential for the alloy systems
In this section, we examine the performance of the developed potential in describing the physical properties of the binary and ternary alloy systems. The settings for molecular statics and MD simulations (e.g., code, timestep, thermostat, barostat, and cell size) described in the previous section ("Performance of the developed unary potentials") were used. By default, results  [27].
The listed quantities correspond to the thermal expansion coefficient ε (10 -6 /K), the heat capacity at constant pressure C P (J/mol K), the melting temperature T m (K), the enthalpy of melting H m (kJ/mol), and the volume change upon melting �V m /V solid (%). a Ref. [ were obtained from molecular statics simulations at 0 K unless they are specifically designated as results from MD simulations. The physical properties calculated in this section can be divided into two groups. The first group comprises properties for intermetallic compounds ("Physical properties of intermetallic compounds" section) that are closely related to the atomic configurations used in the parameter optimization. The comparison of these properties indicates the accuracy of the fitting. The other group comprises properties of solid and liquid solutions ("Physical properties of solid solutions" and "Physical properties of liquid solutions" sections), where most of them were not used directly in the parameter optimization process. The comparison of these properties indicates the transferability of the developed alloy potentials.

Physical properties of intermetallic compounds
As explained above, each binary system exhibits significantly different characteristics in the phase diagram. The binary Ag-Cu phase diagram [34] is relatively simple and shows only a eutectic reaction without any stable intermetallic compounds. The binary Ag-Sn phase diagram [34] shows a stable intermetallic compound phase at the Ag 3 Sn composition, and it also indicates the hcp solid solution (~ 12-23 at.% Sn) in addition to the fcc solid solution on the Ag-rich side. On the Sn-rich side, there is a eutectic reaction (~ 4 at.%) associated with a minimum melting point, and this is one reason why Ag is added to Sn for soldering applications. The binary Cu-Sn system [34] shows a quite complex characteristic of the phase diagram and exhibits several stoichiometric and non-stoichiometric intermetallic compounds. In this system, the present work aims to derive relevant potential parameters that reproduce the properties of stable intermetallic compounds. We specifically focused on the Ag 3 Sn, Cu 3 Sn, and Cu 6 Sn 5 compounds that are often observed at the interface between a solder alloy and Cu through the interdiffusion [7][8][9]. The binary potentials developed in the present work are intended to adequately describe the stability of such intermetallic compounds at finite temperatures. We also focused on the properties related to solid solutions and liquid solutions over a wide composition range.
Tables 3, 4, and 5 list the performances of the developed alloy potentials in reproducing the physical properties of various stable and hypothetical compounds for the Ag-Cu, Ag-Sn, and Cu-Sn systems, respectively. Since the Ag-Cu system does not exhibit any stable intermetallic compound, predictions by the developed potential were thus compared with the present DFT results on several selected hypothetical compounds. As listed in Table 3, the properties of various compounds are reasonably reproduced while the bulk modulus of C1-AgCu 2 and B1-AgCu compounds could not be obtained due to the instability. For the Ag-Sn system, the developed potential also well reproduces the properties of intermetallic compounds as listed in Table 4. In particular, properties of the experimentally reported stable D0 a -Ag 3 Sn compound are reasonably reproduced, showing an error of ~ 2.5% for lattice constants. However, a hypothetical L1 2 -Ag 3 Sn compound is more stably maintained with the same composition, although the energy difference between the stable and hypothetical compounds is very small (~ 0.003 eV/atom). Even in the case of the Cu-Sn system, the properties of various compounds are reasonably reproduced by the present potential as listed in Table 5. In the Cu 6 Sn 5 composition, two phases are presented in the phase diagram [34]: one phase with the space group of C2/c is stable at relatively low temperatures, and the other phase with the space group of B8 1 is stable at relatively high temperatures. In fact, a disordering of Cu and Sn atoms in the Sn-sublattice is expected because the exact stoichiometry of the B8 1 structure is not Cu 6 Sn 5 but equiatomic CuSn. Due to the size limitation of the DFT calculation in realizing atomic disordering, we validated the performance of developed potentials using the equiatomic B8 1 -CuSn compound. The low temperature stable C2/c-Cu 6 Sn 5 compound is of particular importance since it usually appears during the soldering process. In the Cu 3 Sn composition, there are also two compounds in the phase diagram [34]: one phase with a space group of Cmcm is stable at relatively low temperatures and the D0 3 phase with a space group of Fm3m is stable at relatively high temperatures. In fact, the composition range of the D0 3 phase is not confined to the Cu 3 Sn composition. This phase exhibits significant offstoichiometry (~ 16-28 at.% Sn) at high temperatures. The Cmcm-Cu 3 Sn compound is a complex structure with 80 atoms in a unit cell. This structure is a kind of superstructure consisting of multiple unit cells of orthorhombic Cu 3 Sn structure with a space group of Pmmn . The formation of the Cmcm-Cu 3 Sn structure requires multiples of a Pmmn-Cu 3 Sn unit cell with 8 atoms and additional atomic displacements. We also evaluated the phase stability of a hypothetical L1 2 ( Pm3m ) compound in the stoichiometric Cu 3 Sn composition and compared it to that of the stable compounds. The order of formation enthalpy at 0 K predicted by the DFT calculation, i.e., E f (Cmcm-Cu 3 Sn) < E f (Pmmn-Cu 3 Sn) < E f (L1 2 -Cu 3 Sn) < E f (D0 3 -Cu 3 Sn), is correctly reproduced by the developed potential.
As listed in Tables 4 and 5, the experimental reported stable compounds are also stably maintained by the present potential at 0 K, enabling the calculation of physical properties. However, there is a slight reversal of phase stability between the stable and hypothetical compounds (e.g., D0 a -Ag 3 Sn vs. L1 2 -Ag 3 Sn), and thus it is not clear whether such compounds are indeed stably maintained at finite temperatures. Therefore, the examination of the phase stability was extended to finite temperatures to see the possibility of unwanted phase transformations to other structures (i.e., to hypothetical compounds). The stability of intermetallic compounds at finite temperatures was investigated The following quantities are listed: the lattice constants a, b, and c (Å), the enthalpy of formation E f (eV/atom), and the bulk modulus B (10 12 dyne/cm 2 ). a Ref. [ by performing NPT ensemble MD runs. Initial simulation cells with the experimentally reported stable structures (D0 a -Ag 3 Sn, C2/c-Cu 6 Sn 5 , and Cmcm-Cu 3 Sn) were gradually heated from at 0 K to 700 K and cooled again to 0 K with heating and cooling rates of ± 0.5 K/ps. As shown in Fig. S2 (supplementary material), each intermetallic compound retains its original crystal structure during the thermal process, confirming its stability at finite temperatures. For the ternary Ag-Cu-Sn system, there is a lack of detailed experimental information on ternary properties such as the presence or absence of ternary intermetallic compounds. We thus considered hypothetical Cu 2 MnAl-type Heusler compounds (Ag 2 CuSn, AgCu 2 Sn, and AgCuSn 2 ) for the evaluation of ternary properties. Table 6 lists the calculated lattice constant, the enthalpy of formation, and the bulk modulus of those hypothetical compounds. Overall trends in physical properties expected by the present DFT calculation were properly reproduced by the established ternary potential.

Physical properties of solid solutions
In addition to the properties of intermetallic compounds, we examined the performance of the developed potential in reproducing physical properties related to solid solutions. Table S5 lists the examined properties of Ag-rich fcc, Cu-rich fcc, and Sn-rich bct (body-centered tetragonal) binary solid solutions. Predictions by the developed potential and the present DFT calculation were compared focusing on the heat of solution, the vacancy-solute binding energy, the solute-solute binding energy, and the solute migration energy of solid solutions. Here, the binding energy is defined as the energy difference between a situation where two defects or solute atoms are close to each other and a situation where those defects or solute atoms are separated and non-interacting.
As listed in Table S5, the performance in reproducing solution properties is different for each binary system. The properties of fcc solid solutions are generally better reproduced compared to those of Sn-rich bct solid solutions. In the Ag-Cu system, the reproducibility of properties is fairly good for both Ag-rich and Cu-rich fcc solid solutions. In the Ag-Sn and Cu-Sn systems, the properties of Ag-rich and Cu-rich fcc solid solutions are also well reproduced except for the dilute heat of solution. For the dilute heat of Sn in an Ag-rich fcc solution, the present potential predicts a negative vale (− 0.151 eV), while the DFT calculation predicts a slight positive value (0.062 eV). Moreover, there is a rather quantitative discrepancy in the dilute heat of the Cu-rich fcc solid solution. Overall trends predicted by the DFT calculation are well reproduced for the binding energies for solute-vacancy and solute-solute interactions and the solute migration energy in fcc solid solutions. To further support this statement, we further examined the composition dependence of the lattice constant in fcc solid solutions as shown in Fig. 2. Previous experiments [45,46] on Cu-rich fcc solid solutions reported an increase in the lattice constant with increasing amounts of Ag or Sn. The experimental data [46] on the Agrich fcc solid solution exhibit a decrease in the lattice constant with the increase in the amount of Cu. The composition dependence of the lattice constant calculated by the developed potential faithfully follows these experimental trends.
We also examine the performance of developed potential in reproducing the physical properties associated with ternary solid solutions. Table S6 compares the predictions for the binding energies between different kinds of solute atoms. The attractive (+) and repulsive (−) binding tendencies between different solute elements are successfully reproduced for the binding energies defined in Ag-rich and Cu-rich fcc solid solutions. For the binding energies in Sn-rich bct structures, the present potential does not reproduce the trend in DFT data. We therefore stress that the present potential should be carefully utilized in future simulations of Sn-rich bct solid solutions. In contrast, the performance of the present potential is significantly better in reproducing properties of fcc solid solutions on Ag-rich and Cu-rich sides.

Physical properties of liquid solutions
Since the main application of the Ag-Cu-Sn alloy system is in soldering, the reproducibility of liquid phase properties is particularly important. As a transferability test for such a purpose, we first examined the composition dependence of the enthalpy of mixing of the liquid phase. The examination was performed by performing NPT ensemble MD runs on the liquid phase at zero external pressure. An initial liquid structure of a certain composition was prepared by randomly distributing atoms in a simulation cell. After the relaxation of the prepared simulation cell for 60 ps at the designated temperature, the enthalpy of a system was obtained by averaging values for additional simulation steps of 30 ps. The calculated mixing enthalpies of the liquid solution for each binary system are shown in Fig. 3. There were previous experimental works [45][46][47] on the enthalpy of mixing over the full compositional space of all binary alloys. As shown in Fig. 3, the Ag-Cu system exhibits a positive deviation from the ideal mixing, while the Ag-Sn and Cu-Sn systems exhibits a negative deviation. This is consistent with the characteristics of binary phase diagrams, where the Ag-Cu system does not show any intermetallic compound for the entire composition range while the Ag-Sn and Cu-Sn systems show several intermetallic compounds. The present potential satisfactorily reproduces such trends while absolute values for the Ag-Cu system deviate somewhat from the experimental values. In the Ag-Sn and Cu-Sn systems, the experimental trend exhibits the minimum mixing enthalpy at a composition of around 25% Sn. The developed potential successfully reproduces this trend, demonstrating an ability to investigate phenomena associated with the liquid solution over a wide composition range. Finally, we further examined the diffusivity of liquid phase considering a possible application of developed potential for the soldering process. Meanwhile, this example provides additional information about the transferability of development potential because the diffusivity of the liquid phase was not directly considered as a fitting target during the optimization of potential parameters. The diffusivities of different elements in the liquid phase were examined by performing a series of MD simulations. At the atomic scale, the diffusivity can be calculated from the time dependence of the mean square displacement (MSD), which can be readily obtained by the MD simulation. The diffusivity was calculated using the Einstein formula which relates the MSD with the diffusivity ( D): where t is the time, < R 2 (t) > is the MSD of atoms in a simulation cell, and 6 is the dimensional factor. The MSD of atoms in the liquid phase was obtained by performing NPT ensemble MD runs at zero external pressure and designated temperatures. First, an initial simulation cell with 16,224 atoms was maintained at a sufficiently high temperature (2000 K) for enough time (40 ps) to obtain the liquid structure. For the simulation on the diffusivity of alloying elements, we used a random distribution of solutes in a simulation cell with a dilute composition Figure 3: Enthalpy of mixing of the (a) Ag-Cu, (b) Ag-Sn, and (c) Cu-Sn binary liquid phase at 1500 K calculated using the present 2NN MEAM potential in comparison with experimental data (Ref. [45][46][47]) and CALPHAD based modeling (Ref. [63]). Experimental data points represent the results measured at various temperatures.

Figure 4:
Diffusivities of elements in the liquid phase at different temperatures calculated using the present 2NN MEAM potential in comparison with experimental data (Refs. [48,49]). For the calculation of impurity diffusivities (Ag or Cu), the alloying element was randomly distributed in each simulation cell with a concentration of 1 at.%.
(1 at.%). Then, the temperature decreased to a target temperature using the rescaling of average velocity of atoms, and a relaxation at the designated temperature was conducted for 40 ps. A final MD run was performed at the same temperature for a long time (2000 ps) to obtain the time dependence of the MSD. Figure 4 shows calculated diffusivities of atoms in the liquid phase. The results at various temperatures (500-800 K) are presented in the Arrhenius plot and compared with the reported experimental data [48,49]. Despite deviations between different experimental works [48,49], the developed potential reproduces general trends in the activation energy of diffusion (i.e., a slope of the Arrhenius curve fitting). The calculated self-diffusivities of pure Sn-based on the developed potential are closer to the experimental data by Careri et al. [48], although an underestimation at low temperatures is clearly presented. Compared to the experimental data by Ma et al. [49], the developed potential exhibits an overall underestimation of the self-diffusivity of Sn. For the diffusivities of solute elements in Sn-based alloys, the present potential correctly reproduces overall trends in the activation energy and the slightly enhanced diffusivities of solute elements compared to the diffusivity of Sn in the liquid phase, e.g., D Ag > D Sn and D Cu > D Sn .

Conclusion
We provide a robust interatomic potential for the ternary Ag-Cu-Sn system on the basis of the 2NN MEAM formalism. The potential parameters for related alloy systems are decided in relation to the recently improved unary description for pure Sn. We improved the unary descriptions for pure Ag and Cu, and further extended the determined unary potentials to binary and ternary alloy systems. To ensure sufficient performance in realizing atomistic simulations on various applications, the potential parameters were optimized based on the force-matching method using the DFT database of energies and forces. The developed potential exhibits enough accuracy and transferability to the physical properties of pure metals, intermetallic compounds, solid solutions, and liquid solutions. As a possible example, the diffusivities of alloying elements in the liquid solution were examined by performing MD simulations. The proposed interatomic potential can be further utilized in future studies to investigate various phenomena associated with soldering applications at the atomic scale. The outlined strategy for obtaining a robust interatomic potential for the present target alloy system can be applied to other alloy systems in a straightforward manner.

Methods
A DFT database of atomic forces and energies related to various atomic configurations was prepared before the optimization of potential parameters. The established database includes a wide variety of configurations resulting from different atomic situations, e.g., configurations with defects, applied deformation, and at finite temperatures to ensure sufficient transferability of the development potential. Detailed information on the selection of such atomic configurations is presented in a previous study [33].
The DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP) [50][51][52] and the projector augmented wave (PAW) method [53] within the Perdew-Burke-Ernzerhof generalized-gradient approximation (GGA) [54] for the exchange-correlation functional. The recommended pseudopotentials from the VASP library (' Ag, ' 'Cu, ' and 'Sn_d') were used for all calculations. A cutoff energy of 400 eV for the plane-wave basis set and the Methfessel-Paxton smearing method with a width of 0.1 eV were used. A k-point mesh of 21 × 21 × 21 was selected for the primitive unit cell with a single atom, and the similar k-point density was employed for supercells of other unary, binary, and ternary alloys.
The lattice constants and bulk modulus at 0 K were calculated by employing the Birch-Murnaghan equation of state [55,56] fitted to the volume change of a supercell range from 0.99 to 1.01 V 0 , where V 0 is the equilibrium volume. The single-crystal elastic constants were obtained by applying a particular strain range from − 1 to + 1% to an obtained equilibrium supercell. The obtained total energies as a function of strains were then fitted to third-order polynomials to calculate corresponding elastic constants. In the calculations of defect configurations, individual atomic positions were relaxed at a constant volume and cell shape, obtained from corresponding configurations without defects. Ionic relaxations were performed using the conjugate gradient algorithm with convergence criteria for energy and forces set to 10 -6 eV and 10 -2 eV/Å, respectively. For the calculation of the migration energy of an atom to a nearby vacancy, the nudged elastic band (NEB) method with the climbing-image extension [57,58] was employed to obtain a proper saddlepoint configuration. Phonon properties were obtained based on the direct force constant approach [59] as implemented in the 'Phonopy' code [60][61][62]. To ensure the reliability of phonon calculations, a larger supercell (256 atoms) and more precise convergence criteria for energy (10 -8 eV) and forces (10 -4 eV/Å) were used.
For the preparation of atomic configurations at finite temperatures, two-step ab initio MD simulations [50] were conducted. In the first step, the ab initio MD simulations were performed based on relatively less accurate DFT convergence parameters (e.g., k-point density and cutoff energy) to ensure an efficient sampling of the configuration space with a longer simulation time. In this step, the MD runs were conducted for a total of 2000 steps with a timestep of 1.5 fs using a single k-point and the default cutoff energy of the PAW potential. In the second step, uncorrelated atomic snapshots were extracted from the previous step, and they were recalculated using a denser k-point mesh and higher cutoff energy to obtain accurate energies and forces for the force-matching process.