Hydrogen induced vacancy clustering and void formation mechanisms at grain boundaries in palladium

Hydrogen has a significant impact on the formation of vacancies, clusters and voids in palladium and other metals. The formation of vacancy-hydrogen complexes in bulk Pd and at 3 and 5 grain boundaries was investigated using first-principles calculations and thermodynamic models. Equilibrium vacancy and cluster concentrations were evaluated as a function of temperature and hydrogen partial pressure based on the Gibbs energy of formation including vibrational and configurational entropies. Vacancies were found to be significantly stabilized by association with interstitial hydrogen, leading to enhanced concentrations by several orders of magnitude. Vacancy clusters were further stabilized at grain boundaries, with equilibrium concentrations reaching site saturation for clusters comprising up to three vacancies. Nanovoids were investigated based on Wulff constructions from calculated surface energies of the (0 0 1) and (1 1 1) terminations as a function of temperature and coverage of hydrogen adsorbates. The most stable termination changed from (1 1 1) in vacuum to (0 0 1) in H 2 , and the surface energies were lowered due to hydrogen adsorbates. Consequently, voids were also stabilized in the presence of hydrogen. Coalescing of vacancies into nanovoids was found to be thermodynamically unfavorable due to the loss of configurational entropy. It was therefore concluded that enhanced concentrations of vacancies and clusters does not directly cause the formation of voids. The formation of voids in Pd-based membranes was discussed in terms of microstructural characteristics, and strain due to chemical expansion and plastic deformation. © 2020 Acta Materialia Inc. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
The role of hydrogen on the structural and functional properties of metallic materials has long been a subject of interest for a broad range of applications. Metals are in many cases directly exposed to hydrogen-containing environments, for instance when used as sensors [1] , hydrogen-selective membranes [2] , catalysts [3] and hydrogen plasma facing materials in nuclear fusion reactors [4] . Moreover, hydrogen is a typical product of corrosion and otherwise omnipresent as the most abundant element in the known universe.
Exposure to hydrogen can lead to structural damage and deterioration of the mechanical integrity of the material. Depending on the specific metal and environmental conditions, the degradation can manifest itself macroscopically in the form of embrittlement [5 , 6] , cracking [7 , 8] , blistering and pore formation [9 , 10] . Combined with a range of microstructural changes, one of the primary mech-anisms behind these degradation processes is the enhanced concentration of metal vacancies in the form of complexes with interstitial atomic hydrogen. This phenomenon -referred to as superabundant vacancies -was first reported by Fukai et al. for Pd and Ni [11 , 12] , and later for a range of other metals and related alloys [13] .
The vacancy complexes can cluster into nanovoids or cavities, eventually becoming large enough to contain molecular and gaseous H 2 [14 , 15] . Formation of larger clusters by coalescence of m monovacancies and n interstitial hydrogen can be described by the reaction In some cases, the H 2 pressure can build-up inside the voids and act as a driving force for further growth by plastic deformation of the surrounding lattice [16] . In relation to the possibility of forming such pressurized bubbles, there is a crucial difference between the behavior of different metals, which can be distinguished by the Gibbs energy for dissolution of interstitial hydrogen relative to H 2 gas. For a system in equilibrium, the chemical potential of hydrogen should be constant throughout, and the pressure of con-fined H 2 can therefore only increase up to a certain level determined by the equilibrium with hydrogen dissolved in the lattice. Given the logarithmic dependence of the chemical potential on pressure, i.e., μ H ∝ kT ln ( p H 2 ) , bubbles with pressures in the MPa range and higher can be achieved. Still, hydrogen must in these cases be introduced under non-equilibrium conditions, such as by electrochemical loading [17] , ion implantation [18] or exposure to hydrogen plasma [19] .
In the case of palladium and other metals that display high solubility and diffusivity of atomic hydrogen, the pressure of H 2 inside voids cannot exceed the hydrogen partial pressure in the surrounding atmosphere except when the system is pushed far out of equilibrium at very short timescales, e.g., nanosecond laser pulse heating of NbD 0. 7 and VH x [20 , 21] . Otherwise, the chemical potential of hydrogen is allowed to equilibrate and this process is not kinetically hindered in palladium due to the fast kinetics of adsorption/desorption. In this context, the formation of bubbles in palladium containing the tritium isotope of hydrogen may serve as an instructive example. Tritium is highly soluble in palladium due to its chemical similarity with hydrogen, but decays naturally over time to 3 He, which exhibits very low solubility and diffusivity. Association between 3 He and vacancies can be anticipated and further nucleation to bubbles with increasing He pressure eventually results in cracks through which the gas can escape [22] . This process thereby represents a similar bubble formation mechanism as in supersaturated metals that exhibit relatively low diffusivity of hydrogen. Another related bubble formation process has been reported for N-doped ZnO and ZnO-GaN semiconductors due to evolution of N 2 bubbles after ion implantation [23] and annealing [24] , respectively. In these cases, the difference in solubility arises from the distinct oxidation state of nitride anions in the host lattice and molecular nitrogen.
In particular for Pd-based H 2 separation membranes, the formation of voids and pinholes during operation over long time periods can result in reduced selectivity and lifetime [25][26][27][28] . For instance, cavities and pinholes of up to tens of nanometers have been observed in conjunction with increased trans-membrane N 2 leakage under 1-5 bar hydrogen at 300-450 °C [29 , 30] . Notably, the voids mainly formed along grain boundaries and at triple junctions in the films. In understanding the mechanisms of void formation, it is important to also consider the role of chemical expansion due to dissolution of hydrogen and transition to hydride phases; The evolution of internal stress in Pd films subjected to hydriding cycles can result in increased dislocation densities and irreversible microstructural changes [31 , 32] . Higher dimensional defects such as grain boundaries therefore appear to play an important role in void formation, possibly as nucleation sites for vacancies to coalesce.
Ab initio calculations have been utilized to gain atomistic insight into the formation of vacancy-hydrogen complexes of different size in several systems. Nazarov et al. performed a comprehensive study of hydrogen interacting with mono-and divacancies in twelve different fcc metals including Pd; the results consistently showed energetically favorable trapping of multiple hydrogen interstitials in the vacancies [33] . Vekilova et al. [34] reported similar findings for vacancy-hydrogen complexes in Pd. You et al. considered six fcc and six bcc metals and concluded that vacancy formation was particularly favorable in V, Nb, Ta, Pd and Ni due to low formation energies of the complexes which further decreased with the number of hydrogen [35] . Several similar studies have been undertaken on the formation of mono-and divacancy clusters and the interaction with hydrogen in W [36 , 37] , Mo [38] , Ni [39 , 40] , and Cu [41] . Geng et al. investigated nanovoids of up to 27 vacancies in α-Fe and modelled the pressure of molecular hydrogen after the inner surfaces were covered with hydrogen adatoms [42] . Hou et al. recently investigated bubble formation in W due to self-clustering of interstitial hydrogen at a critical concentration and developed predictive models for hydrogen trapping including H-H interaction parameters on nanovoid surfaces [43 , 44] .
Comparatively few studies have considered the role of higher dimensional defects such as grain boundaries and dislocations on the formation and clustering of vacancy-hydrogen complexes. Zhou et al. reported favorable segregation of vacancy-hydrogen complexes to the highly coherent 3 twin boundary in Ni alloys, and the complexes were further stabilized by macroscopic plastic deformation [45] . Favorable segregation of vacancy-hydrogen complexes was also reported for the 3 boundary in α-Fe [46] . In the presence of edge dislocations in α-Fe, annihilation of vacancies and complexes was found to be favored over clustering that could lead to nanovoids [47] .
In the present contribution, we further investigate vacancy and void formation mechanisms in the bulk and along grain boundaries of palladium in the presence of hydrogen by means of density functional theory (DFT) calculations. The Gibbs energy of formation of vacancy-hydrogen complexes and clusters with up to four vacancies ( m = 1-4), were considered in the bulk as well as in the proximity of the 3 (1 1 1) and 5 (2 1 0) grain boundaries, which are among the most abundant small-angle boundaries in Pd [48] . Equilibrium concentrations were thereby calculated as a function of temperature and hydrogen partial pressure. However, a vacancy-hydrogen complex and a nanovoid can be considered to be quite different from a thermodynamic perspective; while formation of a certain amount of vacancies at finite temperatures results from the gain in configurational entropy of the crystal, a nanovoid represents an increase in surface area and is thermodynamically unstable relative to the single crystal. Therefore, the energetics of void formation was also evaluated based on the surface energies of the inner surfaces of nanovoid in the presence of hydrogen.

