First-Principles Calculations for the Energetics of ... - ACS Publications

Feb 7, 2017 - Hiroshige Matsumoto,. ⊥. Hitoshi Takamura,*,† and Shu Yamaguchi*,#. †. Department of Materials Science, Graduate School of Enginee...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/cm

First-Principles Calculations for the Energetics of the Hydration Reaction of Acceptor-Doped BaZrO3 Hiroki Takahashi,*,†,‡ Isamu Yashima,‡ Koji Amezawa,§ Koichi Eguchi,∥ Hiroshige Matsumoto,⊥ Hitoshi Takamura,*,† and Shu Yamaguchi*,# †

Department of Materials Science, Graduate School of Engineering, Tohoku University, 6-6-02 Aramaki Aza Aoba, Sendai, Miyagi 980-8579, Japan ‡ Mitsui Mining & Smelting Co., Ltd., 1333-2 Haraichi, Ageo, Saitama 362-0021, Japan § Institute of Multidisciplinary Research for Advanced Materials, Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai, Miyagi 980-8577, Japan ∥ Department of Energy and Hydrocarbon Chemistry, Graduate School of Engineering, Kyoto University, Nishikyo-ku, Kyoto, 615-8510, Japan ⊥ International Institute for Carbon-Neutral Energy Research (I2CNER), Kyushu University, 744 Motooka, Nishi-ku, Fukuoka, 819-0395, Japan # Department of Materials Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan S Supporting Information *

ABSTRACT: The fundamental factors that influence the hydration of BaZrO3 (BZO) doped with trivalent cation M3+ (Al, Sc, Ga, Y, In, and Lu) for proton conductors were investigated by means of density functional theory calculations which take the configuration of complex defects into account. The creation of oxygen vacancies is favored for Aland Ga-doped BZOs and leads to small hydration energies with stable proton sites at the nearest neighbor (1NN). Meanwhile, Y-, In-, and Ludoped BZOs prefer protons at the second nearest neighbor (2NN). The stability of those defects can be formulated in the context of the energies of oxygen vacancy formation and hydration. BZOs with larger dopants gain more hydration energy by structural relaxation with protons located at 2NN. By isolating the associated complex defects, it is possible to increase the negative hydration energy, which in effect improves the degree of hydration of BZOs.



INTRODUCTION Proton conductors have been attracting much attention as potential electrolytes in intermediate temperature solid oxide fuel cells (IT-SOFCs). A perovskite-type oxide is one of the candidate materials for proton conductors.1−3 In this oxide, represented by the general formula ABO3, the substitution of an M3+ trivalent cation for the B4+ tetravalent cation leads to the creation of oxygen vacancies. Not only is the M3+ cation on the B4+ site, which is regarded as an acceptor (M′B), compensated by the formation of an oxygen vacancy, but also compensation for the electronic holes occurs depending on atmosphere and temperature. This reaction is expressed by the following equation using Kröger-Vink notation: × OO + 2h• = V O•• + 1/2O2 (g)

In this reaction, the H2O molecule reacts with an oxygen vacancy and a host oxide ion, resulting in the generation of two OH groups, generally referred to as protonic defects. These protonic defects migrate from one oxide ion to the adjacent oxide ion in a process known as protonic conduction. Among a number of perovskite-type oxides, BaZrO3 (BZO) exhibits high proton conductivity when doped with Y and has chemical stability against water vapor and CO2.4−6 BZO is one of the promising proton conductors for use in IT-SOFCs.3 To achieve higher proton conductivity, various trivalent dopants have been studied.4,7,8 The proton conductivity σ can be simply written as a function of the proton concentration, c: σ = q×μ×c

where μ denotes proton mobility and q is the electrical charge. Here, proton concentration, c, relies on the thermodynamics of hydration. One of the solutions to improve the proton

(1)

If the formation of an oxygen vacancy is energetically favored, the oxygen vacancies formed can be used for a subsequent hydration reaction according to the following: H 2O +

× OO

+

V •• O

=

2OH•O © XXXX American Chemical Society

(3)

Received: September 14, 2016 Revised: February 6, 2017 Published: February 7, 2017

(2) A

DOI: 10.1021/acs.chemmater.6b03907 Chem. Mater. XXXX, XXX, XXX−XXX

Article

Chemistry of Materials

Figure 1. Schematic representation of (2M′Zr, V•• O ) configuration for the calculation of total energy. M1 and M2 denote two doped trivalent cations. The V•• O site is equidistant and the nearest from the two M cations in each M−M configuration (1−5NN). 2NN and 3NN sites locate on a plane with different z-coordinates of 1/8 and 1/4, respectively. hydration energy (Ehydr). Those energies can be evaluated by the following equations, respectively:

