# First Principles Modeling of the Temperature Dependent Ternary Phase Diagram for the Cu-Pd-S System

###### Abstract

As an aid to the development of hydrogen separation membranes, we predict the temperature dependent phase diagrams using first principles calculations combined with thermodynamic principles. Our method models the phase diagram without empirical fitting parameters. By applying thermodynamic principles and solid solution models, temperature-dependent features of the Cu-Pd-S system can be explained, specifically solubility ranges for substitutions in select crystalline phases. Electronic densities of states calculations explain the relative favorability of certain chemical substitutions. In addition, we calculate sulfidization thresholds for the Pd-S system and activities for the Cu-Pd binary in temperature regimes where the phase diagram contains multiple solid phases.

## 1 Introduction

Sulfidization impedes the use of palladium-based membranes for hydrogen sequester and capture. Due to the presence of hydrogen sulfide (HS) in any practical feed gas, sulfur rapidly covers the surface of the membrane, rendering its catalytic properties inert. Pure palladium membranes have limited application for hydrogen sequester and capture due to irreversible formation of bulk sulfides Chen and Ma (2010) and embrittlement due to lattice distortion caused by hydrogen atoms Amandusson et al. (2001) diffusing through the bulk. Alloying palladium with other metals, such as silver Uemiya et al. (1991) or gold Chen and Ma (2010), fixes many of the problems associated with palladium membranes, but these materials still suffer from sulfidization and are expensive. In the 2000s it was discovered that alloying palladium with copper Iyoha et al. (2007); O’Brien et al. (2010) dramatically reduces the rate of sulfidization and further decreases the price of these membranes, though still not eliminating sulfidization. In order to realize hydrogen capture and sequester on an industrial scale, new membranes with different chemical compositions that resist sulfidization need to be devised, and their reactivity to sulfur needs to be examined. Computational methods offer an attractive way to rapidly prototype different chemical families.

In this paper we computationally model the temperature-dependent phase diagram for the Cu-Pd-S system using first principles methods. Prior experimental work Karup-Møller and Makovicky (1999) determined the ternary phase diagram for the Cu-Pd-S system. We seek a simple and semi-quantitative model to gain insight into the essential features of the system, in contrast to complicated models with many fitting parameters that can yield close agreement with experiment at the cost of obscuring the underlying physics.

First, we describe our method for calculating ternary phase diagrams and apply it to Cu-Pd-S. Our model begins with first principles total energies, moves to focus on substitutional entropy, then examines vibrations. Comparison with experimental results shows general agreement. Certain details of energetics are explored through examination of the electronic density of states. We then apply our model to predict sulfidization thresholds for Pd interacting with S dimers and activities in Cu-Pd binaries.

## 2 Methods and Calculations

In order to calculate a phase diagram using first principles methods, all plausible crystal structures must be considered, in principle. In the case of Cu-Pd-S, a wealth of experimental crystallographic data gives an excellent starting point, especially for the binary compounds. We generalize these binaries to ternary structures formed from substitutions of the ternary species into either lattice sites or interstitial sites. Once the T=0K total energies for all likely crystalline structures have been calculated, it is straightforward to generate a T=0K phase diagram by creating the convex hull of total energies.

Subramanian and Laughlin Subramanian and Laughlin (1991) determined the phase diagram for Cu-Pd experimentally, and Li et al. created a CALPHAD model Li et al. (2008). Above T=871K, the entire composition range of solid Cu-Pd exists in fcc solid solution. Below this temperature other solid phases are present. -CuPd has a CsCl structure and exhibits a Cu-rich composition range extending to a composition of 58% Cu. Trimarchi and Zunger Trimarchi and Zunger (2008) showed using an evolutionary algorithm that first principles calculations predict cP2 to be the expected ground state at equiatomic composition. At copper concentrations near 75%, phases based on 1D and 2D long period superstructures (LPS) compete. Studies using ab initio methods to examine this region of the phase diagram have already been performed Baärthlein et al. (2009); Ceder et al. (1990). As we are primarily concerned with the palladium-rich side of the phase diagram, we only calculate the ideal 1D LPS, CuPd.tP28, (using a notation of [chemical formula].Pearson symbol) and ignore the 2D LPS and phase solubility ranges in the 1D LPS.

Determination of the phase diagram for Cu-S Lee et al. (2007) using first principles calculations is problematic, as it contains an abundance of phases, many of which are non-stoichiometric or have large unit cells. We include only Cu-S structures with well-defined crystallographic refinements, as we do not expect this discrepancy to affect our results on the palladium-rich side of the phase diagram. The Pd-S phase diagram Okamoto (1992) is relatively simple, containing PdS, PdS, PdS, and PdS phases as line compounds at low temperature, with PdS stabilized at T=829K. We summarize the phase diagrams here only for comparison; our method does not require pre-existing knowledge of the phase diagrams.

In order to add temperature dependence to the phase diagrams, we introduce free energy models for select phases of interest, all of which are based on substitutional solid solution models with sublattice filling (where for some structures, the ”sublattice” is the entire lattice). When calculating the ternary phase diagram and the activities, we work in the Gibbs ensemble, where we specify temperature, pressure, and chemical concentration. We neglect pressure dependence in the solid phases’ free energies, as solid compressibilities are low. Phonon free energies are included, in the harmonic approximation, giving free energies of form

