Abstract

Quantum Sutton-Chen (Q-SC) potentials for molecular dynamic (MD) simulation were derived for the Pd-Si system, which were then used to obtain an atomistic description of melting and transport properties for palladium metal, metallic silicon and their alloys. Melting and structural properties were investigated by analysing the radial distribution function, enthalpy, density, and diffusion coefficient as a function of temperature. The agreement between the MD/Q-SC results and experimental values for the estimated melting points and structural properties was excellent for both pure elements: Pd and metallic Si, while melting of alloys was shown to be complicated by chemical association between the interacting constituents, which led to difficulty in the dissociation of long-range order and thus significant overshoot in calculated melting points owing to high heating rate for achievable MD execution.

1. Introduction

Computational modelling has become extensively used in many areas of condensed matter physics and materials science over recent decades ¹. The first ab initio MD (AIMD) effort was made by Car and Parrinello, who used the so-called plane wave MD method to study condense matter ². Since then, the investigation of liquid metals and their alloys has turned into a foremost area in MD simulation. It is well understood that ab initio MD is exceptionally effective in ascertaining properties of materials with high precision, though it has the shortcoming of being untractable for systems bigger than a few hundred atoms due to heavy computing cost. On the other hand, classical MD are based on analytical interatomic potentials, with the advantage of much lower computing cost for practical simulation of systems of reasonable sizes for the investigations of liquid properties. The interatomic potentials can be either fitted to first-principles calculations or optimised to reproduce physical properties such as lattice constants, cohesive energy, vacancy formation energy and surface energy ^{3, 4}.

MD interest in liquid metals covers various properties such as melting point, statistic structures (i.e. pair correlation function) and transport properties. Until the 19^th century, the metalloid was referred to originally as nonmetals or semi-metals. The term metalloid was initially reported by Pinkerton during 1800, to depict a mineral group of pyroxene with metallic magnificence ⁵. Extension MD simulation into metal-metalloid systems are more challenging due to the complexity of describing the interaction potentials associated with directionality. In addition to successful MD simulation using Q-SC potential of various metals and metallic alloys, e.g. face-centred-cubic (FCC) metals palladium, nickel, copper, silver and their binary alloys ^{6, 7, 8} more recently efforts have also been made to study the physical and structural properties of metalloid elements and their alloys with metals; for example, e.g. Silicon ⁹ and Pd-Si alloys ^{10, 11, 12}. In 1985, Tersoff ¹³ proposed an empirical potential model for structural properties of silicon, which provided an accurate description of bonding and geometry for silicon. Later, the above-mentioned model was upgraded to calculate the aspects related to the structural properties and energy of the entire covalent systems ⁹. Simultaneously, to prepare an amorphous silicon through MD simulation, the new potential was developed by Stillinger and Weber (SW) ¹⁴. To develop the amorphous system, the metal was first melted at high temperature, followed by gradual cooling to low temperature in a short span of time. The outcome had revealed the structural and thermodynamic properties of amorphous silicon, which were in agreement with experimental data. To obtain an in-depth analysis of the disordered silicon phases along with its structural, dynamical and electronic properties, “ab initio” MD study was carried out by Car and Parrinello ¹⁵, wherein amorphous silicon was prepared by rapid cooling from the melt. The results predicted during the simulation studies were in line with the then available experimental data. From the combination of both these simulations, “the tight binding” approach was developed for the simulation related to silicon by fitting the parameters obtained from first principle calculations.The newly developed “tight binding” model successfully reproduced the physical and structural properties of silicon, which was in line with the data obtained from “first calculations” ¹⁶.

In line with the above simulations, Ding and Anderson ¹⁷ reparameterized the SW potential for MD studies of amorphous germanium. They effectively used the above-altered SW potential for MD studies of structural properties of germanium. The outcome of the simulation provided excellent structural properties for both the phases i.e. amorphous and crystalline of germanium.