conductivity of BZOs is to increase the concentration of protonic defects, which in effect increases the degree of hydration. Even though various techniques, including thermogravimetry and nuclear magnetic resonance (NMR), have been developed to clarify the degree of hydration in the proton conductors,9−11 the slow kinetics of hydration make it difficult to analyze the degree of hydration with high accuracy. This is especially true at low temperatures below 300 °C. As such, computational simulations are useful to investigate the effect of structural and chemical factors on the hydration reaction of BZOs under ideal conditions. To precisely evaluate the hydration reaction of acceptordoped BZOs, the energetics of at least two distinct reactions needs to be systematically studied; i.e., the first is the formation of oxygen vacancies, and the second is their hydration.12 In addition to these reactions, the local structural configurations of the dopants, oxygen vacancies, and protons, which potentially affect hydration properties, have to be taken into consideration. Hiraiwa et al. have reported that the local configurations of dopants and oxygen vacancies have an effect on structural changes of Y-doped BZO, but the effect on the hydration properties was not discussed.13 Dawson et al. have reported the effect of the configurations of protons and some dopants on the binding energy of protons. The relationship between the formation of oxygen vacancies and hydration, however, has not been thoroughly elucidated.14,15 Thus, while the hydration properties of acceptor-doped BZOs have indeed been studied by computational simulation, the focus has been on basically one of the three factors mentioned above.16−21 To improve the hydration properties of acceptor-doped BZOs, however, all the three factors need to be taken into consideration, simultaneously. In this study, we calculate both the formation energy of oxygen vacancies and the hydration energy for acceptordoped BZOs and take the local structural configurations of complex defects, which consist of dopants, oxygen vacancies, and protons, into account using first-principles calculations to clarify the factors which control the hydration property of BZOs.



× Evac = Et(2M′Zr , V •• O ) + 1/2μO − Et (2M′Zr , OO) 2

(4)

× E hydr = Et(2M′Zr , 2OH•O) − Et(2M′Zr , OO , V •• O ) − Et (H 2O)