(1) |

Electronic free energy differences between competing phases of interest (PdS versus PdS and well-ordered fcc-CuPd versus -CuPd at 50% Cu), were found to be less than 1.5 meV/atom at temperatures as high as 1000 K. These differences are within the margin of error expected for our fits for analytic free energy models, and to simplify our modeling we have chosen to neglect .

For lattice models, two types of defects are common, substitutions and interstitials. Due to atomic size differences, the only interstitials expected to be relevant are sulfur interstitials in Cu-Pd phases. We calculate interstitial defect costs using the well-ordered structures Cu.cF4, Pd.cF4, CuPd.cF4, and CuPd.cP2, with sulfur occupying a single interstitial site for various supercell sizes. The least unfavorable substitution, octahedral in CuPd.cF4, lies approximately 2 eV above the convex hull with significant distortion, and thus we expect contributions to sulfur solubility by occupation of interstitials to be negligible.

We here summarize the configurational free energy G of a simple substitutional solid solution model, valid in the dilute limit. In a substitutional solid solution, we have a set of site classes i, each with N sites, an associated enthalpic cost of performing a chemical substitution of species into site class i, E, and a substitution concentration of x. In the substitutional solid solution model, it is assumed that each substitution is uncorrelated with all other chemical substitutions, and all spatial configurations have identical energies. This implies a linear dependence of with respect to composition, giving an enthalpic contribution of

(2) |

where denotes the enthalpy of the base structure with no substitutions. Including an ideal entropy of mixing

(3) |

leads to a Gibbs free energy . Here, H and E can be determined from first-principles calculations. However, we control only the total impurity concentration for a species , x, not the number of substitutions in a given site class x. Since we work with a fixed number of sites, we fix the total concentration of substitutions, and thus minimize the free energy over all x subject to the constraint of species number conservation.

To better elucidate the model, we restrict to the special case where only one type of substitution (with an enthalpic cost of ) is allowed in only one sublattice containing sites. A case with two sublattices will be considered later. Other special cases include the binary regular solution model (, where is the total number of atoms) and the well-known binary ideal solution model (, ). The configurational free energy takes the form

(4) |

where is the number of substitutions performed and is the concentration of substitutions in the sublattice. The appropriate quantity to consider for thermodynamics is the free energy per atom, , expressed in terms of the , which can be shown to be

(5) |

where , and is a measure of the number density of sublattice sites. Of particular importance is the entropic contribution to the free energy,

(6) |

which we use for all phases modeled throughout this paper.

In order for a given substitution type to appreciably affect the free energy, the two necessary features are large values, corresponding to a sublattice dense in the total lattice and thus a significant entropic contribution from the sublattice, and small positive (or any negative) enthalpic cost for the substitutions, such that at a given temperature the entropic portion of the free energy can overcome the enthalpic cost (or reinforce the enthalpic gain).

We use VASP (the Vienna Ab-Initio Simulation Package) Kresse and Hafner (1993); Kresse and Furthmüller (1996) a plane wave ab-initio package implementing PAW potentials Blöchl (1994) to determine total energies. Previous work by Hu et al. Hu et al. (2010) compared the accuracy of different pseudopotentials in the Pd-S system, and recommended PBEsol Perdew et al. (2008) as most suitable for calculations involving sulfur. The total energy calculations were performed by fully relaxing atomic positions and lattice parameters until energies converged to within 0.1 meV/atom. Subsequent convergence in k-point density was performed. A common energy cutoff of 273 eV was used. In order to model solubility ranges, calculations were performed both in unit cells and larger supercells. All phonon calculations were done using density functional perturbation theory, and were performed in at least a 2x2x2 supercell of the unit cell, with the exception of PdS.cI46, where due to its large size only a 2x2x2 supercell of the primitive cell was used. Though a supercell method with chemical ordering is being used to model phonon free energies for certain disordered phases, Wu et al Wu et al. (2003) found that the high-temperature vibrational entropy difference between L1 and SQS8 configurations is only 0.03, within their margin of error and corresponding to an upper bound on the free energy difference of 2.6 meV/atom at 1000K. Table 1 shows a list of structures and their total energies.

Phase | Pearson Symbol | Prototype | Group # | Space Group | H | dE |

Cu Owen and Yates (1933) | cF4 | Cu | 225 | Fmm | 0 | 0 |

Pd Harris et al. (1972) | cF4 | Cu | 225 | Fmm | 0 | 0 |

S Dimer | N/A | N/A | N/A | N/A | 0 | 0 |

CuPd Okamura (1970) | tP28 | CuPd | 99 | P4mm | -103.0 | 0 |

CuPd | cP2 | ClCs | 221 | Pmm | -116.9 | 0 |

CuS Jr. (1979); Leon et al. (1981) | mP144 | Cu2S | 14 | P2/c | -496.4 | 0 |

CuS Fjellvåg et al. (1988) | oC24 | CuS | 63 | Cmcm | -715.4 | 0 |

CuS King and Prewitt (1979) | cP12 | FeS | 205 | Pa | -821.5 | 0 |

PdS Grønvold and Røst (1962) | tP10 | PdSe | 114 | P2C | -418.9 | 0 |

