Ab initio study and thermodynamic modeling of the Pd-Si-C system

Computational Materials Science 171 (2020) 109238

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Ab initio study and thermodynamic modeling of the Pd-Si-C system a,⁎

b

Chao Jiang , Isabella J. van Rooyen , Subhashish Meher a b c

T

c

Fuels Modeling and Simulation Department, Idaho National Laboratory, Idaho Falls, ID 83415, USA Fuel Design and Development Department, Idaho National Laboratory, Idaho Falls, ID 83415, USA Materials Science and Engineering Department, Idaho National Laboratory, Idaho Falls, ID 83415, USA

A R T I C LE I N FO

A B S T R A C T

Keywords: Ab initio calculations Evolutionary algorithm CALPHAD TRISO

Understanding phase stability in the Pd-Si-C ternary system is important for predicting the high-temperature behavior of tristructural isotropic (TRISO) nuclear fuels for high-temperature gas-cooled reactors. Here, we employ ab initio evolutionary search to predict the structure and energetics of ordered compounds in Pd-Si-C and its binary sub-systems, Pd-Si, Si-C and Pd-C. In addition to confirming the experimentally-known crystal structures for Pd2Si, PdSi, and SiC, our evolutionary search uncovers a previously unreported orthorhombic Pd5Si structure with Cmcm symmetry that is thermodynamically stable at low temperatures. The vibrational formation entropies of Pd3Si, Pd2Si, and PdSi are calculated using the direct force-constant method within the harmonic approximation. Furthermore, ab initio molecular dynamics simulations have been performed to obtain the atomic structures and energetics of liquid alloys. With the incorporation of the present theoretical results, an ab initio informed thermodynamic description of the Pd-Si-C ternary system has been obtained in this study.

1. Introduction Silicon carbide currently finds application as one of the layers in tristructural isotropic (TRISO) nuclear fuels for high-temperature gascooled reactors by providing high-temperature mechanical strength and, also, a diffusion barrier for fission-product containment [1]. Palladium, a major fission product of uranium and plutonium, has been observed in the TRISO fuel SiC layers, both intergranularly and intragranularly, with movement along grain boundaries alongside various other fission or transmutation products, to form palladium silicides and other complex multi-element compounds [2–4]. The only identified SiC layer failure mechanism in fuel from the U.S. Department of Energy Advanced Gas Reactor program is one where formation of a crack in the TRISO particle inner pyrolytic carbon layer allows fission products (primarily Pd) to concentration at and interact with the SiC layer [5]. As these fission products are examined in post-irradiation examination to gain knowledge on the location, compound compositions, and/or phases, it is unclear what compounds will be present during the hightemperature irradiation cycle to eliminate the effects of segregation during cooling after irradiation. Therefore, knowledge of the phase equilibria in the Pd-Si-C ternary system is useful for understanding the safety and performance of the TRISO fuels for advanced nuclear reactors, particularly at the elevated temperatures (≳1600 °C) characteristic of accidents.

⁎

In this paper, we report a theoretical investigation of the zerotemperature phase stability of ordered phases in the Pd-Si-C ternary system using ab initio evolutionary methodology [6,7]. The thermodynamic properties of the high-temperature liquid phases are further obtained using ab initio molecular dynamics (AIMD) simulations. Using the present ab initio results for both solid and liquid phases, combined with experimental phase equilibrium data for Pd-Si [8–11], Si-C [12] and Pd-C [13–15] systems, an ab initio informed thermodynamic description of the Pd-Si-C system has been constructed in this work using the calculation of phase diagrams (CALPHAD) approach. 2. Computational details We use the ab initio evolutionary methodology [6,7] to identify hitherto unsuspected low-energy ground-state structures in the Pd-Si-C system. Evolutionary structure searches are performed in the variablecomposition mode using the Universal Structure Predictor: Evolutionary Xtallography (USPEX) code [16]. The first generation of candidate structures is created randomly, and the subsequent generations are produced through heredity, permutation, mutation, and atomtransmutation operations. The energetically most-favorable structures of the previous generation always survive to the next generation. Due to limited computing resources, only simulation cells containing less than 16 atoms per cell are employed in our search. The underlying geometric

Corresponding author. E-mail address: [email protected] (C. Jiang).

https://doi.org/10.1016/j.commatsci.2019.109238 Received 18 June 2019; Received in revised form 7 August 2019; Accepted 27 August 2019 0927-0256/ © 2019 Published by Elsevier B.V.

Computational Materials Science 171 (2020) 109238

C. Jiang, et al.

240-atom supercells give the formation energy of Pd5Si with P21 structure to be −0.362 eV/atom and −0.367 eV/atom, respectively, with respect to fcc Pd and diamond-structured Si. The good agreement between the two calculations indicate that even a 120-atom supercell is sufficiently large to give fully converged result. Importantly, while the existence of configurational disorder due to the 40% structural vacancies in the Si sublattice of P21 Pd3Si0.6 may lower its free energy by 0.016 eV/atom at 1000 K, assuming an ideal configurational entropy for mixing of ΔSideal = −kB (0.6ln0.6 + 0.4ln0.4) per Si site, this configurational free energy is still insufficient to reverse the relative stability between Cmcm and P21 structures. Future experimental investigation of the crystal structure of Pd5Si compound would be of great interest.