Recently, metalloid alloys, such as silicon-germanium (SiGe) ¹⁸ were studied using AIMD to investigate structural and electronic properties. The results indicated that the total pair correlation function for the amorphous phase of silicon germanium (a-SiGe) was in line with the experimental data. Furthermore, they demonstrated the metallic character of the liquid phase of silicon germanium (l-SiGe) by evaluating its conductivity and derived semiconductor nature for the amorphous phase (a-SiGe).To extend MD systems for the metal and metalloid alloy such as Pd₈₀Si₂₀ , studies were conducted using the Yukawa pair potentials by Takagi et al. ¹⁰. These revealed that a similar structural factor was found in the experimental data obtained by Takagi et al. and the same number of the nearest neighbours of Pd around Si as the Gaskell proposed a structural model ¹². MD simulation for various compositions of Pd-Si alloys using the Lennard–Jones pair potentials ¹¹ were investigated the local atomic structure based on the analysis of Voronoi polyhedra. The result showed that the medium range order (MRO) increased rapidly, as the composition of Si decreased. Similarly, the various composition of liquid Fe-Si systems using AIMD ¹⁹ were studied to investigate the structural, and bonding properties at 1823K. The calculated “density of states” showed that the Fe-Si was bonded due to hybridization of Fe (3d) and Si (3p). These studies investigated the pair-correlation functions in the liquid and obtained reasonably good agreement with experimental data ²⁰.

Since the 1980’s, great achievements have been made in the formulation interatomic potential with the consideration of some many-body effect with the inclusion of volume or density dependence to describe metallic bonding; such as, for example, the effective medium theory, ²¹ the embedded atom method (EAM), ²² Finnis-Sinclair (FS) ²³ and Sutton-Chen (SC) ²⁴ type potentials. Sutton-Chen potentials have been widely used due to their simple mathematical formulation and fairly long-range characteristics; in particular, they were successfully used in MD studies of thermal and mechanical properties of some transition metals and metal alloys. In addition, Quantum Sutton–Chen potentials (Q-SC) ⁶, which involves parametrization from fitting to first-principles properties, were used successfully in various investigations of metallic systems; glass formation, crystallisation, surface phenomena, clusters, nanowires, single crystal plasticity and transport properties.

Owing to relatively low computing cost, Q-SC potentials are particularly convenient for simulations materials phenomena with a thermal contribution. For instance, they have been used by Ozdemir et al. ²⁵ in MD prediction of the cohesive energy, density, melting point, heat capacity, thermal expansion coefficients, elastic constants, and as well as bulk modulus for the ordered Cu–Pd alloys. Similarly, Kart et al. ²⁶ studied melting, cooling, and crystallisation of binary Pd_xAg_1-x alloys and they obtained a wide range of accurate values for physical properties of the liquid, crystal, and glass phases. The resultant Q-SC potentials provided a satisfactory prediction of temperatures for melting, glass transition, and crystallisation during heating and cooling processes. Various efforts were made to study the metalloid silicon materials, using classical potentials such as the simple Lennard-Jones pair-wise potentials ¹², three-body potentials, ²⁷ Stillinger and Weber, ²⁸ Tersoff, ¹³ and EAM potentials ²⁹.

While many researchers have used the SC potentials to study FCC metals such as Pd, Cu, Ni, Ag, Au, Ir, Pb, Al, Pt, Rh and their alloys, to the best of our knowledge, none has attempted to use (Q-SC) to study metallic alloys containing Si. In this work, we have conducted a comprehensive study of the Pd-Si system; including physical, thermal and mechanical properties. Firstly, we derive the Q-SC many-body potentials, and then apply them to simulate important materials phenomena associated with cooling and heating, in the temperature range from 0 to 2000K.

2. Computational Details

First-principles calculations, within the framework of the density functional theory (DFT), were carried out using the well tested VASP code, ³⁰ with the generalized gradient approximation (GGA) for the exchange-correlation energy. The cohesive energy and lattice parameters of palladium, metallic silicon and their alloys were calculated using a 7×7×7 Monkhorst– Pack grid for the Brillouin zone sampling and a 300 eV energy cut-off. The energy cut-off for the plane-wave basis set and k-point grids was selected to ensure energy convergence to less than 1 meV per atom.

An FCC closed packed silicon structure (Figure 1) was built for fitting the Sutton Chen potential which based on the FCC closed packed metals. In order to find the lowest energy and stable structure of FCC crystal, the volume scan optimisation run and relaxed with the assistance of DFT.

