Effect of uniaxial strain on the site occupancy of hydrogen in vanadium from density-functional calculations

We investigate the influence of uniaxial strain on the site occupancy of hydrogen in vanadium, using density functional theory. The site occupancy is found to be strongly influenced by the strain state of the lattice. The results provide the conceptual framework for the atomistic description of the observed hysteresis in the to phase transition in bulk, as well as the preferred octahedral occupancy of hydrogen in strained V layers.


I. INTRODUCTION
Vanadium is a transition metal with electronic configuration [Ar]3d 3 4s 2 that forms a body-centered cubic structure (bcc).The bcc structure contains tetrahedral and octahedral interstitial sites that can accommodate hydrogen 1,2 .In bulk vanadium the hydrogen is found to reside in tetrahedral sites at low concentrations (α-phase), while in the β-phase H occupies octahedral sites 1,2 .The symmetry of the hydrogen induced local strain field is strongly depending on the site occupancy, which is reflected in the hydrogen induced expansion of the lattice 1,2 .When hydrogen resides in a tetrahedral site, the local strain field is close to spherical, while it is almost uniaxial when hydrogen resides in the octahedral sites.
In a body-centered tetragonal structure there are three types of tetrahedral sites and three types of octahedral sites and there are in total three octahedral and six tetrahedral sites per metal atom (Figure 1).The six tetrahedral sites comprise of four T z and two T xy sites (T xy refers to either one of the equivalent T x or T y sites).The three octahedral sites comprise of one O z and two O xy sites.The local elastic response of the lattice arising from the presence of hydrogen in these interstitial sites [8][9][10][11][12][13] gives rise to a local strain field, which is the cause of the hydrogen induced expansion 9 .The expansion can be viewed as the sum of the hydrogen induced local strain fields and can therefore, in principle, be used to determine the preferred site occupancy of hydrogen.
The volume changes depend strongly on the boundary conditions and site occupancy, enabling the polarization of the local strain fields.This can, for example, be experimentally accomplished by the use of clamping 3,9 as the preferred site occupancy of hydrogen in vanadium is linked to the strain state of the structure [4][5][6][7] .The hydrogen induced volume changes of a clamped epitaxial film is restricted to the direction perpendicular to the surface.For example, a single crystal vanadium (001) film on a MgO(001) substrate will exhibit a lattice expansion (or contraction) in the (001) direction, independent of site occupancy 3,4 .Clamping of V can therefore be used to change the volume expansion, as only 1/3 of the strain field will propagate to the surface, when hydrogen is residing in tetrahedral sites.By the same token, the uniaxial component of the local strain field arising from a O z site occupancy will reach the free boundaries and therefore not restrict the expansion.Furthermore, when hydrogen resides in O x or O y sites, the local strain field can not give rise to any expansion due to the constraint implemented through elastic boundary conditions (clamping).Clamping and straining a V layer will therefore strongly affect the site occupancy, and thereby the observed volume changes 3,4,7,9 .
The goal of the present computational study is to provide a conceptual framework on the influence of strain on site occupancy in Vanadium.The results are generalised and can be used as a guidance when discussing hydrogen occupancy in any transition metal.Thus, we use computational methods to perform "virtual experiments" using an atomistic framework and first-principles methodology to explore the influence of polarization of the hydrogen induced local strain fields generated in bulk vanadium.The results provide a plausible atomistic description of the O z site occupancy at low concentrations in strained V layers 3 , change of site occupancy with hydrogen concentration in biaxially strained lattices as well as providing insight on the origin of the observed hysteresis 2 in hydrogen absorption and desorption in bulk V.The calculations were performed using the Vienna Ab initio Simulation Package (VASP) [14][15][16][17] .The interactions between the electrons and the nuclei was obtained using the projector-augmented-wave method 18,19 .The Perdew-Burke-Ernzerhof (PBE) 20,21 approach was employed to approximate the exchange and correlation terms in the density functional theory (DFT) 22,23 method.A conjugate gradient algorithm was used to relax the atomic nuclei positions to a local minimum in the total energy landscape.
In order to reduce H-H interactions resulting from the imposed periodic boundary conditions while still studying a system that is small enough to be computationally manageable for a large number of calculations, a supercell consisting of 128 vanadium atoms (4 × 4 × 4 bcc unit cells) was constructed to mimic bulk vanadium in which the lowest possible ratio of hydrogen to vanadium [H/V] is 1/128 (corresponding to 0.775 at.% of hydrogen).Due to the rather large dimensions (11.9 Å × 11.9 Å × 11.9 Å) of the supercell, only the Γ point was used in sampling the Brillouin zone.Comparisons of the total energy, using a 3 × 3 × 3 k-point mesh further established that Γ point sampling of the Brillouin zone is sufficient for a quantitative investigation.Zero-point energy calculations were deliberately excluded, not because they were unimportant, but rather because their correct treatment requires and deserves a separate dedicated study.
The x and y lattice vectors were fixed for all calculations, only allowing lattice relaxation in the z-direction.The motivation for this approach is to mimic the experimental conditions of hydrogen uptake in a superlattice where the bottom layer of a thin film of vanadium is held in place through strong bonds to a substrate.This constraint results in a one dimensional lattice expansion, perpendicular to the substrate's plane.This approach also mimics the bct lattice of the β phase when c/a 1.1 (see Figure 1).
It is sufficient to consider the strain in one direction, for example the z-direction, to capture the effect of strain on site occupancy.This implies also that we can treat the O x and O y sites to be equivalent when the hydrogen is occupying O z sites in the β-phase.These are identical in the sense that a rotation by 90 • around the z-axis will map the O x site onto the O y site and vice versa.