optimizations and total-energy calculations are performed using the projector augmented wave (PAW) method [17] for describing the electron-ion interactions and generalized gradient approximation of Perdew-Burke-Ernzerhof (GGA-PBE) for exchange-correlation functional [18], as implemented in plane-wave code, the Vienna ab initio Simulation Package (VASP) [19]. The cutoff energy for plane-wave basis sets is set at 400 eV. Dense Monkhorst-Pack k-point meshes are used to ensure high numerical accuracy. By computing the quantummechanical forces and stress tensor, the unit cell volume and shape, as well as the internal atomic positions of all structures, are fully relaxed using a conjugate-gradient scheme. For compounds with an ordered crystal structure, the most important contribution to entropy is lattice vibrations. To obtain the vibrational formation entropies of palladium silicides relative to a mechanical mixture of pure fcc Pd and diamond-structured Si, we have calculated their phonon spectra using the direct force-constant approach [20], as implemented in phonopy code [21]. Large supercells are used in our phonon calculations to allow the inclusion of long-range force constants in the fitting. The magnitude of atomic displacements is chosen to be 0.02 Å. From the fitted force constants, one can construct a dynamical matrix for arbitrary wave vector in the Brillouin zone and obtain the phonon frequencies via a solution of the eigenvalue problem for the dynamical matrix. Integration over the Brillouin zone gives the phonon density of states (DOS), from which the high-temperature limit of the harmonic-vibrational formation entropy of a compound can be calculated as ΔSvib = −kB ∫ ln(ν )Δg (ν ) dν [22,23], where kB is the Boltzmann constant, ν is the phonon frequency, and Δg (ν ) is the phonon DOS difference between a compound and its constituent pure elements. Finally, to obtain the mixing energies of liquid alloys, we have performed AIMD simulations [24] in the canonical (NVT) ensemble with the cell volume adjusted to give zero total pressure at a given temperature. Due to constraint of computing resources, we use periodic 128-atom cubic supercells to model the liquid state. A single Γ point is used to sample the Brillouin zone. The time step is 2 fs. All liquid alloys are first equilibrated at their respective temperature for 10,000 timesteps (20 ps). For temperature control during the equilibration stage, velocity rescaling is performed every 50 time steps. To obtain the thermodynamic properties of liquid phases, the AIMD run is then performed for a total duration of 5000 timesteps (10 ps) using a Nose thermostat for temperature control [25]. The total energy (potential + kinetic) of a liquid can then be calculated from the time average of the AIMD trajectory.

3.1.2. Pd-C In the Pd-C system, our evolutionary search identifies no stable compounds. This finding is in agreement with a recent ab initio study by Seriani et al. [26] reporting positive formation energies for all palladium carbides. Interestingly, in addition to the energetically competitive metastable Pd6C structure with C-2/m symmetry found by Seriani et al. [26], our search also finds a low-energy trigonal Pd7C structure with R-3 symmetry. Here, Pd atoms occupy the Wyckoff 18f (0.0543, 0.8123, 0.1680) and 3b (0, 0, 1/2) positions, with C atoms occupying the 3a (0, 0, 0) position. After full structural relaxations, we obtain the lattice constants a and c to be 7.561 Å and 7.026 Å, respectively. While a Pd6C unit cell (Fig. 2a) contains six Pd layers, a Pd7C unit cell (Fig. 2b) contains three closely packed layers of Pd atoms. The Pd layers in both structures are stacked along the 〈1 1 1〉 direction in a ABCABC sequence to form a fcc lattice. Carbon atoms then occupy the octahedral interstitial sites of the fcc Pd lattice in such a way that each Pd atom in Pd6C is coordinated to exactly one C atom. In Pd7C, only Pd atoms at a Wyckoff 18f position is coordinated to one C atom, while a Pd atom at Wyckoff 3b position has no nearest neighbor C atoms. Using fcc Pd and diamond as reference states, the formation energies of Pd6C and Pd7C compounds are calculated to be only 0.037 eV/atom (0.257 eV per C atom) and 0.029 eV/atom (0.231 eV per C atom), respectively. 3.1.3. Si-C In the Si-C system, our evolutionary search finds cubic 3C-SiC (space group F-43m) (Fig. 3a) and hexagonal 2H-SiC (space group P63mc) (Fig. 3b) to be the only structures exhibiting negative formation energies, with the former being 2.9 meV/atom more stable than the latter. In addition, our search finds a metastable tetragonal Si5C2 structure with I-42 m symmetry (Fig. 3c), whose formation energy is only 0.050 eV/atom with respect to pure Si and C in the diamond structure. In all three structures, a C atom is sp3 hybridized and is tetrahedrally coordinated by four nearest-neighbor Si atoms.