The total energy (E) of the crystal was calculated for a set of chosen volumes (V), and the resulting E vs V data was then fitted to an integrated 3^rd-order Birch–Murnaghan equation-of-state (EOS) ³¹ to calculate the bulk modulus of the Pd-Si system. The GGA functional, used throughout this work, is known to provide an accurate description of structural properties of the closed-packed, metallic and intermetallic systems ³². Electronic properties (transport, magnetism) might require a more accurate treatment of non-local effects (correlation, exchange) in some cases, but are not within the scope of the present work.

The DLPOLY MD package ³³ along with Q-SC potentials was used to gain insights into melting and solidification processes at the atomic scale. Çağin et al. ⁶ were the first to introduce quantum corrections to SC-type potentials by fitting to DFT results, to reproduce experimental properties such as density, cohesive energy, bulk modulus, and phonon frequencies. Q-SC potential parameters for palladium and metallic silicon were obtained by fitting to some experimental data such as lattice parameters, cohesive energy and bulk modulus.

The Sutton-Chen interaction potential comprises a pairwise potential V(r_ij) between atoms i and j, responsible for the Pauli repulsion between the core electrons, and a local energy density ρ_i term accounting for cohesive interactions with atom i. The total potential energy has the following form:

(1)

where,

(2)

(3)

In Eqs. (1) – (3), r_ij is the distance between atoms i and j, a_ij is a length parameter leading to dimensionless V(r_ij) and ρ_i; c_i is a dimensionless parameter scaling the attractive vs repulsive terms, sets the overall energy scale, and n and m are positive integers such that n > m. Interatomic terms are derived from one-centre terms and according to:

(4)

a_ij, m_ij and n_ij interatomic terms are set as the arithmetic mean values of the corresponding one-centre terms.

In order to extract the parameter values, Rafii-Tabar and Sutton (RTS) ³⁴ have developed the so-called random binary fcc metal alloy method, in which the two type of atoms occupy randomly the lattice sites, with the required average concentration. Table 1 lists the Q-SC potential parameter values derived here for Pd–Si alloys, using the RTS approach.

We performed constant-pressure and constant-temperature molecular dynamics simulations (NPT), based on the extended Hamiltonian formalism from the work of Andersen, ³⁵ Parrinello and Rahman, ³⁶ Nosé, ³⁷ and Hoover ³⁸. This scheme combines the Nosé canonical variable shape size and the Parrinello–Rahman variable shape size ensembles. NPT MD allows the study of phase transformations while permitting changes to the cell shape and size. The NPT simulations started from a box with periodic boundary conditions for a system of 1000 to 1372 particles. The temperature was increased initially from 300 K to 2000 K in 100 K increments, at different heating rates ranging from 10 K/ps to 1 K/ps. Smaller temperature increment of 10K was used near the melting temperature to enable more accurate exhibition of system properties associated with the solid-liquid transition. The system was run for 50000 time steps at 2000 K, in order to achieve an equilibrium liquid state as the initial configuration for the cooling process, which was conducted using NPT MD. The simulation time step was chosen as 0.002 ps. The cut off selected for our interaction potential at a range of two lattice parameters where the forces are trivial and the half of a lattice distance was also added to this range to consider the temperature effect. The macroscopic features of the system, such as volume, enthalpy, cohesive energy and transition temperature were obtained from the NPT MD heating. Finally, additional steps of simulation were performed to investigate transport properties such as diffusion.

3. Results and Discussion

The Q-SC potentials for the Pd-Si system derived in this work are verified through MD calculation of structural properties for both elemental phases and binary alloys. The results for lattice constant, cohesive energy and bulk modulus are presented in Table 2 & Table 3, in which they are compared with our DFT data and available theoretical and experimental data from the literature. The calculated lattice parameters and cohesive energies using the Q-SC are in excellent agreement with the experimental data. The cohesive energy for the metastable Si structures i,e. face-centred cubic (FCC) Si, is in good agreement with literature. On the other hand, the cohesive energy for the stable diamond cubic (DC) structure is slightly lower than the reported data; possibly due to the Q-SC potentials ⁶ were based on FCC Si .The MD (Q-SC) calculated bulk modulus is also in good agreement with available experimental data. Interestingly for Pd the MD (Q-SC) modulus even agrees with experimentally determined modulus better than the DFT result. For the metastable FCC Si structure, the current DFT and MD (Q-SC) values are well below previously reported DFT value using the FASTSTRUCTURE code with the Harris functionals for reduced computing cost ⁵. Overall, the Q-SC potentials derived from this work is seen to produce reliable predictions for the structural properties of the Pd-Si system.