PdS Matkovič et al. (1976) | cI46 | PdS | 217 | I3m | -589.7 | 0 |

PdS Gaskell (1937) | tP16 | PdS | 84 | P4/m | -892.8 | 0 |

PdS Grønvold and Røst (1957) | oP12 | PdSe | 61 | Pbca | -952.4 | 0 |

CuS Fjellvåg et al. (1988) | hP12 | CuS | 120 | Ic2 | -712.9 | 2.5 |

CuS Janosi (1964) | tP12 | CuS | 96 | P422 | -479.3 | 3.3 |

CuPd Guymont and Gratias (1976) | cP4 | AuCu | 221 | Pmm | -99.7 | 3.4 |

PdS Røst (1968) | oC16 | PdS | 40 | Ama2 | -492.8 | 8 |

CuPd Soutter et al. (1971) | cF4 | Cu | 225 | Fmm | -106.9 | 10 |

CuPd Geisler and Newkirk (1954) | tP20 | Cu4Pd | 84 | P4/m | -69.6 | 12.8 |

CuPdS | mP18 | CuPdSe | 14 | P2/c | -743.0 | 13.2 |

CuS Yamamoto and Kashida (1991) | cF12 | CaF | 225 | Fmm | -503.1 | 20.6 |

Cu | cI2 | W | 229 | Imm | 36.3 | 36.3 |

## 3 Results and Discussion

### 3.1 Phase Diagrams

Given total energies, we calculate the enthalpy of formation for our structures at T=0K,

(7) |

where h is the calculated T=0K total energy for the base structure, is the composition of a given species in the structure, and is the T=0K total energy for pure species in a chosen reference structure. By convention, this reference structure is chosen to be the equilibrium structure at T=0K, so that at T=0K all enthalpies of formation for the pure species vanish. In the Cu-Pd-S system, the standard reference structures would be Cu.cF4 for copper, Pd.cF4 for palladium, and S.mP48 for sulfur.

As we are interested in interaction of gaseous sulfur with copper-palladium membranes, we have chosen to use gaseous sulfur dimers rather than the standard S.mP48 when calculating phase diagrams. In the case of the ternary phase diagram, all structures of interest have relatively low sulfur composition, so the mismatch in reference phase compared to experimental phase diagrams is negligible.

The quantity of importance for constructing temperature dependent phase diagrams is

(8) |

where [h] is the enthalpy of formation of the convex hull at a structure’s stoichiometry. By definition, this quantity is non-negative, and it vanishes only when the structure lies on the convex hull. This quantity is a measure of how far above the convex hull a particular structure lies at T=0K. Structures with large dE’s require large entropic contributions to free energy in order to overcome enthalpic costs.

Figure 1 shows the T=0K binary enthalpies of formation. Observe that palladium and sulfur react strongly, as can be seen by the maximum enthalpy of formation of 0.95 eV/atom. This strong enthalpy of formation reflects the tendency of pure palladium membranes to sulfidize. Next, copper and sulfur react almost as strongly, with a maximum enthalpy of formation of 0.82 eV/atom. Finally, copper and palladium have a maximum enthalpy of formation of only 0.12 eV/atom, indicating that the interspecies binding between copper and palladium is relatively weak. This is expected from the Cu-Pd phase diagram, as all known phases are variants of fcc (the T=0, p=0 structure for both Cu and Pd) or bcc (which is reached from fcc by a martensic transition Bruno et al. (2001)). We find that the known structure CuS.oC24 is mechanically unstable, with an imaginary phonon mode that stabilizes a slight monoclinic distortion, and this distorted monoclinic structure, CuS.mC24, is the true ground state structure according to DFT. The stable phases predicted correspond to known low temperature stable phases, with the exception of CuS.tP12, which is not observed stable at temperatures as low as 300K (its stability is mostly likely due to the use of S as a reference structure), and CuS.mC24, the experimentally observed stable structure being CuS.hP12.

Experimental results Karup-Møller and Makovicky (1999) show PdS.cI46 and PdS.tP10 having solubility ranges for copper, which implies the existence of substitutional or interstitial defects. PdS.cI46 contains four Wyckoff site classes: 8c and 24g sites containing palladium and 8c and 6b sites containing sulfur (see Fig. 2). For PdS.tP10 there are only two site classes, an 8g containing palladium and 2a containing sulfur.

Shown in Figure 3 is our calculated T=0K phase diagram for the Cu-Pd-S system (supporting data is in Table 2). We find the 8c site class in PdS is favorable to copper substitution. At T=0K it is enthalpically favorable to occupy 6 of the 8 possible lattice sites, with a sudden upturn for further filling of the 8c site class. The only other substitutions that are found favorable are S and Cu in -CuPd, Cu in Pd.cF4, and Pd in CuS. Substitutions of Cu in the -CuPd phase have a solubility range from 50% Cu to 55.5% at T=0K. We find the -CuPd phase is only favorable to Cu substitutions on the Pd sublattice, and not Pd substitutions on the Cu sublattice, which is experimentally supported by observations that the palladium content of -CuPd never exceeds 50%. As none of the Cu- and Pd- rich structures are in competition with CuS.tP12, and the entropic contribution is small, we ignore substitutions in CuS.tP12. We also find a small sulfur solubility range in the Cu sublattice of -CuPd at T=0K, which is supported by experimental phase diagrams at T=673K and T=823K.

