Structural properties of Y-doped BaZrO3 as a function of dopant concentration and position: A density functional study

Solid State Ionics 232 (2013) 1–12

Contents lists available at SciVerse ScienceDirect

Solid State Ionics journal homepage: www.elsevier.com/locate/ssi

Structural properties of Y-doped BaZrO3 as a function of dopant concentration and position: A density functional study D. Zeudmi Sahraoui, Tzonka Mineva ⁎ UMR 5253 CNRS/ENSCM/UM2/UM1, Institut Charles Gerhardt Montpellier, 8 rue de l'Ecole Normale, 34296 Montpellier cedex 5, France

a r t i c l e

i n f o

Article history: Received 24 July 2012 Received in revised form 6 November 2012 Accepted 6 November 2012 Available online 23 December 2012 Keywords: Y-doped BaZrO3 Hydrogen insertion Structural properties Quantum chemical DFT calculations

a b s t r a c t Relative stabilities, structural details and hydrogen binding sites in models of Y-doped barium zirconates with 12.5, 25 and 37.5% yttrium have been studied using Density Functional Theory (DFT) based periodic approach. Two stable crystal structures were obtained for the possible Zr substitutions in the 25 and 37.5% Y-containing BaZrO3. Tetragonal space group symmetries were found for Y ≥ 25% and a volume reduction is obtained with increasing the acceptor dopant percentage. Structural and electronic properties remain very similar in all the considered models. A well established charge difference is noticed only for oxygen sites in the three possible Zr\O\Zr, Zr\O\Y and Y\O\Y configurations. Examination of H binding strength in various OH local environments indicated a general tendency to stabilize the O\H bonds in the Zr\OH\Y configuration(s) compared to the Zr\OH\Zr ones. A clear tendency of increasing the O\H stabilization with increasing Y content from 12.5 to 25% is not obtained. The strongest O\H bonds are formed in the 37.5% Y-doped BaZrO3 model. Formation of O\H bonds in Y\OH\Y configurations was not found that suggest a smaller probability for H-uptake in Y-doped zirconates with high Y\O\Y concentrations. © 2012 Elsevier B.V. All rights reserved.

1. Introduction The potential use of Y-doped barium zirconate as an electrolyte in Protonic Ceramic Fuel Cells, (PCFC) operating in the temperature range of 400–700 °C [1–4] has prompted the research activities towards detailed description of its structural characteristics, found determinative for the proton conduction rate [5]. The lattice cell constant, the spatial and local symmetry, and the proton defect concentrations are among the most important factors that have been so far pointed out to affect the hydrogen mobility in acceptor-doped BaZrO3 (BZ) [5]. Despite the progress of going inside in the structure–proton conductivity relationship [3–11], there are still several open problems to be overcome in order to meet the requirements for technological applications. The major issues, related with the material preparations, discussed in the literature are associated with the poor sinterability and grain boundary conductivities of acceptor-doped zirconates, including Y-doped BaZrO3 (BZY) [12–18]. Along with the research efforts in developing new synthesis routes, it emerged also the necessity to go deeply insight into the basic issues related with the mechanism of proton conduction in solid-state materials. To this aim structural studies using Fourier Transform Infrared Spectroscopy (FTIR), powder X-ray diffraction (XRD), High-Resolution Transmission Electron Microscopy

⁎ Corresponding author. Tel.: +33 4 67 16 34 68; fax: +33 4 67 16 34 70. E-mail address: [email protected] (T. Mineva). 0167-2738/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.ssi.2012.11.011

(HRTEM) and electron diffraction analysis [19] focused on the establishment of the grain morphology and its relation to the proton conductivity in annealed 20% Y-doped BaZrO3 nanograins. The small grain size (of about 10 nm) shows low proton conducting capacity, whereas the increasing grain sizes to well-crystallized structure enhanced the sample proton conductivity. In a more recent studies [17,18] it was demonstrated that after extended annealing at temperatures of ~ 1700 °C the grain boundary conductivity is enhanced along with increasing dopant concentration in the grain boundary region. Local structural order/distortion of doped zirconates was studied by means of IR and Raman spectroscopes [9,20], inelastic neutron scattering and X-ray absorption spectroscopy methods [21]. In these studies highly distorted structures are associated with lowering of H-transfer capabilities. The Extended X-Ray Adsorption Fine Structure (EXAFS) technique was applied to obtain more precise results about the local environment of cations and lattice dynamics in Y-doped barium zirconates as a part of more systemic investigations of doped barium zirconates and cerates [7,8,10]. Based on these EXAF spectra, the Zr\O distances could be calculated in zirconates with different dopant ions, such as In, Ga, Gd, Y and Sc. The general conclusion was that Y 3+ and Gd 3+ produced the highest local disorder, which seems to be the reason for the high bulk conductivity, compared to Sc 3+ and In 3+. The quantum chemical computational methods, which are able to discriminate local structures at different sites and the associated charges, relative stabilities and proton transfer activation energy

2

D.Z. Sahraoui, T. Mineva / Solid State Ionics 232 (2013) 1–12

BaZrO3 consist of 2 × 2 × 2 supercells. The Y percentages of 12.5 (BZY12.5), 25 (BZY25) and 37.5% (BZY37.5) were modeled by replacing one, two or three Zr atoms in the 2 × 2 × 2 supercell, respectively. In the BZY25 and BZY37.5 structures two possible supercells were considered •• for different Y-positions (Fig. 1). Oxygen vacancies (V O ), compensating the negative charge defects, created upon replacement of Zr 4+ cations by Y3+ acceptor-dopant, have been modeled as well. To this aim BZY25 (2× 2 × 2) supercell was considered, because the substitution of two Y 3+ ions compensates electrostatically the removal of one oxygen (Scheme 1). The composition of the supercell, which contains oxygen vacancy is thus Ba8Zr6Y2O23. Because the effect of oxygen vacancies was obtained negligible on the ionic charges, electronic structures and the inter-atomic distances, except for the first shell inter-atomic dis•• •• tances in the Zr\V O \Zr and Zr\V O \Y configurations, the oxygen defects were not considered for the other two dopant percentages. Further on, insertion of hydrogen ions has been modeled for the three Y 3+ concentrations. The charge balance of the supercells is maintained by inserting one, two and three hydrogen ions in the BZY12.5, BZY25 and BZY37.5 models, respectively. The composition of the supercells representing 12.5, 25 and 37.5% of Y dopants thus becomes Ba8Zr7YHO24; Ba8Zr6Y2H2O24 and Ba8Zr5Y3H3O24, respectively. For the exchange-correlation (XC) functional, the hybrid B3LYP approximation [32,33] was used. The atomic orbitals were described with Gaussian basis sets as follows: Ba – HAY&WADT – small core effective core potential [34]; all electron basis sets for Zr [35], Y [36], O (8-411d11G) [37] and H (H_3-1p1G) [38]. The input parameters (basis sets and approximation for the exchange-correlation functional) have been chosen based on the benchmark calculations for crystal cell constant; interatomic distances, vibrational frequencies and electronic band gap energy of the ideal BaZrO3 cubic crystal. Test calculations have been carried out with LDA (VWN) [39], GGA (PBE) [40] and hybrid (B3LYP [32,33] and B3PW [41,42]) functionals. These preliminary results demonstrated that the best comparison with experimental cell constant and band gap energy values are achieved with the hybrid functionals. Whereas, the LDA and GGA produced reasonable cell constant parameters when coupled with effective core potential bases, the band gap energies are strongly underestimated by about 2–3 eV. This finding agrees very well with previous studies predominantly