Simulating alloys with realistically accurate structural configurations, particularly disordered solution phases without particular preferences for lattice site occupation, is a more difficult task.

A first approach is to generate an FCC based supercell and apply relaxation to find the lowest energy structure. Similarly, the alternative approach is to use the real structure from the available (ICSD database) literature. For both approaches, the atomistic calculations allowed us to determine the equilibrium lattice parameters, cohesive energies and bulk modulus for the Pd-Si metallic alloys by fitting the energy variation with respect to the cell volume to the Murnaghan equation of state, ⁵¹ as shown in Figure 2. As we have seen from the Figure 2, each curve in the cohesive energy exhibits a minimum that corresponds to the equilibrium volume, V0, of the condensed phase ⁴⁰. Atoms repel each other when the volume is smaller than V0, and they attract each other at larger volumes. The measure of the interatomic force corresponds to the tangent of the energy curve. It is seen from Figure 1 that the Q-SC calculated attractive energies are in much better agreement with that from first principles DFT calculation than that from classical EAM potentials. The excellent agreement in the calculated V0 between MD (Q-SC) and DFT is consistent with the agreement in calculated equilibrium lattice parameters (Table 2). On the other hand, Q-SC or EAM potentials overestimate the repulsive energy, making the Pd-Si alloys harder under pressure.

For metallic alloys, such overestimation in repelling energy is not expected to cause a considerable problem for the structure-energy landscape, as such condensed systems are known to be impressible.

For the binary Pd-Si_FCC alloys, the MD(Q-SC) calculated cohesive energies were found to follow a decreasing trend from −4.339 to −4.073 eV/atom with increasing Pd content, in agreement with our first principles calculations that yield −4.297 to −4.118 eV/atom (Table 2). The decreasing cohesive energy with increasing Pd content is also in line with enlargement in lattice parameters, for which there is an excellent agreement between MD(Q-SC) values and DFT results (3.8998 to 3.9153 vs 3.8838 to 3.9082 Å). The MD calculated bulk modulus also follows the change in DFT calculated data (MD 143.91 to 174.63 GPa v.s. DFT 144.44 to 165.74 GPa). Such change in bulk modulus and lattice parameter with alloy composition is consistent with the values for pure FCC Si and Pd, i.e. they vary approximately linearly correlated to the Pd concentration. It is worth noting that the evolution of cohesive energy with alloy composition shows negative deviation from the weighted sum of constituent elements, indicating preferred chemical interaction between Pd and Si. This is supported by the existence of various stoichiometric compounds between the two elements, with preferred Pd-Si interaction leading to ordered occupancies of different sublattice sites.

As indicated in Table 3 the MD (Q-SC) cohesive energies obtained for Pd-Si_ICSD alloys using structures from “the inorganic crystal structure database” followed a similar trend to that of Pd-Si_FCC alloys. One notes that the DFT calculated conhesion energies are in excellent agreement with available experimental data which corresponded to the ICSP structures (ground state stable phases), though the MD(Q-SC) results are seen to underestimate the cohesion evidently. In addition, the discrepancies between the DFT and MD results for both cohesion and bulk moduli are the smallest for the Pd₅Si compound, which corresponds the highest Pd content. Looking back at data in Table 2, one notices excellent agreement between the MD and DFT results. This offers reasonably good indication that while the Q-SC model can provide fairly accurate description for metallic faces based on closepact structures (e.g. FCC), it does fall short in precise picture the energy landscape for compound phases associated with significant directionality in bonding, which is typical in the ground state stable structures of metal-Si alloys ⁴⁸.

