A tale of two phase diagrams: Interplay of ordering and hydrogen uptake in Pd-Au-H

Due to their ability to reversibly absorb/desorb hydrogen without hysteresis, Pd--Au nanoalloys have been proposed as materials for hydrogen sensing. For sensing, it is important that absorption/desorption isotherms are reproducible and stable over time. A few studies have pointed to the influence of short and long range chemical order on these isotherms, but many aspects of the impact of chemical order have remained unexplored. Here, we use alloy cluster expansions to describe the thermodynamics of hydrogen in Pd--Au in a wide concentration range. We investigate how different chemical orderings, corresponding to annealing at different temperatures as well as different external pressures of hydrogen, impact the behavior of the material with focus on its hydrogen absorption/desorption isotherms. In particular, we find that a long-range ordered L1$_2$ phase is expected to form if the \ce{H2} pressure is sufficiently high. Furthermore, we construct the phase diagram at temperatures from \unit[250]{K} to \unit[500]{K}, showing that if full equilibrium is reached in the presence of hydrogen, phase separation can often be expected to occur, in stark contrast to the phase diagram in para-equilibrium. Our results explain the experimental observation that absorption/desorption isotherms in Pd--Au are often stable over time, but also reveal pitfalls for when this may not be the case.


I. INTRODUCTION
The prospect of hydrogen as a replacement for fossil fuels continues to generate interest in metallic hydrides. Some metals are of particular interest in this context as they can be reversibly loaded and unloaded with hydrogen by tuning the partial pressure of H 2 gas in the environment; Pd is one such example. When loaded with hydrogen, the Pd hydride has a hydrogen density that is orders of magnitudes larger than H 2 gas under the same conditions, and it is therefore a candidate for storing hydrogen [1]. Another potential application of Pd hydrides is sensing [2,3]. Pd nanoparticles that are exposed to hydrogen quickly form a hydride [4], and when doing so, their optical properties change. This change can be easily detected, and Pd nanoparticles can thus be used as reversible hydrogen sensors [5,6]. The use of Pd for hydrogen sensing applications has, however, a major disadvantage; Pd hydride formation at room temperature is associated with a first-order phase transition from a hydrogen-poor α phase to a hydrogen-rich β phase, which causes the sensor response to be a highly nonlinear function of H 2 pressure. Since the phase transition is also associated with hysteresis [7][8][9], the response moreover differs between the loading and unloading half-cycles. These aspects are generally unfavorable for sensing applications.
Fortunately, these drawbacks can be alleviated by alloying Pd with Au. With around 20 mol-% of Au in Pd, the system can be loaded continuously with hydrogen [10,11], eventually making the sensor readout an almost linear function of H 2 pressure [12][13][14][15]. The introduction of another chemical species thus fundamentally changes the thermodynamics of the material and improves its properties. From a scientific as well as practical standpoint, however, an alloy is significantly more complex than a pure metal, and multiple important questions arise: what is the optimal composition, how are the chemical elements ordered under different circumstances, and how does this ordering influence the properties of the ma- * erhart@chalmers.se terial? The influence of overall composition can usually be relatively easily optimized by experimental screening. The chemical ordering, on the other hand, is difficult to both control and measure experimentally, and can be expected to be influenced by conditions during and after fabrication. Moreover, although Pd-Au nanoalloy sensors tend to be stable over time [16], changes may occur over long time-scales due to slow kinetics, which makes the effects tedious to experimentally assess, while they may still be detrimental to the material. Although the study of hydrogen in Pd and its alloys has a long and rich history [7,[17][18][19][20][21][22][23], many fundamental aspects of the Pd-Au-H system are still unknown, and to fully exploit this material it is of paramount importance that its thermodynamics are understood, as it determines the driving forces that underpin the evolution of the atomic scale structure over time.
A few previous studies have successfully used combinations of density-functional theory (DFT) calculations and statistical methods to predict properties of hydrogen in Pd-Au. For example, Mamatkulov and Zhdanov [24] recently developed a model based on DFT calculations that successfully predicted absorption isotherms, highlighting that once octahedral sites without surrounding Au atoms become few, the phase transition and its associated hysteresis disappears. These calculations assumed a Pd-Au lattice with random configuration, which remained unaffected by insertion of hydrogen. This metastable equilibrium, with equilibrium only on the hydrogen sublattice, is commonly referred to as para-equilibrium [25][26][27]. If the metal atoms are allowed to rearrange in response to hydrogen such that full equilibrium is reached, a significantly more complex phase diagram can be expected [28]. This was indicated experimentally by Lee et al. [29], who found that long-range order (LRO) was formed in Pd-Au with 19% Au when subjected to a high pressure of H 2 at high temperatures. Importantly, alloys with this LRO had strikingly different pressure-composition isotherms than those that had not been treated in a high-pressure H 2 gas and lacked LRO. A similar conclusion was reached by Chandrasekhar and Sholl et al. [30], who computed hydrogen solubility as a function of chemical order at Au contents 4% and 15%. These results show that the different degrees of equilibration matter arXiv:2103.14340v1 [cond-mat.mtrl-sci] 26 Mar 2021 for the properties of the Pd-Au-H system, and since many applications depend on reproducible pressure-composition isotherms, chemical order should not be ignored. In this light, it is obvious that a full accounting of the impact of chemical order on the thermodynamics of Pd-Au hydride is desirable.
Here, we make a detailed atomic scale study of the thermodynamics of hydrogen in Pd-Au. To this end, we use alloy cluster expansions (CEs) fitted to DFT to describe the energetics of the system ranging from 0% to 50% Au in Pd and in the full range from no hydrogen to 100% hydrogen (i.e., 1 H atom per 1 Au or Pd atom). This approach, which has been successfully applied to similar materials in the past [30][31][32][33][34], allows for very fast evaluation of the energy of a system with hundreds of atoms with an accuracy approaching that of DFT. Using Monte Carlo (MC) simulations, we can elucidate the thermodynamics and, in particular, study how chemical order evolves under different circumstances.

A. Structure generation
In Pd-Au hydride, hydrogen primarily occupies octahedral interstitial sites in the face-centered cubic (FCC) lattice formed by Pd-Au (although in some configurations, tetrahedral site may be preferred [24,35]). The full system can thus be viewed as a rock salt structure, i.e., two intertwined FCC sublattices, one occupied by Pd and Au, and the other populated with H and vacancies. In the following, we will state concentrations per sublattice, such that, for example, the H concentration is 100% if all sites on the H/vacancy sublattice are occupied by H, whereas energies are stated per formula unit (f.u.) with a f.u. being two sites, one per sublattice.
Structures were generated by exhaustive enumeration [36,37] of supercells of the primitive structure with up to 12 atoms in the cell (6 per sublattice). The limited supercell size precludes structures with Au concentration between 0 and 16.67% but as this is a composition interval of particular interest, we generated an additional 200 structures in this interval having up to 32 atoms in the cell, heuristically chosen so as to properly span the space of cluster vectors. Finally, in order to properly describe the important limit of pure Pd hydride, we generated 367 structures without Au, having up to 24 atoms in the cell.

B. DFT calculations
For generation of training data, we used the projector augmented wave formalism as implemented in the Vienna ab initio simulation package (version 5.4.1, PAW 2015) [38,39] with the vdW-DF-cx exchange-correlation functional [40,41]. Wave functions were expanded in a plane wave basis set with a cutoff of 400 eV, the Brillouin zone (BZ) was sampled with a Γ-centered grid with a k point spacing of 0.2Å −1 , and occupations were set using Gaussian smearing with a width of 0.1 eV. Atomic positions and cell shape were relaxed until residual forces were below 25 meV/Å and stresses below 1 kbar. Once converged, we ran an additional single-point calculation with a k point spacing of 0.1Å −1 and the BZ was integrated using the tetrahedron method with Blöchl corrections.

C. Alloy cluster expansions
In an alloy CE, the energy of a configuration of atoms σ is decomposed into contributions from single sites ("singlets"), pairs of sites, triplets of sites, etc, commonly referred to as clusters. Denoting such clusters α, the energy can be expressed as where the effective cluster interaction (ECI) J α measures the energy associated with the cluster α (with J 0 being a constant, the "zerolet"), m α is the number of symmetrically equivalent clusters α, and Π α (σ) α quantifies the chemical ordering of the clusters α. Specifically, Π α (σ) α is an average over all symmetrically equivalent clusters α, where each term is a product of so-called point functions, meaning that the product runs over all lattice sites i in the cluster α. Here, we explicitly treat unoccupied interstitial sites as vacancies, which makes it possible to view the interstitial sublattice as an alloy. In our case, the point functions then take the value Θ i = −1 if site i is occupied by Au or H, or Θ i = 1 if site i is occupied by Pd or a vacancy (but the sublattices are distinct, meaning that, e.g., hydrogen cannot occupy a site on the Pd-Au sublattice).
Although it can be shown that this basis set forms a complete orthonormal basis and thus can be used to exactly express any function of the configuration σ [42], in practice it has to be truncated. Physical intuition tells us that clusters with few sites and short interatomic distances should be more important, and by analyzing the Akaike information criterion (AIC) and Bayesian Information Criterion (BIC) [43] (Fig. S3), we found that it was sufficient to include pairs no more than 1.75a 0 apart (where a 0 is the lattice parameter of the conventional (cubic) cell) and triplets that had all interatomic distances shorter than 1.25a 0 . With these cutoffs, we end up with 18 and 42 symmetrically distinct pairs and triplets, respectively.
The ECIs J α were deduced by fitting the linear equation (1) to the energies calculated with DFT using a Bayesian approach [44]. We assumed that the errors were normally distributed with variance σ 2 1 , and used an inverse gamma function as prior for σ 1 . Furthermore, we used Gaussian priors for the ECIs (having mean µ = 0 and variance σ 2 2 ), and used an inverse gamma function as prior for σ 2 . (We note that the fitting process was largely insensitive to the choice of priors due to the large training data set.) We then sampled the posterior distribution using Markov chain Monte Carlo (MCMC) simulations as implemented in the EMCEE software [45] with 300 walkers and 417,000 MC steps (100 times the auto-correlation length). From the resulting posterior distributions, we constructed our final model as the average over the posterior distribution for each ECI (we used the average rather than the peak of the posterior since the posteriors were almost perfectly symmetric and the two measures therefore coincided to within less than 0.5 meV/f.u., and the average is insensitive to noise and choice of bins). All of the results presented below were obtained with this CE unless otherwise stated. In addition to this CE, we randomly selected 10 independent CEs from the MCMC trajectory, which were used to investigate the impact of variations in the ECIs on thermodynamic properties of the model. For details on priors and posterior distributions, see It is usually practical to fit CEs to mixing energies, defined such that the pure phases have zero mixing energy. In our system, we have four pure phases: Au and Pd with and without hydrogen. Only linear transformations are permissible from a thermodynamic perspective since only first-order terms cancel. This implies that only three (out of four) boundary points can be set strictly to zero. Here, we make the (intuitively logical) choice to set Pd, PdH and Au to zero, leaving AuH using the following definition of the mixing energy: where N is the number of f.u.s in the structure (i.e., the total number of Au and Pd atoms). Generation of structures with exhaustive enumeration always yields disproportionately many structures close to 50% concentrations, as it corresponds to the largest number of symmetrically distinct structures. To counter this bias, we constrained the CE to always reproduce the mixing energy of Pd, PdH, and Au to their exact mixing energy values (0 meV/f.u. per definition).

D. Monte Carlo simulations
MC simulations were carried out according to a standard Metropolis scheme, in which the chemical identity of one or more atoms is changed in each step, and the change is accepted with probability P = min {1, exp (−∆ψ/k B T )}. The change in the potential ∆ψ can be modified to account for different physical situations. Here, we used three different approaches.
The first approach is to sample the canonical ensemble by letting ψ be the change in potential energy, ψ = E, and always swap the identities of two unlike atoms on the same sublattice. This has the advantage of preserving the overall concentration of the chemical species and is thus reminiscent of typical experimental conditions for the Pd-Au lattice. In the case of hydrogen, however, one experimentally controls the chemical potential rather than the concentration, and it becomes impossible to relate the conditions to a partial pressure of H 2 since the chemical potential is not an observable of the canonical ensemble.
Second, if we let ψ = E + µ H N H , where µ H and N H are the chemical potential of hydrogen and the number of hydrogen atoms, respectively, we sample the grand canonical ensemble with respect to the H sublattice. In this case, the chemical potential of hydrogen is specified rather than the concentration, similar to experiment, and it is possible to relate the conditions to a partial pressure of H 2 (see Sect. II E). Moreover, one has access to a free energy derivative since If we view vacancies as a chemical species, we may also refer to this as sampling in the semi-grand canonical (SGC) ensemble (with the difference in chemical potential simply being the chemical potential of hydrogen as the chemical potential of vacancies is zero) and this term will be used below. Finally, we also used the variance-constrained semi-grand canonical (VCSGC) ensemble [46,47], in which ψ = E + N k B Tκ(c +φ/2) 2 , whereκ andφ control the variance and mean of the concentration c, respectively. With this approach we may, unlike the SGC, stabilize the system inside multiphase regions and thereby integrate free energies across phase boundaries using the free energy derivative where c is the average observed concentration. By equating Eq. (4) and (5), we may also relate an observed state to an external partial pressure of H 2 .
It should be noted that the two sublattices may be sampled in different ensembles in the same MC simulation. In this work, we have used several combinations of the three abovementioned ensembles as will be detailed below. For each set of independent parameters, we ran 1,000 MC sweeps (1 sweep = N steps, where N is the number of sites in the cell), of which the first 50 were discarded to allow for equilibration. The simulations were commenced from a random H-vacancy sublattice with equimolar concentrations, whereas the starting configuration of the Pd-Au sublattice was either random, or, in the case of random equilibrium and para-equilibrium (see Sect. III D) as well as for the simulation of isotherms, fixed in a predetermined configuration. The simulations were done with a supercell corresponding to 6 × 6 × 6 times the conventional cell, for a total of 1,728 sites (864 per sublattice). All MC simulations were carried out with the MCHAMMER module of the ICET package [36].

E. Conversion between chemical potential and partial pressure of hydrogen
Experimentally, the chemical potential of hydrogen, µ H = 1 2 µ H2 , is controlled via the partial pressure of H 2 , p H2 . Treating H 2 as an ideal gas, we can write the following relation between its chemical potential and pressure, Here, p • H2 is a reference pressure and µ • H2 (T ) is the temperature dependent chemical potential of H 2 at this reference pressure. The latter can in principle be evaluated with a combination of DFT calculations and thermochemical tables, but since small errors in µ • H2 (T ) have a large impact on p H2 , we chose to make a two-parameter fit. Specifically, we extracted experimental data of the partial pressures of the α → β transition in pure Pd at different temperatures, and fitted a linear function to these (ln(p), T ) data points. We then fitted a linear function to the points (2µ H , T ) at which we observed the phase transition in our own simulations. The conversion is then given by where ∆P (T ) is the difference between the two linear fits, specifically ∆P (T ) = 1.2 meV/K · T − 0.3547 eV (for more details, see Fig. S1). It should be noted that breakdown of the ideal gas approximation at very high H 2 pressures will make our conversion between chemical potential and pressure increasingly inaccurate.

F. Short and long range order quantification
We quantify short-range order (SRO) with the Warren-Cowley order parameter, where A and B refer to the two species involved (Au and Pd or H and vacancy, depending on sublattice), n A-B is the total number of unlike bonds A-B in a certain shell (such as nearest neighbors), Z is the total number of bonds (regardless of decoration), and c A and c B is the overall concentration of A and B, respectively. To quantify SRO between the two sublattices, we generalize the above expression, where n Pd-H,Au-vac is the total number of pairs that are either Pd-H or Au-vacancy.
To quantify LRO on the Pd-Au sublattice, we use a standard expression for the structure factor, where the sum runs over all pairs i, j of sites in the sublattice (with a total of N sites), R i is the position of atom i, and f i = 1 if atom i is Au and −1 if it is Pd.

G. Phase diagram construction
Once mixing free energy curves F mix (c H , c Au , T ) have been calculated from MC simulations, phase diagrams can be constructed by identifying the regions where the free energy curve lies above its convex hull. (This construction relies on the assumption that the interface between the two phases is incoherent; hysteresis stemming from coherency will in general increase the apparent solubility [48,49].) The boundaries of these regions constitute the phase boundary, and to obtain a continuous phase boundary, we fitted, primarily as a guide to the eye, the observed boundary points to the function where a is a fitting parameter and L(c) is a fourth-order Redlich-Kister polynomial with five fitting parameters.
The phase diagram at full equilibrium is more intricate to construct since it involves integration of data with noise over a path in two dimensions. To this end, we averaged the integration over multiple paths and constructed the phase diagram based on several convex hull constructions. For details, see Supplementary Note 1.

A. Cluster expansion construction
A CE operates under the assumption that the atoms reside on a fixed lattice. In practice the effect of relaxation from the ideal lattice sites is, however, incorporated effectively through the fitting of the ECIs J α . For this treatment of relaxation to be feasible, the relaxed structure must not deviate too strongly from the ideal lattice [50]. An analysis of our DFT calculations revealed that in many structures, in particular those with a high Au content, relaxations from the ideal lattice sites were large, sometimes more than half the distance to the nearest neighbor on the other sublattice (Fig. S2a). To avoid training the CE on data that is not well represented by a lattice model, we excluded all structure in which at least one atom had relaxed more than 0.5Å from its ideal site, which is equivalent to approximately 25% of the distance to its nearest neighbor. Furthermore, the relaxed cell shape sometimes deviated from the cubic cell. In particular, pure AuH with 100% H is mechanically unstable and relaxes into a monoclinic structure (it should be noted that pure Au hydrides are considered unstable even under very high H 2 pressures [51,52]). To further investigate the crystal structure, we therefore also analyzed the shearing of the cell. Here, we define the shear strain tensor ∆A by subtracting the volumetric strain from the Biot strain tensor [53] (defined as U − I, where U is in turn defined by the polar decomposition of the deformation tensor, F = UP), and a scalar measure of the strain is obtained by taking the Frobenius norm of ∆A. On average, higher concentrations of Au and H tend to increase the shear strain (Fig. S2b).
Since structures with a very high Au content are not experimentally relevant except under fairly exotic conditions, we chose to exclude all structures with Au concentration above 50%. With the limitation to structures with no displacements larger than 0.5Å and no Au content above 50%, all structures with a shear strain ||∆A|| F > 0.4 were effectively excluded. In the end, we were left with 1,305 structures, out of which 1,255 were used for training the CE, with the remaining 50 used as test set to evaluate the performance of the final model.
The mixing energies span an energy range of approximately −100 meV/f.u. to 500 meV/f.u. (Fig. 1a). This wide range of energies is a challenge for the CE model, since some of the important phenomena occur on a much more narrow energy scale. Consider, for example, the two-phase region of pure Pd hydride. This two-phase region is the result of a concave mixing energy curve between c H = 0 and 0.67 (Fig. 1b), and the distance between the lowest energy structures and the convex hull is at most a few meV/f.u. (with the limited cell size in our training data). Meanwhile, if we increase the Au content, relevant variations occur on a scale that is at least one order of magnitude larger. To ensure that the important variations at low Au content are properly captured, we multiplied the cluster vectors and target energies with a weight w(c Au ) = 1 + 10(0.5 − c Au ) prior to fitting, chosen to make the variations in mixing energy occur within the same order of magnitude over the full range of Au concentrations. We note that non-uniform weighting of structures as well as the imposing of constraints to exactly reproduce some structures is done at the cost of a larger error for other structures, and inevitably a larger root mean square error over the full test set. Yet, it is more important to capture the relevant aspects of the energy landscape properly rather than having the lowest possible root mean square error.
The CE constructed in this fashion yields a root mean square error for the mixing energies of 8.3 meV/f.u. over the training set and of 15.3 meV/f.u. over the test set. To put these number in perspective, it is instructive to compare them to the wide span of mixing energies in this system (Fig. 1). The largest errors are frequently found among structures with a large content of Au, which partially explains the large difference in error between training and test set, as the latter contained a higher proportion of structures with high Au content. The predicted energy of all structures (including structures that relaxed too far relative to the reference structure (see Fig. S2) but excluding those with c Au > 0.5) has a good albeit not perfect correlation with the corresponding DFT energies (Fig. 1a, R 2 = 0.934). The energies of structures with long relaxation distances are sometimes severely overestimated. This is unsurprising since they were not included in the training data and cannot be expected to be well described by the CE in any case. Many of the poorly described structures are those with higher mixing energies in a series of configurations with the same composition ( Fig. 1b-g), and hence are in general less relevant for the thermodynamics of the system.
The ECIs of the constructed CE exhibit physically sound behavior in the sense that the strongest interactions are those that occur at a short distance, and triplet interactions are generally smaller than pair interactions (Fig. 2). Based on the ECIs, we can already at this point anticipate roughly how the system will behave. The strongest interaction is the one between Au/Pd and its nearest neighbor H/vacancy. We can anticipate that H will be attracted to Pd and vacancies to Au. H-H interactions, on the other hand, are comparatively weak, whereas the positive, short-ranged Au/Pd-Au/Pd ECIs indicate that the formation of Pd-Au pairs is energetically favorable.

B. Ordered phases in hydrogen-free Pd-Au
Although the focus of our work is on the thermodynamics of hydrogen in Pd-Au, we begin our analysis with the Pd-Au alloy in absence of hydrogen as ordering turns out to be also relevant with regard to the hydrogenated system, as discussed below. As our CE is only fitted to structures with 50% Au or less, we restrict our analysis to this composition interval. The mixing energies of structures with 12 Pd/Au atoms or less are negative, meaning mixing is favorable. It is, however, still possible that phase separation occurs between ordered (intermetallic) phases. LRO in Pd-Au has been reported in experiments on thin films [54,55] and nanoparticles [56], but has generally not been found in the bulk, suggesting that the energetics causing ordering is either enhanced by surface stress or has critical temperatures too low for ordering to occur readily in bulk samples. By enumerating structures with up to 12 atoms in the cell, our CE identifies three structures on the convex hull between 0 and 50% Au: the trivial case of pure Pd, AuPd 3 in the L1 2 configuration (Fig. 3a) and AuPd in the L1 0 configuration. These findings are largely consistent with other computational studies, although other structures have occasionally been found on the convex hull [57,58].
To investigate at below which temperature the L1 2 phase will appear, we ran MC simulations in the canonical ensemble, i.e., with the composition fixed at 25% Au in Pd, and tracked the LRO parameter (see Eq. 10) as a function of temperature. The results revealed a critical temperatures at 180 K (Fig. S6), with a span among the sampled CEs from less than 150 K to 260 K.

C. Hydrogen in pure Pd
We are now ready to add hydrogen to the system. Since the pure Pd hydride (no Au) has been more thoroughly studied than the alloy hydride, we begin in this limit. We carried out VCSGC-MC simulations on the H/vacancy lattice with temperatures between 200 K and 700 K in steps of 50 K. The thus obtained isotherms (Fig. 4a) indicate a low solubility of H in Pd until a "threshold pressure" has been reached. For sufficiently low temperatures, there is then a phase transition from a hydrogen-poor α phase to a hydrogen-rich β phase, which can be detected by integrating the free energy derivative and identify the regions where the free energy curve lies above its convex hull (as outlined in Sect. II G). The corresponding phase diagram is largely consistent with the experimental phase diagram constructed by Wicke and Blaurock [59] (Fig. 4b). The critical temperature, which is challenging to estimate accurately both experimentally and by means of MC simulations, differs approximately 140 K (560 K experimentally and 420 K based on our CE). Although this is a fairly large discrepancy, it should be noted that small variations in the ECIs have a large impact on the critical temperature; at temperatures close to the critical temperature, uncertainties in the phase boundary become large (quantified by ten sampled CEs; the spread from these ten CEs is indicated with red horizontal bars in Fig. 4b). Discrepancies between simulation and experiment are thus not unexpected, especially given that it is also challenging to construct the phase diagram experimentally.
D. Hydrogen loading of fixed Pd-Au lattice: random and para-equilibrium We are now ready to include Au in the MC simulations. To this end, we note that the two sublattices are markedly different from a kinetic perspective. Hydrogen diffuses easily even at low temperatures, while Pd and Au diffuse over a much longer timescale. To reach full equilibrium, in which the Pd and Au atoms rearrange in response to a change in hydrogen environment, a very slow experiment with H 2 pressure kept constant over extended periods of time, is required. If the H 2 pressure is more quickly increased/decreased, we may instead assume that the Pd-Au sublattice is frozen and that a metastable equilibrium is reached by rearrangement of the hydrogen atoms only. The chemical order on the Pd-Au sublattice would then typically be dictated by the conditions during fabrication of the alloy, which usually involve a form of thermal annealing. Here, we distinguish two extremes of such a metastable equilibrium: para-equilibrium, in which the Pd-Au ordering is equilibrated at 300 K in the absence of H 2 , and random equilibrium, in which the Pd-Au lattice is randomized, corresponding to quenching of the system after equilibration at a very high temperature in absence of H 2 .
We can model these situations by only carrying out VCSGC-MC flips on the H/vacancy sublattice. For the random equilibrium, we commence from a simulation cell in which the Pd and Au atoms have been randomly distributed, whereas for para-equilibrium, we run a canonical MC simulation with only Pd and Au and randomly pick a configuration as the fixed Pd-Au lattice. To suppress spurious effects from the particular choice of Pd-Au lattice, we averaged our results over five instances of both cases.
These MC simulations reveal weak SRO (Fig. 5). We first note that the equilibration of the Pd-Au lattice in absence of hydrogen yields a weak propensity for the system to form unlike (Pd-Au) nearest neighbor bonds, whereas bonds that are alike (Au-Au and Pd-Pd) are weakly favored among nextnearest neighbors ( Fig. 5a; in random Au-Pd the corresponding order parameters are zero by construction). For nearest neighbor Pd/Au-H/vacancy pairs (Fig. 5b-c), the SRO parameter is negative, indicating that H is more likely to occupy sites next to Pd than next to Au. This confirms what has been observed both with Mössbauer spectroscopy [60] and firstprinciples calculations [24], namely that it is energetically unfavorable for hydrogen to occupy the sites closest to a Au atom. For next-nearest Pd/Au-H/vacancy neighbors (Fig. 5de), on the other hand, the effect is reversed, albeit weaker. The change of sign between first and second nearest-neighbor was observed also by Sonwane et al. [61] in a model based on DFT calculations. It is furthermore consistent with the common argument that chemistry causes repulsion between neighboring Au and H, whereas dilation of the lattice by a Au atom gives rise to a more long-ranged elastic attraction [20].
The most pronounced difference between the two types of equilibria is found for the most short-ranged H/vacancy pair ( Fig. 5f-g), for which the system in para-equilibrium has a stronger tendency to form H-H pairs. For next-nearest neighbor H/vacancy pairs (Fig. 5h-i)   Compared to experimental results [59] (orange line in b), the critical temperature is underestimated, but the uncertainty from the choice of CE is significant, as indicated by the spread of the miscibility gap at 300 and 400 K from ten sampled CEs (red lines indicate difference between minimum and maximum value predicted by these ten CEs Given the quantitative difference in chemical ordering between random Au-Pd and para-equilibrium, we may now look for consequences for their respective thermodynamics. It turns out, however, that their H 2 pressure-composition isotherms (Fig. 6a-b) are very similar. It is only at fairly high Au concentrations, approaching 25% Au, that the isotherms of random and para-equilibrium exhibit discernible differences. At high Au concentrations and high H 2 pressures, the hydrogen uptake is slightly higher in random equilibrium than in is negative, which, with our definitions (see Eq. (9)), means that H tends to sit close to Pd and vacancies close to Au. In the second Au/Pd shell (c-d), a weak reverse trend is seen throughout most of the H concentration range. In the first H/vacancy shell (f-g), H-H pairs are favored, only weakly in random equilibrium but stronger in para-equilibrium. The behavior is more complex in the second H/vacancy shell, with ordering tendencies being dependent on both H and Au concentrations.
para-equilibrium (Fig. 6d). Although this difference is small, there seems to be a weak but consistent trend; when SRO emerges due to lower annealing temperature, the hydrogen uptake at 300 K goes down at pressures above approximately 10 mbar (Fig. S7). The difference seems to be the largest at a hydrogen pressure around 1 bar, where the fully random alloy with 25% Au absorbs almost 4 percentage points more hydrogen than the one equilibrated in 300 K. These results seem to be in agreement with Chandrasekhar and Sholl [30]. To summarize, the SRO that emerges at low temperatures in Pd-Au makes the material absorb slightly less hydrogen at high H 2 pressures. The discontinuity in the Au-poor isotherms, which is the hallmark of the phase transition from α to β, disappears quickly when the Au content is increased and it does so in almost exactly the same way for random and para-equilibrium. Consequently, the respective phase diagrams are essentially identical (Fig. 6c). Our results predict a critical Au concentration for the α + β two-phase region around 8-9%. This is significantly lower than experimental measurements, which yield estimates that vary from 10-15% (at 303 K) [11] to around 17% (at 298 K) [10]. The prediction of the solvus line is, however, sensitive to small variations in the ECIs (quantified by the spread obtained by sampling ten CEs and indicated by red horizontal bars in Fig. 6c). Thus both our prediction and the experimental value are associated with significant uncertainties.

E. Hydrogen loading in full equilibrium
We may now ask what will happen if we allow the Pd-Au lattice to rearrange as we expose it to hydrogen. This situation is commonly referred to as full (or complete) equilibrium. It should be noted that this is an idealization that is very time-consuming to achieve in practice, as in most experiments the H 2 pressure is not maintained long enough to allow for Pd and Au to diffuse to a sufficient extent. Yet, it is important as it provides the thermodynamic driving force for the changes that do occur, although they may only rarely take the system to full equilibrium. To investigate the full equilibrium, we carried out MC simulations in the VCSGC ensemble on both sublattices. Subsequently we obtained the free energy by integration across the concentration plane using many different integration paths, and sampled the convex hull based on these different integration paths (for details, see Supplementary Note 1). Thereby we obtained a heat map of the probability that certain compositions are on the convex hull, which can be interpreted as (an isothermal cut of) the phase diagram ( Fig. 7; see Fig. S10 for data at 250 and 500 K). The two-phase region at low Au content observed in para and random equilibrium is present in the full phase diagram at 300 K as well, but is accompanied by a much larger multi-phase region at higher Au concentrations (Fig. 7a). At 400 K, the former two-phase region is almost gone, but the multi-phase region at higher Au contents is still present. This new multi-phase region that appears in the full phase diagram, as opposed to random and para-equilibrium, is largely driven FIG. 6. Pressure-composition isotherms at 300 K for Pd-Au hydride with (a) a random Pd-Au sublattice and (b) in para-equilibrium, i.e., with the Pd-Au lattice equilibrated at 300 K in absence of H. Phase transitions from the H-poor α phase to the H-rich β phase are indicated with a dashed line. The corresponding two phase diagrams (c) are almost identical, with the two-phase region closing at approximately 8-9% Au. Red lines at 6.3% Au indicate the range of values obtained from ten sampled CEs. By studying the hydrogen uptake at fixed H2 pressure as a function of Au concentration (d), it is clear that differences between random and para-equilibrium are found primarily at high H2 pressures, where at the same pressure the random system absorbs a slightly larger amount of hydrogen than the system in para-equilibrium. Larger variations are also seen at 10 −2 bar (orange line), since this it close to the plateau pressure. The phase diagram in para-equilibrium is indicated with transparent grey in (d).
by a particularly stable composition interval around 25% Au and approximately 10-30% H, where at 300 K isobars ranging from 10 −3 bar to 10 2 bar all converge. A closer look at this composition interval reveals that LRO develops here (Fig. S9). Specifically, the Pd-Au sublattice orders in the L1 2 phase, in which the atoms are arranged such that 25% of the hydrogen sites are octahedral "cages" with only Pd nearest neighbors (Fig. 3b). As has been previously reported [24] and can be expected from the ECIs (Fig. 2), occupation of hydrogen at such sites is energetically particularly favorable, and renders this long-range ordered phase particularly stable. At 500 K (Fig. S10b), LRO no longer forms and the phase diagram is largely featureless, i.e., a solid solution can be expected at any concentration (except possibly at high Au content and extremely high H 2 pressures).
The phase diagram in Fig. 7a thus predicts that if a Pd-Au alloy at 300 K with, say, approximately 15% Au is subjected to a H 2 pressure between approximately 10 −3 bar and 10 2 bar, phase separation will occur. Of the resulting two phases, one will be the L1 2 phase with between 10 and 30% H. The character of the other phase depends on the H 2 pressure; if the pressure is below 10 −2 bar, the second phase will be the dilute α phase (with Au content strongly dependent on H 2 pressure) whereas if the pressure is above approximately 10 −2 bar, the second phase will be the hydrogen-rich β phase with approximately 5% Au.
The existence of a long-range ordered phase in hydrogentreated Pd-Au has been reported experimentally by Lee et al. [29], who observed a super-lattice phase in Pd-Au with 19% Au. The structure of this phase could not be determined but the authors reported it to be more complex than L1 2 . It should perhaps not be ruled out, however, that the off-stoichiometric composition as well as defects such as anti-phase boundaries, which occurred frequently in some of our simulations, may have played a role in these experiments. Fig. 7 is the single-phase region that the 10 3 bar isobar traces out (pink line in Fig. 7). Here, the SRO parameter for nearest neighbor Pd-Au pairs switches sign, indicating a transition from an excess of Pd-Au pairs to an excess of Au-Au and Pd-Pd pairs (Fig. S8). Inspection of the MC trajectories at these composition reveals that the system has started to cluster into hydrogen-poor Au clusters and hydrogen-rich Pd clusters, both only a few atoms large. Our model consider this phase thermodynamically stable and not merely a first step towards full phase separation. It seems plausible that a phase like this may provide a favorable balance between the short-ranged repulsion between Au and H and the long-ranged elastic attraction caused by dilation of the Pd lattice by adjacent Au clusters. It is not surprising that this phase has not been reported experimentally given the very high H 2 pressure required for it to form ( 10 3 bar). It should be noted that the pressure specified in this region is highly approximate, as this is outside the applicability range of the ideal gas law.

F. Impact of annealing conditions on hydrogen solubility
Although full equilibrium will usually not be reached in the time frame of a typical hydrogen loading experiment, there may be circumstances in which full equilibrium is approached. For example, annealing of Pd-Au is sometimes done in the presence of hydrogen in order to prevent oxidation. Although the pressure of H 2 is then typically too low to have a significant impact, the phase diagram Fig. 7 shows that complex behavior may emerge in the presence of hydrogen, especially if the temperature is not too high. We may consider a fairly typical situation in which the alloy is annealed at a particular temperature and H 2 pressure, after which hydrogen absorption/desorption isotherms are measured at another temperature. Since the latter absorption/desorption is usually carried out at a much lower temperature than the annealing, and during a much shorter period of time, it is reasonable to assume that the Pd-Au sublattice gets frozen in during anneal-ing and does not change when measuring the isotherm. The chemical ordering on the Pd-Au sublattice would then be determined entirely by the conditions during annealing.
We now mimic such an experiment by simulating isotherms at 300 K in Pd-Au with 25% of Au annealed in different conditions. To this end, we first run MC simulations at a specified annealing temperature with canonical MC swaps on the Pd-Au sublattice, and SGC MC flips on the H-vacancy sublattice, using a chemical potential corresponding to a fixed H 2 pressure. We then pick five random snapshots from the resulting trajectory, remove the hydrogen, fix the Pd-Au sublattice, and run MC simulations at 300 K with SGC flips on the H-vacancy sublattice only, using a wide range of hydrogen chemical potentials.
Inspection of 300 K isotherms with Pd-Au annealed in 400 K (Fig. 8a) reveals that the hydrogen concentration depends strongly on the conditions during annealing. When annealed in pressures of 10 −1 bar or lower (blue and orange lines), the isotherms behave as in the random or paraequilibrium case, with an almost linear isotherm (when plotted on a logarithmic scale). For higher pressures, however, the isotherms are markedly different. After annealing in 1, 10 or 100 bar, the uptake of hydrogen at low pressures is much higher, after which the concentration of hydrogen stays at about 25% up to very high pressures, meaning that at sufficiently high pressures, the hydrogen content is in fact higher in Pd-Au samples that were annealed in the absence of hydrogen. The difference between these two kinds of isotherms is that annealing in 1-100 bar induces L1 2 LRO. If the H 2 pressure is raised to 1,000 bar, the ordered phase no longer forms, and the isotherm (brown line in Fig. 8a) becomes more similar to the low-pressure isotherms, although the hydrogen uptake is significantly higher if the H 2 pressure is above a few millibar.
The ordered L1 2 phase is thus clearly distinguishable from the ones lacking LRO already from the isotherms. These findings are consistent with Lee et al. [29]. When annealing in 600 K, on the other hand, the temperature is too high for any LRO to emerge, and the 600 K isotherms are essentially identical regardless of annealing pressure (Fig. 8b). It is expected that much higher pressures are required to impact the system at 600 K compared to 400 K, because at constant pressure when the temperature goes up, the hydrogen content in the material goes down. Nevertheless, when the annealing pressure is 100 bar, the hydrogen content in the system is about 14% during annealing, but the 300 K isotherm is still virtually unaffected. Only when the annealing pressure reaches 1,000 bar (leading to approximately 26% hydrogen in the system during annealing), is the 300 K isotherm clearly distinguishable from the isotherm of Pd-Au annealed in vacuum, and even then the difference is small.
This picture emerges more comprehensively if we study the hydrogen absorbed at 300 K and specific partial pressures (Fig. S11): the content of hydrogen is virtually independent of annealing conditions as long as the phase transition to the L1 2 phase does not occur. Very high annealing pressures, on the order of 10-1000 bar, are required to achieve a significant impact on the content of absorbed hydrogen unless the ordered L1 2 phase is formed.
It is worth stressing that while hydrogen uptake at moderate pressure is enhanced by LRO formation, we observe a different trend in Sect. III D; SRO formation decreases the hydrogen uptake unless the pressure is very low. It is thus advisable not to speak in too general terms about the impact of chemical order on hydrogen solubility, because it depends on the details of the chemical order as well as the H 2 pressure at which the solubility is assessed.
Given that the chemical ordering of the L1 2 phase thus has a significant impact on the nature of absorption of H 2 in Pd-Au, it may have a notable impact on any utilization thereof. We therefore estimate the conditions during which this long-range ordered phase will form (Fig. 8c). Our CE indicates that the ordered structure starts to form below approximately 500 K with a H 2 pressure between approximately 50 and 100 bar. The H 2 pressure range where LRO forms then widens quickly as temperature is decreased, but slower kinetics will of course inhibit order formation at too low temperatures. Lee et al. [29], who studied Pd-Au with 19% Au, assumed that distinct isotherms are a fingerprint of LRO, and while the purpose of their study was not to map out the conditions under which order emerges, their observations are qualitatively consistent with our CE (red crosses in Fig. 8c), but with a higher critical temperature (with LRO persisting up to at least 598 K). It thus seems likely that our CE underestimates the critical temperature by at least 100 K. By repeating the same calculations with 10 different CEs (grey areas in Fig. 8c), we observe qualitatively the same behavior but find that both the critical temperature and the pressure range are very sensitive to small variations in the ECIs. Although none of the sampled CEs exhibit a critical temperature above 600 K, it should be clear that small errors not captured by the CE approach (such as neglect of vibrations or the choice of exchange-correlation functional used to calculate the training data) may cause an error of this magnitude.

IV. CONCLUSIONS
We have comprehensively investigated the impact of chemical order on hydrogenation of Pd-Au alloys using a computational approach based on DFT calculations, CE models, and MC simulations. Although relatively large relaxations from octahedral hydrogen sites hamper the ability of the CE models to exactly reproduce formation energies calculated with DFT, we found that our CE reproduced thermodynamic properties of Pd-Au hydrides that are well-established experimentally, with a quantitative agreement on par with what can be expected from a CE approach. This applies especially given the large impact of small changes in the ECIs, which is likely always inherent in this approach or indeed any other approach to derive phase diagrams from interatomic potentials or firstprinciples data. Also, we acknowledge that zero-point energy as well as phonons can sometimes play an important role in hydrides due to the low mass of hydrogen, but inclusion of these effects would have been computationally prohibitive. The general impact of these effects may warrant future studies with 25% Au at 300 K. After annealing in 400 K (a), the isotherms exhibit a markedly different behavior if the H2 pressure during annealing was between 1 and 100 bar (the red 10 bar isotherm is hidden under the green 1 bar isotherm). After annealing in 600 K (b), the isotherms are almost identical regardless of H2 exposure during annealing. The origin of the change in isotherm in (a) is the formation of LRO under certain conditions as quantified in (c). The CE predicts that LRO will form in full equilibrium when Pd-Au with 25% Au is subjected to H2 pressure and temperature corresponding to the area between the black lines. The corresponding area predicted with ten sampled CEs is indicated with transparent blue, one per CE. Darker color thus means more CEs predict order formation at that point. Out of these ten CEs, two have a critical temperature too low to be visible in this figure. The conditions investigated by Lee et al. [29] with signs of order formation in Pd-Au with 19% Au are indicated with red crosses. Colored squares and triangles indicate the annealing conditions for isotherms in (a) and (b), respectively. and may explain some of the quantitative disagreement with experiment.
Our results provide a rationale for the experimental observation that absorption/desorption isotherms are relatively stable over time; as long as the L1 2 phase does not form, isotherms will stay similar even if chemical order changes. This is manifested by the similarity between isotherms in random and para-equilibrium. Our results predict that a small reduction in the ability to absorb hydrogen may be observed if Pd-Au annealed at a very high temperature is allowed to reach equilibrium at a much lower temperature.
Under long-term exposure to hydrogen, however, the situation changes substantially. The emergence of the ordered L1 2 phase introduces complexity in the phase diagram, and Pd-Au stored at room temperature and, say, 1 mbar H 2 may over time exhibit very different absorption isotherms, as the result of emerging LRO. The hydrogen content in air is of course orders of magnitude lower than what we predict is required for this phase to form, but repeated exposure to high H 2 pressures may result in formation of LRO and an altered isotherm.
Does the LRO vanish when the hydrogen is removed? This question is difficult to answer, because from a kinetic point of view, our simulation approach is not able to describe the stability of the ordered phase in absence of hydrogen (or vice versa). We do predict that the L1 2 phase is the ground state at this composition also in absence of hydrogen, but its critical temperature is well below 300 K. It may, however, be noted that there are (hydrogen-free) Pd-Au phase diagrams in the literature where this phase is expected to form above room temperature [54,55,58], that is, even without exposure to hydrogen. Although the existence of this phase, let alone its critical temperature, remains debated, it seems reasonable to assume that the L1 2 phase might be fairly stable at room temperature even after hydrogen is removed.
The emergence of a different isotherm upon ordering might constitute a challenge as the sensor would have a different readout at the same H 2 pressure. It may, however, also present an opportunity. The ordered phase absorbs significantly more hydrogen at low pressures. A sensor consisting of the ordered phase is thus likely (depending on the nature of the readout) to be significantly more sensitive to small changes in H 2 pressure. If the phase is sufficiently stable, this may provide a relatively simple way to improve a Pd-Au hydrogen sensorannealing in hydrogen to boost sensitivity.