barriers, while accounting for the electronic structures of each constituent cation in the solid, have been applied as well. Although these computational methods are based on various approximations, for instance consideration of idealized models, temperature equal to 0 K, assumption of well defined long range ordering under the periodic boundary conditions, etc., they have so far provided sound results about thermodynamics and energetics of defect formation in BaZrO3 [22], proton mobility, assessed by extensive mapping of migration pathways, as a function of dopant ionic radius for dopant concentrations below 12.5% [23], and the proton conduction pathways [24], and activation energies of the proton transfer between two oxygen sites in 12.5% of Y-doped BaZrO3 crystals [25]. Energies of water adsorption on bared and Y-doped (001)-BaZrO3 have been also computed [26,27]. More recent studies on the Y-doped-BaZrO3 [28] and on the dry or hydrated BaZrO3 [29] grain boundary structures, pointed out the increased energy barrier for the proton migration from one grain to an adjacent grain in BZY [28] and the oxygen vacancy segregation in the BZ structures [29] as possible reasons for the experimentally established low proton conductivity in the grain boundary regions. All these works employed methods based on the Density Functional Theory (DFT) in its periodic schemes. With the use of ab-initio path-integral molecular-dynamics simulations it was demonstrated that accounting for nuclear quantum effects has an important impact on the computed transition state barriers and that the reorientation step becomes rate limiting below 600 K [30]. To address the questions of how precisely the local structures around the dopants and the hydrogen binding strengths vary upon increasing Y-dopant concentrations, we provide in the present study structural and energetic details of models representing three different concentrations of Y dopant (12.5, 25 and 37.5%) in BaZrO3. The results showed only a little influence of the long-range ordering and crystal symmetries on the local properties. Instead, a well pronounced dependence of atomic charges and H-binding energies on the local environment around oxygen sites is obtained. 2. Models and computational details Periodic Density Functional Theory based computations were carried out with Crystal06 computer program [31]. Models of Y-doped

O2

O1 O3

BZY12.5

BZY37.5-M1 1.36

BZY25-M1 0.78

BZY25-M2 0.0

BZY37.5-M2 0.0

Fig. 1. Optimized structures of BZY with 12.5 (BZY12.5), 25 (BZY25) and 37.5% (BZY37.5) of Y. Oxygens in Zr\O\Zr, Zr\O\Y and Y\O\Y configurations are labeled as O1, O2 and O3, respectively. Relative energies in eV of the BZY25 and BZY37.5 models M1 and M2 are reported below the structure names. The green, red, gray and cyan colored balls denote yttrium, oxygen, zirconium and barium ions, respectively.

D.Z. Sahraoui, T. Mineva / Solid State Ionics 232 (2013) 1–12

a)

3

b)

Scheme 1. Schematic presentation of BZY25 supercell with oxygen vacancy: a) oxygen vacancy in Zr\V O•• \Zr configuration and b) oxygen vacancy in Zr\V O•• \Y. Characteristic bond distances (in Å) and bond angles (in degrees) as obtained from the full geometry optimization are reported as well.

based on the DF plane wave methods using LDA – or GGA – types for the exchange-correlation approximation and reporting a band gap being underestimated by ~ 2 eV [23,43,44]. This is a problem inherent to the method, which is solved by employing hybrid approximations for the XC-functionals within the DFT. Our test calculations showed that either B3LYP or B3PW hybrid XC-functional performed very comparably. The use of already optimized and available small-core effective core potentials bases [35] for Zr in the Y-doped structures led to strong underestimation of the cell constant parameters compared to the experimental values. A strong improvement was obtained when employing all-electron bases for Zr. Therefore we choose the B3LYP and all-electron bases for Zr and Y. At this level of theory, the fully optimized cell constant of the ideal BaZrO3 with a space group symmetry Pm-3m, is 4.25 Å agreeing with previous DFT-based calculations [22,24,26,27,44]. Experimentally, 4.19 Å have been reported [45]. The band gap energy computed by us is 5.59 eV, in a significantly better agreement with the experimental value of 5.3 eV [46] when compared to the available LDA [43] and GGA [23,44] results. Bilić and Gale [47] found a minimum energy cell of the ideal BaZrO3 with orthorhombic space group symmetry. The normal mode analysis within the vibrational frequency calculations, we have carried out for the optimized ideal cubic barium zirconate, showed only positive values. This clearly indicates that the cubic cell obtained with the present approximations of the computational method, is also a local minimum structure on the potential energy surface. In addition, we note that the computed vibrational frequencies of 99 cm−1 (Ba\ZrO3 stretching), 196 cm−1 (O\Zr\O bending) and 505 cm−1 (Zr\O stretching) match very well the experimental frequencies of 115 cm−1 (Ba\ZrO3 stretching), 210 cm−1 (O\Zr\O bending) and 505 cm−1 (Zr\O stretching) [48]. Based on the obtained agreement between our benchmark calculations and the experimental data concerning the cell constant, band gap energy and vibrational frequencies of the ideal cubic BaZrO3, we have chosen to employ the DFT-B3LYP method for the present study of Y-doped BaZrO3 structural properties. For all the calculations, computational tolerances for controlling the accuracy of Coulomb and exchange, and correlation series were set equal to: 10 −7, 10−7, 10−7, 10−7 and 10−24. A shrinking factor of 4 × 4 was used to generate the commensurate k-points grid in the reciprocal space within the Pack–Monkhorst method [49]. Atomic positions and crystal cells were fully optimized by employing the analytical geometry optimization based on a quasi-Newton optimization scheme with Broyden–Fletcher–Goldfarb–Shano (BFGS) algorithm for the Hessian matrix updates as it is implemented in Crystal06 program.