3. Results and discussion 3.1. Ab initio evolutionary search at zero temperature

3.1.4. Pd-Si-C In the Pd-Si-C system, no stable ternary compound has been identified in our evolutionary search. This conclusion is in agreement with the experimental study of the Pd-Si-C system by Bhanumurthy and Schmid-Fetzer using X-ray diffraction and scanning electron microscopy [27].

3.1.1. Pd-Si Experimental phase-diagram data for the Pd-Si binary system have been extensively reviewed by Baxi and Massalski [11]. In this system, there exist five palladium silicides: monoclinic Pd5Si (space group P21), orthorhombic Pd9Si2 (space group Pnma), orthorhombic Pd3Si (space group Pnma), hexagonal Pd2Si (space group P-62 m), and orthorhombic PdSi (space group Pnma). For both Pd2Si and PdSi, experimentallyknown stable structures have been successfully identified by our evolutionary search. Importantly, our search finds a previously unreported orthorhombic Cmcm structure to be the lowest-energy structure for Pd5Si, which is 0.049 eV/atom more stable than the experimentally reported monoclinic P21 structure [10]. In this new Cmcm structure (Fig. 1a), Pd atoms occupy the Wyckoff 8f (0, 0.9463, 0.6248), 8f (0, 0.6534, 0.5476), and 4c (0, 0.7564, 1/4) positions, while Si atoms occupy the 4c (0, 0.4574, 1/4) position. The lattice constants a, b, and c are calculated to be 3.774 Å, 8.835 Å, and 10.850 Å, respectively. We model the P21 structure by randomly removing 40% of all Si atoms in a 120-atom (Fig. 1b) and 240-atom (Fig. 1c) supercell, giving a final Pd3Si0.6 composition. Our ab initio calculations using the 120-atom and

3.2. Phonon calculations Fig. 4 shows the calculated phonon-dispersion relations for Pd3Si, Pd2Si, and PdSi along high-symmetry directions in the Brillouin zone. Any imaginary phonon mode signals instability towards phase transformation into lower-energy structures via spontaneous atomic movements. It is evident that all three palladium silicides are mechanically stable at zero pressure because all phonon frequencies are positive throughout the Brillouin zone. For PdSi, this finding is in agreement with a previous ab initio study of the Pd-Si system by Turchi and Ivashchenko [28]. We further calculate the high-temperature limit of the harmonic-vibrational formation entropies of Pd3Si, Pd2Si, and PdSi 2

Computational Materials Science 171 (2020) 109238

C. Jiang, et al.

Fig. 1. Crystal structure of (a) orthorhombic Pd5Si with Cmcm symmetry, (b) monoclinic Pd3Si0.6 with P21 symmetry modeled using a 120-atom supercell, and (c) monoclinic Pd3Si0.6 modeled by a 240-atom supercell. Green, blue and red spheres represent Pd atoms, Si atoms, and Si vacancies, respectively.

higher first peak of gC–C(r) in Si-C liquids (Fig. 5d). Our analysis thus indicates that the interaction between C and Pd (Si) is not strong compared to those between C atoms themselves. Fig. 6a further shows our AIMD-calculated mixing energies of Pd-Si alloys at various compositions and temperatures to be in good agreement with the experimental data from Castanet et al. [29]. AIMD results for Pd-C and Si-C alloys at 3000 K are shown in Fig. 6b and c, respectively. Due to the extremely high melting point of carbon, we use diamond instead of liquid as the reference state for C when calculating the formation energies of liquid Pd-C and Si-C alloys.

to be 0.125 kB/atom (1.039 J/mol/K), 0.339 kB/atom (2.821 J/mol/K), and 0.244 kB/atom (2.030 J/mol/K), respectively. The positive vibrational-formation entropies can be understood because the phonon DOS of a mechanical mixture of Pd and Si has more weight in the highfrequency region compared with the palladium silicides (Fig. 4d). 3.3. Ab initio molecular dynamics simulations Fig. 5 shows the partial-pair correlation functions (PCFs) of liquid Pd-Si, Pd-C and Si-C alloys at various temperatures. The definition of PCF can be found in Ref. [24]. For liquid Pd-Si alloys, the first peak of gPd-Si(r) is clearly higher than those of both gPd-Pd(r) and gSi-Si(r), which indicates that Pd–Si bonds are energetically favored over Pd–Pd and Si–Si bonds in Pd-Si alloys. Consequently, Pd-Si alloys should be classified as an ordering-type system. With increasing Pd concentration, the height of the first peak of Si-Si partial PCF decreases rapidly. Interestingly, when Pd content is greater than 75%, the first peak even completely disappears, suggesting a high degree of short-range ordering since a Si atom is surrounded by Pd atoms only and sees no other Si atoms in its first coordination shell (Fig. 5c). In contrast to Pd-Si alloys, the C–C partial PCFs of liquid Pd1−xCx and Si1−xCx (x = 1/8 and 1/4) alloys exhibit a prominent first peak despite the low-C contents, suggesting a strong clustering tendency of C atoms. The position of the first peak is 1.35 Å, which is rather similar to the C]C double bond length (1.34 Å). A visual inspection of the atomic structures of liquid Pd-C and Si-C alloys shows that C atoms tend to aggregate to form dimers and chainlike complexes, with such a tendency being clearly stronger in the latter (see insets in Fig. 5d). Such a conclusion is consistent with the