While these are the only stable T=0K substitutions, we expect that at higher temperatures, substitutions into other sites that had suitably low enthalpic cost will be stabilized by entropy: PdS.tP16 (2c), PdS.cI46 (24g), and PdS.tP10 (8e) sites being possible candidates (see Table 2). However, we chose to exclude modeling the PdS.tP16 2c substitution range, as due to the small concentration of 2c sites in PdS.tP16, combined with competition with other phases with strong entropic effects, we expect the concentration range to be negligible.

To determine solubility ranges more precisely, we model phases analytically (see Table 3), using first principles calculation data to suggest suitable enthalpic models. All models are based off the regular solution model, though with modifications to better reflect the behavior of the phase. We choose the models to smooth out fluctuations in the first principles calculated energies while still keeping the essential trends seen intact. The phases modeled are PdS, PdS, -CuPd, and the CuPd-FCC phase. In all cases we use the ideal entropy (Eq. 6) for sublattice filling. Our model for PdS differs from a regular solution model in that we use a piecewise linear function (see Fig. 6) for the enthalpic part of the free energy. The enthalpic contribution to the free energy of PdS is increasingly favorable up to 6 substitutions (x=0.13) on the 8c sublattice, but from 6 to 8 substitutions (x=0.174) the favorability is reduced. Previous crystallographic work by Matkovic Matkovič et al. (1976) observed a crystalline phase of PdS with 14% Cu concentration, the Cu occupying 8c sites. However this work was performed at 853K and did not suggest any mechanism by which these substitutions are stabilized. For PdS, we use the dilute limit form of the regular solution model (), though it was necessary to add a small quadratic term to the free energy to fit our first principle results. Both CuPd-FCC and -CuPd were fit using a first-order regular solution model, with CuPd-FCC having the entire lattice available for substitution, allowing usage of the standard form of a regular solution model. -CuPd has only the Pd sublattice available for substitution by Cu, which may be modeled using a regular solution model by replacing terms involving with terms involving , where denotes the Cu concentration. We also compute a SQS configuration for fcc at 50 Pd%Wolverton (2001) with purely random correlations until the 8th nearest neighbor and find it lies above the convex hull by 40 meV/atom.

Experimental phase diagrams show a sulfur solubility range near CuPd trending in the palladium-rich direction. We find dilute S substitutions onto the Cu sublattice in -CuPd, creating a phase trending in the palladium-rich direction, are stable at T=0K. Subsequent substitutions become rapidly unfavorable. From this, we expect limited entropic stabilization of additional substitutions of S in -CuPd at higher temperatures, leading to a roughly constant solubility range.

To include phonon free energy in solid solution phases, we calculate phonon free energies in the harmonic approximation at well-ordered configurations at select compositions and linearly interpolate between these compositions across the phase region. To model the free energy of gaseous sulfur, we use empirical data taken from the NIST-JANAF tables, as explained later in the paper.

Base Structure | Site Class | Chemical Species | Enthalpic Cost | Maximum Impurity Concentration |

Change | (meV/atom) | in Site Class | ||

-CuPd | 2a | Pd to Cu | 0 | 12.5 |

-CuPd | 2a | Cu to S | 0 | 3.8 |

-CuPd | 2a | Pd to S | 0 | 3.8 |

Pd.cF4 | 4a | Cu to Pd | 0 | 21.9 |

PdS.cI46 | 8c | Pd to Cu | 0 | 12.5 |

PdS.cI46 | 24g | Pd to Cu | 12.9 | 8.33 |

PdS.tP10 | 8e | Pd to Cu | 4.7 | 3.13 |

CuS.cP12 | 4a | Cu to Pd | 0 | 6.25 |

PdS.tP16 | 2c | Pd to Cu | 5.3 | 50.00 |

Phase | Substitution Type | Phonon Concentrations | ||

(Cu Composition Range) | (eV/atom) | |||

PdS | Cu@Pd8e (0-0.8) | 8/10 | 0, 0.8 | |

PdS | Cu@Pd8c (0-0.175) | 8/46 | 0, 0.175 | |

-CuPd | Cu@Pd2a (0.5 - 1) | 1/2 | 0.5, 1 | |

CuPd-FCC | Cu@Pd4a (0-1) | 1 | 0, 0.50, 0.75, 1 |

Figures 4 and 5 compare our calculated phase diagrams with experiment for T=673K and 823K, respectively. Our T=673K phase diagram matches the experimental results qualitatively. Note, however, that the ”CuS” region of the experimental phase diagram contains multiple structures, many of which are non-stoichiometric, whereas ours contains only one structure. The phase labeled ”CuPd” (or 1D LPS) has a slight solubility range experimentally whereas we model it as a well-ordered line compound. Our T=823K diagrams show more disagreement, as there is only one Cu-S structure that is stable and a PdS structure is stabilized at this temperature. More recent phase diagrams have PdS being stabilized at T=829K, however. This is likely due to entropic effects other than configurational and phononic (in the harmonic approximation) that were not included. As can be seen, PdS is part of our collection of calculated structures but is not a stable low-temperature phase. At both temperatures we attempt to model sulfur solubility in the -CuPd phase but find sulfur solubility to be unstable at temperatures of interest. This is surprising, as our T=0K phase diagram has a sulfur solubility range in -CuPd similar in shape to what is experimentally observed at higher temperatures, and the enthalpic behavior of our sulfur substitutions supports the observed temperature independence of the sulfur solubility range.