Density functional theory calculations
DFT calculations were performed using the projector augmented wave method implemented in the Vienna ab initio simulation package (VASP) and the generalized gradient approximation due to Perdew, Burke and Ernzerhof (GGA-PBE) [49][50][51] . The calculations were carried out with a plane wave cut-off energy of 500 eV and an electronic convergence criterion of 10 -5 eV. Acentered 8 × 8 × 8 Monkhorst-Pack k-point grid was used for the cubic F m 3 m unit cell [52] . The atomic positions were relaxed until the residual forces were within 0.01 eV Å -1 (0.02 eV Å -1 for cells containing defect clusters at grain boundaries). The optimized bulk lattice parameter a 0 was 3.944 Å .
Defect clusters in the bulk were studied using 4 × 4 × 4 supercells (256 atoms) and a corresponding 2 × 2 × 2 k-point grid. In addition to isolated metal monovacancies, vacancy clusters v m with m = 2-6 nearest neighbor vacancies were introduced in linear, trigonal planar, tetrahedral and octahedral configurations. Hydrogen was successively introduced on tetrahedral or octahedral sites as far apart as possible in the clusters, forming v m H n clusters ( Fig. 1 ). Adjacent to metal vacancies, these sites become symmetrically equivalent to 3-fold hcp sites on the (1 1 1) surface and 4-fold sites on the (0 0 1) surface, respectively. The concentration dependent dissolution enthalpy of hydrogen into bulk Pd and the associated chemical expansion was calculated for filling of 1-16 of the 32 octahedral sites in a 2 × 2 × 2 supercell including relaxation of the lattice parameter.
The vibrational frequencies of Pd and hydrogen interstitials in bulk and associated with a vacancy were calculated using the finite displacement method. Before these calculations, the vibrating species were further relaxed to an accuracy of 10 -4 eV Å -1 . The 3 (1 1 1) and 5 (2 1 0) coincidence site lattice boundaries were constructed with the boundary planes separated by 5 (1 1 1) layers (13.7 Å ) and 14 (2 1 0) layers (13.2 Å ), respectively. Due to the periodic boundary conditions, the cells contained two equivalent grain boundaries. The super cell size was relaxed in steps of 0.1 Å in the direction perpendicular to the grain boundary ( zaxis), which resulted in an optimized grain separation that was unchanged for 3 and extended by 0.4 Å per boundary for 5. Defect clusters were considered in supercells with an interface area of 16.7 Å × 19.3 Å for 3 (576 atoms) and 15.8 Å × 17.6 Å for 5 (480 atoms) and the corresponding k-point grids were 2 × 2 × 1.
Surface energies were calculated for the (1 1 1) and (0 0 1) terminations and the contribution from atomic hydrogen adsorbates on 3-fold and 4-fold sites, respectively, were included. The calculations were performed using slabs with a surface unit cell of √ 2 a 0 × √ 2 a 0 and thickness of 7 atomic layers (28 atoms) and 6 × 6 × 1 k-point sampling.