3.4. Thermodynamic modeling The Pd-Si-C ternary system and its binary subsystems (Pd-Si, Pd-C, Si-C) have been critically reviewed and modeled using the CALPHAD approach in various previous studies [30–33]. The present work represents an effort to improve upon these previous thermodynamic assessments by extensively using ab initio data to add constraints during the model parameter-optimization process, which may otherwise lead to many different combinations of parameters that can represent the phase diagrams equally well [34]. In this work, the liquid phase in the Pd-Si-C system is described by an (Pd, Pd2Si, Si, C) association model [35]. The fcc phase is described using a (Pd, Si)1(C, Va)1 two-sublattice model, in which Si and C atoms occupy the substitutional and octahedral interstitial sites of the fcc Pd lattice, respectively. The diamond phase in the Si-C binary system is modeled as a simple random solution. Finally, the homogeneity ranges of Pd5Si, Pd9Si2, Pd3Si, Pd2Si, PdSi, and SiC are neglected in this work and they are simply treated as line

Fig. 2. Crystal structure of (a) monoclinic Pd6C with C2/m symmetry and (b) trigonal Pd7C with R-3 symmetry. Green and red spheres represent Pd and C atoms, respectively. 3

Computational Materials Science 171 (2020) 109238

C. Jiang, et al.

Fig. 3. Crystal structure of (a) 3C-SiC, (b) 2H-SiC, and (c) tetragonal Si5C2 with I-42m symmetry. Blue and red spheres represent Si and C atoms, respectively. In Si5C2, Si atoms occupy Wyckoff 8i (0.7979, 0.7979, 0.3512) and 2a (0, 0, 0) positions, and C atoms occupy the 4d (0, 1/2, 1/4) position. The lattice constants a and c are calculated to be 4.266 Å and 10.897 Å, respectively.