3. Results and discussion The effect of dopant concentration on cell parameters and lattice symmetry, interatomic distances, relative energies and Mulliken atomic charges have been first examined by studying model structures containing 12.5, 25 and 37.5% of Y without hydrogen in the cell. For 25 and 37.5% of Y dopаnts two possible positions of Y in different environments were found to produce minimum energy structures. The influence of oxygen defects on the geometrical and electronic properties has been studied as well for BZY25 models. Hydrogen atoms were subsequently inserted and the geometry optimizations carried out without any constrain allowed us to estimate the structural changes induced by the presence of hydrogen in the BZY unit cells. Geometrical properties and binding energies of hydrogens bound to three different types of oxygens in the Zr\O\Zr (O1), Zr\O\Y (O2) and Y\O\Y (O3) configurations have been compared in order to understand the role of the local and long range structural arrangements on the defect sites and H-binding sites and energies. 3.1. Geometrical and energetic properties of BZY models Optimized structures of the considered BZY models are presented in Fig. 1. Neutral atomic centered wave functions were used to describe the atomic wave functions in order to keep the BZY supercells neutral. The substitution of one, two and three Zr by Y atoms in the cells representing 12.5, 25 and 37.5% BZY, respectively, gave rise to minimum energy structures with spin polarized density equal to one in BZY12.5, to two in BZY25 and to three in BZY37.5. The energetically favored structures in the BZY25 and BZY37.5 are those where two neighbor Zr atoms are substituted. The configurations with the same Y percentages, but alternating Zr\O\Y\O\Zr are obtained to be less stable by 0.78 and 1.36 eV (Fig. 1), for the BZY25 and BZY37.5 crystals. Most of the experimental studies [3,6,11,14] reported on cubic cells of Y-doped BaZrO3 ceramics with Y concentrations varying from 10 to 40%. The XRD measurements were predominantly performed for powder samples obtained from different synthesis conditions, and especially prepared at different sintering temperatures. Yamazaki et al. [6] investigated samples with Y-dopant concentration in the 20 to 40% range and established a linear increase of the crystal cell constants with the Y concentration. These results do not agree with the previous ones of Kreuer et al. [5] showing that lattice cell constants, respectively cell volumes, increases up to 20% Y concentration in the BaZrO3 and significantly decreases for the 25% of Y dopants. Moreover, the space

4

D.Z. Sahraoui, T. Mineva / Solid State Ionics 232 (2013) 1–12

symmetries are established to vary as well as with the Y concentration. The ceramics with 2 and 5% Y are characterized with cubic lattices, which became tetragonal with increasing Y content up to 20% and again cubic for 25% Y. A dependence of the space symmetry on the Y dopant percentage and on the Y position in the BZY systems has been also obtained from our calculations, as follows from the results reported in Table 1. The symmetry of the optimized BZY12.5 is cubic. It breaks down to P4/mmm and Pmmm in the first (M1) and second (M2) model configurations in BZY25 and becomes again cubic for BZY37.5 (M1). The relative energy calculations favor however, the lower symmetry BZY25 and BZY37.5 structures (structures M2 in Fig. 1) with Pmmm (BZY25) and P4/mmm (BZY37.5) space group symmetries. Comparison with the results in Ref. [5], also summarized in Table 1, shows a very reasonable agreement between the optimized cell parameters of BZY25 for the two M1 and M2 structures. In the case of 12.5% Y model, the optimized cell parameters/cell volume are better agreeing with the experimental data for 15% Y concentration than for 10% Y. To the best of our knowledge, previous theoretical calculations are available for models of barium zirconates with Y ≤ 12.5 [22,24], giving no details about the optimized cell constant parameters or symmetries of acceptor-doped BaZrO3. The largest local structural distortion is found in the YO6 octahedrons with Y\O distances differing by 0.02 Å in BZY25-M2 models. The Zr\O distances differ by a maximum of 0.005 Å. The averaged inter-atomic distances, presented in Table 2, show a slight decrease

Table 1 Cell parameters in Å, volume in Å3 and symmetry for Y-doped BaZrO3 structures, presented in Fig. 1, with 12.5, 25 and 37.5% of Y at different positions in the (2 × 2 × 2) supercells and with oxygen vacancies in the BZY25, labeled as BZY25\V O•• . Comparison with available experimental data for Y-doped BaZrO3 is also provided. Structure

Cell constants

Volume

Symmetry

BZY12.5 BZY25-M1

a = b = c = 4.233 a = 4.214 b = c = 4.207 a = 4.220 b = 4.196 c = 4.226 a = 4.185 b = c = 4.210 a = 4.142 b = 4.229 c = 4.184 a = b = c = 4.183 a = 4.203 b = c = 4.175 a = b = 4.215 c = 4.204 a = 4.204 a = 4.196 a = b = 4.231 c = 4.213 a = 4.216 a = b = 4.241 c = 4.226 a = 4.210 a = 4.210 a = 4.216 a = 4.220 a = 4.230

75.835 74.583

Pm-3m P4/mmm

74.830

Pmmm

74.175

P4/mmm

73.289

Pmmm

73.192 73.261

Pm-3m P4/mmm

74.706

P4mm

74.300 73.877 75.446

Pm-3m Pm-3m P4mm

74.938 76.017

Pm-3m P4mm

76.618 74.618 74.995 75.151 75.687

Pm-3m Pm-3m Pm-3m Pm-3m Pm-3m

BZY25-M2

BZY25-V O•• a BZY25-V O•• b

BZY37.5-M1 BZY37.5-M2 BZY10c BZY10d BZY10e BZY15c BZY16f BZY20c BZY20g BZY20h BZY25c BZY30g BZY40g

Oxygen vacancy in Zr\V O•• \Zr configuration. Oxygen vacancy in Zr\V O•• \Y configuration. c From Ref. [5]. d From Ref. [3]. e From Ref. [11]. f From Ref. [14]. g From Ref. [6]. h From Ref. [19] for samples annealed at 800 °C, whereas 4.219 and 4.214 Å are reported for the samples annealed at 1250 °C and at 1500 °C, respectively. a

b

Table 2 Averaged inter-atomic distances in Å in Y-doped BaZrO3 with different Y concentrations and oxygen vacancies in the BZY25 supercell, labeled as BZY25-V O•• . Structure

Zr\O

BZa BZY15a BZ BZY12.5 BZY25-M1 BZY25-M2 BZY25-V O•• b BZY25-V O••c BZY37.5-M1 BZY37.5-M2

2.10 2.11 2.12 2.12 2.10 2.11 2.11 2.13 2.09 2.09

Y\O 2.23 2.13 2.12 2.11 2.09 2.08 2.10 2.10

Zr\Ba 3.63 3.65 3.68 3.67 3.68 3.73 3.66 3.61 3.63 3.63

Y\Ba 3.65 3.61 3.56 3.44 3.49 3.59 3.61 3.42

a From Ref. [10] for the anhydrated samples. Note that the distances in the hydrated samples increases by a maximum of 0.01 Å. b Oxygen vacancy in Zr\V••O\Zr configuration. c Oxygen vacancy in Zr\V••O\Y configuration.

