Thermodynamics of the Atomic Distribution in Pt3Pd2, Pt2Pd3 and their Corresponding (111) Surfaces

In this study, we have developed solid-state models of platinum and palladium bimetallic catalysts, Pt3Pd2 and Pt2Pd3, which are rapidly thermally annealed at 800 °C. These models were constructed by determining all the unique atomic configurations in a 2 × 2 × 1 supercell, using the program Site-Occupation Disorder (SOD), and optimized with the General Utility Lattice Program (GULP) using Sutton-Chen interatomic potentials. Each catalyst had 101 unique bulk models that were developed into surface models, which were constructed using the two-region surface technique before the surface energies were determined. The planes and compositions with lowest surface energies were chosen as the representative models for the surface structure of the bimetallic catalysts. These representative models will now be used in a computational study of the HyS process for the production of hydrogen.


Introduction
As research into alternative energy becomes increasingly important, hydrogen has been proposed as a potential replacement for hydrocarbon fuels. 1 Hydrogen has the advantages of being environmentally friendly, readily available and possessing a high energy density. 2 However, hydrogen as replacement fuel is held back by many challenges, including storage and production technology shortfalls. A current industrial-scale hydrogen production method is electrolysis, where water is split into hydrogen and oxygen. Unfortunately, water splitting using electrolysis is an energy-intensive process, and the goal is therefore to create more efficient and environmentally sound routes to produce hydrogen. One such route could be the hybrid sulphur (HyS) process, which was first proposed by the Westinghouse Corporation. 3 The HyS process produces hydrogen via the electro-oxidation of aqueous SO 2 at the anode of an SO 2 depolarized electrolyzer (SDE). 4 Various metal catalysts (platinum, palladium and gold), serving as the anode in the SDE, have been reviewed by Díaz-Abad et al. 5 Platinum was shown to be the most effective catalyst, with a high catalyst activity and high resistance to corrosion in the acidic environment present at the anode. However, scarcity and cost are major fundamental limitations associated with Pt, and research into lowering the amount of platinum used in catalysts is therefore vital if hydrogen is to become a viable replacement energy source.
One study into alternative platinum-based catalysts was carried out by Falch et al. 6 who found potential replacement catalysts by depositing multi-metal thin films through physical vapour deposition onto glassy carbon electrodes as the substrate. The study showed that thin-film combinations of Pt 3 Pd 2 and Pt 2 Pd 3 were two of the best performing catalyst as anode materials for SO 2 oxidation, achieving onset potentials of 0.587 V and 0.590 V, respectively. 6 The onset potentials of these catalysts compare well with the value reported for standard platinum as the catalyst (0.598 V), 7 whilst using less expensive alloys. However, one problem related to these catalysts is the tendency of the thin metal films to delaminate, to certain degrees, from the glassy carbon surface during testing. This problem was addressed by rapid thermal annealing at 800°C, which reduced the delamination phenomenon and increased stability. 4,8 It was shown that the annealing drastically changes the surface structure of the catalysts, which results in the improved durability of the catalyst. The XRD analysis of the catalysts indicates that the (111) Miller plane peak increases with annealing temperature, while the (200) and (311) planes, as identified by Mahapatra and Datta, 9 remain mostly unchanged. The (111) peak also correlates well with the catalyst activity, 10 indicating that this (111) surface enables active site formation. These potential replacement catalysts perform to a similar degree as the standard platinum catalysts in terms of activity, but they have a longer life-time compared to the standard platinum catalysts. 8a Given the promise of Pd/Pt alloys as catalysts in the HyS process, this study aims to develop solid-state models of platinum and palladium bimetallic catalysts, Pt 3 Pd 2 and Pt 2 Pd 3 , annealed at 800°C to be used in the modelling of the HyS process for the production of hydrogen. To this end, we will identify and construct the most stable Pt 3 Pd 2 and Pt 2 Pd 3 bi-metallic bulk material, followed by the creation of (111) surface models.