Thermo-Calc [37], which works by minimizing the weighted sum of differences between model-predicted and theoretical/experimental values using nonlinear least-square regression. To greatly reduce the degrees of freedom for the optimization problem, both the a and b parameters for compounds Pd3Si, Pd2Si, and PdSi are fixed at their ab initio calculated formation enthalpy and vibrational formation entropy values, respectively. The binary interaction parameters for fcc Pd-Si 1 1 0 0 (LPd , Si : Va and LPd, Si : Va ) and diamond Si-C (LC , Si and LC , Si ) solution phases, which describe the deviation of a real alloy from the ideal solution behavior, are directly calculated using the Special quasirandom structure (SQS) approach [38,39]. For the two end members of the (Pd,

compounds. Here, the Gibbs free energy of a AxBy compound can be A B generally written as Gm x y = xGAΦA + yGBΦB + a + bT , where the coefficients a and b are the enthalpy and entropy of formation of the compound, respectively. GiΦ denotes the Gibbs energy of pure element i in crystal structure Φ, which can be taken from the Scientific Group Thermodata Europe (SGTE) pure element database [36], with the reference state being the enthalpies of the pure elements in their stable states at 298.15 K, commonly referred to as Standard Element Reference (SER). More details of the thermodynamic models can be found in Table 1. Model optimization is performed using the Parrot module in

Fig. 4. Calculated phonon dispersion relations for (a) orthorhombic Pd3Si, (b) hexagonal Pd2Si, and (c) orthorhombic PdSi along various high-symmetry directions in the Brillouin zone. In (d), the phonon DOS of hexagonal Pd2Si and a mechanical mixture of fcc Pd and diamond-structured Si with the same composition are shown. 4

Computational Materials Science 171 (2020) 109238

C. Jiang, et al.

Fig. 5. Partial-pair correlation functions of liquid Pd1−xSix (x = 0, 1/8, 1/4, 3/8, 1/2, 5/8, 3/4, 7/8, and 1) alloys at 1600 K from the present AIMD simulations. Results for gPd-Pd(r), gPd-Si(r), and gSi-Si(r) are shown in (a), (b), and (c), respectively. In (d), the C–C partial-pair correlation functions of liquid Pd1−xCx and Si1−xCx (x = 1/8 and 1/4) at 3000 K from present AIMD simulations are shown. The atomic structures of liquid Si3C and Pd3C alloy at 3000 K are shown as insets in (d). Green, blue, and red spheres represent Pd, Si, and C atoms, respectively.

describing the Pd-Si-C system are given in Table 1. As shown in Figs. 7–9, our model-predicted Pd-Si, Pd-C, and Si-C phase diagrams are in overall good agreement with the experimental phase-equilibrium data. The AIMD results for liquids are also well reproduced (Fig. 6). Table 2 further reports detailed comparisons between model-calculated and experimentally measured temperatures and compositions for invariant reactions. Fig. 10 shows the calculated isothermal sections of Pd-Si-C ternary system at 800 and 1000 °C compared to the X-ray diffraction experiments by Bhanumurthy and Schmid-Fetzer [27]. As shown, most of the experimental phase-identification results are consistent with our calculated phase diagrams. One large inconsistency exists concerning the phase stability of the PdSi compound. Pd0.45Si0.45C0.1 alloy annealed at 800 °C did not show the presence of PdSi compound. However, PdSi was observed in Pd0.45Si0.45C0.1 annealed at a higher temperature of 850 °C [27]. While this experimental observation agrees with the proposed eutectoid decomposition of PdSi into a mixture of Pd2Si and Si at 824 °C [9], it is in direct contradiction to present and previous [28,32] ab initio calculations showing that PdSi is thermodynamically stable at zero temperature. For ordered compounds such as PdSi and Pd2Si, configurational entropies are negligibly small. We can thus calculate the free-energy change for the eutectoid reaction 2PdSi = Pd2Si + Si at finite temperatures as ΔFtotal (T ) = ΔEstatic + ΔFvib (T ) , where ΔEstatic is static-lattice total energy difference readily accessible from ab initio calculations at 0 K. ΔFvib (T ) denotes non-configurational contribution to the free-energy change due to lattice vibrations, which can be calculated from phonon DOS within the harmonic approximation [22,23]. Using the compound formation energies reported by Turchi and Ivashchenko [28] and Zhou et al. [32], we have calculated ΔEstatic to be 0.026 and

Si)1(C, Va)1 two-sublattice model describing fcc phase, rocksalt PdC and rocksalt SiC, their formation energies are further obtained from present ab initio calculations. Interestingly, our phonon calculations reveal that rocksalt PdC is actually mechanically unstable at zero pressure, which renders its thermodynamic properties undefined. In this work, the formation energy of rocksalt PdC with respect to fcc Pd and diamond is set at a value of 0.518 eV/atom (50,000 J/mol), which is reasonably large to prevent its appearance in phase-diagram calculations. The Gibbs energy of the SiC compound is taken from the study of Grobner et al. [33]. The remaining unknown model parameters, mostly related to the liquid phase, are then simultaneously fitted to both experimental phase-equilibrium data, such as liquidus, solidus, and invariant reaction temperatures, and ab initio energetics for both solid and liquid phases. For the Pd-Si system, the liquidus data from Langer and Wachtel [8,9] significantly differ from those from Wysocki and Duwez [10] for the Liquid/Pd2Si and Liquid/Pd3Si boundaries. Because the results from Langer and Wachtel [8,9] are consistent with the recent differentialscanning colorimetry measurements by Zhou et al. [32], they are adopted in the present study. As a consequence, the Pd3Si compound melts congruently, rather than via a peritectic reaction as was proposed by Baxi and Massalski [11]. For the Pd-C system, the solvus data have been determined by Selman et al. [13] based on detailed lattice-parameter measurements as a function of C content using X-ray diffraction. It is worth noting that the C solubility values from Selman et al. [13] are well below those reported by Siller et al. [15]. The former experimental data are adopted in this work. The binary interaction parameter between C and Pd2Si associate in the liquid phase, LC0, Pd2 Si , is fitted to the AIMD-calculated formation energies for ternary Pd0.875−xSixC0.125 liquid alloys at 3000 K (Fig. 6d). The final optimized model parameters 5

Computational Materials Science 171 (2020) 109238

C. Jiang, et al.

Fig. 6. AIMD calculated energies of (a) Pd1−xSix, (b) Pd1−xCx, (c) Si1−xCx, and (d) Pd0.875−xSixC0.125 liquid alloys at various temperatures. For Pd-Si system, results from experimental measurements by Castanet et al. [29] are also show. The solid lines are predicted by the present thermodynamic modeling at 1800 K for Pd1−xSix and 3000 K for Pd1−xCx, Si1−xCx and Pd0.875−xSixC0.125 alloys. Liquid Pd, liquid Si, and diamond are used as the reference states in calculating the liquid formation energies.

neutron irradiated TRISO particle [4]. Our ab initio calculations using PBE suggest that this structure is 0.197 eV/atom higher in energy than the ground-state orthorhombic Pd3Si structure. Clearly, phase stability in the Pd-Si-C system can be modified by other factors such as irradiation and/or internal stresses in a sample.