(5) where Et denotes the total energy taking defects and species in parentheses into account and μO2 denotes the chemical potential of the O2 molecule. Here, Et(2M′Zr, O×O) implies the coexistence of two holes in the system, as shown in eq 1. These values were calculated using the CASTEP code, which is based on the density functional theory (DFT).22 DFT calculations were performed with plane-wave expansions using the generalized gradient approximation in the Perdew-Burke-Emzerhof form (GGA-PBE). The ionic cores were represented by ultrasoft pseudopotentials. The plane-wave cutoff energy was set to be 380 eV, and a 4 × 4 × 4 supercell of BaZrO3 unit cell was used with a single k-point sampling at Γ-point. Geometry optimizations were performed by relaxing all atom positions and lattice volumes under the following conditions: all the residual forces were smaller than 0.03 eV/Å, and the stress was smaller than 0.05 GPa with an energy convergence of 10−9 eV/atom for self-consistent calculations. Density of states (DOS) was calculated with 0.1 eV smearing width and 3 × 3 × 3 k-point sampling. In this study, the chemical potential of the O2 molecule (μO2) was calculated with the oxygen gas partial pressure and temperature taken into account.23,24 Here, the oxygen gas partial pressure was set to be 0.21 atm on the assumption that BZOs are equilibrated in air. Complex Defect Modeling. In addition to substituting two trivalent cations for two Zr cations, one oxygen vacancy was introduced into each supercell to maintain overall electroneutrality. Al, Sc, Ga, Y, In, and Lu were selected as the trivalent acceptor dopants. Because these dopants have a closed-shell electronic structure, it is possible to calculate the systems as nonmagnetic and it is not necessary to take the spin degrees of freedom into account. To ′ ) and the consider the interaction between the trivalent cation (MZr oxygen vacancy (VO••), geometry optimization calculations were performed by changing the initial distance between two M cations from the nearest neighbor (1NN) to the fifth nearest neighbor (5NN). For the sake of simplicity, V•• O was placed so that it was equidistant and the nearest from the two M cations in each M−M configuration (1− 5NN) as illustrated in Figure 1. Additionally, geometry optimization calculations for other configurations with longer M′Zr−V•• O distance as shown in Figure S1 were also performed. The formation energy of oxygen vacancies was calculated for the most stable configuration of ′ , V•• the complex defect (2MZr O ).

METHODS

Energy Calculations. According to the above eqs 1 and 2, the energetics of hydration of acceptor-doped BZOs can be divided into two energies: the formation energy of oxygen vacancies (Evac) and the B

DOI: 10.1021/acs.chemmater.6b03907 Chem. Mater. XXXX, XXX, XXX−XXX

Article

Chemistry of Materials In the case of the hydrated structure, as shown in eq 2, water vapor reacts with an oxygen atom and an oxygen vacancy, and then, two protons are incorporated as a form of OH. Here, there appear to be several possible configurations of the proton site. It has been reported that the proton occupies the nearest neighboring (1NN) oxygen to Sc in Sc-doped BZO,25 whereas in the case of Y-doped BZO, it occupies the second nearest neighboring (2NN) oxygen to Y.26,27 For this reason, the sites of two protons were explored within the range of 1NN and 2NN oxygen for each dopant. The total energy of each configuration was then calculated to evaluate the most stable configuration as a whole system. The same calculation was also performed on cubic yttrium stabilized zirconia (YSZ). Even though YSZ is comprised of similar elements to Y-doped BZO, YSZ shows no proton conductivity due to negligible proton concentration.28−30 This comparison with YSZ clarifies the important factors responsible for the increase in the proton concentration of BZO.



RESULTS Optimization of Defect Configuration. Prior to calculating the energies related to the hydration reaction of acceptor-doped BZOs, their total energy was calculated with respect to various defect configurations. The computationally derived lattice constant of a nondoped BZO unit cell was calculated as a = 4.229 Å, while the experimentally obtained lattice constant was reported as a = 4.194 Å.31 The calculated value is in good agreement with the experimental value within 1% error, indicating the quality of this calculation condition is such that the structure of BZO can be reproduced. Moreover, the convergence test of the hydration energy with respect to cutoff energy was performed in the case of Sc-doped BZO as an example. The hydration energy was converged within 0.01 eV when a 380 eV cutoff energy was used (Figure S2). This cutoff energy ensures sufficient energy convergence to compare the hydration energy of various doped BZOs. The total energy of BZO with complex defects depends on the local configurations of two trivalent dopants and one oxygen vacancy (2M′Zr, V•• O ). The relative total energy with respect to the distance between dopants M is shown in Figures 2a and S3. The most stable configuration is obviously the case of the shortest M′Zr−M′Zr distance in all trivalent dopants: this configuration is expressed as (MZr ′ −V•• O −M′Zr), where the dopant M and oxygen vacancy are adjacent, as shown in Figure 3a. In other words, this associated complex defect is the most stable form among the representative configurations. The weak interaction between complex defects results in the convergence of the relative total energy when the M′Zr−M′Zr distance exceeds 2NN. The defect is considered isolated when the configuration is 5NN. It should be noted that this distance is the longest of all configurations investigated in this study. The association energy for the dopant and oxygen vacancy is defined as the difference in the relative total energy between the isolated configuration (5NN) and the associated configuration (1NN). The association energy for the dopants calculated is shown in Figure 2b. Y-doped BZO has the smallest association energy, and those of Al- and Ga-doped BZOs are larger. The association energy shown in Figure 2b can be compared with the results reported by Sundell et al.18 For that purpose, ΔEsingle = −(2ΔEM‑Vo + ΔEM‑M) can be calculated using the data in ref 18, since their interaction energy (ΔEM‑M and ΔEM‑Vo) is for single defect pairs. The results are summarized in Table 1 together with our association energy (ΔEcomplex) for the (2M′Zr, V•• O ) complex defects. The trend of ΔEcomplex and ΔEsingle with respect to dopant is the same, and the two energies for Scdoped BZOs are in good agreement with each other.

Figure 2. (a) Relative total energy of M-doped BZOs (M = Al, Sc, Ga, Y, In, and Lu) with (2M′Zr, V•• O ) complex defect with respect to M−M configurations (1−5NN). Each configuration corresponds to Figure 1. (b) The association energy of M-doped BZOs with the complex defect.

Figure 3. Relaxed structures of (a) Y-doped BZO and (b) YSZ with complex defects. Both (a) and (b) correspond to the most stable configuration. Ba atoms are omitted for clarity. Crystal structure is visualized using VESTA.38

Meanwhile, as for In- and Y-doped BZOs, ΔEcomplex is smaller than ΔEsingle. This implies that the association between dopants and an oxygen vacancy in the complex defect are weaker than that in single defect pairs. This weaker association can be attributed to the repulsive interaction between cations. The repulsive effect becomes more significant in the complex defect containing dopants when the ionic radius is larger. The opposite is also true for the case of Ga-doped BZO with smaller dopant of Ga. The association between dopants and an oxygen vacancy in the complex defect becomes stronger than that in C

DOI: 10.1021/acs.chemmater.6b03907 Chem. Mater. XXXX, XXX, XXX−XXX

Article

Chemistry of Materials Table 1. Association Energy of (2MZr ′ , V•• O ) Complex Defect (ΔEcomplex) and Corresponding Association Energy Calculated from Single Defect Pair (ΔEsingle)a

a

dopant M

ΔEcomplex

ΔEsingle

Sc Ga Y In

0.73 1.52 0.31 0.78

0.72 1.33 0.64 1.05

of doped BZOs to reduce the computational costs. In the case of Sc-, Ga-, and In-doped BZOs, no imaginary frequencies emerged for the vibrational calculations. This implies that the calculated free energies were accurate for the Sc-, Ga-, and Indoped BZOs.34,35 Meanwhile, imaginary frequencies emerged for Al-, Lu-, and Y-doped BZOs under the same conditions. Bilić and Gale indicated that ultrasoft pseudopotentials, especially semicore states of Ba, affect the stability of BZO.36 This suggests that the stability of doped BZOs could be evaluated more precisely if the pseudopotentials used are improved. Even though the calculation of vibrational contribution to the total energy should be further improved and taken into account for comparison purposes, the vibrational contribution, i.e., difference between the solid and dashed line in Figure 4, is almost identical regardless of the doping elements. In this study, therefore, the formation energy of oxygen vacancies at 0 K (without vibrational contribution) will be used for comparison. The formation energy of oxygen vacancies computed by Sundell et al. was reported to be 1.21 18 eV for V•• Because our results O as an acceptor-doped BZO. range between 0.02 and 0.76 eV, which is smaller than those reported by Sundell et al., it can be concluded that the oxygen vacancies are readily formed by acceptor doping. In an earlier study of acceptor-doped BZOs, Bévillon et al. reported that the oxygen vacancies can be easily formed in the order Sc (most difficult case), Y, In, and Ga (the easiest case).37 With the exception of Y-doped BZO, this trend is in good agreement with our results. The discrepancy in Y-doped BZO may be because the concentration of the dopants in Bévillon’s model (25 at %) was considerably higher than ours (3.1 at %). The large structural relaxation, resulting from high concentrations of large dopants such as Y, has a significant effect on the total energy of the relaxed structure. In the case of YSZ with a fluorite-type structure, the most stable configuration of complex defects was one where they are adjacent to each other, in the form (Y′Zr−V•• O −Y′Zr). This was also the case with acceptor-doped BZOs. The significantly negatively large VO•• formation energy of YSZ for this configuration, as shown in Figure 4, indicates that the oxygen vacancies in YSZ are quite stable. The optimized structure of YSZ is also shown in Figure 3b. During the structure relaxation process, the oxygen atom adjacent to V•• O moves very close to the original V•• O site. This relaxation accompanied by the large diffusion of the oxygen atom is distinctive to YSZ. Hydration Energy. As noted in the section of complex defect modeling, there are several configurations of protons in hydrated BZOs. Figure 5a shows the schematic representation of the hydration process that can be distinguished by proton (H) sites relative to two M atoms in the associated defects (2 types). In addition to the stable associated (2MZr ′ , V•• O) configuration, the hydration process for the isolated configurations (3 types) is also taken into account, as shown in Figure 5b. Table 2 shows the hydration energy of each hydrated type for (a) associated and (b) isolated configurations. In the case of the associated configuration, type A-1 (1NN, 1NN) is more stable for Sc-doped BZO than type A-2 (1NN, 2NN), whereas type A-2 is more stable for Y-doped BZO. Figure 6 shows the geometry optimized hydrated structures of Sc-doped (type A1) and Y-doped (type A-2) BZOs. For Y-doped BZO (Figure 6b), the distortions of YO6 and ZrO6 octahedra with protons is larger than those for Sc-doped BZO (Figure 6a). Meanwhile, in the case of the isolated configurations, type B-3 (2NN, 2NN) is the most stable for Y-doped BZO among all series of type B,

Unit: eV.

single defect pairs with negligibly small repulsive interaction between the dopants (0.03 eV for Ga in ref 18). Formation Energy of Oxygen Vacancies. The formation energy of oxygen vacancies in acceptor-doped BZOs was then evaluated using the most stable (M′Zr−V•• O −M′Zr) configuration as a function of temperature, as shown in Figure 4. The

Figure 4. Formation energy of oxygen vacancies of M-doped BZOs and YSZ (M = Al, Sc, Ga, Y, In, and Lu). The dashed and solid lines correspond to the calculation condition with and without phonon vibrational contribution, respectively.

temperature dependence of the formation energy is determined by the oxygen gas chemical potential shown in eq 4. The formation energy negatively increases as the temperature increases. All dopants have negative formation energy above 800 K (see solid lines in Figure 4). This implies that the formation of oxygen vacancies is energetically favored above this temperature. It should also be noted that Al-, Ga-, and Ydoped BZOs can readily form oxygen vacancies. To evaluate the formation energy more precisely, the phonon vibrational contribution to the free energy of BZOs is taken into account as shown below:32,33 •• Evac = Et(2M′Zr , V •• O ) + E vib(2M′Zr , V O ) × − TSvib(2M′Zr , V •• O ) + 1/2μO − {Et (2M′Zr , OO) 2

+

× Evib(2M′Zr , OO )

× − TSvib(2M′Zr , OO )}

(6)

where Evib denotes the sum of the zero point vibrational energy and the temperature dependence of enthalpy and Svib denotes the vibrational contribution to entropy. The phonon vibrational calculations were performed with a smaller 2 × 2 × 2 supercell D

DOI: 10.1021/acs.chemmater.6b03907 Chem. Mater. XXXX, XXX, XXX−XXX

Article

Chemistry of Materials

Figure 6. Relaxed hydrated structures of (a) Sc-doped BZO with type A-1 and (b) Y-doped BZO with type A-2.

doped BZO using DFT calculations with 2 × 2 × 2 supercell (see ref 39). The most stable configuration reported in the previous study is in good agreement with the result in this study (type A-2), even though the dopant concentration of type A-2 is different with that in the previous study. This implies that it is not the dopant concentration but the local structural configuration around the dopant that determines the hydration energy. Meanwhile, in the case of the surface reaction, the reported dissociative adsorption energy of a water molecule on the BZO surface, at −1.09 eV,40 is negatively larger than our result. This can be attributed to the difference in the stability between protons adsorbed on a surface and protons in a bulk. The hydration energy of more stable configurations in the associated and isolated conditions of the complex defect is shown in Figure 7. In the case of the associated configuration of

Figure 5. Schematic representation of the hydration process by the reaction of H2O molecule. (a) Associated (type A) and (b) isolated ′ , V•• (type B) configurations of complex defects (2MZr O ) before hydration. 1NN and 2NN denote H atom site at the first and second nearest neighbor from dopant M, respectively.

Table 2. Hydration Energy of Each Hydrated Type in the (a) a Associated and (b) Isolated (2M′Zr, V•• O ) Configurations (a) Associated Type A-2 (1NN, 2NN)

ΔEhydr (=EA‑1 − EA‑2)

dopant M

Type A-1 (1NN, 1NN)

Al Sc Ga Y In Lu

−0.34 −0.69 −0.19 0.23 −0.23 −0.24

dopant M

Type B-1 (1NN, 2NN)

Type B-2 (1NN, 2NN)

Type B-3 (2NN, 2NN)

Al Sc Ga Y In Lu

−0.97 −0.78 −0.90 −0.35 −0.64 −0.66

−1.12 −0.93 −1.04 −0.52 −0.81 −0.83

−0.65 −0.59 −0.63 −0.61 −0.62 −0.88

0.24 −0.53 0.16 −0.22 −0.31 −0.38 (b) Isolated

−0.58 −0.16 −0.35 0.45 0.08 0.14

Figure 7. Hydration energy of M-doped BZO and YSZ with the associated and isolated configurations of complex defect (M = Al, Sc, Ga, Y, In, and Lu).

a Unit: eV. Configuration types correspond to the representation in Figure 5.

BZOs, Sc-doped BZO has the largest negative hydration energy, whereas Y- and Ga-doped BZO have sufficient but smaller negative hydration energies. It should be noted that the hydration energy of the isolated configuration is negatively larger than that of the associated configuration because of the difference in the stability of (2MZr ′ , V•• O ) configuration without hydration, as shown in Figure 2a. This implies that isolated configurations, which are not energetically favored (as shown in Figure 2a), are more desirable from the perspective that protons are more stable with larger negative hydration energy.

whereas the type B-2 (1NN, 2NN) is stable for Sc-doped BZO. This difference in stable proton sites between Sc- and Y-doped BZOs is in good agreement with reported experimental results.11,25,26 Furthermore, the trend of the stable configuration of a proton and dopant reported by Dawson and Tanaka, based on the difference in the total energy,15 is in good agreement with our result. In the case of plural protons, the stability of a system with two protons has been studied for YE

DOI: 10.1021/acs.chemmater.6b03907 Chem. Mater. XXXX, XXX, XXX−XXX

Chemistry of Materials

Article



DISCUSSION The role of dopants on the oxygen vacancy formation and hydration and their relationships are discussed in terms of the energies calculated. First of all, it is worth noting that the oxygen vacancy formation among Sc-, Ga-, and In-doped BZOs varies considerably (Figure 4). The large negative formation energy of oxygen vacancies for Ga-doped BZO predisposes this BZO to the creation of oxygen vacancies. Figure 9 shows the

That is, these results suggest that the higher negative hydration energy of BZOs makes them better hydrated and, therefore, more suitable for high temperature applications by the increase of the dehydration temperature. On the other hand, it has been pointed out that a higher negative hydration energy causes proton binding/trapping (see in eq 5) and results in a decrease in mobility.15,27,41 However, it is possible to compensate for the decrease in mobility by reducing the thickness of BZO specimens. It needs to be noted that, since hydration energy, i.e., proton concentration, is a thermodynamic property, it is necessary for BZOs to have a sufficient hydration energy to maintain sufficient proton concentration at elevated temperatures. In fact, the decrease in the transport number of protons is the focus of considerable attention recently as a key factor in the degradation of the performance of proton-conductor-based SOFCs.42 One of the reasons for this is the desorption of protons from the surface. This suggests that the increase in the hydration energy leads to an improvement in the total performance of the SOFC by keeping protons in a bulk. Meanwhile, even at low temperature applications, BZOs with larger negative hydration energy can remain well-hydrated. It should also be noted that, in the case of YSZ, the positive hydration energy in any configuration suggests that hydration does not occur at all (Figure 7). It is because oxygen vacancies in YSZ are too stable to cause the hydration of YSZ (Figure 4). The hydration energy shown above corresponds to the full hydration state, as can be seen in eq 5. We also calculated the hydration energy as a function of the degree of hydration for Sc-doped BZO. A 4 × 4 × 4 supercell with 8 Sc and 4 V•• O (12.5 at % Sc) was used to deal with the intermediate hydration state, with Sc atoms not located within 2NN of each other, i.e., the isolated configuration. Using this model, the degree of hydration was set to be 25%, 50%, 75%, and 100%, which correspond to 1, 2, 3, and 4 H2O molecules incorporated with oxygen vacancies, respectively. As shown in Figure 8, the hydration energy negatively increases in proportion to the hydration level. This suggests that the interaction among the defects is not significant for this dopant concentration (12.5 at %).

Figure 9. Local DOS of the O atom sandwiched between two M atoms in M-doped BZO (M = Sc-, Ga-, In-, and nondoped). The Fermi level of each system is set at 0 eV.

local DOS of the oxygen atom sandwiched between two M atoms for Sc-, Ga-, In-, and nondoped BZOs before the formation of an oxygen vacancy. Ga- and In-doped BZOs have more states at shallow energy levels (−0.3 to 0 eV from the Fermi level) than Sc- and nondoped BZOs. This nature of the DOS is consistent with the energetically favored oxygen vacancy formation shown in Figure 4. That is, the electronic states affect the formation energy of oxygen vacancies: this is particularly true of electronic states near the Fermi level with an M−O bond. We next examine the relationship between the formation energy of oxygen vacancies and the hydration energy. This relationship governs the stability of defects introduced by acceptor doping. Figure 10 shows the relationship between the formation energy of oxygen vacancies and the hydration energy of the associated (2MZr ′ , V•• O ) configuration. For Al- and Gadoped BZOs with a stable hydrated configuration (black squares), the oxygen vacancies are stable (Evac ≈ 0 eV); however, the hydration energy is negatively small (Ehydr ≈ −0.2 to −0.35 eV). Meanwhile, oxygen vacancies are not readily formed in the case of Sc-doped BZO due to its higher positive energy (Evac ≈ 0.75 eV). It should also be noted that its hydration energy is negatively large (Ehydr ≈ −0.7 eV). The same trade-off relationship can be seen for the unstable hydrated configurations (red triangles). As shown below, this trade-off relationship can be expressed as a slope of minus unity. To derive the relationship, the entire reaction from the formation of oxygen vacancies through the hydration is expressed by adding eqs 1 and 2:

Figure 8. Hydration energy as a function of the degree of the hydration for 12.5 at % Sc-doped BZO with 4 × 4 × 4 supercell.

× 2OO + 2h• + H 2O(g) = 2OH•O + 1/2O2 (g)

F

(7)

DOI: 10.1021/acs.chemmater.6b03907 Chem. Mater. XXXX, XXX, XXX−XXX

Article

Chemistry of Materials

proton tends to be trapped at 1NN due to both the Coulombic interaction between positively charged proton and negatively charged acceptor and the local distortion (size effect) around the dopant.15,17,26 Meanwhile, larger dopants, i.e., Y-, In-, and Lu-doped BZOs, prefer the type A-2 configuration, in which one of the two protons occupies the 2NN site. This change in the configuration may be attributed to an increase in the local distortion. To quantify the local distortion with respect to dopants, the differences in the M−O−Zr angle (Δθ) shown in Figure 5 and the hydration energy (ΔEhydr) between type A-1 (1NN, 1NN) and A-2 (1NN, 2NN) are introduced as follows: Δθassociated = θA‐2 − θA‐1 , ΔE hydr_associated = E hydr_A‐2 − E hydr_A‐1

(12)

The difference in the hydration energy (ΔEhydr) is plotted as a function of differences in the M−O−Zr angle (Δθ) in Figure 11. As can be seen, Y-, In-, and Lu-doped BZOs with the type Figure 10. Relationship between the formation energy of oxygen vacancies and the hydration energy with the associated configuration of complex defect. The strings in parentheses mean the type of hydrated process that corresponds to Figure 5. Iso-E′h•→H+ lines with a slope of minus unity (see eq 11) are plotted for three cases: E′h•→H+ < 0, E′h•→H+ = 0, and E′h•→H+ > 0.

Eq 7 includes the decomposition of water into hydrogen and oxygen as follows: H 2O(g) = 1/2O2 (g) + H 2(g)

(8)

Subtracting eq 8 from eq 7 yields × 2OO + 2h• + H 2(g) = 2OH•O

(9)

Here, this reaction energy Eh•→H+ is given by Eh•→ H+ = Evac + E hydr − E H2O

(10)

where EH2O denotes the reaction energy of eq 8. Given that EH2O is constant (≈ + 2.5 eV), eq 10 can be simplified as follows: E hydr = −Evac + Eh′•→ H+

Figure 11. Relationship between the difference in M−O−Zr angle and the difference in hydration energy when the proton configuration is changed from type A-1 to A-2 or from type B-2 to B-3.

(11)

where E′h•→H+ is equal to Eh →H + EH2O. Eq 11 implies two important points: (1) the stability of protons (Ehydr) and oxygen vacancies (Evac) shows a trade-off relationship with a slope of minus unity, and (2) E′h•→H+ is a good index for the prediction of proton conduction; a negatively large E′h•→H+ is essential for good oxide-based proton conductors. The relationship among Ehydr, Evac, and E′h•→H+ expressed by eq 11 is also plotted in Figure 10. Unstable hydrated types (red triangles) lie on a line with E′h•→H+ > 0; meanwhile, stable hydrated types (black squares) are on a line with E′h•→H+ < 0. In short, for a proton conductor with negatively large E′h•→H+, ionic defects are energetically favored. In addition, higher Ehydr (lower right direction in Figure 10) is more favored on the iso-E′h•→H+ lines to be a better proton conductor. We then examine how the hydration reaction is affected by the doping elements, with special attention to the hydrated configuration. The stable and unstable hydration configuration in Figure 10 can be classified into two types. For smaller dopants, i.e., Al-, Ga-, and Sc-doped BZOs, protons at 1NN (type A-1) are favored. This trend is not surprising given that •

A-2 configuration show negatively larger Δθ and ΔEhydr than those for Al-, Ga-, and Sc-doped BZOs. This implies that local stress induced by the presence of the large dopant and proton at 1NN can be relaxed by a change in the proton site (1NN to 2NN) accompanied by the rotation of the polyhedron with the dopant included. In other words, Y-, In-, and Lu-doped BZOs gain more hydration energy by structural relaxation. The same trend can be observed for isolated configurations, even though the magnitude of Δθ and ΔEhydr differs. In such cases, differences in the M−O−Zr angle (Δθ) and the hydration energy (ΔEhydr) between type B-2 (1NN, 2NN) and B-3 (2NN, 2NN) are defined as follows:

+

Δθisolated = θB‐3 − θB‐2 , ΔE hydr_isolated = E hydr_B‐3 − E hydr_B‐2

(13)

The distorted local structure and its dynamics, which should affect proton migration, are issues which undoubtedly will require much attention from researchers in this field in the future. G

DOI: 10.1021/acs.chemmater.6b03907 Chem. Mater. XXXX, XXX, XXX−XXX

Chemistry of Materials



To improve proton concentration for practical applications, the separation of the associated defects into isolated ones would be clearly beneficial from the viewpoint of hydration. As shown in Figure 7, negatively larger hydration energy is obtained if the associated configuration is isolated. The association energy of Y-doped BZO is 0.31 eV, which is the smallest of all trivalent dopants, as shown in Figure 2b. This energy is equal to 1800 K per one Y atom, suggesting that the density of protonic defects through hydration can be increased if the isolated defect state can be kept by quenching the sample from a high temperature. This would appear possible given that BZOs are generally prepared by sintering at temperatures above 1800 K.3,5,7 On the other hand, in the case of Sc- and Ga-doped BZOs with association energies of 0.73 and 1.52 eV for the isolation of dopants, both would require extremely high temperatures in excess of 4000 K. As such, isolating the complex defects must be considered impossible in the case of Sc- and Ga-doped BZOs. In this work, we focused on the formation of oxygen vacancies and the hydration of BZO at the static condition and investigated the key factors associated with the local structural configurations of defects. It should be acknowledged that other factors, including protonic mobility and the kinetic effect, were not taken into consideration. Further analysis involving a dynamic simulation by evaluating the migration barrier energy needs to be conducted.



CONCLUSION



ASSOCIATED CONTENT

Article

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. *E-mail: [email protected]. ORCID

Hitoshi Takamura: 0000-0002-4841-4582 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS



REFERENCES

This work was supported by the 69th committee on Materials Processing and Applications, University-Industry Cooperative Research Committees, Japan Society for the Promotion of Science (JSPS).

(1) Kreuer, K. D. Proton conductivity: Materials and applications. Chem. Mater. 1996, 8, 610−641. (2) Norby, T. Solid-state protonic conductors: principles, properties, progress and prospects. Solid State Ionics 1999, 125, 1−11. (3) Iwahara, H.; Yajima, T.; Hibino, T.; Ozaki, K.; Suzuki, H. Protonic conduction in calcium, strontium and barium zirconates. Solid State Ionics 1993, 61, 65−69. (4) Imashuku, S.; Uda, T.; Nose, Y.; Taniguchi, G.; Ito, Y.; Awakura, Y. Dependence of dopant cations on microstructure and proton conductivity of barium zirconate. J. Electrochem. Soc. 2009, 156, B1− B8. (5) Kreuer, K. D. Aspects of the formation and mobility of protonic charge carriers and the stability of perovskite-type oxides. Solid State Ionics 1999, 125, 285−302. (6) Ryu, K. H.; Haile, S. M. Chemical stability and proton conductivity of doped BaCeO3-BaZrO3 solid solutions. Solid State Ionics 1999, 125, 355−367. (7) Han, D.; Shinoda, K.; Sato, S.; Majima, M.; Uda, T. Correlation between electroconductive and structual properties of proton conductive acceptor-doped barium zirconate. J. Mater. Chem. A 2015, 3, 1243−1250. (8) Babilo, P.; Uda, T.; Haile, S. M. Processing of yttrium-doped barium zirconate for high proton conductivity. J. Mater. Res. 2007, 22, 1322−1330. (9) Kjølseth, C.; Wang, L. Y.; Haugsrud, R.; Norby, T. Determination of the enthalpy of hydration of oxygen vacancies in Y-doped BaZrO3 and BaCeO3 by TG-DSC. Solid State Ionics 2010, 181, 1740−1745. (10) Yamazaki, Y.; Babilo, P.; Haile, S. M. Defect Chemistry of Yttrium-Doped Barium Zirconate: A Thermodynamic Analysis of Water Uptake. Chem. Mater. 2008, 20, 6352−6357. (11) Oikawa, I.; Takamura, H. Correlation among Oxygen Vacancies, Protonic Defects, and the Acceptor Dopant in Sc-Doped BaZrO3 Studied by 45Sc Nuclear Magnetic Resonance. Chem. Mater. 2015, 27, 6660−6667. (12) Amezawa, K.; Takahashi, H.; Kuwabara, A.; Unemoto, A.; Kawada, T. Local structural arrangements around oxygen and hydrogen-related defects in proton conducting LaP3O9 investigated by first principles calculations. Int. J. Hydrogen Energy 2012, 37, 7995− 8003. (13) Hiraiwa, C.; Han, D.; Kuramitsu, A.; Kuwabara, A.; Takeuchi, H.; Majima, M.; Uda, T. Chemical Expansion and Change in Lattice Constant of Y-Doped BaZrO3 by Hydration/Dehydration Reaction and Final Heat-Treating Temperature. J. Am. Ceram. Soc. 2013, 96, 879−884. (14) Dawson, J. A.; Miller, J. A.; Tanaka, I. First-Principles Insight into the Hydration Ability and Proton Conduction of the Solid State Proton Conductor, Y and Sn Co-Doped BaZrO3. Chem. Mater. 2015, 27, 901−908.

We investigated the formation energy of oxygen vacancies and hydration energy with taking the local structural configurations of defects into account, as factors influencing the hydration property of acceptor-doped BZOs. This computational analysis has shown that the key factors influencing proton concentration are the stability of oxygen vacancies and the availability of a stable protonic site with structural distortions in hydration and that the degree of hydration is increased when the associated defects are isolated. It was shown that the creation of oxygen vacancies is favored for Al- and Ga-doped BZOs, leading to small hydration energies with stable proton sites at 1NN. On the other hand, the formation of oxygen vacancies is difficult in the case of Sc-doped BZO, but negatively large hydration energy is obtained. This trade-off relationship can be formulated as an equation with a slope of minus unity. Y-, In-, and Lu-doped BZOs gain more hydration energy by structural relaxation with protons located at 2NN in the hydrated state. This study suggests that increasing the isolation of the associated defects in BZOs wherever possible will result in increased hydration energies, making them more suitable for use in SOFCs.

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.chemmater.6b03907. Schematic representations of the configuration of Mdoped BZOs with complex defect (Figure S1), the convergence test of the hydration energy for Sc-doped BZO (Figure S2), and relative total energy of M-doped BZOs with complex defect with respect to M−M configurations (1−5NN) (Figure S3) (PDF) H

DOI: 10.1021/acs.chemmater.6b03907 Chem. Mater. XXXX, XXX, XXX−XXX

Article

Chemistry of Materials (15) Dawson, J. A.; Tanaka, I. Proton trapping in Y and Sn Co-doped BaZrO3. J. Mater. Chem. A 2015, 3, 10045−10051. (16) Stokes, S. J.; Islam, M. S. Defect chemistry and proton-dopant association in BaZrO3 and BaPrO3. J. Mater. Chem. 2010, 20, 6258− 6264. (17) Bjørheim, T. S.; Kuwabara, A.; Ahmed, I.; Haugsrud, R.; Stølen, S.; Norby, T. A combined conductivity and DFT study of protons in PbZrO3 and alkaline earth zirconate perovskites. Solid State Ionics 2010, 181, 130−137. (18) Sundell, P. G.; Björketun, M. E.; Wahnström, G. Thermodynamics of doping and vacancy formation in BaZrO3 perovskite oxide from density functional calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 104112. (19) Björketun, M. E.; Sundell, P. G.; Wahnström, G. Structure and thermodynamic stability of hydrogen interstitials in BaZrO3 erovskite oxide from density functional calculations. Faraday Discuss. 2007, 134, 247−265. (20) Bjørheim, T. S.; Kotomin, E. A.; Maier, J. Hydration entropy of BaZrO3 from first principles phonon calculations. J. Mater. Chem. A 2015, 3, 7639−7648. (21) Lindman, A.; Erhart, P.; Wahnström, G. Implications of the band gap problem on oxidation and hydration in acceptor-doped barium zirconate. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 245114. (22) Clark, S. J.; Segall, M. D.; Pickard, C. J.; Hasnip, P. J.; Probert, M. J.; Refson, K.; Payne, M. C. First principles methods using CASTEP. Z. Kristallogr. - Cryst. Mater. 2005, 220, 567−570. (23) Moriwake, H.; Fisher, C. A. J.; Kuwabara, A. First-Principles Calculations of Rare-Earth Dopants in BaTiO3. Jpn. J. Appl. Phys. 2009, 48, 09KC03. (24) Chase, M. W. NIST-JANAF Thermochemical Tables, Fourth ed.; American Institute of Physics: New York, 1998. (25) Oikawa, I.; Ando, M.; Kiyono, H.; Tansho, M.; Shimizu, T.; Maekawa, H. On the symmetry of defects in perovskite-type protonic conductors: A Sc-45 NMR study. Solid State Ionics 2012, 213, 14−17. (26) Blanc, F.; Sperrin, L.; Lee, D.; Dervişoğlu, R.; Yamazaki, Y.; Haile, S. M.; Paëpe, G. D.; Grey, C. P. Dynamic Nuclear Polarization NMR of Low-γ Nuclei: Structural Insights into Hydrated YttriumDoped BaZrO3. J. Phys. Chem. Lett. 2014, 5, 2431−2436. (27) Yamazaki, Y.; Blanc, F.; Okuyama, Y.; Buannic, L.; Lucio-Vega, J. C.; Grey, C. P.; Haile, S. M. Proton trapping in yttrium-doped barium zirconate. Nat. Mater. 2013, 12, 647−651. (28) Wagner, C. Die Löslichkeit von Wasserdampf in ZrO2-Y2O3Mischkristallen. Ber. Bunsenges. Phys. Chem. 1968, 72, 778−781. (29) Tandé, C.; Pérez-Coll, D.; Mather, G. C. Surface proton conductivity of dense nanocrystalline YSZ. J. Mater. Chem. 2012, 22, 11208−11213. (30) Jiang, J.; Hertz, J. L. Intermediate temperature surface proton conduction on dense YSZ thin films. J. Mater. Chem. A 2014, 2, 19550−19555. (31) Levin, I.; Amos, T. G.; Bell, S. M.; Farber, L.; Vanderah, T. A.; Roth, R. S.; Toby, B. H. Phase equilibria, crystal structures, and dielectric anomaly in the BaZrO3−CaZrO3 system. J. Solid State Chem. 2003, 175, 170−181. (32) Refson, K.; Clark, S. J.; Tulip, P. R. Variational density functional perturbation theory for dielectrics and lattice dynamics. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 155114. (33) Baroni, S.; Gironcoli, S. D.; Corso, A. D.; Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory. Rev. Mod. Phys. 2001, 73, 515−562. (34) Sternik, M.; Parlinski, K. Free-energy calculations for the cubic ZrO2 crystal as an example of a system with a soft mode. J. Chem. Phys. 2005, 123, 204708. (35) Sternik, M.; Parlinski, K. Lattice vibrations in cubic, tetragonal, and monoclinic phases of ZrO2. J. Chem. Phys. 2005, 122, 064707. (36) Bilić, A.; Gale, J. D. Ground state structure of BaZrO3: A comparative first-principles study. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 174107.

(37) Bévillon, E.; Dezanneau, G.; Geneste, G. Oxygen incorporation in acceptor-doped perovskites. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 174101. (38) Momma, K.; Izumi, F. VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data. J. Appl. Crystallogr. 2011, 44, 1272−1276. (39) Gomez, M. A.; Fry, D. L.; Sweet, M. E. Effects on the Proton Conduction Limiting Barriers and Trajectories in BaZr0.875Y0.125O3 Due to the Presence of Other Protons. Han'guk Seramik Hakhoechi 2016, 53, 521−528. (40) Polfus, J. M.; Bjørheim, T. S.; Norby, T.; Bredesen, R. Surface defect chemistry of Y-substituted and hydrated BaZrO3 with subsurface space-charge regions. J. Mater. Chem. A 2016, 4, 7437− 7444. (41) Björketun, M. E.; Sundell, P. G.; Wahnström, G.; Engberg, D. A kinetic Monte Carlo study of proton diffusion in disordered perovskite structured lattices based on first-principles calculations. Solid State Ionics 2005, 176, 3035−3040. (42) Nakamura, T.; Mizunuma, S.; Kimura, Y.; Yamauchi, K.; Mikami, Y.; Kuroha, T.; Amezawa, K. Design of proton ceramics fuel cells considering mixed conduction in solid electrolyte. In 18th International Conference on Solid State Protonic Conductors, Oslo, Norway, 2016.

I

DOI: 10.1021/acs.chemmater.6b03907 Chem. Mater. XXXX, XXX, XXX−XXX