Melting temperatures of Pd, Si and their alloys were determined typically from the discontinuity in density and enthalpy. Such significant change in structural configuration from a crystalline solid to disordered liquid can also be demonstrated in other properties such as the pair distribution function, mean square displacements and associated diffusion coefficient. Figure 3 shows the variation in density and enthalpy during heating materials at 5K/ps, using the Q-SC for MD simulation. Both density and enthalpy experienced a significant change in the slope at an elevated temperature for each alloy. Take the Pd₄Si alloy, for example, a corresponding kink appears in both the enthalpy and density curves at the same temperature. As melting is a first order phase transition, significant atomic activation is involved in the dissociation of crystallinity and long-range ordering, and thus fast heating will lead to delayed crystal-liquid transition, leading to higher melting temperature. It is assumed that a homogeneous system without a free surface (WFS) is one of the reasons for the higher melting temperature ⁷. Moreover, as a result of the rapid heating rate, the system may not have had sufficient time to generate distribution of defects in equilibrium; thereby leading to a slightly higher melting temperature. The stronger binding or cohesion in a strongly bond condensed matter such as metals therefore makes melting a rather sluggish process. This makes it rather difficult to determine equilibrium melting temperatures of metallic alloys using feasible heating rate within the framework of MD simulation, due to the very short simulation time achievable.

Recently, Z. Du et al. ⁵⁵ were accessed the Pd-Si system and were not reported in regards with Pd₄Si phase; conversely, it was mentioned earlier by Baxi and Massalski ⁵². Besides, most of the accessible binary phases were studied excluding a stable Pd₂Si phase in the Pd-Si phase diagram as shown Figure 4.

The melting temperatures derived from Figure 3 during heating at 5K/ps (data listed in Table 4), is superimposed on the experimentally determined melting temperature in Figure 4. In spite of some significant overshoot in the melting temperatures of alloy phases, the variation in calculated melting temperatures with compositional change is strikingly consistent with the shape of the phase diagram. In particular, minima in melting temperature are shown to correspond to valleys corresponding to eutectic compositions. It is shocking to see that the calculated melting temperatures of elemental materials Pd (1810±10 K) and FCC Si (1680±10 K) are shown to be in excellent agreement with the phase diagram values, given the heating rate in the MD simulation is over a billion times of that used in metallurgical practice. This indicates that for the pure element, due to the spontaneous nature in self-wetting, the kinetic barrier for solid-liquid transition is rather trivial ⁵⁶. At the melting temperature, we can see that the densities of Pd and FCC Si are 10.27 and 2.56 g/cm³, respectively, which are consistent with experimental values (10.49 and 2.65 (DC Si) g/cm³).

On the other hand, self-wetting is less likely in alloys due to additional freedom in free energy associated with compositional changes associated with the configurational energy landscape in the alloy system. Complexity associated with dissociation of chemical ordering in an alloy solid makes melting of an alloy kinetically more sluggish than in pure elements, and thus the high heating rate used in the MD simulation leads to a significant overshoot in the calculated melting temperature. The lower melting temperatures for the ground state crystal structure correspond to the underestimated cohesion for these phases.

For the Pd-Si system, we were further approached to estimate the melting temperature for selected real alloys with free surface using the slab method ⁵⁷ via MD (Q-SC). In this method, the supercell was constructed with exposing (100) surface by an orthorhombic cell of dimension [3a, b, c] as shown in Figure 5. A real crystal is positioned at the centre of supercell in order to remain equal mass distribution. As a result, the central one-third of volume of the supercell was filled with the real crystal and the volume of remaining two-third unoccupied. The system undergoes a first–order transition at the melting point, when a solid system is heated to the liquid state. The melting point of the system was calculated by identifying the abrupt changes of the cohesive energy and bond orientational order ⁵⁸ while heating Pd₃Si alloy as shown in Figure 6. In this calculation, we used one of the parameter of bond order called Q₆. The calculated melting (solidus) temperature using the slab method, for the PdSi, Pd₃Si, Pd₄Si and Pd₅Si alloys are found to be 1035 K, 1275 K, 1280 K and 1280 K, respectively. This method provides better result for Pd₃Si and Pd₄Si in comparison with WFS method; however, PdSi and Pd₅Si display enhanced result for WFS method. Consequently, the error percentage for the calculated melting using Q-SC potential is slightly higher compared to the experimental melting point for Pd-Si system. Although, the Q-SC potential does not agrees fully to the Pd-Si system, conversely, this potential still provides the useful description for the liquid Pd-Si system.