The PdS.cI46 copper solubility range is not appreciably affected by the temperature, as it is dominated by enthalpy rather than entropy. The 24g site is unfavorable for copper substitution. The high multiplicity of the 24g site class gives a large entropic contribution to the free energy at high temperature relative to the 8c site class. However, we find that the enthalpic cost (see Figure 6) associated with occupying 24g sites is too high relative to the entropic contributions at the given temperatures and only affects the solubility range by 1% to 2%. Thus the solubility range consists almost completely of 8c sites and is inappreciably affected by temperature. For this reason, our model (Table 3) only includes 8c sites and not 24g sites. For PdS.cI46, at both measured temperatures the experimental copper solubility limit is 20%, and our calculated value is 17.5%, the mismatch likely due to Cu occupation of 24g sites. Our results explain why this solubility range appears to be independent of temperature; it is not stabilized due to entropic contributions but due to enthalpy arising from details of the electronic structure (see subsection 3.2).

In contrast, Cu substitutions in PdS.tP10 are not enthalpically favorable at T=0K but are entropically stabilized at higher temperature. As can be seen in the above phase diagrams, the copper solubility for PdS.tP10 is temperature-dependent in both the calculated and experimental work. While our results agree qualitatively with experimental results, there is still a discrepancy in the copper solubility range. For PdS.tP10, at T=673K the experimental value is 12.7% and ours is 2.5%, and for T=823K the values are 18.6% and 5% respectively. Solid solution models overestimate the configurational entropy, as including correlations between substitutions will alter the energies of distinct configurations with the same composition and reduce the multiplicity. We have thus overestimated the entropy of the PdS phase, whose solubility range extends to almost complete filling of the 8c site, relative to the PdS phase, whose filling in the 8e sublattice is sparse. Thus we would expect more accurate free energy models will lead to decreased entropy in the PdS phase with comparatively small changes in PdS, increasing the solubility range of PdS. This will not affect the stability of PdS as it is primarily enthalpically stabilized.

The large solubility range in the palladium rich side of the CuPd-FCC phase at T=0K is also a factor in the incorrect solubility ranges for PdS.tP10, as we have found that altering the fitting parameters for the CuPd-FCC phase appreciably affects the solubility range for PdS.tP10. As we also encountered an issue with sulfur solubility in -CuPd, with the sulfur solubility range destabilizing relative to CuPd-FCC at temperature below temperatures of interest, and we have a large solubility range on the Pd-rich side at T=0K for the CuPd-FCC phase, it appears that one major source of error is inaccurate modeling of the CuPd-FCC phase. We will return to this issue in the activities portion of this paper.

For our Cu-Pd binary at high temperatures, we find that CuPd.tP28 becomes unstable relative to the CuPd-FCC phase at T=827K, and -CuPd becomes unstable relative to CuPd-FCC and CuPd.tP28 at T=536K and a composition of 59.5%, significantly below the experimentally observed temperature of T=871K but close to the experimental composition of 58%. This relative instability of -CuPd at temperatures of interest is responsible for the observed instability of any sulfur solubility in -CuPd in our T=673K and T=823K phase diagrams. We will show later that first-principles inspired corrections to the enthalpic portion of the free energy for CuPd-FCC will bring the phase transition closer in line with experimental results, though as CuPd-FCC phase has maximal entropy and -CuPd phase (with only Cu substitutions in Pd sites) has minimal entropy in near the equiatomic composition region, it is likely that the assumption of ideal entropy in the FCC phase also plays a role.

But even with these discrepancies, it is clear that we have predicted major details of the Cu-Pd-S ternary phase diagram, specifically the mechanism for stabilization of the solubility ranges of PdS.tP10 and PdS as well as details about the stability of the -CuPd system, using only first principles calculations with no empirical fitting.

### 3.2 Electronic Structure

The question remains why the PdS.cI46 8c site class favors copper in the first 6 substitutions, but substitution in the final two 8c sites are unfavorable. Consider the constituent binaries: copper binds less strongly with sulfur than palladium, and copper binding with palladium is relatively weak, so we would expect copper substitution into Pd-S binaries to be slightly unfavorable. If it were to turn out to be favored, then we would expect Cu to fully occupy the 8c class. First-principles methods provide electronic structure information that can resolve this issue.

Figure 7 shows the electronic densities of states for substitutions in PdS.cI46. Each colored solid line corresponds to the electronic density of states for a different number of substitutions in the 8c site class of PdS.cI46. It can be seen that the densities of states follow a rigid band model. Here copper is being substituted for palladium, and copper has one more electron in its valence shell than palladium. Under a rigid band model, performing a substitution adds one electron per unit cell without appreciably changing the valence electronic eigenstates, increasing the Fermi level relative to the fixed band structure.