(up to 0.02 Å) for the Zr\O and Y\O separation with increasing Y % in the cells. Y\O distances are always slightly longer than the Zr\O in agreement with recent EXAFS results for 15% Y [10]. These EXAFS data show however a larger difference between the Y\O and Zr\O distances compared to our calculations. A very good agreement is found for the other inter-atomic distances. The more pronounced distance variation with yttrium increase is obtained for r(Y\Ba): Y\Ba shortens from 3.61 Å in BZY12.5, to 3.44 Å in BZY25, to 3.42 Å in BZY37.5. The O\O separation between the oxygen atoms in the three different octahedrons, i.e. ZrO6\Zr, ZrO6\Y and YO6\Y, present in the doped structures have been pointed out among the most important factors affecting the H-transfer barriers in the proton transfer reactions [25]. We found that these oxygen–oxygen separations are however only a little dependent on the Y percentage and position, especially for the ZrO6\Zr and ZrO6\Y octahedrons (Table 3). A more pronounced O\O shortening is found for the YO6\Y configurations. The O\O distances vary between: 2.99 (BZY12.5)–2.93 Å (BZY37.5) in the ZrO6\Zr octahedrons; 2.98 (BZY12.5)–2.94 Å (BZY37.5) in ZrO6\Y octahedrons; and 3.01 (BZY25)–2.90 Å (BZY37.5) in YO6\Y octahedrons. Note that O⋯O separation larger than 2.9 Å is known for many perovskites [50]. The presence of oxygen defects lead to further structural disorder in the oxygen vacancy environment and shrinking of the cell parameters as follows from the results reported also in Table 1. A higher local disorder is observed for the model with Y\V •• O \Zr configuration (Scheme 1b). However, full lattice relaxation of structures with protonic defects (BZHY models discussed in Section 3.3) yielded again an expansion of the cell constants, thus compensating the obtained shrinking due to the oxygen defects prior the formation of hydroxyl ions. Latter result correlates well with the established cell constant expansions in the as-synthesized samples at T = 130 °C BZY20, attributed to fully incorporated protonic defects and absorbed water [19]. Significant changes of inter-atomic distances and angles are obtained in the local environment of V •• O , which are shown in Scheme 1. Interestingly, we found that on the averaged scale, the increased local structural disorder in the near V •• O region does not influence remarkably the Y\O, Zr\O, Ba\Zr and Ba\Y inter-atomic separations (Table 2) when comparing them to those in the BZY25 model without oxygen vacancies. 3.2. Electronic structures of BZY models The influence of dopant variations on the total density of states (TDOS) spectra, atomic charges as obtained within Mulliken scheme

D.Z. Sahraoui, T. Mineva / Solid State Ionics 232 (2013) 1–12 Table 3 O⋯O separations in Å within a ZrO6 octahedron adjacent to a ZrO6 octahedorn (ZrO6intra-Zr); within a ZrO6 octahedron adjacent to a YO6 octahedron (ZrO6-intra-Y); and within a YO6 octahedron (YO6-intra) in the BZY12.5, BZY25 and BZY37.5 supercells.

5

Table 4 Mulliken atomic charges in e− and band gap energy in eV in Y-doped BaZrO3 structures, presented in Fig. 1, with 12.5, 25 and 37.5% of Y at different positions in the (2×2×2) supercell models and with oxygen vacancies in the BZY25, labeled as BZY25-VO•• .

Structure

ZrO6-intra-Zr

ZrO6-intra-Y

YO6-intra

Structure

Zr

O1

O2/O3

Y

Ba

Band gap

BZY12.5 BZY25-M1

2.99 ± 0.01 2.98 ± 0.01

2.98 ± 0.01 2.96 ± 0.01

BZY25-M2

2.98 ± 0.01

2.98 ± 0.01

BZY37.5-M1

2.971

BZ BZY12.5 BZY25-M1 BZY25-M2 BZY25-V O•• a BZY25-V O•• b BZY37.5-M1 BZY37.5-M2

2.91 2.92 2.90 2.90 2.89 2.90 2.90 2.90

−1.58 −1.55 −1.51 −1.50 −1.52 −1.53 −1.57 −1.47

–/– −1.33/– −1.28/– −1.38/−1.01 −1.37/– −1.32/– −1.19/– −1.24/−0.85

– 0.90 0.75 0.88 0.74 0.80c 0.58 1.00; 0.54

1.82 1.82 1.82 1.81 1.81 1.81 1.82 1.81

5.59 5.53 4.78 4.61 4.65 4.65 4.74 4.62

BZY37.5-M2

2.944 2.954 2.926

2.941 2.951 2.937 2.983 2.965 2.970

3.014 3.004 2.998 2.971 2.981 2.973 2.979 2.972 3.030 2.946 2.853 2.961

and the band gap energies has been analyzed as well. The yttrium uptake did not lead to a meaningful variation in the total DOSS structures, shown in Fig. 2. The predominant atomic orbital contributions to different bands of TDOSS are also indicated in Fig. 2. As it follows, the DOS structures of BZY closely resemble that of the ideal BZ cubic cell. Atomic charges, q, are presented in Table 4 together with the band gap values. The presence of Y dopants does not influence Zr and Ba atomic charges, remaining in all structures independent on their symmetry and cell parameters, equal to 2.9 e − for Zr and 1.8 e − for Ba. Similarly, only a little variation is obtained for the q(O1) ~− 1.5 e − in the Zr\O\Zr configurations (see Fig. 1), whatever is the Y percentage and position in the cell. A more significant decrease of oxygen charge down to ~− 1.2 to − 1.3 e − and ~− 1.0 to − 0.9 e −is found for the O2 (Zr\O\Y) and O3 (Y\O\Y) centers. This values are again little affected by the precise cell symmetry and cell parameters. Moreover, there is a reasonable agreement between our atomic charges and those reported in Ref. [5], also computed within Mulliken scheme, although using a semi-empirical Density Functional Tight Binding Method for the electronic structure calculations and averaging over twenty seven configurations of supercells with three

a b c

Oxygen vacancy in Zr\V O•• \Zr configuration. Oxygen vacancy in Zr\V O•• \Y configuration. The charge of yttrium ion in the Zr\V O•• \Y configuration is 0.71 e−.

Y-dopants and three hydrogens, produced from the molecular dynamics simulations at T = 1000 K. These atomic charges are 1.796, 1.283, 0.679, − 1.085 and − 1.092 e − for Ba, Zr, Y, O(Zr) and O(Y) ions, respectively [5]. The decrease of the negative charge at oxygen ions next to Y 3+ dopant, found by us, could be related to the partially covalent antibonding character of Y\O bonds. As follows from the overlap population matrix calculations, a predominant anti-bonding for the Y\O bonds is indicated by a negative value of the related element in the overlap population matrix of about − 1.1 for the O1\Y and O2\Y bonds and of about − 1.3 to − 1.5 for the O3\Y ones. The other matrix elements are zero, thus indicating only ionic bindings. It is worth noting that the atomic charges (Table 4) and the bonding character are found practically unchanged in the models with oxygen vacancies. Only the charge of the oxygen in front of V •• O decreases down to − 1.01 e − (Scheme 1a) and − 1.18 e − (Scheme 1b). The Y\O bonding character and the atomic charges are preserved also for the H-containing structures (vide infra). This confirms again that there is an important charge redistribution and structural distortions in the nearest surrounding of the defect (acceptor-dopant, oxygen vacancy, or protonic defect), whereas the ions and geometry properties in the defect-free octahedrons in the host oxide appear to be nearly unaffected. 3.3. Structural parameters of hydrogen containing BZHY models