The pair correlation function g(r) is defined as:

(5)

Where g(r) is the probability of finding an atom at a distance r from the reference atom located at the origin, considering an average density . This function can be extracted from X-ray and neutron diffraction experiments ⁵⁹. g(r) can also be computed from atomic trajectories in MD. As shown in Figure 7, it shows sharp peaks for crystalline phases (black line) and broad hillocks in liquids (red line); secondary hillock splitting characteristic of an amorphous structure (blue line).

Figure 8 shows the total and partial pair correlation functions g(r) for Pd₃Si_FCC at high temperatures. The first peak of g(r) is found around 2.70 Å as shown in Figure 8(a). The peak position is not dependent on temperature. However, the peak heights decrease with increasing temperature. The total pair correlation is largely typical of liquid at a temperature of 1600K, though there are some remnant crystalline characteristics in the material at 1500K, as is shown by the existence of minor hillocks around 4 and 6.5 Å. The evolution of structural changes can also be explicitly depicted through trajectory records of structural configurations, such as these shown in Figure 9 in snapshots of simulated PdSi structures: regular configurational pattern for a crystalline structure (Figure 9a) and liquid disorder structure with clustering tendency of either Pd or Si (Figure 9b).

The partial pair correlation functions g_Pd-Si(r), g_Pd-Pd(r) and g_Si-Si(r) can also be estimated, provided that the density in Eq. (5) is replaced by the partial density where is the density of the system and c_i and c_j are concentrations of elements i and j in the alloy. g_Pd-Pd(r) values for the Pd₃Si alloys at various temperatures are plotted in Figure 8(b). The first peak in g_Pd-Pd(r) is found around 2.73 Å and exhibits very little shift with the temperature. This peak position is very close to the experimental value of 2.71 Å in pure liquid Pd; the height of the principal peak increases with decreasing temperature. g_Si-Si(r) and g_Pd-Si(r) values at different temperatures are shown in Figure 8(c) and Figure 8(d). The statistics for g_Si-Si(r) are not regular but the intensity of the first peak can be clearly seen increasing as the temperature decreases. The shapes of g_Pd-Si(r) are similar to that of total g(r) at the corresponding temperature. The position of the first peak in g_Pd-Si(r) is 2.73 Å and there is a little shift with temperature. The minor peak at 6.5Å for 1500K in (a) is characteristic of Pd-Pd distances (Figure 8.b), which disappears at a higher temperature. It is an evidence that Pd atom moves slightly longer distances to form the second closest neighbours ⁶⁰.

On the other hand, for the Si-Si distances (Figure 8.c), the large spectral peak at 5.5Å likewise vanishes at high temperature because of change in atomic structure of silicon in the liquid phase. At the melting point, the bonding turns out to be progressively more metallic while holding the trivial level of covalence ⁶⁰. In any case, the Pd-Si (Figure 8.d) peak at 5Å drastically weakens and broadens at a high temperature above 1600K. Overall partial pair correlations offer further insight that the melting temperature of the Pd₃Si_FCC alloy is above 1700K, when the atomic association between the Si-Si and Pd-Si nearly disappear. Figure 10 shows the structural differences of various Pd-Si alloys at melting temperature, from g(r) values in the liquid state. The change of peak height and position shows the structural differences within alloys at their melting temperature. The shifting of peaks with temperature in pair correlation functions depends on the chemical associations of Pd-Pd, Si-Si and Pd-Si bonds. The first peak heights increase and positions shift to the right with increasing palladium content in the Pd-Si alloys. This is consistent with the decreased cohesion and enlarged equilibrium volume in the Pd-Si alloys (Table 3).

In MD simulations, self-diffusion coefficients can be determined from the temporal evolution of means square displacement (MSD) by the Einstein relation:

(6)