In all cases, a pseudogap lies near the Fermi level. In the case of no substitutions, the pseudogap lies to the right of the Fermi level. As the number of copper substitutions increases, the increased electron count drive the pseudogap to the left, towards the Fermi level. At six substitutions the Fermi level lies almost directly at the minimum of the pseudogap. This is a well-known stabilization mechanism for alloys Friedel (1988). As the number of substitutions increases above six, the Fermi level is driven away from the pseudogap, destabilizing the structure and leading to the sudden increase in the enthalpy of formation (see Fig. 6). This pseudogap stabilization mechanism at T=0K, combined with large enthalpic cost of 24g substitutions, explains the temperature independence of the PdS solubility range.

### 3.3 Predominance Diagrams

We now calculate the predominance diagram for the Pd-S system. We work with a fixed quantity of palladium interacting with a reservoir of gaseous S, with temperature and pressure being free parameters. For an ideal gas, this implies that we control the chemical potential of the gas, and thus we work in the semi-grand canonical ensemble:

(9) |

where G is the previously defined Gibbs free energy, is the chemical potential for the sulfur dimers, and is the half the number of sulfur atoms in a given phase (as the chemical potential is for sulfur dimers). To determine sulfidization thresholds, we calculate the semi-grand canonical free energy for all phases of interest, Pd.cF4, PdS.tP10, PdS.cI46, PdS.oC16, and PdS.tP16. The phase with the lowest free energy will be stable. All Pd-containing phases of interest are line compounds, containing only total energy and phonon contributions.

The chemical potential for S complicates the calculation, as it is a gaseous phase. We incorporate empirical results into our method, in a manner consistent with the results obtained from first principles calculations. For the free energy of the gas, we have

(10) |

Here is the chemical potential per molecule, is a reference pressure (by convention taken to be 1 bar, as will be done here), and is the chemical potential per molecule measured at the reference pressure,

(11) |

where H and S are obtained from the NIST-JANAF tables M. W. Chase (1998). obtained from the JANAF tables has a physically unrealistic infinite value at T=0K due to the presence of 1/T terms in both the H and TS contributions. As shown in Figure 8, is nearly linear for temperatures of interest, so we fit to a linear equation using the region to eliminate the T=0K pole. There is still an issue of constant offsets to energy; the JANAF tabulated values are relative to the enthalpy at T=298K, and whereas all our calculated values are computed at T=0K. In order to incorporate the empirical data into our method, we must shift its zero of energy. In particular, at T=0K, we set the empirical data to match what first principles calculations would have predicted,

(12) |

As we have measured all energies relative to the tie-plane between pure elements, , and so the zero-point energy is the only contribution at T=0K taken from first principles calculations.

Shown in Table 4 are our predicted sulfidization thresholds relative to experimental predominance diagrams obtained by Taylor Taylor (1985). Predictions with and without phonon free energies are shown. It is necessary to include phonons to obtain reasonable agreement with experimental results, as when phonon free energies are excluded, the relative errors in pressure are as large as four orders of magnitude. When including phonon free energies, in all cases our pressure estimates are within a factor of 6 of experimental results. It should also be be noted that Taylor’s phase boundaries are themselves based on fits from experimental enthalpy data. One difficulty is that Taylor has a liquid region in the high-temperature, low pressure region of the predominance diagram, whereas we are only interested in solid phases, so experimental thresholds from the solid regions were extrapolated to liquid regions to give numbers suitable for comparison. This extrapolation is valid as Taylor’s thresholds in the solid regions were found to be extremely linear. Another difficulty is that PdS is stable for a small pressure range in between temperatures of 830K and 920K, and it would likely be stable for higher temperatures if not for the liquid phase, but our calculations do not predict it stable in this region. This is most likely due to errors in the phonon free energy, which as already shown appreciably affects results.

T (K) | Transition | Calculated | Calculated | Experimental Taylor (1985) |

Without Phonons | With Phonons | |||

1000 | PdS to PdS | -1.1 | -1.8 | -2.2* |

1000 | PdS to PdS | -1.4 | -2.6 | -2.4* |

1000 | PdS to Pd | -3.3 | -5.3 | -5.2 |

800 | PdS to PdS | -3.1 | -3.5 | -3.7 |

800 | PdS to PdS | -3.5 | -4.5 | -4.7 |

800 | PdS to Pd | -5.9 | -7.6 | ** |

### 3.4 Activity Calculations

In this last section, the activities of species in binary Cu-Pd are calculated, as these can be used to predict sulfidization thresholds for Cu-Pd in the presence of S. We wish here to find the activities across the entire composition range with multiple phases.

The definition for the activity of a species A as a function of composition , pressure P, and temperature T is DeHoff (2006)

where is the activity for species A, is the chemical potential for species A, and is a reference chemical potential for species A. As we are specifying , P, and T, the natural choice of thermodynamic potential is the Gibbs free energy,

where and are the total number of A and B atoms, respectively, and the thermodynamic derivative giving the chemical potential is

However, the quantity of importance for phase transitions is the Gibbs free energy per atom,

We may express the chemical potential in terms of intensive quantities as

where is the composition of species i, and we have chosen to parameterize the free energy by the composition for A. Similarly,