Fig. 2. Density of states spectra in BZ (a), BZY12 (b), BZY25 (c) and BZY37.5 (d). Spectra in dashed lines are those of the M2 models in BZY25 and BZY37.5. Atomic orbitals, predominantly contributing to a particular band, are indicated in the top of panels (a) and (b).

Hydrogen insertion leads to further reduction of space symmetry of the fully optimized structures and induces distortions up to various extents depending on the H-positions and mutual orientations. In the BZY12.5, BZY25 and BZY37.5 models, the insertion of one, two and three H atom(s) to the O1, O2 and O3 sites was considered. All the optimized minimum energy structures are collected in Appendix A, Figs. A1–A5, in order to illustrate the induced structural deformation due to the H insertion into the host oxide lattices. The Y\O distances in the hydrogen free and hydrogen containing octahedrons are also reported in Figs. A1–A5. The resulted minimum energy structures in BZHY12.5 are schematically drawn in Scheme 2 (top). Local minima structures were obtained for hydrogen bound to either O1 or O2 sites with bond distances of 0.985 Å (O2-Pm), 0.978 Å (O1-Pm) and 1.22 Å (O1-Pmmm). For the three H-containing structures, there is a volume increase of ~ 1.3% compared to the hydrogen-free BZY12.5 models. The Zr\OH and Y\OH distances increase in the interval of 0.15– 0.25 Å, whereas nearly negligible distance variation is found for the Zr\O/Y\O bonds in the H-free octahedrons. The Y\O2H distance of 2.25 Å is very close to the value of 2.29 Å, reported for the Y\OH bond distance in a similar configuration of protonic defect and yttrium dopant concentration [23]. Due to the H-bond formation,

6

D.Z. Sahraoui, T. Mineva / Solid State Ionics 232 (2013) 1–12

O1

Zr

O1

Zr

O1

Zr

H O1 H

O1

Zr

O1-Pm

O1

O1

Zr

O2

H

O2

O1-Pmmm

Y

O2-Pm H O1

H

H O1

O1

Zr Zr

Zr

O1 O1

H

H

O1 H

O1-O1-a

O1-O1-p

O1-O1-perp

Scheme 2. Schematic presentation of the considered hydrogen binding sites and H-mutual orientations: (top) BZHY12.5 structures and (bottom) BZHY25 structures. Analogously, same O\H mutual orientations were considered for the O2 and O3 sites.

the respective O⋯O distances reduce down to 2.896 Å in O1-Pm and to 2.400 Å in O1-Pmmm structures in Scheme 2. In the contrary, the distance between HO2⋯O elongates slightly up to 3.03 Å. The elongation of the HO\Zr distances, very comparable to our results, was also reported in Ref. [51] for the insertion of hydrogen proton and atom into an undoped BaZrO3. Several possible local minima structures were found for the insertion of two H atoms in the cells with two Y 3+ (BZHY25). The H-insertion caused 1–2% volume variation. The volume is obtained to shrink or expand depending on the H-mutual orientations and precise oxygen binding site (s) as follows from the results in Table 5. The space symmetries in most of the structures are P1 and P2. Similarly to the BZHY12.5 cases, the Zr\OH and/or Y\OH distances in the H-containing octahedrons increase by 0.12–0.15 Å. The structural variations depend certainly on the mutual orientations of the two H-atoms, which are schematically illustrated in Scheme 2 (bottom) for the Zr\OH\Zr configurations. Same H-binding sites and orientations were considered for the Zr\OH\Y and Y\OH\Y possible configurations. It matters also whether the two H atoms bound to oxygens in two adjacent octahedrons as presented in Scheme 2 (bottom), or the two binding oxygens are in two different planes. All the relaxed structures are displaced in Appendix A section, Figs. A2–A4. A negligible effect of the different long-range ordering in the two possible M1 and M2 BZHY25 supercell models is found on the local geometrical properties around the OH sites (Table 5), except that hydrogen binding to O3 centers in M2 cells was not obtained either in BZHY25 or in BZHY37.5. During the optimization procedure, the Y\OH\Y angle bends strongly, thus leading to a very short Y 3+\Y 3+ distances causing a strong repulsion between the two Y 3+ cites followed up by an energy divergence. This effect is obtained for several mutual H-orientations and combinations of H-insertion

involving at least one O3 site. Most likely, the low charge density at the O3 site of ~− 1.0 e − (Table 4) is not sufficient to form a stable O\H bond and we observe a continuous elongation of the H\O3 bond length, followed by a strong bending of Y\OH\Y angle. Configurations with Y\O\Y defect regions have not been considered previously and our results represent indeed the first evidences for the limitation of protonic defect formations at O3 sites in Y-doped barium zirconates. H-insertion in BZY37.5, also led to more pronounced structural deformations that did not allowed reaching minimum energy states for the most of the initially started combinations of H-binding to O1 and O2 centers. The fully optimized minimum energy structures were obtained only for one O1\O1\O1, one O2\O1\O1 and one O2\O2\O2 H bindings in M1 configuration (Fig. A5). Interestingly, several initial HO1\HO1\HO1 configurations resulted in an H-transfer from O1 to O2 site during the optimization, thus leading in all cases to HO1\HO1\HO2 combination. In these cases, strong decrease of Zr\O\Zr angle was obtained, followed by a HO1⋯O2 distance shortening down to 2.35–2.40 Å, allowing thus to transfer the hydrogen to O2 site. This O1⋯O2 distances resemble the optimal O⋯O separation distance that was proposed (from the DFT calculations) necessary to achieve an activation energy of nearly 0.4–0.5 eV for the H-transfer in BZHY12.5 model [25]. A rough estimation of the energy difference between the energy at the geometry point, where an H-atom is shared between the O1 and O2 centers, and the final minimum energy, corresponds indeed to 0.4–0.5 eV. The obtained shrinking of the cell volume for the BZH37.5 compared to BZH25, which is related with a slight decrease of the interatomic separations (Table 2) and especially those of the O⋯O intra-octahedral distances (Table 3) is most probably the reason for the observed H-transfer to the next oxygen ion, thus reducing the number of possible minimum energy states.

D.Z. Sahraoui, T. Mineva / Solid State Ionics 232 (2013) 1–12

In the three considered Y percentages, the O\H bond distances are between 0.96 and 0.99 Å in the most of the stable structures. Latter result is in good agreement with previous DFT data concerning the O\H distances in 12.5% Y-doped barium zirconate [23]. A very little variation of the atomic charges is found (by less than 0.05 e −), except for the binding oxygens, whose charge became less negative by 0.1–0.2 e −. In all BZHY structures, the hydrogen has a charge of 0.3–0.4 e −.

Table 5 Averaged inter-atomic distances in Å, volume (V) in Å3 in BZHY with Y=12.5, 25 and 37.5% for a particular H-binding site(s). In the structures with more than one H atoms the mutual orientation between H-atoms is labeled with “a”, “p”, “perp” for anti-parallel, parallel or perpendicular O\H-orientations for H-insertions in the same plane (Scheme 2). The prime symbol is used to distinguish the H-insertion in different planes. The HO…O are the shortest distances between the hydrogen bound oxygen and the neighbor oxygens in the same octahedron. Structure