A. Strain and site occupancy
Before presenting the results from ab initio calculations we will provide a conceptual framework for the effect of strain on site occupancy, using a simple hardsphere model.We will use this approach to estimate the relative energetics of hydrogen occupation in octahedral and tetrahedral sites, solely based on the available interstitial space in a V single crystal.The maximum sphere radius that can be accommodated in the interstitial space formed by metal atom spheres arranged in a bcc lattice is 0.155 for octahedral and 0.291 for tetrahedral sites, in units of the metal atom sphere radius.In the atomistic model used here the vanadium has a Wigner-Seitz radius of 1.217 Å and the corresponding value for hydrogen is 0.370 Å.An octahedral site in V has therefore 0.155 • 1.217 Å= 0.189 Å of spherical radius available, which is much smaller than the hydrogen radius.A tetrahedral site provides 0.291•1.217Å= 0.354 Å, which is close to the hydrogen radius.From this consideration, one can see directly that it is not energetically favourable for hydrogen to occupy octahedral sites in an unstrained lattice, because of a large overlap between H and V electrons, raising the total energy through a Born-Mayer repulsion 24 .In the tetrahedral sites, the corresponding density overlap is much smaller, favouring occupation of tetrahedral sites.When the lattice is under uniaxial tensile strain (i.e., c/a > 1) the maximum sphere radius that can be accommodated in the T z and O z sites is more and more shifted in favour of the O z sites.The maximum sphere radius that can be accommodated in T z and O z sites becomes equal for c/a = 1.118.
When the lattice is expanded in the z-direction the T z and O z sites are energetically favored in comparison to their x and y oriented counterparts.This is rather obvious for the octahedral sites but not as clear for the tetrahedral sites.For the O z sites, a tensile strain in the zdirection will increase the spacing between the vanadium atoms that sit closest to hydrogen, while for O xy sites, the closest vanadium atoms lie in the xy-plane, which are geometrically unaffected by the uniaxial strain in the z-direction.In the experiments by Pálsson et al. 3 no distinction is made between hydrogen occupation of T z and T xy sites.
Figure 2 shows the ab initio calculated local strain fields in vanadium caused by hydrogen occupying either a T z or an O z site in a 128 vanadium atoms supercell.The arrows indicate the direction and the magnitude of the displacement of the vanadium atoms.Only the strain on the vanadium atoms in the O z and T z sites are shown (i.e., 4 atoms for a tetrahedral site and 6 atoms for an octahedral site).The isotropic strain field from hydrogen occupying a tetrahedral site and the strongly anisotropic strain field from occupying an octahedral site is clearly seen in Figure 2. The "top" and "bottom" vanadium atoms in the octahedron (i.e., the two vanadium atoms that possess the same x and y coordinates) are much closer to the hydrogen than the vanadium atoms in the tetrahedron; hence, the former are pushed farther away.In the absence of hydrogen the calculated lattice parameter is 2.99 Å .When hydrogen is placed in the O z site, the "top" and "bottom" V-atoms are displaced, increasing their mutual distance to 3.35 Å.This local strain corresponds to an increase of 12.1% in spacing between the V atoms which is in excellent agreement with the experimental results of 12.7% for β-phase VH To quantify the qualitative ideas obtained from the hard-sphere model, we used ab initio total energy calculations to determine the preferred hydrogen occupancy.Figure 3 shows a plot of the energy for a single hydrogen in a supercell occupying a T z , T xy , O z or an O xy site as a function of the c/a ratio ([H/V] = 1/128, corresponding to 0.775 at.% of hydrogen).The vertical axis shows [E(V + H) − E(V )] where E(V + H) is the total energy of the metal-hydrogen system and E(V ) is the total energy of the hydrogen free vanadium supercell, both calculated at the same c/a ratio.When the lattice is uniaxially strained ( c/a > 1), the O z sites will "open up" as described above and becomes more energetically favoured.This is easily inferred from the results in Figure 3 since the slope of the O z line is larger than than that of T z , thus at some c/a ratio the site occupancy of O z will be favoured as compared to the T sites.As seen in Figure 2, the strain is very large in the z-direction for hydrogen occupying an O z site.A comparison of the strain fields from hydrogen occupation of O z and T z sites shows a larger increase of available spherical radius for the hydrogen occupying an O z site.This, together with the favourable effect of the uniaxial tensile strain for the occupation of the O z sites makes the O z sites energetically favorable already when c/a > 1.051.The hard sphere model yielded a transition at c/a = 1.118 which can be considered as satisfying when considering the simplicity of the model.