where is the mean square displacement (MSD), D is the diffusion coefficient, and C is a constant. MSD as a function of time, calculated from Q-SC potentials for Pd, metallic Si and their alloys are plotted in Figure 11(a-c). Changes in MSD slopes allow obtaining the melting points as well.

According to Figure 11 (a, b & c), MSD curves show a linear dependence on time for the elemental matter Si ,Pd, their alloy PdSi also exhibits a linear regime with a change of slope at 0.2ps. The radical difference in the slopes for Pd can be attributed to phase equilibrium involving the solid phase at or below the melting point, when stronger atomic association due to lowered thermal excitation leads to remarkably lowered diffusion coefficient. This can be envisaged with the so-called cage effect ⁶¹ created by the interactive nearest neighbour atoms. The microscopic reason for the low diffusivity at lower temperatures is that the associated atoms are trapped in the cage and it takes a long time to escape from it when thermal excitation is inadequate. The cage becomes more and more rigid with decreasing temperature, thus the time needed to escape increases accordingly.

It is interesting to monitor the variation of the diffusive behaviours in the Pd-Si system during heating. Table 5 lists the diffusion coefficients of metallic Pd-Si alloys, calculated near the melting region using Q-SC potentials and Eq. 6. Diffusion coefficients decrease with increasing Pd content (Figure 12), with a drop of more than one order of magnitude from PdSi to Pd₅Si away from melting point.

As is shown in Table 5 diffusion coefficient of Si is evidently larger than that of Pd., though the cohesive energy of Si is larger than that of Pd (Table 2). This may be attributed to the evidently lower bulk modulus for Si than for that for Pd, so that it is easier to break the less rigid atomic association cages of the former than that for the latter. The diffusion coefficient of Pd is found to be 4.58 (nm2/ns) in agreement with the result of Kart et. al ¹⁵ (4.98 (nm2/ns)) who carried out simulation using the HPN (constant enthalpy and pressure) ensemble in a cubic box of 864 atoms at 1853 K from Green-Kubo relations. Due to lack of diffusion data for Pd-Si alloys in the literature, hence it was not possible to compare our results with literature.

In order to calculate the activation energy for Pd, we have calculated diffusion coefficient from the Einstein equation at various temperatures, as shown in Figure 13, and fitted the data to an Arrhenius equation:

(6)

where D_o is the self-diffusion coefficient prefactor and E_a is the diffusional activation energy.

The estimated values for D_o and E_a are given in Table 6. The activation energy increases with increasing amount of palladium in the alloys, which is consistent with the decreasing diffusion coefficient (Figure 13). Our calculated D_o and E_a values are also in agreement with the available experimental data for Pd. On the other hand, D_o and E_a values for FCC silicon are significantly different from the literature data based on diamond Si. Such discrepancy may be attributed the remarkable difference in the nature of bonding between diamond and FCC Si phases, with the former being dictated by covalent bonding but the latter being largely of metallic nature. Indeed, the Q-SC potentials only account for some volumetric contribution for many body interactions, which is inadequate to account for directionality in covalently bonded structures ³³.

4. Conclusion

In this work, the Quantum Sutton-Chen (Q-SC) potentials for molecular dynamic (MD) simulation were derived for the Pd-Si system for the first time. They were revealed to enable the reliable atomistic description of the metallic or intermetallic alloys in the binary system; including solid-liquid transition and associated structural and transport properties.

The kinetic barrier for melting a Pd-Si alloy is significantly larger than that for melting a pure element, so that significant overestimation of melting temperatures for Pd-Si alloys cannot be avoided owing to MD heating rates being billions of orders higher than practical heating rates in metallurgical practice.

Diffusional process in Pd-Si alloys is dictated by the slower diffusion of Pd atoms, so that higher Si concentration leads to larger inter-diffusion coefficients.

While Q-SC potentials are suitable for modelling systems largely of metallic characteristics in bonding, such FCC based models are not considered adequate for materials phenomena involving significant directionality in bonding due to valence orbital hybridisation. Thus, it is not shocking that the Q-SC potential for Si was inaccurate in describing the diamond cubic bonding in Si. The reasonably good description of Pd-Si alloys using Q-SC is thus, attributed to their dominant metallic nature, with silicide compounds in this system being classified as intermetallic compounds.

References