H-binding site(s)

O\H

Zr\O

HO…O

V

BZHY12.5

O1-Pm

0.978

76.21

O1-Pmmm

1.226

(×18) (×2) (×19)

2.896

BZHY12.5

2.408

76.19

BZHY12.5

O2-Pm

0.985

(×20)

O1\O1-a

1.051

3.027 4.312 2.930

76.3

BZHY25-M1

74.15

BZHY25-M2

O1\O1-a

0.996

(×15) (×2)

2.983

74.45

BZHY25-M1 BZHY25-M1

O1\O1-p O1\O1-perp

0.967 0.997 0.982

2.124 ± 0.02 2.244 ± 0.01 2.125 ± 0.04 2.242 (×2) 2.126 ± 0.05 2.303 2.114 ± 0.05 2.248 2.236 2.110 ± 0.04 2.218 ± 0.02 2.362 2.112 ± 0.03 1.996 2.115 ± 0.04 2.280 2.209 2.107 ± 0.04 2.298 1.999 2.116 ± 0.06 2.344 2.118 ± 0.05

(×23)

2.800 2.867 2.680

75.61 75.32

2.667 2.866 3.100

74.84

74.92

(×18)

2.967 2.892 2.953 2.865

BZHY25-M1

O2\O1-a′

BZHY25-M1

O2\O1-p′

O1: 0.978 O2: 0.986 O1: 0.969 O2: 0.987

BZHY25-M1

O2\O2-a

0.988

BZHY25-M2 BZHY25-M2

O2\O2-a O2\O2-p

0.990 0.974

BZHY25-M2

O2\O2-p′

0.995 0.997

BZH37.5-M1

O1\O1\O1

0.987 1.042 0.994

2.124 ± 0.04 2.053 2.121 ± 0.07 2.371 2.106 ± 0.04 2.308 2.239 ± 0.02 2.086 ± 0.08 2.280 2.224 ± 0.02

(×17)

(×10)

(×17)

3.4. Hydrogen binding energies in BZYH models The O\H bond strength is related to the hydrogen transfer in the sense that a sufficiently strong O\H bond would lead to a fast O\H rotation towards the neighbor oxygen [50]. Therefore, our aim was to compare the binding strength of different O centers in the three possible configurations and to examine the role of the Y uptake in the BZ on the O\H bonds discriminating between oxygen centers in different local and long-range environments. In Fig. 3 are presented the binding energies, ΔE (H) per H ion, computed as ΔE(H) = 1 / n[E(BZHY) − E(BZY) − n.E(H)] with n being the number of H ions. The H-binding energies have been computed for the fully relaxed structures. These results show that the preferred cites for H-binding are O2, independently on the Y content. The stronger stabilization of hydrogen at O2 sites correlates well with previous DFT-GGA calculations reporting the Zr\OH\Y structures to be more stable than the Zr\OH\Zr ones in BZHY12.5 supercell models [23,25]. It suffice that one of the hydrogen binds to O2 site in BZHY25 or BZHY37.5 to increase ΔE(H). As follows, ΔE(H) is spread in the interval of − 1.5 ÷ − 2.17 eV for the H\O1\bonds; in the interval of − 2.40 ÷ − 2.68 eV for the HO1\HO2 or HO1\HO1\HO2 configurations; and in the interval of − 2.1 ÷ − 3.12 eV for the HO2-bindings. In most of the structures, a stabilization is obtained due to the formation of linear or nearly linear (OH⋯O angle > 160°) short H-bonds (OH⋯O length of ~ 1.7 ÷ 1.9 Å). Indeed, increasing ΔE(H) in BZHY37.5 in comparison to the lower Y percentages can be associated with the formation of three H-bonds due to the insertion of three hydrogens. As expected, there is an influence of the O\H distance separations if more than one H is inserted. For instance, insertion of two hydrogens to the O2 centers in different planes in BZHY25 produces stronger H-binding energy compared to that in structures with H\O2 bonds in the same plane. The higher repulsion between closely spaced H atoms decreases ΔE(H). The mutual H-orientation (in the BZHY25) matters only if the two O\H bonds are in the same octahedron: the anti-parallel H-orientations led to more stable structures than the parallel ones (Scheme 2, bottom).

74.90

(×16) (×18)

74.59 75.06

(×17) (×15) (×2) (×11) (×3)

BZH37.5-M1

O2\O1\O1

O1: 0.993 O1: 0.997 O2: 0.988

2.072 ± 0.05 (×10) 2.327 2.147 ± 0.03 (×4)

BZH37.5-M1

O2\O2\O2

0.994 1.003 0.989

2.062 ± 0.03 (×4) 2.156 ± 0.03 (×5) 2.243 ± 0.02 (×3) 2.380

2.747 3.067

73.23

2.736 2.673 2.505 2.825 2.628 2.879 3.051 2.619 2.686 2.614 2.772 2.670 2.688 2.433

74.32

7

73.25

71.89

ΔE (H) (eV)

-3,2 -2,7 -2,2 -1,7

P1 -O2-

-P1-p

-P2-p

'

2-O2 M1-O

2-O2 M2-O

-P1-a 2-O2

2-O2 M2-O

-P1-a

m O2-P

2-O2

M2-O

M1-O

P1 -O1-

-P1-p

'

2-O1 M1-O

-P1-a

'

2-O1 M1-O

P1 -O1-

2-O1 M1-O

erp -P1-p

1-O1 M1-O

1-O1

-P1-a

-P2-p

M1-O

1-O1 M1-O

-P1-a

1-O1

mmm

1-O1

M2-O

M1-O

O1-P

O1-P

m

-1,2

Fig. 3. O\H binding energies in BZHY12.5 (triangles), BZHY25 (squares) and BZHY37.5 (circles) minimum energy structures. Binding to only O1 and only O2 sites are denoted with empty and filled symbols in dark blue, respectively. Oxygen binding sites, cell symmetries and model type, M1 or M2, for Y≥25% are used to label the minimum energy BZHY structures.

8

D.Z. Sahraoui, T. Mineva / Solid State Ionics 232 (2013) 1–12

As follows from the results in Table 5 and in Fig. 3, very similar OH binding geometries and strengths have been obtained in the two possible BZHY25 models M1 and M2. This shows that there is a negligible influence of the Y arrangements and symmetry of the supercell on the H-binding. Despite the increased H binding with increasing of yttrium, our results demonstrate that the determinative factors for the hydrogen binding strengths are the local arrangements around the O centers. For oxygens between two Y ions, formation of stable H\O3 bonds was not obtained, meaning that very low probability to stabilize protonic defects between two Y-dopants is expected. One can speculate that the presently obtained limits of hydroxyl ion formations at the oxygen sites in the Y\O\Y could be also a factor related to the small hydrogen transfer in the grain-boundary regions, together with other reasons so far proposed, such as an increase of the proton migration barrier [28], charge depletion [18] or segregation of oxygen defects [29]. For comparison, we have also considered an ideal cubic BaZrO3 with one to three H atoms and protons. In all minimum energy structures, a negative charge of ~− 0.25 e − is found at H sites with a volume expansion of the cell up to ~ 25% compared to the ideal barium zirconate. The computed H interaction energies are positive, thus an O\H bond is not formed. Test calculations of inserting H + were also carried out. The geometry optimization led to very comparable structures as those where H atoms were inserted. Latter results indicate that, a hydrogen insertion into the ideal cubic BZ cell is not found with the present theoretical approach. We suggest that some structural deformation or symmetry reduction would be necessary for the O\H stabilization in the absence of acceptor-dopants in barium zirconates.

