Rational design of Nb-based alloys for hydrogen separation : A first principles study

We have investigated the effect of alloying metal elements on hydrogen solubility and mechanical integrity of Nb-based alloys, Nb15M1 (where M = Ca–Zn, Ge), using first principles-based calculations. In general, the chemical interaction between the interstitial H and metal is weakened as the alloying element is changed from an early to a late transition metal, leading to lower H solubility and higher resistance to H embrittlement. This effect becomes more pronounced when a smaller alloying element is used due to stronger elastic interaction between interstitial H and metal atoms. These finding may provide scientific basis for rational design of Nb-based hydrogen separation membranes with tailored H solubility to effectively suppress H embrittlement while maintaining excellent hydrogen permeation rate.


I. INTRODUCTION
With increasing demand for clean and sustainable energy, hydrogen separation and purification becomes increasingly important, especially in reforming of carbon-containing fuels (fossil and renewable fuels such as biomass) and CO 2 capturing.It is also an important step in hydrogen production from other renewable sources such as photo-electrolysis of water, where efficient collection of hydrogen may dramatically enhance the efficiency and production rate.In an integrated gasification combined cycle (IGCC) or a steam reforming process, for example, hydrogen must be separated from CO 2 and other byproduct gases efficiently at low cost. 1 In some cases, the purity of hydrogen is vital to subsequent applications (e.g., as a fuel for PEM fuel cells or feedstock for higher-valued chemicals).However, the conventional gas separation technologies, such as pressure swing adsorption (PSA) and cryogenic distillation, are energy intensive and inefficient, more so for higher purity of hydrogen.2][3] In particular, Pd based alloys (e.g., Pd-Ag, Pd-Au, and Pd-Cu) have been extensively studied because of their high selectivity, excellent catalytic activity, and moderate permeability of hydrogen. 1,2 nfortunately, broader commercialization of Pd-based membranes is hindered mainly by the cost of the materials.5][6] With a BCC structure, the group-V metals can dissolve larger amount of hydrogen and allow faster hydrogen diffusion (and hence higher permeability) than Pd and Pd-based alloys of an FCC structure.However, the high solubility of hydrogen often results in H embrittlement fracture under service conditions.a Electronic mail: byungki.ryu@samsung.com.b Electronic mail: hc0303.park@samsung.com.For example, a Nb membrane suffers from embrittlement fracture under typical IGCC operating conditions (400 • C, 7 atm), where the atomic ratio of H to Nb varies from 0.7 to 1.0. 3,6 hen the H concentration (H/M) is greater than 0.7, the structure of the metal-hydrogen solid solution system changes from α to β (or γ ) phase, which is susceptible to brittle fracture.This mechanical failure has severely hindered commercial application of these membranes.It is thus very important to suppress the formation of the latter phase in the membranes.
][10][11][12] Computational studies of Ti-and Pd-alloys suggest that the interactions between H and metal alloying atoms may be elucidated by the atomic size effect and/or Miedema's "reverse stability" rule. 8,9 owever, the origin of H stabilization/destabilization mechanism in metals and alloys are yet to be revealed.Although Aboud and Wilcox found a correlation between the H solution energy and H charge state in pure Nb and Pd metals, 11 they are unable to quantitatively explain the interactions between H and metal atoms in alloys.
In this paper, we report our prediction of hydrogen solution energy in BCC Nb metal and Nb based alloys using density functional theory (DFT) calculations. 13We systematically examined the effects of alloying elements on hydrogen stabilization, which may offer the scientific basis for rational design of alloy membranes with desired H solubility by controlling the chemical and elastic interaction energy between H and metal ions.