0.014 eV/atom, respectively. These values are in reasonable agreement with the present PBE calculations giving ΔEstatic = 0.021 eV/atom. Here, a positive ΔEstatic indicates the thermodynamic stability of PdSi at 0 K, and vice versa. To confirm the prediction of PBE, we have also calculated ΔEstatic using local density approximation (LDA) and new PBE for solids (PBEsol) [40] functionals. Due to its overbinding tendency, LDA predicts more negative formation energies than PBE and a more positive value of 0.061 eV/atom for ΔEstatic . Calculations with PBEsol also yield a positive value of ΔEstatic = 0.056 eV/atom. Thus, all three exchange-correlation functionals predict the T = 0 K stability of PdSi. As shown in Fig. 11, after taking into account the vibrational free-energy contributions, ΔFtotal (T ) remains positive over a wide temperature range, suggesting that the decomposition of PdSi is thermodynamically unfavorable. Another notable inconsistency is that experiments by Bhanumurthy and Schmid-Fetzer [27] could not confirm the existence of Pd5Si compound in Pd0.8Si0.1C0.1 alloy annealed at 800 °C for 500 h. Instead, the Pd9Si2 compound was observed. Presumably, the annealing time may not be sufficient to reach the true equilibrium state in this alloy due to either sluggish kinetics or small thermodynamic driving force for the formation of Pd5Si. Finally, it is worth mentioning that Pd3Si with a cubic L12 structure (space group Pm-3m) has been recently observed in the SiC layer of a

4. Summary In the present study, we have performed an extensive ab initio study of the ternary Pd-Si-C system. While our ab initio evolutionary search confirms many experimentally observed crystal structures in this system, it finds a new orthorhombic Cmcm structure for Pd5Si that is 0.049 eV/atom more stable than the experimentally reported monoclinic P21 structure. In agreement with experimental studies, no stable ternary compound has been identified by our evolutionary search. Our ab initio phonon calculations show that palladium silicides Pd3Si, Pd2Si and PdSi are all mechanically stable at zero pressure since they exhibit no imaginary phonon frequencies (soft modes) throughout the Brillouin zone. AIMD simulations have also been performed to model the liquid state. Analysis of partial PCFs indicates a strong short-range ordering tendency in Pd-Si liquids and a strong clustering tendency of C atoms in Pd-C and Si-C liquids. Combining the calculated results for both solid and liquid phases with experimental phase-equilibrium data, an ab 6

7

=

i = Pd, Si

(∑ +

' 1 LPd , Si : Va (yPd

)

)

' −ySi ))

fraction of i in the first and second sublattice, respectively

+

'' ' 0 yVa yPd ySi' (LPd , Si : Va

+

Line compound: (Pd)2(Si)1

Line compound: (Pd)3(Si)1

Line compound: (Pd)1(Si)1

Line compound: (Pd)9(Si)2

Line compound: (Pd)5(Si)1

Line compound: (Si)1(C)1

Pd3Si

PdSi

Pd9Si2

Pd5Si

SiC

dia Gexcess = x C xSi (LC0, Si + LC1 , Si (x C −xSi))

dia dia Gm = ∑i = C , Si xi Gidia + RT ∑i = C , Si xi lnxi + Gexcess

Random solution model: (C, Si)

fcc Gexcess

∑i = C , Va yi'' lnyi'' )

fcc fcc ∑j = C , Va yi yj Gi :fcc j + Gideal + Gexcess