of Y content. The results showed a symmetry and volume dependence on acceptor dopant concentration and on the dopant position in the cell. Interatomic distances, atomic charges at Zr and Ba sites, crystal bands structures, and DOS spectra undergoes small variations. A more pronounced dependence on the local structural environment has been found for the atomic charges at different oxygen sites, which does not vary with increasing Y-concentrations. Oxygen defects led to strong structural distortions only in the defect near region without causing a significant modification on the averaged interatomic distances or ionic charges. Various possibilities for H insertions in the BZHY12.5, BZHY25 and BZHY37.5 models have been examined as well. Hydrogen binding strength is strongly dependent on the OH local environment. The most favored binding is found for the O2 sites, slightly increasing with yttrium uptake. In the structures with two or three acceptor dopants, BZHY25 and BZHY37.5, containing respectively two or three hydrogens, the O\H binding is determined also by the mutual H-orientations for the O\H bonds in the same plane. There is an overall tendency of increasing the H-binding energy from HO1 to HO2\HO1 to HO2\HO2, which is independent on the Y percentage in the barium zirconate. Interestingly, hydrogen insertion to the O3 site (Y\OH\Y configuration) was not obtained. This suggests a possible decrease of protonic defect concentration in doped barium zirconates with higher concentration of Y\O\Y configurations. As follows from the relative energy values, the most stable BZY25 and BZY37.5 models are those with Y\O\Y configurations. Therefore, it could be expected that increasing of Y content in the host oxide would increase the probability for formation of Y\O\Y regions.

Acknowledgments 4. Conclusion A systematic study of local and long range structural properties variation of BZY (2 × 2 × 2) supercell models and hydrogen binding parameters and energies have been provided for 12.5, 25 and 37.5%

This work was supported by the European Union (EU-NMP-IndiaCP-FP 233482 HYPOMAP). Part of the calculations was performed at the HPC resources of CINES under the allocation x2011086396 made by GENCI (Grand Equipement National de Calcul Intensif).

Appendix A. Optimized structures of hydrogen containing Y-doped BaZrO3 The optimized minimum energy structures of Y-doped BaZrO3 (2 × 2 × 2) supercell models containing protonic defects are presented in Figs. A1–A5. Models with 12.5, 25 and 37.5% of Y contents are shown in Fig. A1, Figs. A2–A4 and Fig. A5, respectively. The Y\O distances are also reported in the Figures. For the other details see Section 3.3.

O1-Pm

O1-Pmmm

O2-Pm

Fig. A1. Optimized minimum energy BZHY12.5 structures. The green, red, gray, cyan and blue colored balls denote yttrium, oxygen, zirconium, barium and hydrogen ions. The Y\O distances in Å are also reported.

D.Z. Sahraoui, T. Mineva / Solid State Ionics 232 (2013) 1–12

M1-O1-O1-P1-a

M2-O1-O1-P1-a

M1-O1-O1-P2-p

M1-O1-O1-P1-perp

9

Fig. A2. Optimized minimum energy BZHY25 structures with two hydrogens bound to O1 sites (labeled as O1\O1). The model structure (M1/M2), crystal symmetries (P1/P2) and hydrogen mutual orientations (a = anti-parallel; p = parallel and perp = perpendicular) are also indicated in the structure names. The planes containing OH bonds are shown next to the 3D structure. The green, red, gray, cyan and blue colored balls denote yttrium, oxygen, zirconium, barium and hydrogen ions, respectively. Characteristic Y\O distances in Å are also reported.

M1-O2-O1-P1-a'

M1-O2-O1-P1-p' Fig. A3. Optimized minimum energy BZHY25 structures with one hydrogen bound to O2 and the next one — to O1 sites (labeled as O2\O1). The model structure (M1), crystal symmetries (P1) and hydrogen mutual orientations (a = anti-parallel and p = parallel) are also indicated in the structure names. The planes containing OH bonds are shown next to a 3D structure. The prime symbol is used to distinguish for the OH bonds in two different planes. The green, red, gray, cyan and blue colored balls denote yttrium, oxygen, zirconium, barium and hydrogen ions, respectively. Characteristic Y\O distances in Å are also reported in blue for the YO6\H octahedrons and in black for the YO6 octahedrons.

10

D.Z. Sahraoui, T. Mineva / Solid State Ionics 232 (2013) 1–12

M1-O2-O2-P1-a

M2-O2-O2-P1-a

M2-O2-O2-P2-p

M2-O2-O2-P1-p' Fig. A4. Optimized minimum energy BZHY25 structures with two hydrogens bound to O2 sites (labeled as O2–O2). The model structure (M1/M2), crystal symmetries (P1/P2) and hydrogen mutual orientations (a = antiparallel and p = parallel) are also indicated in the structure names. The planes containing OH bonds are shown next to the 3D structure. The prime symbol is used to distinguish for the OH bonds in two different planes. The green, red, gray, cyan and blue colored balls denote yttrium, oxygen, zirconium, barium and hydrogen ions, respectively. Characteristic Y\O distances in Å are also reported in blue for the YO6\H octahedrons and in black for the YO6 octahedrons.

D.Z. Sahraoui, T. Mineva / Solid State Ionics 232 (2013) 1–12

11

M1-O1-O1-O1-P1

M1-O2-O1-O1-P1

M1-O2-O2-O2-P1 Fig. A5. Optimized minimum energy BZHY37.5 structures with three hydrogens. The model structure (M1), crystal symmetries (P1/P2) and type of oxygen binding centers are used to form the name of the structures. Two of the O\H bonds are displaced in the same plane and the third one is in a perpendicular plane in M1\O1\O1\O1\P1 and M1\O2\O2\O2\P1, and in a parallel plane in M1\O2\O1\O1\P1. The green, red, gray, cyan and blue colored balls denote yttrium, oxygen, zirconium, barium and hydrogen ions, respectively. Characteristic Y\O distances in Å are also reported in blue for the YO6\H octahedrons and in black for the YO6 octahedrons.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