Computational Methods
Various bulk configurations, at the desired ratios of Pt 3 Pd 2 and Pt 2 Pd 3 , were constructed by substitution of Pd into all available sites of the pure Pt bulk and vice versa, using the program Site-Occupation Disorder (SOD) 11 in a similar manner to the recent study conducted by Botha et al. 12 The unit cells of Pt and Pd were constructed with F-type lattice, in a face-centred cubic (fcc) cell with space group Fm-3m (number 225) by entering the crystallographic properties into SOD. The cell size was limited to a 2 × 2 × 1 supercell as larger supercells were unfeasible due to the large numbers of configurations generated for the different alloy compositions that would need to be considered, which would be beyond the computational capabilities of the SOD program. In this instance, SOD can treat a 2 × 2 × 1 cell with 16 atoms within a reasonable timeframe. However, a 2 × 2 × 2 cell with 32 atoms leads to so many different configurations, that its computation cannot be completed within a sensible timeframe. 13 The unit cell of Pt was expanded into a 2 × 2 × 1 rectangular cell with a and b lengths of 7.847 Å and c of 3.924 Å, and cell angles of 90.0°, based on the lattice constant at 25°C obtained from Arblaster. 14 The Pd cell had lengths a, b of 7.780 Å and c with a length of 3.890 Å, and cell angles of 90.0°, again based on lattice constants at 25°C obtained from Arblaster. 15 In addition to generating all inequivalent distribution configurations for each alloy composition, SOD also provided various properties and energies of each bulk configuration, i.e. Boltzmann probability distributions, bulk energies and thermodynamic energies at different temperatures. The temperatures used were in the range of -272°C to 900°C, as 900°C was the highest temperature used in the annealing of the catalysts as reported by Falch et al. 8a The optimized structures were analyzed, and the most probable bulk configurations determined using the Boltzmann probabilities. The list of configurations obtained from SOD was geometry optimized using the molecular mechanics based General Utility Lattice Program (GULP), 16 with the interatomic potentials from the Sutton-Chen library of potentials. 17 The geometry optimization was carried out under constant pressure, until the equilibrium volume was reached.
The surface energies for Pt, Pd and all 101 inequivalent configurations for both Pt 3 Pd 2 and Pt 2 Pd 3 catalysts were calculated. However, before the surface energies could be calculated the cubic symmetry of the bulk was re-obtained through resizing the bulk to a 2 × 2 × 2 cell by doubling in the z-direction, thereby increasing the simulation cells to 32 atoms. The Pt 3 Pd 2 bulk materials now contained 12 Pd and 20 Pt atoms, with cell dimensions of 7.8472 Å and cell angles of 90.0°. The Pt 2 Pd 3 bulks now contained 12 Pt and 20 Pd atoms, with cell dimensions of 7.7804 Å and cell angles of 90.0°. Slab structures without dipole were created using the program METADISE 18 (minimum energy techniques applied to dislocation interface and surface energies), to cut geometrically optimized bulk structures along the (111) Miller plane. Additionally, METADISE was used to determine the symmetrical terminations of each configuration, which is required for the determination of the surface energies. 19 The surface energy was determined by employing a two-region method 18 , as shown in Fig. 1. Region 1 represents the atoms near the surface which are all allowed to relax fully, while Region 2 represents the rest of the bulk material where the atoms are kept fixed in their bulk-optimized positions. 13,20 The addition of Region 2 is necessary to ensure that the potentials acting on the atoms in Region 1 at the interface with Region 2 and close to the surface are modelled correctly. The surface thickness (Region 1) and the bulk thickness (Region 2) must be optimized to model a sufficiently large system to achieve convergence and increase accuracy, whilst keeping the system small enough to minimize computing time. This optimization was done for one of the configuration structures and assumed to be the same for the remainder of the compositional configurations.
The total energy of a surface (U Total ), consisting of the two regions, is calculated using Equation (1), where the subscript denotes the interactions between the regions. However, the energy of Region 2 (U 22 ) is not required as it is cancelled out in the surface energy calculations, as it describes the energy of an infinitely large bulk below the surface, but it is required to obtain convergence of the calculations. To calculate the energy of the relaxed surface (U Surface ), the energy of Region 1 (U 11 ) and half the energy interaction between Region 1 and Region 2 (U 12 ) is required, as shown in Equation (2). The surface energy (ΔU SE ), written in terms of Equation (3), can now be calculated by subtraction of the bulk energy (U Bulk ) from the energy of the U Surface per unit surface area (A). U Bulk is obtained by multiplying the energy value of the original 16-atom cell by the number of times this cell is occurs in Region 1 and Region 2. A lower positive value for the surface energy (ΔU SE ) of a solid indicates a more thermodynamically stable surface. 13 U Total = U 11 + U 12 + U 22 (1) To validate the surface energy calculations, surface energies of the pure Pt and Pd metals were calculated and compared to values reported in the literature. The Pt and Pd bulk metals were constructed using the same process that was used to create the bulk of the Pt 3 Pd 2 and Pt 2 Pd 3 . Both Pt and Pd bulks were optimized and cut along the (111) Miller plane and constructed into surface models with the same Region 1 and 2 thicknesses as stated earlier. Surface energies of Pt (0.90 J m -2 ) and Pd (0.99 J m -2 ), as reported by Todd and Lynden-Bell 21 were obtained using a similar method and also employing the Sutton-Chen interatomic potentials. Similarly, Kimura et al. 22 reported surface energies, using Sutton-Chen potentials, for Pt and Pd at 0.913 J m -2 and 1.003 J m -2 .