RT (∑i = Pd, Si yi' lnyi'

'' 2 − yVa

1

' '' 0 = yPd yC'' yVa LPd : C , Va '' ' yi and yi denote the site

fcc Gideal

Gmfcc =

Two-sublattice model: (Pd, Si)1(C, Va)1

Pd2Si

Diamond

Fcc

(

i, j = C, Pd, Pd2Si, Si yi denotes the site fraction of i.

∑i yi Giliq + RT ∑i yi lnyi + ∑i ∑j > i yi yj Li0, j

T

This work

SiC SER SER GSi − HSi = −88583.96 + 271.1462 ∗ T − 41.27945 ∗ T ∗ ln(T) − 0.00436266 ∗ T2 + 800000/T + 2 ∗ 10−7 ∗ T3 : C − HSi

fcc Pd5 Si diamond GPd = −238830.8 − 5.77951 ∗ T : Si − 5GPd − GSi

fcc Pd9 Si2 diamond GPd = −477193.3 − 8.99104 ∗ T : Si − 9GPd − 2GSi

PdSi − G fcc − G diamond = − 103796.1 − 4.06055 ∗ T GPd Si : Si Pd

fcc Pd3 Si diamond GPd = −233496.5 − 4.15741 ∗ T : Si − 3GPd − GSi

fcc Pd2 Si diamond GPd = −199511.7 − 8.46359 ∗ T : Si − 2GPd − GSi

LC1 , Si = 222448.0

Ref. [33]

This work

This work

This work

This work

This work

This work

This work

This work

1 LPd , Si : Va = − 59312.4

LC0, Si = 474869.7

This work

0 LPd , Si : Va = − 152232.8

This work

This work

This work

This work

This work

This work

This work

This work

This work

Source

0 LPd : C , Va = − 22969.9 − 29.09838 ∗ T

diamond GSifcc: C − GSi − GCdiamond = 90202.4

fcc fcc diamond GPd = 100000.0 : C − GPd − GC

0 LPd , Pd2 Si = − 112203.1 + 30.48502 ∗ 0 L Pd = −1479.3 − 48.25860 ∗ T 2 Si, Si 0 LPd , Si = − 121248.5 LC0, Si = −3163.5 + 2.29617 ∗ T LC0, Pd = −58098.1 + 8.87260 ∗ T LC0, Pd2 Si = 70537.9

liq liq liq G Pd − 2GPd − GSi = −217639.9 + 5.49331 ∗ T 2 Si

1 1 + 2y Pd2 Si

Association model: (Pd, Si, Pd2Si, C)

Liquid

liq Gm =

Parameter

Model

Phase

Table 1 Optimized thermodynamic parameters for the Pd-Si-C system (in SI units).

C. Jiang, et al.

Computational Materials Science 171 (2020) 109238

Computational Materials Science 171 (2020) 109238

C. Jiang, et al.

Fig. 7. (a) Model-predicted Pd-Si binary phase diagram in comparison with liquidus and invariant-reaction measurements from Langer and Wachtel [8,9], Wysocki and Duwez [10], and Zhou et al. [32]. (b) Pd-rich part of the phase diagram.

initio-informed thermodynamic description of the Pd-Si-C system has been developed through a highly constrained parameter-optimization process in order to mitigate the uncertainties associated with the large possible solution space of CALPHAD modeling. It is also worth mentioning that ab initio calculations using PBE, PBEsol, and LDA functionals all indicate that orthorhombic PdSi is thermodynamically stable against decomposition into Pd2Si and Si at low temperatures, which directly contradicts the experimentally observed eutectoid decomposition of PdSi. Future studies to resolve this discrepancy between theory and experiment will be of great interest. CRediT authorship contribution statement Chao Jiang: Conceptualization, Investigation, Formal analysis, Methodology, Visualization, Writing - original draft. Isabella J. van Rooyen: Conceptualization, Investigation, Funding acquisition, Project administration, Writing - review & editing. Subhashish Meher: Investigation, Writing - review & editing.

Fig. 8. Model-predicted Pd-C binary phase diagram in comparison with solvus measurements from Selman et al. [13] and Siller et al. [15] and invariant-reaction temperature from Nadler and Kempter [14].

Conflict of interest We have no conflict of interest to declare. Acknowledgements This work was sponsored by the U.S. Department of Energy, Office of Nuclear Energy, under U.S. Department of Energy Idaho Operations Office Contract DE-AC07-05ID14517, as part of the Advanced Reactor Technology Program. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The authors would like to acknowledge the efforts of John Stempien for the review of this paper. All calculations are performed using the Falcon supercomputing system at Idaho National Laboratory.

Fig. 9. Model-predicted Si-C binary phase diagram in comparison with liquidus and invariant reaction measurements from Scace and Slack [12].

Data availability The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study. 8

Computational Materials Science 171 (2020) 109238

C. Jiang, et al.

Table 2 Invariant reactions in the Pd-Si, Pd-C, and Si-C systems. The values shown in parentheses are from experimental measurements. System

Invariant Reaction

Type

Temperature (K)

Liquid composition

Pd-Si

Liquid = Pd2Si Liquid = Pd3Si + Pd2Si Liquid = Pd3Si Liquid + Pd3Si = Pd9Si2 Liquid = Pd5Si Liquid = Fcc + Pd5Si

Congruent Eutectic Congruent Peritectic Congruent Eutectic

1663 1323 1346 1096 1102 1090

– 26.6 at.% Si (26.0 at.% Si [8]) – – – 15.0 at.% Si (15.5 at.% Si [10])

Pd-C

Liquid = Fcc + Graphite

Eutectic

1777 (1777 [14])

–

Si-C

Liquid + Graphite = SiC

Peritectic

3089 (3103 [12])

18.1 at.% C (19.0 at.% C [12])

(1667 (1323 (1343 (1096 (1108 (1083

[8]) [8]) [8]) [10]) [8]) [10])

Fig. 10. Model-predicted Pd-Si-C isothermal sections at (a) 800 and (b) 100 °C, in comparison with experimental measurements from Bhanumurthy and SchmidFetzer [27]. The symbols indicate the experimentally explored compositions. The green lines are tie-lines, and red triangles indicate three-phase equilibrium. 36281, 2017. [4] S. Meher, I.J. van Rooyen, T.M. Lillo, Sci. Rep. 8 (2018) 98. [5] J.D. Hunn, C.A. Baldwin, T.J. Gerczak, F.C. Montgomery, R.N. Morris, C.M. Silva, P.A. Demkowicz, J.M. Harp, S.A. Ploger, Nucl. Eng. Des. 306 (2016) 36. [6] A.R. Oganov, J.H. Chen, C. Gatti, Y.Z. Ma, Y.M. Ma, C.W. Glass, Z.X. Liu, T. Yu, O.O. Kurakevych, V.L. Solozhenko, Nature 457 (2009) 863. [7] Y.M. Ma, M. Eremets, A.R. Oganov, Y. Xie, I. Trojan, S. Medvedev, A.O. Lyakhov, M. Valle, V. Prakapenka, Nature 458 (2009) 182. [8] H. Langer, E. Wachtel, Z. Metallkd. 74 (1983) 535. [9] H. Langer, E. Wachtel, Z. Metallkd. 72 (1981) 769. [10] J.A. Wysocki, P.E. Duwez, Metall. Trans. A 12 (1981) 1455. [11] H.C. Baxi, T.B. Massalski, J. Phase Equilibria 12 (1991) 349. [12] R.I. Scace, G.A. Slack, J. Chem. Phys. 30 (1959) 1551. [13] G.L. Selman, P.J. Elison, A.S. Darling, Platinum Metals Rev. 14 (1970) 14. [14] M.R. Nadler, C.P. Kempter, J. Phys. Chem. 64 (1960) 1948. [15] R.H. Siller, W.A. Oates, R.B. Mclellan, J. Less-Common Metals 16 (1968) 71. [16] A.R. Oganov, C.W. Glass, J. Chem. Phys. 124 (2006) 244704. [17] G. Kresse, D. Jouber, Phys. Rev. B 59 (1999) 1758. [18] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [19] G. Kresse, J. Furthmuller, Phys. Rev. B 54 (1996) 11169. [20] K. Parlinski, Z.Q. Li, Y. Kawazoe, Phys. Rev. Lett. 78 (1997) 4063. [21] L. Chaput, A. Togo, I. Tanaka, G. Hug, Phys. Rev. B 84 (2011) 094302. [22] A. van de Walle, G. Cedar, Rev. Mod. Phys. 74 (2002) 11. [23] C. Jiang, B. Gleeson, Acta Mater. 54 (2006) 4101. [24] X. Hui, H.Z. Zhang, G.L. Chen, S.L. Shang, Y. Wang, J.Y. Qin, Z.K. Liu, Acta Mater. 57 (2009) 376. [25] S. Nose, J. Chem. Phys. 81 (1984) 511. [26] N. Seriani, F. Mittendorfer, G. Kresse, J. Chem. Phys. 132 (2010) 024711. [27] K. Bhanumurthy, R. Schmid-Fetzer, Z. Metallkd. 87 (1996) 244. [28] P.E.A. Turchi, V.I. Ivashchenko, J. Nucl. Mater. 454 (2014) 308. [29] R. Castanet, R. Chastel, C. Bergman, J. Chem. Thermodyn. 15 (1983) 773. [30] Z. Du, C. Guo, X. Yang, T. Liu, Intermetallics 14 (2006) 560. [31] N. Saunders, Calphad 9 (1985) 297. [32] S.H. Zhou, Y. Guo, R.E. Napolitano, Metall. Trans. A 47 (2006) 194. [33] J. Grobner, H.L. Lukas, F. Aldinger, Calphad 20 (1996) 246. [34] S. Bajaj, A. Landa, P. Soderlind, P.E.A. Turchi, R. Arroyave, J. Nucl. Mater. 419 (2011) 177. [35] F. Sommer, Z. Metallkd. 73 (1982) 72. [36] A.T. Dinsdale, Calphad 15 (1991) 317.

Fig. 11. Calculated free-energy change for the 2PdSi = Pd2Si + Si eutectoid reaction. The red line represents results from phonon calculations. The blue line is calculated using the thermodynamic parameters for PdSi and Pd2Si reported in Table 1. The discrepancy between the two curves at 0 K is due to the neglect of zero-point energies in the present thermodynamic modeling.

References [1] J.B. Malherbe, J. Phys. D: Appl. Phys. 46 (2013) 473001. [2] I.J. van Rooyen, Y.Q. Wu, T.M. Lillo, J. Nucl. Mater. 446 (2014) 178. [3] I.J. van Rooyen, T.M. Lillo, H. Wen, K.E. Wright, J. Madden, J. Aguiar, Advanced Electron Microscopy and Micro-Analytical Technique Development and Application on Irradiated TRISO-Coated Particles from the AGR-1 Experiment, INL/EXT-15-

9

Computational Materials Science 171 (2020) 109238

C. Jiang, et al.

[40] J.P. Perdew, A. Ruzsinszky, G.I. Csonka, O.A. Vydrov, G.E. Scuseria, L.A. Constantin, X. Zhou, K. Burke, Phys. Rev. Lett. 100 (2008) 136406.

[37] J.O. Andersson, T. Helander, L. Hoglund, P.F. Shi, B. Sundman, Calphad 26 (2002) 273. [38] A. Zunger, S.H. Wei, L.G. Ferreira, J.E. Bernard, Phys. Rev. Lett. 65 (1990) 353. [39] C. Jiang, C. Wolverton, J. Sofo, L.Q. Chen, Z.K. Liu, Phys. Rev. B 69 (2004) 214202.

10