K.D. Kreuer, Annu. Rev. Mater. Res. 33 (2003) 333–359. K.D. Kreuer, Solid State Ionics 125 (1999) 285–302. T. Schober, H.G. Bohn, Solid State Ionics 127 (2000) 351–360. H.G. Bohn, T. Schober, J. Am. Ceram. Soc. 83 (2000) 768–772. K.D. Kreuer, St. Adams, W. Munch, A. Fuchs, U. Klock, J. Maier, Solid State Ionics 145 (2001) 295–306. Y. Yamazaki, P. Babilo, S.M. Haile, Chem. Mater. 20 (2008) 6352–6357. F. Giannici, A. Longo, A. Balerna, K.D. Kreuer, A. Martorana, Chem. Mater. 21 (2009) 2641–2649. F. Giannici, A. Longo, K.D. Kreuer, A. Balerna, A. Martorana, Solid State Ionics 181 (2010) 122–125. M. Karlsson, I. Ahmed, A. Matic, S.G. Eriksson, Solid State Ionics 181 (2010) 126–129. F. Giannici, M. Shirpour, A. Longo, A. Martorana, R. Merkle, J. Maier, Chem. Mater. 23 (2011) 2994–3002. A.K. Azad, Ch. Savaniu, S. Tao, S. Duval, P. Holtappels, R.M. Ibberson, J.T.S. Irvine, J. Mater. Chem. 18 (2008) 3414–3418. Y. Jiang, S.V. Bhide, A.V. Vikar, J. Solid State Chem. 127 (2001) 149–159. B. Gross, Ch. Beck, F. Meyer, Th. Krajewski, R. Hempelmann, H. Altgeld, Solid State Ionics 145 (2001) 325–331. U. Anselmi-Tamburini, M.T. Buscaglia, M. Viviani, M. Bassoli, C. Bottino, V. Bascaglia, P. Nanni, Z.A. Munir, J. Eur. Ceram. Soc. 26 (2006) 2313–2318. P.A. Stuart, T. Unno, R. Ayres-Rocha, E. Djurado, S.J. Skinner, J. Eur. Ceram. Soc. 29 (2009) 697–702.

[16] F. Iguchi, T. Tsurui, N. Sata, Y. Nagao, H. Yugami, Solid State Ionics 180 (2009) 563–568. [17] F. Iguchi, N. Sata, H. Yugami, J. Mater. Chem. 20 (2010) 6265–6270. [18] M. Shirpour, B. Rahmati, W. Sigle, P.A. van Aken, R. Merkle, J. Maier, J. Phys. Chem. C 116 (2012) 2453–2461. [19] R.B. Cervera, Y. Oyama, S. Miyoshi, K. Kobayashi, T. Yagi, S. Yamaguchi, Solid State Ionics 179 (2008) 236–242. [20] M. Karlsson, A. Matic, C.S. Knee, I. Ahmed, S.G. Eriksson, L. Börjesson, Chem. Mater. 20 (2008) 3480–3486. [21] M. Karlsson, A. Matic, S.F. Parker, I. Ahmed, L. Börjesson, S. Eriksson, Phys. Rev. B 77 (2008) 104302–106308. [22] P.G. Sundell, M.E. Björketun, G. Wahnström, Phys. Rev. B 73 (2006) 104112–104122. [23] M.E. Björketun, P.G. Sundell, G. Wahnström, Phys. Rev. B 76 (2007) 054307–054316. [24] M.A. Gomez, M. Chunduru, L. Chigweshe, L. Foster, S.J. Fensin, K.M. Fletcher, L.E. Fernandez, J. Chem. Phys. 132 (2010) 214709-1–214709-8. [25] B. Merinov, W. Goddard III, J. Chem. Phys. 130 (2009) 194707-1–194707-6. [26] A.V. Bandura, R.A. Evarestov, D.D. Kuruch, Surf. Sci. 604 (2010) 1591–1597. [27] R.A. Evarestov, A.V. Bandura, Solid State Ionics 188 (2011) 25–30. [28] D.-H. Kim, B.-K. Kim, Y.-C. Kim, Solid State Ionics 213 (2012) 18–21. [29] B.J. Nyman, E.E. Helgee, G. Wahnström, Appl. Phys. Lett. 100 (2012) 061903-1–061903-4. [30] Q. Zhang, G. Wahnström, M.E. Björketun, S. Gao, E. Wang, Phys. Rev. Lett. 101 (2008) 215902-1–215902-4. [31] V.R. Saunders, R. Dovesi, C. Roetti, R. Orlando, C.M. Zicovich-Wilson, F.N.M. Pascale, B. Civalleri, K. Doll, N.M. Harrison, I.J. Bush, P. D'Arco, M. Llunell, in: CRYSTAL 06 user's manual, Universita' degli Studi di Torino, Torino, Italy, 2008. [32] A.D. Becke, J. Chem. Phys. 98 (1993) 5648–5652.

12 [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]

D.Z. Sahraoui, T. Mineva / Solid State Ionics 232 (2013) 1–12 C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785–789. S. Piskunov, E. Heifets, R.I. Eglitis, G. Borstel, Comput. Mater. Sci. 29 (2004) 165–178. http://www.crystal.unito.it/Basis_Sets/zirconium.html. A. Buljan, P. Alemany, E. Ruiz, J. Phys. Chem. B 103 (1999) 8060–8066. L. Valenzano, F.J. Torres, K. Doll, F. Pascale, C.M. Zicovich-Wilson, R. Dovesi, Z. Phys. Chem. 220 (2006) 893–912. C. Gatti, V.R. Saunders, C. Roetti, J. Chem. Phys. 101 (1994) 10686–10696. S.H. Vosko, L. Wilk, M. Nusair, Can. J. Phys. 58 (1980) 1200–1211. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. J.P. Perdew, Electronic Structure of Solids 1991, Akademie Verlag, Berlin, 1991. J.P. Perdew, Wang Yue, Phys. Rev. B 33 (1986) 8800–8802. R. Khenata, M. Sahnoun, H. Baltache, M. Rerat, A.H. Rashek, N. Iles, B. Bouhafs, Solid State Commun. 136 (2005) 120–125.

[44] N. Iles, A. Kellou, K. Driss Khodja, B. Amrani, F. Lemoigno, D. Bourbie, H. Aourag, Comput. Mater. Sci. 39 (2007) 896–902. [45] W. Pies, A. Weiss, in: K.-H. Hellweg, A.M. Hellwege (Eds.), Landolt–Börnstein: Numerical Data and Functional Relationships in Science and Technology, New Series, Group III, vols. 7b1, Springer, Berlin, 1975, (and 7e Springer, Berlin, 1975). [46] J. Robertson, J. Vac. Sci. Technol. B 18 (2000) 1785–1791. [47] A. Bilić, J.D. Gale, Phys. Rev. B 79 (2009) 174107-1–174107-9. [48] C.H. Perry, D.J. Mccarthy, G. Rupprech, Phys. Rev. 138 (1965) A1537–A1538. [49] H.J. Monkhorst, J.D. Pack, Phys. Rev. 13 (1976) 5188–5192. [50] K.D. Kreuer, Solid State Ionics 136–137 (2000) 149–160. [51] M.E. Björketun, P.G. Sundell, G. Wahnström, Faraday Discuss. 134 (2007) 247–265.