Results
In the model of the Pt 3 Pd 2 catalyst with the desired concentration of Pt 60.0 % Pd 40.0 %, a total of 6 out of the 16 Pt atoms in the bulk were substituted by Pd using the SOD program, which led to a Pt concentration of 62.5 % and Pd of 37.5 %. Similarly, the Pt 2 Pd 3 model was constructed by replacing 6 Pd atoms with Pt in the Pd bulk giving concentrations of Pd at 62.5 % and Pt 37.5 %. These model concentrations are similar to the composition of the experimental results, as reported by Falch et al. 6 For both models, SOD generated 8008 possible configurations, but of those only 101 were unique. SOD uses Boltzmann probabilities to predict which configurations are likely to occur at specific temperatures and these probability predictions were used to further reduce the number of required bulk calculations by focussing only on the configurations which are likely to occur at the desired temperature; in this study, a probability threshold of 1 % was used. Figure 2 indicates the number of bulk configurations that are probable at temperatures ranging from -272°C to 900°C for both Pt 3 Pd 2 and Pt 2 Pd 3 systems.
At -272°C, only one configuration is probable in each catalyst, which correlates to the lowest-energy configuration. However, as the temperature increases, the number of configurations increases slowly to 12 and 6 configurations at -72°C for Pt 3 Pd 2 and Pt 2 Pd 3 , respectively. Pt 3 Pd 2 shows an exponential increase in the number of configurations between 0°C and 400°C from 28 to 94 configurations, until the maximum number of configurations, namely 101 configurations, is reached at 800°C. In terms of Pt 2 Pd 3 , the rate at which the number of configurations grows is much slower, only reaching 95 at 900°C. This difference in the number of configurations between the two compositions would suggest that the Pt 2 Pd 3 catalyst requires higher annealing temperatures to achieve more random configurations. However, in experimental work 6 these catalysts are rapidly thermally annealed at 800°C, and, as illustrated, all 101 inequivalent bulk configurations of Pt 3 Pd 2 and 94 of Pt 2 Pd 3 are probable at this temperature. For the sake of completeness and comparison, it was decided to investigate all 101 configurations for both systems. Table S1 and Table S2 found in the supplementary information, list the ten most likely configurations at different temperatures of the Pt 3 Pd 2 and Pt 2 Pd 3 systems, respectively. Included in the tables are relative bulk energies and Boltzmann probabilities.
These tables indicate how the distribution of the configurations is spread more evenly as the temperature increases, highlighting the necessity to take into account not only lowest-energy configurations but also higher energy configurations which become probable at the annealing temperatures. For example, at 25°C configuration 94 of the Pt 3 Pd 2 catalyst has a probability of 25.36 % (highest probability configuration), but at 800°C, it has fallen to sixth position with a probability of only 3.55 %, even though it is the bulk configuration with the lowest energy. This behaviour can be explained by the level of the degeneracy of the different configurations; a higher degeneracy will allow a configuration to become more probable at higher temperatures. In the case of the Pt 2 Pd 3 catalyst, the lowest-energy configuration, namely 94, remains the most probable configuration as predicted by the Boltzmann probabilities. This suggests that the other configurations do not have a high enough degeneracy, or their relative energies are too high to compete effectively with configuration 94 to become the most likely configuration in the distribution. However, we note that the differences between the relative bulk energies, of the configurations, are less than 0.1 kJ, supporting the conclusion that at higher temperatures a larger number of configurations becomes accessible. However, the Boltzmann probability of the bulk configurations is complemented by a description of the surface energies to identify the ground state configuration.
First, surface energy calculations were carried out on different sizes of Regions 1 and 2. Figure 3 shows the results of varying the sizes of Region 1 and Region 2, using the bulk of Pt 3 Pd 2 configuration 1 as a representative example. For Region 1 the surface energy is stable at 3 -4 layers, however, it was decided to double the thickness to seven layers for all proceeding surface energy calculations to ensure convergence for the remainder of configurations. Region 2 shows stability at five layers and again it was doubled in thickness to 10 layers for all surface energy calculations, to ensure that the system converges with the remainder of configurations.
Hence all subsequent surfaces were constructed with a Region 1 with 7 layers and Region 2 with 10 layers.  The calculated surface energies of the two catalyst systems are displayed in Fig. 4 and Fig. 5. The surface energies for Pt 3 Pd 2 vary from 0.817 J m -2 to 1.007 J m -2 , with configuration 3 considerably more stable (lower surface energy) than the rest of the configurations, which will, therefore, be the most likely representative surface configuration of the Pt 3 Pd 2 surface. The surface energy results for the Pt 2 Pd 3 vary less between the minimum and maximum values with an energy range of 0.887 J m -2 to 1.007 J m -2 . There is a group of three configurations that are notably more stable than the rest of the 101 configurations, i.e. configurations 17, 81 and 8 with energies calculated at 0.887 J m -2 , 0.888 J m -2 and 0.889 J m -2 , respectively. These configurations will, therefore, most likely occur in the surface structure of the Pt 2 Pd 3 model.
The four best performing surface configurations, as listed above, visualizations of the top layer of the surfaces are shown in Fig. 6. The surface area of each configuration has been quadrupled in size, by doubling the lengths of a and b, of the periodic cell. This was done to easily identify the arrangement of Pt and Pd, and determine the atomic sequences. The sequences are described in the following in notation, in which the square brackets denote the beginning and end of a sequence and the round brackets denote the individual rows, as seen from left to right. Figure 6a shows Pt 3 Pd 2 configuration 3, the surface structure consists of a row with the following sequence [(Pt), (Pd-Pd-Pd-Pt)]. Figure 6b shows Pt 2 Pd 3 configuration 17, the surface structure contains the following sequence [(Pt-Pt-Pd-Pd), (Pd-Pd-Pd-Pt), (Pd-Pt), (Pd-Pd-Pd-Pt)]. Figure 6c shows Pt 2 Pd 3 configuration 81, with the following surface structure [(Pt-Pd-Pd-Pd), (Pt-Pd)]. Figure 6d shows Pt 2 Pd 3 configuration 8 with the surface structure with the following sequence [(Pt-Pd), (Pd), (Pt), (Pd)]. The top layer arrangement was highlighted to see if any correlation between configurations, which would explain the reason for the lower surface energies. However, as can be seen, the configurations show no correlation to one another, and there are no obvious similarities between them.