These activities depend only on the free energy and the derivative of the free energy at a given composition. They are piecewise analytic functions across the composition range, completely determined by the thermodynamically stable phase at that composition. The activities can therefore be viewed as alternative representation of the underlying phase diagram for the system. They are discontinuous only when the derivative of the free energy changes (the free energy itself must always be continuous, due to the requirement of convexity).

As mentioned before, we believe the magnitude of the calculated free energy of our fcc phase is artificially large, destabilizing the phase, and thus any calculated activities would be flawed near the 50% Cu region of the phase diagram. To support this assertion, here we refit the fcc phase using only low-solubility first principles data on either side of the composition range. This reduces the magnitude of the free energy of the fcc phase while still using first principles derived calculations. This should not be viewed as abandoning our pure first principles phase diagram as shown before, but rather a correction inspired by known experimental results presented in an alternative-but-equivalent format. The SQS estimate for fcc at 50% Pd lies almost directly on this free energy curve, further supporting its usage for quantitative estimates. We also neglect phonon free energy for activity calculations, as for ordered fcc and bcc at 50% Pd, the difference in phonon free energy per atom is only 2 meV at T=1000K.

Shown below in Figure 9 in solid lines are our calculated activities for palladium in the Cu-Pd binary, using the free energy models from Table 3, the modified fcc phase, and our calculated data for CuPd. The phase is the sharply increasing region in the center of the graph. This sharp behavior is due to the entropic portion of the free energy, which has a logarithmic singularity at x=0.5 and sharply varying derivative near this singularity. The smooth curves on the left and right sides of the graph for T1000K denote the fcc phase. Linear free energy functions give constant activities as a function of composition, and thus phase coexistence regimes appear as plateaus in activities, at a fixed temperature. Two sets of plateaus are present. The set of plateaus on the left side are coexistence regimes between the bcc and fcc phases. There are 3 types of plateaus on the right side. For T 900K, the plateaus are due to coexistence with the line compound CuPd, where the plateaus in the region indicate coexistence with bcc and the plateaus at indicate coexistence with fcc. The plateau on the right side of the plot for T=900K is a special case. Our free energy models predict that the CuPd phase becomes thermodynamically unstable at 899K and fcc becomes stable, so the plateau at 900K indicates coexistence between bcc and fcc. At a temperature of 960K, our model predicts that the bcc phase becomes unstable relative to the fcc phase, disappearing at a composition of x=0.576. By decreasing the fcc free energy magnitude we have properly stabilized the phase at temperatures of interest and brought the transition temperature to within 90K of the expected value. The 1000K and 1350K plots are pure fcc throughout the entire composition range.

The activities of Pd at various temperatures were also calculated using CALPHAD (acronym of CALculation of PHAse Diagrams) method. The thermodynamic database developed by Li et alLi et al. (2008) was used. Our activities are systematically low compared to the CALPHAD activities (shown in dashed lines on Figure 9), most likely due to our ideal entropy approximation. As the true entropy in solid solutions is reduced relative to the ideal entropy, the true free energy will be greater, leading to the true activities being greater than our calculated activities. This is especially true for fcc, where the refitted free energy still overestimates the magnitude of the free energy and leads to qualitatively different behavior in the high Pd limit. CALPHAD results give activities that exceed Raoult’s Law in the high Pd limit, which would require positive enthalpies of formation for Cu substitution in Pd using a solid substitution model in the dilute limit. This contradicts first principle calculations which unambiguously predict Cu substitutions in bulk Pd having negative enthalpies of formation (see Figure 1). Clearly more precise models are needed to understand the binary in the high Pd limit. However, for the phase and associated coexistence regimes, the first principles and CALPHAD results are in close agreement, and if one is interested in regions outside the single-phase fcc region, first principles calculations gives reasonably accurate results even with simple models.

## 4 Conclusions

Using only first principles calculations and solid solution models, we were able to semi-quantitatively model the ternary phase diagram for the Cu-Pd-S system without any empirical fitting parameters. We obtained the sulfidization thresholds for the Pd-S system and activities that reflect the essential features of the Cu-Pd binary system at temperatures where multiple phases exist across the composition range. More detailed models are necessary to reach better quantitative agreement with experimental results. However, our work reveals the essential physics underlying the phase behavior. For example, we distinguish enthalpic vs. entropic driven substitution, and we identify a pseudogap stabilization mechanism, yielding a temperature-independent solubility range for Cu in PdS.

## Acknowledgements

We are grateful to B. Gleeson, M. Mihalkovič, and J. Kitchin for useful discussion. This technical effort was performed in support of the Fuels Program of Strategic Center for Coal at DOE National Energy Technology Laboratory under the RES contract DE-FE0004000.

## References