Thermodynamic models 2.2.1. Formation of vacancy clusters
The formation energy of a v m H n cluster can be expressed as where E tot v m H n and E tot bulk are the total electronic energies of the defective and bulk supercells as calculated by VASP, respectively, k is Boltzmann's constant, T is the temperature and p H 2 is the partial pressure of hydrogen. The total energies of palladium and hydrogen were defined from the bulk structure (per Pd atom) and H 2 molecule, respectively.
The vibrational entropy of the v m H n clusters and bulk Pd was calculated from the vibrational frequencies of the atomic species according to [53] where S f v was calculated as the difference in vibrational entropy between that of a Pd atom adjacent to avacancy and that of a bulk Pd atom, and l was the number of Pd atoms adjacent to the vacancy cluster (Table S1). Since bulk Pd was used as reference state in Eq. (1) , there was no change in vibrational entropy for the Pd atoms that were removed upon cluster formation. S f H was calculated as the difference between the vibrational entropy of a 3-fold or 4-fold H atom in a v 1 H 1 cluster and the entropy of molecular H 2 which was taken from thermochemical tables [54] . Additional defect interactions and changes in the vibrational frequencies with increased cluster size and hydrogen content were not considered.
The vibrational frequencies were also used to calculate the zero-point energy (ZPE) of v m H n , bulk Pd and molecular H 2 according to The formation enthalpy of v m H n clusters was calculated according to where ZPE was calculated in the same way as for the entropies.
The change in ZPE of Pd upon vacancy formation in the grain boundaries was assumed to be the same as in the bulk, although one could expect some differences for the 5 boundary due to the more open structure. The binding energy of the v m clusters was calculated relative to the isolated vacancies including ZPE contributions according to The Gibbs free energy of formation of the v m H n clusters was calculated as where S conf v m H n = k ln ( n ) is the configurational entropy of n hydrogen atoms distributed over the total number of sites in the vacancy cluster, N , according to For simplicity, and due to the close proximity between hydrogen on adjacent 3-fold and 4-fold sites, the clusters were only filled with hydrogen on either 3-fold or 4-fold. The number of hydrogen sites N are listed in Table S1. The equilibrium concentration of v m H n clusters was calculated according to [55] c where z is the number of configurations of the v m cluster (Table  S1). The concentration (site fraction) of v m H n clusters was thereby limited to the number of Pd sites in the bulk or grain boundaries.