II. CALCULATION METHOD
In the DFT calculations, 13 we used the generalized gradient approximation (GGA-PBE) 14 for the exchange-correlation functional and the projector augmented wave potentials as ionic pseudopotentials, as implemented in VASP code. 15,16 o describe the accurate electronic structure of Nb alloys, we include the inner shell electrons as valence electrons.The 4p electrons are included in the Nb pseudopotentials.For Sc, Ti, V, Cr, Mn, Fe, Ni, and Cu, the ionic potentials include the 3p electrons.In the case of Co, Zn, Ga, and Ge atoms, 3p electrons are considered as core electrons.All integrations over the Brillouin zone were performed using a 6×6×6 grid of Monkhorst-Pack k-point mesh 17 for a conventional 2×2×2 supercell of BCC.The energy difference between the total energy with 6×6×6 k-point grid and 16×16×16 k-point grid was less than 0.5 meV/atom.The wave functions were expanded in plane-waves with a kinetic cutoff energy (E CUT ) of 300 eV.It is also found that the H solution energy using E CUT = 300 eV is converged within 10 meV, as compared to that using E CUT = 600 eV.The alloy structures were modeled by replacing one of the 16 Nb atoms by an alloying metal atom, resulting in a supercell consisting of 15 Nb atoms and 1 alloying metal atom.We considered elements from Ca to Zn and Ge on the 4 th row of the periodic table as the alloying elements.The H concentrations in the metals, c = H/(Nb+M), were varied from 1/16 to 1 2 in the calculations.All atomic positions and the lattice parameters were fully relaxed until the residual forces were smaller than 0.02 eV/Å and the total energy convergence within 1 meV/supercell.
The solution energy of H in a Nb based alloy is defined as  is negligible in our Nb-based alloy systems.From the phonon mode calculation, the zero point energies of 16-atom Nb supercells without and with H atom are calculated to be 454 and 689 meV, respectively.Thus, the H zero point energy predicted by the phonon mode calculations is slightly different by 3 meV/H, as compared to that predicted by the vibrational frequency calculations.After taking into consideration the vibration frequency of the H 2 molecule (4257 cm -1 ), E ZPE [H T ] is predicted to be 106 meV/H.After taking into consideration the vibration frequency of the H 2 molecule (4257 cm -1 ), E ZPE [H T ] is predicted to be 106 meV/H.For Nb-M alloys, the calculated E ZPE ranges from 75 to 125 meV.As the effect of alloying atoms on E ZPE is nearly the same, the trends of H solution energy do not depend on inclusion of E ZPE .Zero-point energy correction is not included although an estimation of the absolute solubility will depend on the inclusion of these effects.
It is found that Fe is non-magnetic in Nb.When H is positioned at a nearby Fe, however, a magnetic moment is induced.Since the magnetization energy is relatively small (less than 10 meV per Fe atom), we neglected the magnetization effect of Nb alloys on H stabilization.

III. RESULTS AND DISCUSSION
When the H atom is introduced in BCC Nb, it prefers locating at a tetrahedral interstitial site (T) than at an octahedral interstitial site.H at a T-site (H T ) is surrounded by four Nb atoms.The ideal distance from the T-site to the nearest neighbor is 1.86 Å; however, it increases to 1.93 Å after relaxation in this case.On the other hand, H at an octahedral site (H O ) is surrounded by six Nb atoms, which is less stable than H T by 0.293 eV/H (without E ZPE correction) because the size of the octahedral interstitial site is smaller than that of the T-site.Similarly, we found that H T in Nb-based alloys studied, Nb 15 M 1 (M = Ca to Zn and Ge), is also more stable than H O .We further calculated the vibration modes of H T and H O in Nb metal.H O has two imaginary frequencies, as reported by Ouyang and Lee, 12 implying that the O-site is not a local minimum.Note that, the H O in Nb metal is a saddle point configuration and the force on H O is exactly zero due to symmetry.Thus, most H in BCC Nb metal and Nb-based alloys prefers the T-sites at low H concentrations, consistent with a previous study. 18he H T atoms in Nb 15 M 1 alloys can be classified into two groups: H T near-M and H T far-M , depending on the local atomic structure.H T near-M has three Nb and one M as the nearest neighbors.The H 1s orbital interacts with the orbitals of the nearest neighbor metals, forming a H-induced bonding state (BS H ), as shown in [Figs.1(a)-1(d)].Since Ca has empty p and d bands, the interaction between Ca and H is relatively weak, as compared to the interactions with other alloying elements.For those with partially filled d bands, the H 1s state is well hybridized with the partially filled d electrons of the alloying element.For M = d 10 , there is a repulsive interaction between H 1s and the s orbitals of the alloying element.Because of the difference between Nb and the alloying element, the coordination tetrahedron is not symmetrical: the bond length of Nb-H T differs from that of M-H T , which is proportional to the atomic size of M and ranges from 1.65 to 2.18 Å.