B. Concentration dependence of site occupancy
It is not only the initial strain state which is the source of tetragonal distortion.The hydrogen induced volume changes will also influence the c/a ratio in clamped samples and thereby alter the energy balance between the T z and O z sites.Figure 4 compares the energies of T z and O z site occupancy at optimal c/a ratios to identify the critical hydrogen concentration where change in site occupancy occurs.The average energy of 50 structures with random hydrogen distributions for four different hydrogen concentrations is calculated (disordered phase).Change in site occupancy is approximated to occur between [H/V] of 0.34±0.07as this is where the total energy of T z and O z site occupancy becomes equal.The shift in the site occupancy is found to be driven by energetics rather than entropy.The configurational entropy was determined using Boltzmann's entropy formula and the internal energy was approximated as the number of hydrogen atoms in O-sites times an energy penalty of 0.2 eV (i.e., we approximate that moving a hydrogen atom from a T site to an O site will raise the energy by 0.2 eV, in accordance with the difference in energy between T and O site occupancy, cf. Figure 3).For all tested hydrogen concentrations and for a broad temperature range, the internal energy is always found to dominate the entropy part, so that coexistence of T and O sites is concluded unlikely to occur in the low concentration region.The hydrogen induced lattice expansion can however give rise to change in site occupancy.This can take place in both ordered and disordered phases, thus a O z occupancy does not need to imply an ordered β-phase and a change of site does therefore not by necessity imply a disorder-ordered phase transition.

C. Volume expansion and hysteresis effects
Figure 5 shows the resulting uniaxial lattice expansion (quantified by the c/a ratio) as a function of hydrogen concentration [H/V].The relationship between calculated c/a ratio and hydrogen concentration [H/V] is to a very good approximation linear with a slope of 0.120 for T z sites and 0.236 for O z sites, see Figure 5.These results are in excellent agreement with the experimental results by Pálsson et al. 3 which determined the expansion to be 0.1189 (7) for hydrogen occupation in tetrahedral sites and 0.19(1) for octahedral sites.The calculated change in total volume due to hydrogen occupation of T z and O z sites are 1.61 Å3 and 3.14 Å3 , respectively per added hydrogen atom.O z occupancy gives rise to larger increase in volume, as compared to T z occupancy, due to the anisotropy of the local strain field as seen in Figure 2. The strain component in the z-direction is larger for O z than for T z sites, implying that a shift from T z to O z occupancy is accompanied by an increased c/a ratio, which favours O z occupancy.The shift in site occupancy from T z to O z can therefore be viewed as a self-amplified process.The shift in site occupancy resembles therefore in many ways a first order phase transition.Now we will discuss the difference in the lattice response when the hydrogen concentration is increased or decreased.In an unstrained or nearly unstrained lattice hydrogen is exclusively found in T z sites .When increasing the hydrogen concentration from low concentrations, the expansion will open up the O z sites, which become energetically favored above the critical (c/a) crit value of 1.051.The uniaxial lattice expansion will therefore result in a shift in site occupancy from T z to O z at that concentration.(c/a) crit for the change from T z to O z occupancy is marked by a vertical line in Figure 3 and a horizontal dashed line in Figure 5.
When starting at a high concentration all the hydrogen will reside in O z sites.When decreasing the concentration, (c/a) crit =1.051 will be reached at [H/V] = 0.212, resulting in a shift in site occupancy from O z to T z .Thus, when increasing the concentration the shift from T z to O z is reached at a different concentration as compared to the change of sites from O z to T z sites when decreasing the hydrogen concentration.Thus, a hysteresis with respect to lattice expansion is expected when loading and unloading H under the specified conditions and the thermal excitations are smaller than the energy difference between the two states.These effects do resemble the α to β phase transition in bulk V, with respect to both change of sites as well as observed hysteresis 2 .Furthermore, these results clearly illustrate the effect of clamping on the site occupancy, which can be changed without entering the β-phase in V.When the initial strain of the sample is changed, these boundaries will move as illustrated in Figure 5: With a biaxial compressive strain in the x − y plane, the boundaries will move to lower concentrations and the hysteresis gap will decrease.When c/a will be larger than a critical value, hydrogen will solely reside in O z sites, as inferred from experiments 3 .O z n o i n i t i a l s t r a i n T z n o i n i t i a l s t r a i n O z 3 % i n i t i a l s t r a i n T z 3 % i n i t i a l s t r a i n FIG.5: (Color online) The uniaxial lattice strain c/a resulting from varying the concentration of hydrogen occupying exclusively either Tz sites (blue data points and lines) or Oz sites (red data points and lines) in bcc vanadium.The horizontal black dashed line at c/a = 1.051 marks the critical uniaxial lattice strain for which the hydrogen occupancy of Tz and Oz site becomes equal in energy, as seen in Figure 3.The vertical colored full lines indicate at which hydrogen concentration the critical c/a ratio of 1.051 is reached for occupancy of Tz ([H/V] = 0.424) and Oz ([H/V] = 0.212) sites, respectively, when there is no initial strain (i.e.c/a = 1.00).The dotted lines represent the case of an initial strain of c/a = 1.03 before any hydrogen has entered the system.The critical c/a ratio is then reached at [H/V] = 0.173 for Tz occupancy and 0.084 for Oz occupancy.