- Chen and Ma (2010) C.-H. Chen, Y. H. Ma, J. Membr. Sci 362 (2010) 535–544.
- Amandusson et al. (2001) H. Amandusson, L.-G. Ekedahl, H. Dannetun, J. Membr. Sci 193 (2001) 35–47.
- Uemiya et al. (1991) S. Uemiya, T. Matsuda, E. Kikuchi, J. Membr. Sci 56 (1991) 315–325.
- Iyoha et al. (2007) O. Iyoha, R. Enick, R. Killmeyer, B. Morreale, J. Membr. Sci. 305 (2007) 77–92.
- O’Brien et al. (2010) C. P. O’Brien, B. H. Howard, J. B. Miller, B. D. Morreale, A. J. Gellman, J. Membr. Sci. 349 (2010) 380–384.
- Karup-Møller and Makovicky (1999) S. Karup-Møller, E. Makovicky, N. Jb. Miner. Mh. 12 (1999) 551–567.
- Subramanian and Laughlin (1991) P. R. Subramanian, D. E. Laughlin, J. Phase Equilib. 12 (1991) 231–243.
- Li et al. (2008) M. Li, Z. Du, C. Guo, C. Li, CALPHAD 32 (2008) 439–446.
- Trimarchi and Zunger (2008) G. Trimarchi, A. Zunger, J. Phys.-Condens. Mat. 20 (2008) 295212.
- Baärthlein et al. (2009) S. Baärthlein, E. Winning, G. L. W. Hart, S. Müller, Acta Mater. 57 (2009) 1660–1665.
- Ceder et al. (1990) G. Ceder, D. de Fontaine, H. Dreysse, D. M. Nicholson, G. M. Stocks, B. L. Gyorffy, Acta Metall. Mater. 38 (1990) 2299–2308.
- Lee et al. (2007) B. J. Lee, B. Sundman, S. I. Kim, K. G. Chin, ISIJ Int. 47 (2007) 163–171.
- Okamoto (1992) H. Okamoto, J. Phase Equilib 13 (1992) 106–107.
- Kresse and Hafner (1993) G. Kresse, J. Hafner, Phys. Rev. B 43 (1993) 558–561.
- Kresse and Furthmüller (1996) G. Kresse, J. Furthmüller, Phys. Rev. B 54 (1996) 11169–11186.
- Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50 (1994) 17953–17979.
- Hu et al. (2010) R. Hu, M. C. Gao, O. N. Doğan, P. King, M. Widom, CALPHAD 34 (2010) 324–331.
- Perdew et al. (2008) J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, K. Burke, Phys. Rev. Lett. (2008).
- Wu et al. (2003) E. J. Wu, G. Ceder, A. van de Walle, Phys. Rev. B 67 (2003) 134103.
- Owen and Yates (1933) E. A. Owen, E. L. Yates, Philos. Mag 15 (1933) 472–487.
- Harris et al. (1972) I. R. Harris, M. Norman, W. E. Gardner, J. Less-Common Metals 29 (1972) 299–309.
- Okamura (1970) K. Okamura, J. Phys. Soc. Jpn 28 (1970) 1005–1014.
- Jr. (1979) H. T. E. Jr., Science 203 (1979) 356–358.
- Leon et al. (1981) M. Leon, N. Terao, F. Rueda, Phys. Status Solidi A 67 (1981) K11–K14.
- Fjellvåg et al. (1988) H. Fjellvåg, F. GrÃ¸nvold, S. StÃ¸len, A. F. Andresen, R. Müller-Käfer, A. Simon, Z. Krist. 184 (1988) 111–121.
- King and Prewitt (1979) H. E. King, C. T. Prewitt, Am. Mineral. 64 (1979) 1265–1271.
- Grønvold and Røst (1962) F. Grønvold, E. Røst, Acta Cryst. 15 (1962) 11–13.
- Matkovič et al. (1976) P. Matkovič, M. El-Boragy, K. Schubert, J. Less-Common Met. 50 (1976) 165–176.
- Gaskell (1937) T. F. Gaskell, Z. Kristallogr. 96 (1937) 203–213.
- Grønvold and Røst (1957) F. Grønvold, E. Røst, Acta Cryst. 10 (1957) 329–331.
- Janosi (1964) A. Janosi, Acta Cryst. 17 (1964) 311–312.
- Guymont and Gratias (1976) M. Guymont, D. Gratias, Phys. Status Solidi A 36 (1976) 329–334.
- Røst (1968) E. Røst, Acta Chem. Scand. 22 (1968) 819–826.
- Soutter et al. (1971) A. Soutter, A. Colson, J. Hertz, Mem. Etud. Sci. Rev. Metall. 68 (1971) 575–591.
- Geisler and Newkirk (1954) A. H. Geisler, J. B. Newkirk, J. Met. 6 (1954) 1076–1082.
- Yamamoto and Kashida (1991) K. Yamamoto, S. Kashida, J. Solid State Chem. 93 (1991) 202–211.
- Topa et al. (2006) D. Topa, E. Makovicky, T. Balić-Z̆unić, Can. Mineral. 44 (2006) 497–505.
- Bruno et al. (2001) E. Bruno, B. Ginatempo, E. S. Giuliano, Phys. Rev. B 63 (2001) 174107.
- Wolverton (2001) C. Wolverton, Acta Mater. 49 (2001) 3129–3142.
- Friedel (1988) J. Friedel, Helv. Phys. Acta 61 (1988) 538–555.
- M. W. Chase (1998) J. M. W. Chase, J. Phys. Chem. Ref. Data 9 (1998) 1–1951.
- Taylor (1985) J. R. Taylor, Met. Trans. B 16B (1985) 143–148.
- DeHoff (2006) R. DeHoff, Thermodynamics in Materials Science, 2nd ed., CRC Press, 2006.