Surface energies
For adsorption of hydrogen onto a surface site * , i.e., 1 2 H 2 + * = * H , the equilibrium hydrogen coverage, H , can be expressed from the equilibrium constant, K H , according to where H ads H is the coverage dependent adsorption enthalpy including ZPE contributions, calculated by VASP. The vibrational part of the adsorption entropy and ZPE were obtained in the same manner as for the v m H n clusters, while the change in vibrational properties of Pd upon adsorption could be neglected due to the low mass of hydrogen. The configurational entropy of the adsorbates was calculated as Equilibrium hydrogen coverages were obtained by numerical solution of Eqs. (11)(12)(13) where E tot slab is the total energy of the slab unit cell, r is the ratio between the number of Pd atoms in the slab and bulk cells, A is the area of the surface, and hkl are the Miller indices of the specific surface. The calculated surface energies can be used to estimate the energy required to generate a void of a certain size and shape. Since the inner walls of a void are equivalent to the surfaces of a particle, the equilibrium shape of the void can be found by generating particles based on Wulff constructions [56] . A minimum energy polyhedron satisfies the following criterion where d hkl is the distance to a specific surface from the polyhedron center, and the c W parameter determines the size of the polyhedron. The pressure of molecular hydrogen in the voids was the same as in the surrounding atmosphere as defined from the constant chemical potential of hydrogen throughout the system.