IV. SUMMARY
The preferred interstitial site occupancy in vanadium with constrained boundaries has been studied using calculations based on density functional theory.The energetics of hydrogen atoms in a bcc-bct supercell were investigated to provide a conceptual understanding of the experimentally observed shifts in site occupancy 3 .In the investigated range of c/a from 1.00 up to 1.07, the tetra-hedral (T z ) sites are energetically favored for hydrogen occupation in comparison to the octahedral (O z ) sites in the c/a range from 1.00 to 1.051.The octahedral sites are energetically favored for hydrogen occupation when c/a > 1.051.The forces exerted on the vanadium lattice by hydrogen atoms occupying interstitial sites will alter the global strain state which in turn triggers a shift in site occupancy above the critical value of c/a = 1.051.This self-amplified process can be understood by the obtained strain field from octahedral (O z ) site occupancy which has a larger z-component than that obtained for a tetrahedral (T z ) site occupancy.The increase in c/a as a function of hydrogen concentration [H/V] is linear and in good agreement with previously obtained experimental results.
The different slopes of c/a as a function of [H/V] for tetrahedral (T z ) and octahedral (O z ) site occupancy has the consequence that the condition for shift in site occupancy is met at different hydrogen concentrations [H/V] when starting from high (O z ) or low concentration (T ).This leads to the theoretical prediction of a hysteresis in the hydrogen loading-unloading process, in which the switch from T z to O z site occupancy and the reverse switch from O z to T z occur at different hydrogen concentrations.The experimentally observed 3 coexistence of tetrahedral and octahedral hydrogen occupation in the [H/V] concentration range of 0.065-0.068might be an indication that such a hysteresis behavior could indeed be found in vanadium.

V. ACKNOWLEDGMENTS
Financial support from the Swedish Research Council is gratefully acknowledged.The project is part of the COST Action MP1103.O.E. also acknowledges the KAW foundation, eSSENCE, STandUP for Energy, and the ERC (project 247062 -ASD).The calculations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC).

FIG. 1 :
FIG. 1: (Color online) The different types of (a) tetrahedral and (b) octahedral interstitial sites in bcc vanadium are illustrated here.Large dark spheres represent vanadium atoms and the small red, blue, light blue, and light red spheres represent (according to their respective labels) different interstitial positions that hydrogen can occupy.The z-axis is aligned along the vertical direction, while the x-and y-axes lie in the horizontal plane.

FIG. 2 :
FIG. 2: (Color online) The strain on surrounding vanadium atoms by hydrogen occupying a (a) Tz site or (b) Oz site.The arrows represent the displacement vectors, i.e., how much the V atoms are "pushed" away by the H atom.The length of the arrows has been scaled by a factor of 30.

FIG. 3 :
FIG. 3: (Color online) Energy difference as a function of an externally applied global uniaxial lattice strain c/a where ∆E = E(V + H) − E(V ).The dashed vertical line at c/a = 1.051 marks the critical uniaxial lattice strain for which hydrogen occupancy of Tz and Oz sites becomes energetically equivalent.

FIG. 4 :
FIG.4: Energy difference between Tz and Oz occupancy at optimal c/a ratios as a function of hydrogen concentration.∆E = (ET z − EO z )/NH where NH is the number of hydrogen atoms in the simulation.Data points are for average values and the bars indicate ± the standard deviation (defined as the square root of the variance).Lines are second order polynomial fits.