Conclusions
This study aimed to construct solid-state models of two bi-metallic catalysts, namely Pt 3 Pd 2 and Pt 2 Pd 3 , which could then be used in the modelling of the HyS process for the production of hydrogen. The bulk materials for both catalyst compositions were constructed by using the program SOD to create   unique bulk configurations by replacing single atoms from the pure metal by the alloying metal. The unit cell of the two compositions produced a total of 8008 configurations with 101 unique configurations for both Pt 3 Pd 2 and Pt 2 Pd 3 systems. The results obtained from the bulk calculations indicated that all 101 configurations had to be considered, as the high-temperature annealing process of the catalysts would allow most of the configurations to occur. Hence, another method of choosing a representative surface model was needed. Thus, we created surface models from the different bulk configurations, which were used to calculate the surface energies of each of the catalyst systems, with the configurations with the lowest surface energies selected as being representative of the catalyst surface. These surfaces were not constructed from the bulk configuration with the lowest energy, justifying our decision to select the representative surface configurations based on their surface energies. This method of calculating surface energies using Sutton-Chen interatomic potentials, was validated by comparing the surface energies of the pure metal systems against those reported in the literature, with good agreement between our results and the literature values. The calculated surface energies of the Pt 3 Pd 2 catalyst indicate that one configuration had notably lower surface energy and is hence a significantly more stable surface than the rest, with a surface energy of 0.817 J m -2 ; visualizations and Cartesian coordinates can be found in Fig. S1 in the supplementary information. The Pt 2 Pd 3 system had three configurations with comparable low energies, calculated at 0.887 J m -2 , 0.888 J m -2 and 0.889 J m -2 , respectively, and shown in Figs. S2-S4 in the supplementary information.
The stacking sequences of the four lowest energy surfaces show no evident correlation or similarities, which indicates that the minor differences in surface energy that occur are caused by nearest-neighbour interactions. As the compositions of each of the surfaces remain the same, only these configurational differences can be the reason for the variations in surface energies. However, due to the classical nature of interatomic potentialbased simulations, electronic properties of the surfaces cannot be obtained. Future work should therefore include investigations of these surface models using quantum mechanical techniques, such as calculations based on the Density Functional Theory.

Supplementary material
A supplementary document is available with this article that has additional information of the results obtained from the SOD program. Additionally, the supplementary document gives the Cartesian coordinates and visualizations of the four lowest surface energy configurations, namely Pt 3 Pd 2 configuration 3 and Pt 2 Pd 3 configurations 8, 17 Figure S1: Visualisation of (a) across and (b) onto views Pt3Pd2 configuration 3, including the (c) list of Cartesian coordinates of the 32 atom bulk 4. Figure S2: Visualisation of (a) across and (b) onto views Pt2Pd3 configuration 17, including the (c) list of Cartesian coordinates of the 32 atom bulk 5. Figure S3: Visualisation of (a) across and (b) onto views Pt2Pd3 configuration 81, including the (c) list of Cartesian coordinates of the 32 atom bulk 6. Figure S4: Visualisation of (a) across and (b) onto views Pt2Pd3 configuration 8, including the (c) list of Cartesian coordinates of the 32 atom bulk