Energy (eV)
Total DOS (states eV  On the other hand, H T far-M refers to H located far away from the alloying metal element M and is surrounded by four Nb atoms as the nearest neighbors.The coordination tetrahedron is symmetrical, with the average bond length of Nb-H T far-M ranging from 1.90 to 1.94 Å, whereas the length of Nb-H T in pure Nb metal is 1.93 Å.We calculated the effective H volume per atom [vol(H)] in Nb 15 M 1 by subtracting the volume of Nb 15 M 1 alloy lattice from that of Nb 15 M 1 -H hydride lattice.The calculated average vol(H T ) is 2.84 and 2.86 Å 3 /H in 16-atom alloy supercells for c = 1/16 and 1/2, respectively, consistent with the experimental results of 3 Å 3 /H. 3e then investigated the electronic structure of Nb metal, Nb 16 -H T hydride, and H atom in vacuum and drew the total-density of states (DOS) [Fig.2(a)].The vacuum level of Nb metal and H atom were estimated using the reference potential method. 19The work function of the Nb(100) surface is predicted to be 3.587 eV in DFT.The calculated 1s level of an H atom is -7.540 eV below the vacuum level, as expected from the spin-polarized calculations.The average band energy of Nb metal is positioned at -5.075 eV, which is higher than the position of the H 1s orbital level.Thus, it is anticipated that H in Nb metal will obtain excess electrons from Nb atoms.According to Bader charge analysis, 20,21 each H in Nb-H hydride is negatively charged, having 1.68e at around the H, consistent with Ref. 11.However, integration of the charge densities at around T-site in Nb metal shows there is no significant charge transfer from nearest neighboring metal atoms to the H atoms.For a sphere of 1.0 Å, there are 1.50e at around H T in ideal Nb 16 -H T .However, there are already 0.85e at around T-site in Nb metal before H incorporation. Similar argument can be found for H in V metal in Ref. 12.The bonding character between H and metal atoms is not ionic.Rather, the H atom seems to interact covalently with the nearest neighboring Nb metal atoms.The H 1s orbital is hybridized with s, p, and d bands of the nearest neighbor metals, similar to H multicenter bonds in oxides, 22 forming a deep BS H at -10.304 eV below the vacuum level.We expect that the antibonding state may appear above E F , resonant with metal bands [Fig.2(b)].
Fig. 3 shows the projected-DOSs (PDOSs) which are projected on H T in Nb 16 -H T , T-site in Nb, M in Nb 15 -H T and M in Nb 15 M with the projection sphere radii of 1.5 Å.As shown in Fig. 3(a), the PDOS of T-site is decreased at energy ranging from -5 to 0 eV due to the formation of BS H .Meanwhile, the antibonding state between H and metal atoms is appeared in the metal PDOSs of Nb 15 M 1 -H T at above the E F .For M = Nb, Ca, V, Ni, Ge, the integrated charges of M within 1.5 Å at energy range of 0 to 5 eV are 4.36, 2.10, 5.89, 1.50, 0.32 eV.And they are increased to 4.68, 2.36, 6.30, 2.04, 1.40 eV, respectively.Note that the electrons that would occupy the antibonding state are transferred to the Fermi-level, as the substitutional H at the O site in MgO or ZnO donates electrons to the conduction band minimum state. 22For Nb, the E F of the metal-H solid solution is slightly increased by about 70 meV, as compared to that of pure metal.It appears that H atoms in Nb alloys are stabilized by forming BS H with the nearest neighbor metal atoms.As shown in Fig. 3(a), there are already abundant delocalized electrons at T-sites, consistent with Ouyang and Lee's report. 12After H incorporation, however, the electrons of the metal atoms adjacent to the H are transferred to the localized BS H , [Fig.3(b)] resulting in an electronic energy gain, which is defined as where E M and E H denote, respectively, the average band energy of the M nearest neighbor within a sphere radius of 1.5 Å in Nb 15 M 1 alloy and the band energy of BS H state in the Nb 15 M 1 -H T hydride.With the positive electronic energy gain, the H atom overcomes the strain energy when H is accommodated at the T-site.Note that there is a negligible correlation between the H solution energy and the H Bader charge state.In Nb-based alloys, the H Bader charge ranges from -0.52 e to -0.68e, and is not proportional to the H solution energy.Rather, the E SOL [H T near-M ] is strongly correlated with the elec .As drawn in Figs.4(a) and 4(b), E SOL [H T near-M ] decreases with increasing elec .Note that the H solution energy of the ideal structure (i.e., the lattice without structure relaxation) can be determined by the electronic energy gain and can be plotted on the single curve as a function of the elec .Like this, the H is stabilized with the electronic energy gain.After structure and lattice relaxations, the H solution energy is lowered by about 0.096 to 0.367 eV/H.The calculated H solution energy in Nb 15 M 1, E SOL [H T near-M ], ranges from -0.350 to 0.439 eV/H, whereas that in Nb is -0.357 eV/H.For Ca (d 0 ) with empty d bands below the E F , there is little interaction between Ca and H, resulting in a large E SOL in the ideal structure.However, E SOL [H T near-Ca ] is lowered by 0.367 eV/H after structure relaxation.In the case of M = d n (n = 1 to 7), the H 1s state is well hybridized with the partially filled d electrons of M and elec ranges from 5.576 to 5.975 eV, as compared to elec of Nb (5.575 eV).In essence, the transition metal atoms attract H atoms, acting as H traps. Approaching to the late transition metals, E M decreases and so does the electronic energy gain, especially for M = d around -1.988 and -3.241 eV below E F , respectively.Their corresponding H solution energies are -0.114 and 0.003 eV for the ideal structure.For the relaxed structure, E SOL [H T near-M ]'s are expected to be -0.290 and -0.173 eV/H, higher than that of H T in pure Nb metal.For d 10 elements (M = Zn and Ge), the solution energy becomes positive due to the repulsion between H 1s and metal s states (Fig. 1(d)).Similarly, in BCC Ti-SM alloys (where SM = Al, Si, Ga, Ge), the repulsive interactions between H and SM were reported. 9he lattice constant of Nb 15 M 1 differs from that of BCC Nb (a = 3.32 Å in the calculation and a = 3.30 Å in the experiment), depending on the alloying element (See Table I).The volume of Nb 15 M 1 alloy lattice is larger for M = Ca and Sc than that of Nb lattice.As the alloying element is changed from V to Ni, the volume change goes from −3.74 to −12.81 Å 3 , reaching the smallest value for Nb 15 Ni 1 due to the small atomic size of Ni.We note that the volume change is nearly negligible for Ti and d 10 elements.The H solution energy depends on the unit volume of the lattice because of the elastic interactions between H atoms and the lattice.For H T in BCC Nb with c = 1/16, the H solution energy decreases linearly with lattice parameter due to easier strain relaxation in a larger lattice.We define the H solution energy at a given lattice parameter of a, E (a) SOL [H], defined as E (a) SOL = (1/n){E (a) [ , where E (a) [Nb 15 M 1 -nH] and E (a) [Nb 15 M 1 ] are the total energy of the relaxed structure with fixed lattice parameter of a.The calculated E (a) SOL s with c = 1/16 are -0.135,-0.340, and -0.516 eV/H for a = 3.26, 3.32, and 3.38 Å, respectively [see Fig. 4(c)].Similarly, the formation of a metal hydride is accelerated by enlarging the lattice volume.Due to strain relaxation around the H, the H solution energy of H T for c = 8/16 is lowered in the large lattice by -0.069 eV/H, as compared to that for c = 1/16 (we used the same metal-hydride configuration for Nb-H and Nb alloy-H).We would like to point out that the H-H repulsion energy does not change severely, when the lattice is enlarged.For two H atoms at a nearest T-site distance in 54-atom 3×3×3 Nb BCC supercell, the H-H repulsion energy is slightly decreased only by -15 meV, with increasing lattice parameter from 3.32 to 3.41 Å.We also note that the solution energy of H in α-Fe is positive because of the small lattice size and the large strain energy.By performing additional DFT calculations on H T in α-Fe, the E (a) SOL [H T in Fe] is predicted.The E (2.84Å) SOL [H T in Fe] is 200 meV/H, implying that H is hardly incorporated in the Fe metal.However, when the lattice parameter of Fe is larger than 2.9 Å, the E (a) SOL [H T in Fe] becomes negative, resulting in the attractive interaction between H and Fe.We also investigated E SOL [H T far-M ] in Nb 15 M 1 alloys with high H concentration of c = 8/16 (i.e., 8 interstitial H atoms in a 16-atom metal supecell) and confirmed that the large lattice volume enhances H stabilization.The interaction between the alloying element M and H atoms is minimized when the H atoms are placed far away from the M atom (M-H bond length > 2.8 Å).Under this condition, the distance between H atoms should be larger than 2.5 Å to avoid H repulsion.Similar to H in Nb, E SOL [H T far-M ] is increased as the lattice parameter of the alloy and the bond length between H T far-M and Nb nearest neighbors are decreased [Fig.4(b)].For M = Fe and Ni, the lattice parameters of the alloys are expected to be 3.29 and 3.27 Å and E SOL [H T far-M ]'s are -0.337 and -0.275 eV/H, respectively.For M = Zn and Ge, however, the H solution energies are smaller than that for M = Ni because the atomic size of Zn and Ge are larger than Ni.
These results seem to suggest that a late transition metal with smaller atomic size will be a good alloying element to suppress the formation of metal hydrides.For example, the H solubility in  According to a H site-occupying-model in metals, 23 H solubility of BCC Nb metal is about 0.8.Thus, the H solubility in a Nb-based alloy is lower than this limit.Summarized in Table II are the calculated H solubilities of various Nb 15 M 1 alloys.As expected, alloying with early transition metals (e.g., Ti and V) is ineffective in limiting the H solubility.For alloying elements M = Ca to Co, the elastic interactions between M and H play an important role in reducing H solubility, although there are attractive interactions between M and H atoms.In the case of M = Ni, Cu, Zn, and Ge, the H solubility is lowered at both H T near-M and H T far-M .Thus, it is expected that M = Mn-Zn and Ge are good candidates as alloying elements for Nb alloys to suppress the formation of metal hydrides if the BCC phase is stable.

IV. CONCLUSION
In conclusion, we have investigated the H solution energy in Nb-based alloys with several elements on the 4 th row of the periodic table, providing critical insights into the mechanism of H stabilization through chemical and elastic interactions between H and the alloying elements.In general, the chemical interaction between H and metal is weakened as the alloying elements are changed from early to late transition metals, leading to lower H solubility and greater tolerance to H embrittlement.This effect is further enhanced when an alloying element of smaller geometric size is used due to stronger elastic interactions.Accordingly, Nb-based alloys with a late transition metal element of small atomic size may be ideally suited to limit H solubility and suppress H embrittlement.
FIG. 2. (a) Total DOSs of pure Nb metal, Nb 16 -H T solid solution, and H in vacuum.The zero energy is at vacuum level.Fermi energy is denoted by the vertical line.(b) Schematic illustration of the between the T-site electrons and the H 1s state to form the H-induced bonding state.The H 1s orbital interacts with the T-site projected metal bands and forms a bonding state below metal bands and an antibonding state in the conduction band.

FIG. 3 .
FIG. 3. (a) Projected DOSs of T-site in Nb (dotted line) and H T in Nb 16 -H T (solid line).(b) Projected DOSs of alloying metal (M = Nb, Ca, V, Ni, Ge) in Nb 15 M 1 (dotted line) and in Nb 15 M 1 -H T .For the projected DOSs, sphere radii of 1.5 Å are used.The zero energy is at E F .

FIG. 4 .
FIG. 4. (a) Plot of E SOL [H T near-M ] of Nb 15 M 1 -H T (c = 1/16) versus elec = E M -E H .The dashed and solid lines stand for the E SOL [H T near-M ] of ideal and relaxed structures, respectively.(b) Plot of E SOL [H T near-M ] of Nb 15 M 1 -H T (c = 1/16) versus M are shown for the ideal structure (unfilled square) and the relaxed structure (filled square) and compared with the electronic energy gain of elec = E M -E H . (c) E (a) SOL [H T ]'s with given lattice parameter of a are drawn for Nb 16 -H T and Fe 16 -H T (c = 1/16).(d) E SOL [H T far-M ]'s of Nb 15 M 1 -8H T (c = 8/16) solution are compared with the average lattice parameters of the Nb 15 M 1 alloy ( a alloy ), where a alloy is defined as (1/2)[Vol.(Nb15 M 1 )] 1/3 .

8 +
C and 7 atm will be lowered by a factor of exp[-E SOL /k B T] = 0.21 for H T far-M in Nb 15 Fe alloy.By taking into consideration both H T near-M and H T far-M in Nb alloy for c = 1/16 and 8/16, the H solubility (N) may be roughly estimated (using the Boltzmann factor with respect to the H solubility in Nb at 400 • C and 7 atm) as follows, N = N[H near-M T ] + N[H far-M T ] = min.{0.25, (24/96) × exp[− near-M /k B T] × 0.min.{0.75, (72/96) × exp[− far-M /k B T] × 0.8 where near-M and far-M are E SOL [H T near-M ] -E SOL [H T in Nb] for c = 1/16 and E SOL [H T far-M ] -E SOL [H T in Nb] for c = 8/16, respectively.

TABLE I .
The calculated average lattice parameters ( a alloy ) and volumes of (2×2×2) Nb metal and Nb 15 M 1 alloy lattice supercells for various alloying elements.The average lattice parameter a alloy is defined as (1/2)[Vol.(Nb15 M 1 )] 1/3 .The volume change of alloy ( Vol.) is estimated by subtracting the volume of Nb 16 metal lattice from that of Nb 15 M 1 alloy lattice.

TABLE II .
The H solubilities of Nb and Nb 15 M 1 alloys are estimated with respect to that of Nb (80 % at 400 • C and 7 atm) for various alloying elements Nb 15 M 1 (where M = Ni-Ge) would be much smaller than that in a pure Nb metal (about 80 % at 400 • C and 7 atm), because it is energetically unfavorable to position the H atoms immediately next to the M atom (larger H solution energy).Accordingly, the presence of the M ion may effectively block 24 of the 96 T-site nearest neighbors in a (2×2×2) BCC supercell, reducing the effective number of empty T-sites to 75 %.In addition, as the atomic size of the alloying element is reduced, the reduction in lattice volume may further lower the probability for the T-sites to be occupied by H, leading to greater resistance to H embrittlement.In BCC Nb alloys with M = Fe and Ni, E SOL [H T near-M ]'s are smaller by 0.090 and 0.152 eV/H than that in BCC Nb metal.Thus, the H solubility at 400 •