Vacancy clusters in bulk
The different configurations of vacancy clusters are shown in Fig. 2 . The corresponding binding energies showed a favorable −0.04 eV for v 2 and up to −0.14 eV per vacancy for the v 6 cluster.
The v 1 H n , v 2 H n , and v 4 H n clusters were considered further due to their relatively high number of 4-fold sites per vacancy (Table S1).
These clusters are shown in Fig. 3 with the 3-fold and 4-fold sites saturated with hydrogen. The 3-fold hydrogen sites are more numerous and more closely spaced within the cluster than the 4-fold ones. The 4-fold hydrogen sites reside within the atomic planes of Pd surrounding the vacancy cluster and the extent of unfavorable H-H interactions at higher concentrations can thus be expected to be lower. These considerations are supported by the calculated formation enthalpies of the clusters. For the v m H n clusters with 3-fold hydrogen, there was a notable shift in the enthalpy curves' slope around n = 4-6, which is where more closely spaced sites are filled ( Fig. 4 a). In contrast, the formation enthalpy of the clusters with 4-fold hydrogen exhibited an essentially linear behavior as hydrogen was added. Overall, the formation enthalpies increased substantially with cluster size due to the cost of forming Pd-vacancies, while this was partly alleviated by the increased number of hydrogen sites due to the exothermic contribution of adding hydrogen. The formation enthalpies were generally more exothermic for the 4-fold configurations, especially for the higher hydrogen contents. The formation enthalpies of the v m H n clusters not explicitly calculated were taken by linear interpolation between adjacent points, i.e., lines in Fig. 4 . The ZPE contributions to the formation enthalpies of selected v m H n clusters are listed in Table 1 . Notably, ZPE was exothermic for 4-fold hydrogen and endothermic for 3-fold hydrogen, which demonstrates the importance of including ZPE contributions in determining dissolution and adsorption enthalpies for hydrogen. As with the formation enthalpies, the change in ZPE with hydrogen content was quite significant for the monovacancy clusters with 3-fold hydrogen. In the case of clusters with more vacancies, we assume that ZPE does not change noticeably, since it is essentially the same for v 2 H 1 and v 1 H 1 . Furthermore, a linear dependence on coverage according to the v 1 H n values has been used in H f v m H n for all clusters. The calculated vibrational frequencies and ZPE are provided in Table S3. Fig. 4 c shows the contributions to the formation entropy of the individual species in a v 1 H 1 cluster: hydrogen in 3-fold and 4-fold configurations, and a Pd atom adjacent to the vacancy. The formation entropy of the clusters was predominated by the contribution from hydrogen, primarily due to the loss of the translational and rotational degrees of freedom of the H 2 molecule. The vibrational entropy of hydrogen in the cluster was highest for the 4fold site, and the entropic penalty was therefore lower by as much as 0.2 eV per hydrogen at 800 °C. In comparison, the contribution from the increased vibrational entropy of the Pd atoms surround-  ing the vacancy remained a fraction of the contribution of a single hydrogen.
The calculated concentration of v m H n clusters is shown as a function of inverse temperature in Fig. 5 . For the v 1 H n cluster, the configuration fully saturated with 4-fold hydrogen predominated up to about 15 °C ( Fig. 5 a). The hydrogen content decreased with increasing temperature as the entropic contribution, nT S f v m H n , became increasingly important. Eventually, v 1 H 1 , v 1 H 2 and the pure monovacancy exhibited similar concentrations at 800 °C. The concentration of the v 1 H n with 3-fold hydrogen was significantly lower, in particular at lower temperatures. Compared to the single vacancy clusters, the concentrations of the v 2 H n and v 4 H n clusters remained several orders of magnitude lower at all temperatures.
As shown for v 1 H n in Fig. 6 a, the concentration of individual v m H n clusters followed a linear relationship with the hydrogen partial pressure on a logarithmic scale, i.e., log ( . Accordingly, the total concentration of v 1 H n clusters increased substantially with p H 2 as the clusters with successively higher hydrogen content became predominating. Nevertheless, the concentration of the v 2 H n and v 4 H n clusters remained several orders of magnitude lower ( Fig. 6 b). Notably, the concentration of v 1 H n with 3-fold hydrogen showed a significantly lower dependence on p H 2 compared to the 4-fold clusters. This difference reflects the non-linear regime in H f v m H n as the 3-fold sites at closest proximity are filled ( n = 5-8 in Fig. 4 a).
Overall, the concentration of vacancies in the bulk was significantly enhanced in the presence of hydrogen, particularly at temperatures below 800 °C and at H 2 partial pressures above approx. 0.1 bar H 2 . Nevertheless, monovacancy clusters, v 1 H n , predominated larger clusters under all the considered conditions. Thus, the results showed no strong tendency for vacancies to agglomerate into larger clusters or nanovoids in the bulk.

Vacancy clusters at grain boundaries
The optimized structures of the 3 and 5 grain boundaries are shown in Fig. 7 . The 3 twin boundary is simply a stacking fault in the fcc lattice and the Pd site at the interface layer differs from bulk only in the symmetry of the 12-fold coordination. On the other hand, the 5 boundary is more open and distorted, and four different Pd-sites were considered for vacancy segregation and cluster formation. Fig. 8 shows the local structure around monovacancies at the grain boundary sites and the corresponding segregation energy, i.e., total energy difference between a cell with a vacancy residing at the grain boundary and in the bulk. The segregation energies were insignificant or positive for the 3 boundary and sites 1, 3 and 4 in the 5 boundary, and can in general be understood by considering the coordination environment and excess volume of the sites [48] . A substantial segregation energy of −0.67 eV was obtained for site 2 near the 5 boundary, and this site was therefore considered further with addition of hydrogen. A v 2 cluster was constructed by introducing an additional equivalent vacancy on the other side of the boundary plane and a v 3 clus-   ter was made with a third vacancy on site 1 in order to obtain a triangular geometry. The considered v m H n clusters at the 3 and 5 grain boundaries are depicted in Fig. 9 , and the corresponding formation enthalpies are presented in Fig. 10 . Hydrogen was added to 4-fold sites for all clusters, and four additional 3-fold sites were filled for v 3 H n in order to completely saturate the inner surface of this clus-ter. Comparing the formation enthalpy of v 1 H n in the 3 boundary with that of the bulk (dashed line in Fig. 10 ), there was a clear stabilizing effect at the 3 boundary for larger n . The main difference between these clusters is the connectivity between the 4-fold sites, which in the 3 boundary share edges in a pairwise manner, resulting in shorter H-H distances of 2.5 Å along one direction compared to 2.9 Å for v 1 H n in bulk ( Fig. 9 a). The stabilizing effect of the cluster was significantly larger at the 5 boundary, partly due to the segregation energy of the vacancy. Ultimately, the segregation enthalpies amounted to −0.4 eV and −1.4 eV for v 1 H 6 to the 3 and 5 boundary, respectively, and −1.9 eV for v 2 H 8 to the 5 boundary , thus indicating a strong tendency for segregation of the clusters towards the grain boundaries.
The concentration of v 1 H n clusters at the 3 (1 1 1) grain boundary is shown as a function of inverse temperature and compared to the bulk values in Fig. 11 . Despite the similarity in local structure and the negligible segregation energy of pure monovacancies, the concentration of v 1 H n clusters was substantially higher at the 3 boundary compared to bulk, especially towards lower temperatures. Considering the v 1 H 6 cluster in particular, the neg- ative formation enthalpy is clearly evident from the activation energy at lower temperatures.
The concentration of v m H n clusters is further increased at the 5 (2 1 0) grain boundary, in part due to the added segregation energy of the isolated vacancy. As shown in Fig. 12 , the considered site was saturated with v 1 H n clusters over the whole temperature range at 1 bar H 2 . In other words, the vacant site can be considered as a structural part of the grain boundary under these conditions. For the v 2 H n and v 3 H n clusters, the concentrations became saturated below about 175 °C and 20 °C, respectively. This is in remarkable contrast to the corresponding bulk concentrations which was at their highest values of 10 -8 and 10 -16 , respectively, at 800 °C ( Fig. 5 b). Clusters of all three sizes thus compete for the same sites at the 5 boundary at the lowest temperatures. The  predominating cluster at a given temperature may be evaluated by considering the Gibbs free energy of the clusters. From this, the v 3 H n clusters became favored over v 2 H n clusters below −233 °C, as apparent also from the minor gain in H f for v 3 H 15 compared to v 2 H 11 ( Fig. 10 ). Similarly, the v 2 H n clusters were favored over v 1 H n below 37 °C. Fig. 13 shows the equilibrium cluster concentrations as a function of p H 2 at 673 K at the 5 boundary. The saturation concentration was reached at a pressure as low as 0.2 bar for v 1 H n , and at about 5 bar and 30 bar for v 2 H n and v 3 H n , respectively. For the 3 boundary, the saturation limit of v 1 H n was reached at about 530 bar.

Inner surfaces of nanovoids
The thermodynamic parameters of hydrogen adsorption on the (0 0 1) and (1 1 1) surfaces are presented in Fig. 14 . For the (0 0 1) termination, the adsorption enthalpy was essentially unaffected by neighboring adsorbates. On the other hand, interactions between adsorbates were prominent for the (1 1 1) termination, and the surface was considered to be saturated with half of the sites occupied due to the positive adsorption enthalpy. The resulting amount of hydrogen adsorbates per unit area was four hydrogen atoms per 31.1 Å 2 and 26.9 Å 2 for the (0 0 1) and (1 1 1) surface, respectively. In terms of H-H interactions, the overall behavior of the adsorption enthalpies in Fig. 14 a was found to be similar to that observed for filling of the respective 4-fold and 3-fold sites in v m H n ( Fig. 4 b), as well as previous surface studies [57] . ZPE was −0.075 eV and 0.030 eV for the (0 0 1) and (1 1 1) terminations, respectively. The adsorption entropy for the (1 1 1) termination was essentially the same as for the 3-fold sites in v 1 H n and slightly more favorable for (0 0 1) termination in comparison to the 4-fold sites in v 1 H n ( Fig. 4 c).
Fig. 14 c shows the calculated equilibrium hydrogen coverage and surface energy as a function of temperature in 1 bar H 2 for both terminations. The (0 0 1) surface remained essentially saturated with hydrogen over the whole temperature range up to 1100 K. On the other hand, the coverage of the (1 1 1) surface decreased significantly with increasing temperature above 500 K since the adsorption enthalpy and entropy were both less favorable. While the (1 1 1) termination was the most stable in vacuum, the (0 0 1) termination was significantly stabilized in the presence of hydrogen to the extent of becoming the most favorable termination. Thus, restructuring of the palladium surfaces can be expected in hydrogen containing atmospheres. Fig. 15 shows the shape of larger nanovoids based on Wulff constructions using the calculated surface energies in vacuum and 1 bar H 2 at 673 K, as well as corresponding atomistic models of nanovoids of specific sizes. Adsorption of hydrogen clearly leads to a higher fraction of (0 0 1) facets, which further increases with increasing temperature due to the relative decrease in surface energy compared to the (1 1 1) termination ( Fig. 14 c). The atomistic models are discrete in size and in the ratio between facets. The theoretical facet ratio according to the Wulff shapes could therefore not be replicated exactly in the atomistic models, particularly for the smallest sizes.
The calculated surface energies of voids of discrete size are shown in comparison to the Gibbs formation energy of v m H n clusters in Fig. 16 . The energy of the voids increases proportionally with the surface area ( ∝ r 2 ) and becomes orders of magnitude larger than the most stable vacancy clusters. There is significant stabilization of the voids in the presence of hydrogen due the reduced surface energy with hydrogen adsorbates. Based on the large   energetical differences, it is reasonable to expect separate underlying mechanisms for the formation of clusters and nanovoids.

Discussion
The presence of point defects, such as vacancies in palladium, at finite temperatures is a thermodynamic necessity due the gain in configurational entropy of the system. The defect concentrations are determined by their Gibbs energy of formation ( Eq. (10) ), and the monovacancies are significantly stabilized by association with interstitial hydrogen on the internal 3-fold and 4-fold sites of the v 1 H n complexes. Thus, the concentration of monovacancies is enhanced by several orders of magnitude in the presence of hydrogen. However, coalescing of vacancies into larger clusters represents a reduction in the configurational entropy of the system, i.e., the principal cause for the presence of point defects in the first place. This aspect can be appreciated by considering the equilibrium constant of coalescing m monovacancies into a larger cluster, i.e., the reaction in Eq. (1) : where E b v m H n is the binding energy between monovacancies and interstitial hydrogen at a given temperature. From Eq. (16) it is ap-parent that the concentration of v m H n clusters follows a powerlaw dependence of the monovacancy concentration, c m v 1 , which predominates the exponential dependence on binding energy for larger clusters. Accordingly, the concentration of clusters of increasing size is quickly suppressed, as evidenced by the present results ( Fig. 5 b). It can therefore be concluded that the increased concentration of (mono)vacancies in the bulk is a not a direct cause of the formation of nanovoids or pores in hydrogen containing atmospheres.
In the case of grain boundaries, vacancies can be further stabilized on specific sites that exhibit negative segregation energies. Particularly in the presence of hydrogen, the concentration of vacancies may reach full occupancy for such sites, as was calculated for the v 1 H n , v 2 H n and v 3 H n clusters associated with the Pd2 site at the 5 (2 1 0) boundary ( Fig. 12 ). Due to the interconnectivity of the Pd2 sites, these vacancies may form channels along the grain boundary with diameters of around 2.5-4.1 Å (assuming an atomic Pd radius of √ 2 2 a 0 from the close-packed plane). Channels of this size may constitute fast diffusion paths for atomic hydrogen, while possibly also allowing Knudsen diffusion of small molecules such as He used for probing trans-membrane leakage [30] . The maximum size of voids formed by equilibrium concentrations of vacancy clusters will, nevertheless, be limited to the sites within absolute nearest vicinity of the grain boundary where structural distortions and segregation energies may be substantial. Accordingly, the largest diameter of such voids or channels can be expected to be about 0.5 nm in most cases.
With respect to the formation of larger voids, the surface energy of the predominating (0 0 1) and (1 1 1) facets was found to be significantly reduced in the presence of hydrogen due to stabilization by adsorbates ( Fig. 14 c). Nevertheless, surfaces and voids remain thermodynamically unfavorable relative to bulk, meaning that voids will not spontaneously form in palladium single crystals. Voids can take the place of other higher dimensional defects in polycrystalline material and thereby effectively reduce their relative energy. Moreover, formation of voids in the μm-range has been linked to plastic deformation in Pd membranes subjected to pressure gradients [28] . The origin of the formation of nanovoids and pores therefore appears to be closely linked to microstructural features and the presence of strain in the as-deposited films, as well as under operation. Several reports have detailed the formation of voids due to plastic deformation, i.e., in the absence of hydrogen [58][59][60][61] . It is important to consider that palladium undergoes significant chemical expansion due to dissolution of hydrogen, with additional contributions from vacancy clusters. For instance, the 1D chemical expansion reaches 3% for PdH 0.5 ( Figure S1), in good agreement with experimental data from Feenstra et al. [62] . The resulting strain in polycrystalline materials can be relieved by accommodating microstructural changes, and voids may become a relatively low energy alternative in the presence of hydrogen. Thus, small-grained polycrystalline palladium and materials that otherwise exhibit microstrain and low crystallinity shows particular potential for forming and retaining stable hydrogen-covered voids. The Pd-based membranes subject to void formation in previous studies exhibited these microstructural characteristics and significant grain growth occurred during operation [26 , 28-30] . According to these considerations, it is proposed that void formation in Pdbased membranes may be alleviated by annealing in hydrogen-free atmosphere so that crystallization, grain growth and densification may proceed in absence of chemical expansion and under conditions where voids are significantly less stable. While it is challenging to directly characterize v m H n complexes and clusters, the role of grain boundaries and microstrain on the formation of voids may be systematically studied in single-and polycrystalline samples annealed in hydrogen-containing atmospheres using transmission electron microscopy.
In summary, the formation of vacancy clusters and voids have different origins. The role of hydrogen on their stability is, however, similar in that hydrogen terminates the under-coordinated metal atoms of the inner surface of the cluster or void. While the increased concentration of vacancies in the presence of hydrogen may improve the kinetics of microstructural changes, the present work indicates that they do not directly cause the formation of voids. Hydrogen therefore has a confounding effect on the observed vacancy and void formation phenomena in palladium and similar metals. Further studies on hydrogen induced vacancy clustering should consider the role of chemical expansion and structures with ordered metal vacancies [12] .

Conclusions
Palladium vacancies are significantly stabilized in the presence of hydrogen due to favorable association between the vacancies and interstitial atomic hydrogen into v m H n clusters. Accordingly, the concentration of vacancy clusters is enhanced by several orders of magnitude. At grain boundaries, the concentration of clusters can reach site saturation and form structural channels with diameters of up to around 0.5 nm. Even at high hydrogen pressures, larger clusters and nanovoids remain thermodynamically unstable due to the loss of configurational entropy upon coalescing of monovacancy complexes. The formation of nanovoids in Pd-based membranes was ascribed to microstructural characteristics of polycrystalline materials, combined with strain induced by chemical expansion due to dissolution of hydrogen and plastic deformation due to pressure gradients.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.