On the Atomistic Interactions That Direct Ion Conductivity and Defect

Published: November 24, 2011 r 2011 American Chemical Society. 704 dx.doi.org/10.1021/jp2074682 |J. Phys. Chem. C 2012, 116, 704-713. ARTICLE...
3 downloads 0 Views 2MB Size
ARTICLE pubs.acs.org/JPCC

On the Atomistic Interactions That Direct Ion Conductivity and Defect Segregation in the Bulk and Surface of Samarium-Doped Ceria: A Genetic Algorithm Study Arif Ismail, Javier B. Giorgi,* and Tom K. Woo* Centre for Catalysis Research and Innovation, Department of Chemistry, University of Ottawa, 10 Marie Curie, Ottawa K1N 6N5, Canada

bS Supporting Information ABSTRACT: We study the (111) surface of 10.3%, 14.3%, and 18.5% samarium-doped ceria (SDC) using a genetic algorithm (GA) to search for the most energetically stable configurations. In all cases, both Sm ions and oxygen vacancies segregate to the surface, which is similar to experimental findings for 5.3% SDC.1 Importantly, at the optimal doping level of SDC (∼11%), where conductivity is maximal, defect segregation is limited such that vacancies remain 6 Å apart and pairs of vacancies do not form. At higher concentrations, pairs of vacancies are present, which likely contributes to the observed decrease in ionic conductivity. We also investigate the low-energy bulk structure of SDC from 10.3% to 18.5%, at the DFT+U level of theory, which has not been previously reported. The DFT+U energetics allow us to gain further insight on the fundamental interactions that influence ionic conductivity and defect segregation, as well as to confirm the insight reported from classical simulations.2 The low-energy configurations found by our GA search enable future studies of SDC at experimentally relevant concentrations and identify the important interactions at the bulk and surface of the fuel cell electrolyte.

1. INTRODUCTION Ceria (CeO2)-based materials are widely used as the electrolyte in solid oxide fuel cells (SOFCs) due to their high oxygen ion mobility.37 Because the rate of ion mobility affects the power SOFCs can deliver, understanding oxygen ion transport in ceriabased materials is of great interest to build better SOFCs. Cerium’s native oxide is commonly doped with a trivalent ion to improve its ionic conductivity.5 This is because for every two trivalent ions that replace tetravalent Ce ions, an oxygen vacancy is formed to maintain charge neutrality; thus, the resulting lattice has intrinsic vacancies, which are thought to be responsible for the increased ionic conductivity.5 Recently, it was shown that Sm3+ is the best dopant in this regard.4,6 In doped ceria, the dopant ions and oxygen vacancies form associates8 with corresponding association energies which influence vacancy mobility. At temperatures below 500 C, which are the target operating temperatures for SOFCs,911 the relative association energies of different configurations play a key role in diffusion pathways.6,12 It is thus desirable to know what structure factors minimize the dopantvacancy association energy, which will in turn increase vacancy mobility. Toward this goal, a number of computational studies have used simulation to calculate association energies of various configurations of lanthanide-doped ceria. The 3.2% dopant concentration is typically simulated by placing one vacancy and two dopant atoms in a 2  2  2 (96-atom) ceria lattice. Because of the symmetry of the cell, there are only 33 unique configurations, and so it is feasible to evaluate them all. Systematic investigations, r 2011 American Chemical Society

at classical, DFT, and DFT+U levels of theory, have identified the lowest-energy configurations that predominate in 3.2% SDC. Over the years, researchers have progressed from studying 3.2% lanthanide-doped ceria with inexpensive classical methods to more accurate and compute intensive calculations, such as the DFT+U method.6,1315 Dopant concentrations above 3.2% are more relevant to study as the optimal dopant concentration for SDC used in fuel cells is 11.1%.5 However, identifying dopantvacancy association energies at concentrations higher than 3.2% is much more challenging due to the enormous number of configurations possible. For instance, to build 10.3% SDC, one must place three vacancies and six Sm atoms in a 2  2  2 (96-atom) simulation cell, which correspond to over 193 trillion possible configurations. Because the configuration space of 10.3% SDC is so large, it is not feasible to evaluate with a systematic search. In a recent work, we developed a genetic algorithm to effectively recover the lowestenergy structures of doped metal oxides, with even larger search spaces, such as in a 3  3  3 simulation cell.16 This genetic algorithm was used to find the low-energy configurations of 10.2%, 14.9, 17.4%, and 20.0% SDC, at the classical level of theory. Above the optimal doping value of 11.1%, the conductivity of SDC drops as a function of dopant concentration.5 Our GA at the classical level found that vacancies aggregate in SDC at concentrations above 14.9%, which correlates fairly well with the Received: August 4, 2011 Published: November 24, 2011 704

dx.doi.org/10.1021/jp2074682 | J. Phys. Chem. C 2012, 116, 704–713

The Journal of Physical Chemistry C

ARTICLE

study.1 Also, the details of the surface defect interactions were not identified. The vacancyvacancy interactions at the surface could be different from those in the bulk, as a result of segregation, which would affect the surface conductivity. It is also unknown how increasing doping level affects segregation and the structure of the SDC surface. Of particular interest are concentrations near 11% where it is found that SDC has the highest ionic conductivity for use in SOFCs. In this work, we will apply a GA to study the (111) surface structure of 10.3%, 14.3%, and 18.5% SDC. We examine the (111) surface because HRTEM studies showed that particles mainly expose the thermodynamically stable (111) surfaces.1 Also, because of the stability of this surface, the atomic rearrangements are minimal, which is more amenable to investigate with the GA. As with the bulk material, studying the structure of high concentration SDC surfaces is computationally challenging due to the enormous number of configurations that can be made by arranging the defects in a surface slab. As our GA is well equipped to handle large potential energy surfaces,16 a structure search of this magnitude is made feasible. Indeed, to the best of our knowledge, this is the first GA study of the surface of any metal oxide material. The remainder of this Article is organized as follows. Details of the simulations are given in section 2. A thorough DFT+U validation of the GA-retrieved structures of bulk 10.3%, 14.3%, and 18.5% SDC is presented in section 3. In section 4, the surface structures of SDC are discussed and compared to the bulk material at the same concentrations. Throughout, DFT+U calculations are presented to confirm the findings of the GA, which uses classical potentials. Wherever possible, theoretical results are compared to experiment, and any discrepancies are rationalized in detail. Finally, the conclusions are presented in section 5.

drop in conductivity observed above 11.1%. The low-energy configurations demonstrated a preference for the next nearest neighbor (NNN) Smvacancy interaction, although both NN and NNN interactions were present in the favored structures. Several of these “hybrid” NN/NNN configurations recovered by the GA were found to be near-degenerate with classical potentials, which would presumably facilitate vacancy migration, as it involves rapidly changing from one low-energy configuration to the next. Interestingly, this defect distribution is also seen in gadolinium-doped ceria, which is another popular fuel cell electrolyte with conductivity similar to that of SDC.2 In the first part of this Article, we will attempt to corroborate the classical results with first-principles DFT+U calculations, to validate these conclusions. This is important because some studies have shown that DFT+U calculations can yield results different from those of lower level calculations for doped ceria. For example, we reported with DFT+U that Smvacancy interactions do not exclusively favor the NNN site in 6.6% SDC, in contrast to conventional DFT simulations.15 Also, Dholobhai et al. found that Gd3+ prefers the NN site in 3.2% Gd-doped ceria at the DFT+U level14 in contrast to previous findings at the classical level.12 Because DFT+U calculations are very compute demanding, we will employ a 96 atom, 2  2  2 simulation cell, instead of the 284-atom, 3  3  3 simulation cell, used in the classical GA study.2 Thus, the GA using the classical potential will be rerun with the smaller simulation cell on 10.3%, 14.3%, and 18.5% SDC models. Using these results, we will then examine the structures at the DFT+U level of theory. The goal will be to validate the results from the classical potential and possibly further refine our understanding of the influence of defect position on oxygen ion conductivity in doped ceria. In the second part of this Article, we will examine the surface structure of SDC. To date, all computational investigations of SDC have focused on the bulk material. While this provides insight on the conductivity in the bulk, all SOFCs will have a surface interface, which presumably affects the overall conductivity, particularly if the surface structure is very different from the bulk. Understanding the surface structure will provide further insight into the nature of the conductivity in these materials. Moreover, important chemistry takes place at the surface of these materials when used in SOFCs. Key catalytic reactions are believed to take place at the triple-phase boundary, which consists of the anode, the electrolyte surface, and the gaseous fuel. Thus, detailed knowledge of the surface structure of SDC is critical for meaningful computational investigations of the triplephase boundary and fuel oxidation mechanisms. Unfortunately, little is known about the surface structure of SDC. For example, no computational studies have examined SDC surfaces. A limited number of theoretical studies have examined the surface of other fuel cell electrolytes, such as gadoliniumdoped ceria (GDC) and yttrium-stabilized zirconia (YSZ).1719 The simulations based on classical potentials benchmarked to DFT calculations found that the defects segregated to the surface. The degree of segregation was found to be in general agreement with that found in experimental studies.2023 To date, the only experimental study of the SDC surface is a recent Raman investigation of 5.3% SDC, which found some defect segregation at the surface.1 Still, much of the SDC surface structure remains unknown. For instance, unlike the experimental studies of segregation in GDC and YSZ, the degree to which defect segregation occurs was not ascertained by the Raman

2. COMPUTATIONAL DETAILS DFT+U calculations were performed with the Vienna ab initio simulation package (VASP)24,25 and the projector augmented wave (PAW) method.26 The exchange-correlation effects were described with the Perdew Burke Erzenhof (PBE) functional27 within the generalized gradient approximation (GGA). In the DFT+U method, a Hubbard (U-term) correction is added to the conventional KohnSham Hamiltonian to more correctly account for the localization of strongly correlated electrons, and to mitigate the self-interaction error present in conventional DFT.28 A U-value of 5 eV was placed on all cerium ions, which is consistent with previous studies of doped ceria.13,14,2931 For the case of Sm, we used a pseudopotential in which the f electrons are placed in the core, and hence are already localized. Previous studies of lanthanide-doped ceria materials were performed in the same way, where no Hubbard correction is placed on the dopant ion.1315 The cerium 5s2, 5p6, 6s2, 5d1, and 4f1 states, as well as the samarium 5s2, 5p6, 6s2, and 5d1 states, and the oxygen 2s2 and 2p4 states, were treated as valence electrons, whereas the others were placed in the core. Structures were relaxed until the forces on each ion were below 0.02 eV/Å and the total energy was converged within 0.01 meV. A 520 eV plane-wave cutoff energy was used throughout. Studies of similar compounds using VASP have employed a plane-wave cutoff of 400 eV that were found to be converged to approximately 0.01 meV.6,13,14 Gamma point Brillouin zone sampling was used for the energy evaluations, although we confirmed that the use of a finer, 3  3  3 MonkhorstPack grid had negligible influence on the results. 705

dx.doi.org/10.1021/jp2074682 |J. Phys. Chem. C 2012, 116, 704–713

The Journal of Physical Chemistry C

ARTICLE

optimization at the classical level should not affect the recovery of the lowest-energy configurations. Besides, constant volume optimizations with GULP have been used to simulate lanthanidedoped ceria in the past.12 Furthermore, we verify the surface GA results with DFT+U calculations, wherein the cell volume is allowed to relax. In the (111) surface, the atoms are arranged in repeating oxygen metaloxygen (OMO) layers. We varied the slab thickness in our structure search of 10.3% SDC, employing slabs of 4, 6, and 8 OMO layers. For higher concentrations, only the 6 OMO slab was investigated, as it was deemed sufficiently large to observe defect segregation.

Table 1. Buckingham Parameters for Interatomic Potentials35,36 species

A (eV)

F (Å)

C (eV Å6)

O2O2

22 764.3

0.149

27.89

0.3511 0.3034

20.40 0.0

Ce4+O2 Sm3+O2

1986.8 4040.9

Other standard features include a variable cell and Gaussian smearing of 0.20 eV. Spin polarization was allowed, but it was found that all valence electrons were paired following relaxation. The classical simulations were done using the General Utility Lattice Program (GULP).32 These calculations involve shortrange pairwise Buckingham potentials (eq 1) and an Ewald summation33 of long-range interactions. The shell model34 is used to account for the polarizability of the O2 ions. The oxygen shell had a charge of 2.08 e and was tied to the atomic core by a 27.29 eV/Å2 spring constant. The Buckingham parameters, shown in Table 1, were taken from Balducci et al.35 and Senyshyn et al.36 The functions were evaluated between 0.5 and 10 Å for CeO and OO pairs and between 0.5 and 6 Å for SmO pairs. The elevation of the lower cutoff from 0.0 to 0.5 Å was employed to help eliminate unphysical structures from our genetic algorithm; it should in no way affect the energies of relevant configurations. Sij ¼ Aeðrij =FÞ 

C rij6

3. STRUCTURES OF BULK SDC 3.1. Bulk 10.3% SDC. In this section, we intend to compare the DFT+U results for 10.3% SDC to those obtained with a classical potential. 10.3% SDC was modeled with three vacancies and six Sm atoms in a 2  2  2 ceria lattice, which can be arranged in almost 200 trillion possible ways. Using classical potentials, the GA uncovered the lowest-energy structures of this material. (For details on the GA search, see the Supporting Information.) The top 40 configurations were found to exhibit the same structural qualities. Thus, in this material, there exist certain defect interactions that are dominant, which may play a role in the optimal conductivity at this doping level. These interactions are as follows. First, the Sm ions are preferentially located in the next nearest neighbor (NNN) sites around oxygen vacancies. This finding is in contrast to 6.6% SDC where no clear dopantvacancy preference exists.15 Second, the vacancies are all separated by at least 6 Å, as was also seen at lower concentrations.15 These findings are in agreement with classical GAs run on the larger 3  3  3 simulation cell.2 The fact that all lowenergy configurations recovered from the GA have the same distribution of defects suggests that this specific distribution is favored in 10.3% SDC. Yet, the above results are based on classical potentials. To corroborate these results, we have examined the configurations at the DFT+U level. A full DFT+U GA study at this concentration is currently not feasible due to the large search space. A reasonable way to validate the results obtained with the classical potential is to demonstrate that the same defect distribution is favored at the DFT+U level of theory. To assess the correlation between DFT+U and classical calculations, four sets of structurally different configurations were evaluated. The four sets consisted of: (A) the top 28 unique configurations recovered from the classical GA trials, all of which had the same structural characteristics described above; (B) 25 configurations, which had most of their Sm ions located NNN to the vacancies, a trait in common with the true low-energy structures; however, these 25 structures possessed a variety of vacancyvacancy distances, in contrast with the optimal separation found among the low-energy structures of set A; (C) five other configurations with different dopant and vacancy distributions; and (D) high-energy structures selected by a special version of the GA, which reversed the fitness to favor high-energy structures, at the classical level. Two such GA trials with 500 structures were run for 100 generations, and the three optimal structures from these trials were optimized. If the correlation holds, the distribution favored classically should also be favored at the DFT+U level. The results of our structural evaluation are summarized in Table 2. In each set, structures are ordered by increasing DFT+U

ð1Þ

The same force field was used in previous works, wherein we validated its use by comparing it to first-principles DFT and DFT+U calculations.15,16 In this work, we further its validation with DFT+U at higher SDC concentrations and find it is well suited to model this system. It is therefore appropriate to use these potentials in our GA and in future classical studies of SDC. To perform the GA calculations, we have used a code developed in our research group. The details of our genetic algorithm that searches bulk structures are described elsewhere.16 To summarize, it starts by randomly placing dopant ions and oxygen vacancies in a bulk 2  2  2 pure ceria simulation cell, to create the initial population of structures at a given concentration. The fitness of each structure is evaluated on the basis of its energy from a classical geometry optimization. Subsequently, new structures are created through mating and mutation operations starting from the configurations of the previous population; structures with higher fitness are given a greater chance of passing on their structural information to the new population, hence the “survival of the fittest” principle. The population is evolved until the best structures stop changing for a number of generations. A number of these trials are completed to ensure the lowest-energy structures are recovered at each concentration. The GA that searched the surface configurations used the same methods to find low-energy structures. However, instead of a bulk 2  2  2 ceria lattice, the starting cell in the surface GA was a (111) slab of undoped ceria, with 10 Å of vacuum in the z-direction to separate the surfaces. The slab was two unit cells long in the x and y directions and of variable thickness. Another minor difference in the surface GA is that the geometry optimizations were done at constant volume, to prevent the vacuum layer from collapsing, whereas in the bulk GA, the geometries were optimized at a constant pressure of 0 GPa and the cell was allowed to relax. This minor difference in 706

dx.doi.org/10.1021/jp2074682 |J. Phys. Chem. C 2012, 116, 704–713

The Journal of Physical Chemistry C

ARTICLE

Table 2. Relative Energies and Structural Information of Selected 10.3% SDC Structures Smvacancy distances no. NN

no. NNN

no. NNNN

no. NNNNN

ΔEDFT+Ub/eV

ΔEclassicalc/eV

distance range /Å

(∼2.4 Å)

(∼4.5 Å)

(∼5.9 Å)

(∼7.0 Å)

per cell

per cell

A

6.06.6

13

1113

03

03

00.10

00.15

B

6.06.6

4

9

5

0

0.19

0.70

3

9

6

0

0.23

0.81

2

7

8

1

0.46

1.85

2

6

8

2

0.56

2.28

3

5

9

1

0.56

2.41

3 1

7 6

7 10

1 1

0.56 0.62

2.64 2.32

1

5

9

3

0.74

2.96

1

7

7

3

0.85

3.29

0

5

8

5

1.12

4.09

6

0

11

1

0.76

2.46

vacvac set

C

D

a

6.08.1 7.67.6

6

0

12

0

1.30

3.07

2.75.4

10

8

0

0

1.63

4.10

2.75.4 2.75.4

2 0

11 12

5 6

0 0

2.18 2.28

2.13 5.26

2.73.8

0

0

10

8

3.30

7.31

2.73.8

0

0

9

9

3.39

6.13

2.73.8

0

1

7

10

3.41

6.17

a

Range of vacancyvacancy distances, considering the three distances between vacancies in the cell. b Relative total electronic energies from DFT+U calculations. c Relative total energies from classical potentials. Set A: All low-energy configurations recovered from the classical GAs have this distribution. A total of 28 structures of this type were evaluated. Set B: All structures in this set have the same optimal vacancy separation, but different Sm positions. Set C: Other selected distributions. Set D: Highest-energy configurations found at the classical level.

the vacancyvacancy separation distance, and the most favored structures all have vacancies separated by at least 6 Å. The structures in sets A and B of Table 2 show that the NNN position of the dopant is favored. In particular, in the second set, the two structures that were closest to the optimal distribution, with 9 NNN Smvacancy interactions, were lowest in energy out of all structures in this set, at 0.19 and 0.23 eV. All other structures in set B had 7 or fewer NNN interactions, and are higher energy as a result. The structures in set D of Table 2 serve as further justification for the correlation between the classical potentials and DFT+U. These configurations are significantly higher in energy than all other structures, at greater than 3 eV over structures with the favored distribution. This further suggests that the same distributions favored classically are favored by DFT+U, and that those disfavored classically are equally disfavored by DFT+U. It seems as though no other defect distribution is energetically competitive with the favored distribution. Indeed, the thermal energy available at 673 and 773 K, which are targeted SOFC temperatures, is 0.058 and 0.066 eV, respectively. Thus, one can expect only configurations with the favored distribution to be populated significantly at these temperatures. To accurately model the SDC electrolyte, future researchers should base their model on the optimal defect distribution found by our GA. Because the DFT+U energetics parallel those from classical potentials, the following conclusions from our previous classical GA study2 are confirmed. The GA found several structures with the preferred distribution, which are accessible thermally (within kT for operational temperatures) for 10.3% SDC. This contrasts

Figure 1. DFT+U energy of various 10.3% SDC configurations as a function of shortest vacancyvacancy distance. These structures have similar Smvacancy distances.

energy. Importantly, the low-energy configurations recovered from the classical GA (set A) are also the lowest-energy configurations at the DFT+U level. Overall, the results indicate that DFT+U favors precisely the same defect distribution, albeit with a smaller energy spread, as the classical force field. The energies of structures in set B, with similar dopantvacancy distances and different vacancyvacancy distances, help clarify how the vacancy separation affects the DFT+U energy of a configuration. The energies for these 25 configurations as a function of vacancy separation are plotted in Figure 1. As was observed with the classical force field,2 there is a strong dependence on the energy to 707

dx.doi.org/10.1021/jp2074682 |J. Phys. Chem. C 2012, 116, 704–713

The Journal of Physical Chemistry C

ARTICLE

Table 3. Relative Energies and Structural Information of Selected 14.3% SDC Structures Smvacancy distances no. NN

no. NNN

no. NNNN

no. NNNNN

ΔEDFT+Ub/eV

ΔEclassicalc/eV

distance range /Å

(∼2.4 Å)

(∼4.5 Å)

(∼5.9 Å)

(∼7.0 Å)

per cell

per cell

A

6.06.6

36

1620

46

35

00.13

00.29

B

6.06.6

8

11

10

3

0.25

1.65

4

14

10

4

0.50

2.25

1

15

11

5

0.84

3.14

3.89.4

7

11

9

5

0.66

3.25

6.67.6

6

9

12

5

1.52

3.52

3.88.1 5.49.4

4 6

13 10

10 10

5 6

1.58 1.62

4.80 4.55

6.67.6

6

8

14

4

1.65

3.79

6.67.6

4

12

12

4

1.67

3.61

4.78.1

4

11

11

6

1.68

4.05

2.78.1

5

14

10

3

2.13

5.81

2.75.4

16

16

0

0

3.06

7.72

2.75.4

0

0

16

16

5.87

18.13

vacvac set

C

a

a Range of vacancyvacancy distances, considering the four distances between vacancies in the cell. b Relative total electronic energies from DFT+U calculations. c Relative total energies from classical potentials. Set A: All low-energy configurations recovered from the classical GAs have this distribution. A total of 13 structures of this type were evaluated. Set B: All structures in this set have the same optimal vacancy separation, but different Sm positions. Set C: Other selected distributions.

with the results for 6.6% and 3.2% SDC, where much fewer lowenergy structures were found.15 Thus, the observed increase in conductivity as a function of dopant concentration up to 11.1%5 may be due to the fact that more configurations are accessible and that there exist more low-energy pathways for diffusion of oxygen ions in 10.3% SDC than in 6.6% or 3.2% SDC. Our theoretical findings for 10.3% SDC are in agreement with EXAFS measurements of 11.1% SDC. The experiment finds isolated vacancies characterized by a lack of NN Smvacancy interactions.37 This suggests a dominance of another Sm vacancy interaction, such as the NNN. Because the NN dopant vacancy interaction is stronger than the NNN interaction, dopants in NN sites tend to trap vacancies and limit conductivity.37 Thus, the NNN preference of Sm, or rather the lack of NN site preference, allows the vacancies to be more mobile in 10.3% SDC. Furthermore, the experiment suggests that no pairs of vacancies exist in 11.1% SDC, as replicated by the GA result. 3.2. Bulk 14.3% SDC. 14.3% SDC was modeled in a 2  2  2 ceria lattice with eight Sm atoms and four oxygen vacancies. The 40 most fit structures from the classical GA search all had the same defect distribution. This distribution is similar to that favored in 10.3% SDC, which agrees with previous findings in the larger simulation cell.2 The vacancies are separated by at least 6 Å, and the dopant ions prefer the NNN site surrounding the vacancies. To confirm that the same defect distribution is favored with higher level calculations, DFT+U optimizations were run on the top 13 configurations from the GA, as well as other selected configurations. To assess the influence of dopant ion positions on the energy, configurations with the optimal vacancy distribution but different Sm positions were evaluated. Also, various vacancy distributions were evaluated in DFT+U to assess the influence of vacancy separation on the energy. The results are shown in Table 3, where structures are grouped into sets (in the same way as Table 2). No reverse-fitness GA was run on 14.3% SDC,

because the correlation between the two methods has already been proven for 10.3% SDC, which we find to be very similar in structure. It is clear that the DFT+U-favored defect distribution matches the classically favored one, as all top 13 configurations lie within 0.13 eV/cell of the lowest-energy structure. Overall, the relative energies from classical potentials match fairly closely with the DFT+U rankings, again with a wider spread in the classical energies. The optimal vacancy separation remains at 6 Å in 14.3% SDC. This finding is in apparent contrast with experimental results that show ordering of vacancies at concentrations above 11%,5,38 represented by vacancy pairs separated by 4.7 Å or less. However, in experimental samples, there exist nanosized domains with doping levels higher than other parts of the sample,3943 which have been observed by TEM.38 Because vacancies have an association enthalpy with the dopant ions, it is reasonable to expect pairs of oxygen vacancies to form in these domains. Ou et al. suggested that chains of vacancies form in these domains, generating the order observed experimentally.38 In contrast, in our simulations, the entire sample is homogeneously doped at 14.3%. Thus, it is impossible to observe the effects seen experimentally due to local ordering, because our simulation cell is not big enough to house nanosized domains, which are on the order of 10 nm.39 Overall, we find that 14.3% SDC is very similar to the 10.3% SDC material. The same structural qualities are favored at the DFT+U level of theory at both concentrations, which are above and below the 11.1% SDC material used experimentally. Thus, we expect 11.1% SDC to mimic these structural properties, which are a dominance of the NNN Smvacancy interaction and at least a 6 Å separation between vacancies. These structural findings are in agreement with EXAFS measurements.37 3.3. Bulk 18.5% SDC. 18.5% SDC was simulated in a 2  2  2 cell of ceria with 10 Sm atoms and five vacancies. The best 708

dx.doi.org/10.1021/jp2074682 |J. Phys. Chem. C 2012, 116, 704–713

The Journal of Physical Chemistry C

ARTICLE

Table 4. Relative Energies and Structural Information of Selected 18.5% SDC Structures Smvacancy distances vacvac set

a

distance range /Å

no. NN

no. NNN

no. NNNN

no. NNNNN

ΔEDFT+Ub/eV

ΔEclassicalc/eV

(∼2.4 Å)

(∼4.5 Å)

(∼5.9 Å)

(∼7.0 Å)

per cell

per cell

A

6.09.4

48

2227

1014

710

00.07

0.060.13

B

4.76.6

48

2330

713

69

0.170.58

00.15

C

3.89.4

68

2225

1012

69

0.050.16

0.310.41

D

6.09.4

8d

20

17

5

0.86

2.00

10e

11

24

5

1.05

3.82

11f

12

23

4

1.10

2.59

9c 0

20 25

16 16

5 9

1.37 1.69

3.81 4.58

E

F

4.76.6

4.76.6

9

16

17

8

0.17

1.87

10e

15

17

8

0.20

1.94

8

16

17

9

0.35

2.77

10f

14

20

6

0.44

2.30

6d

15

23

6

0.65

4.00

0

24

21

5

1.12

4.34

5c 10e

15 18

22 16

8 6

1.49 0.66

7.49 1.94

8

19

18

5

0.73

1.93

6

20

18

6

0.89

2.55

6

22

18

4

0.91

2.43

Range of vacancyvacancy distances, considering the five distances between vacancies in the cell. b Relative total electronic energies from DFT+U calculations. Sets A and B: All low-energy configurations recovered from the classical GAs have one of these distributions. In total, 34 structures of these types were evaluated. Set C: The third best distribution at the classical level found by a special GA. Three structures of this type were evaluated. Set D: All structures in this set have the same vacancy distribution as set A, but different Sm positions. Set E: All structures in this set have the same vacancy distribution as set B, but different Sm positions. Set F: All structures in this set have two vacancy pairs at 4.7 Å (not found in the GA). c Sm ions clustered in the cell (small Sm-Sm distances). d Sm ions scattered about the cell (large SmSm distances). e Two Sm NN to each vacancy. f Four Sm NN to some vacancies, while others have no Sm NN. a

structures from the GA had one of two vacancy distributions. In all cases, most Sm are NNN to the vacancies, and most vacancies are separated by 6 or 6.6 Å. However, in some structures, there are two vacancies separated by 4.7 Å, whereas in others, all vacancies are separated by at least 6 Å. This result is in contrast to GA results for 10.3% and 14.3% SDC, wherein all favored structures shared the same defect distribution. Our GA study of 17.4% and 20.0% SDC in a 3  3  3 simulation cell found similar 4.7 Å vacancy separations in the most fit structure.2 However, that study did not consider more than the top four structures from the GA, which, in this case, does not provide a full description of the material. In fact, we find upward of the top 80 structures in 18.5% SDC exhibiting one of two vacancy distributions. To accurately determine the relative energies of structures with these distributions, DFT+U optimizations are necessary. Moreover, it is crucial to investigate the energies of structures with alternative defect distributions, to prove that the correlation between classical potentials and DFT+U holds at this concentration. Toward this goal, a special version of the GA was designed to find the third-best defect distribution favored at the classical level of theory. Even with large populations, the previous GA trials found only two different distributions of defects that were favorable, so we needed to alter the fitness to remove these two distributions from consideration. The GA was changed so that if a structure had one of the two favored vacancy distributions, it was assigned a large positive energy penalty to remove it from the

gene pool. From this search, a new vacancy distribution was found, which had most vacancies separated by 6 Å and one pair at 3.8 Å apart, with none 4.7 Å apart. Table 4 contains a summary of the DFT+U energetics of several 18.5% SDC structures. The structures are grouped into sets, in the same way as Table 2. The superscripts c-f denote other structural qualities that are not listed in the table, such as SmSm distance, which serve to distinguish certain configurations that have similar Smvacancy and vacancyvacancy distributions. For instance, superscripts c and d denote relatively short and long Sm-Sm separations, respectively. Sets A and B represent the top 34 configurations found by the classical GAs, which exhibit one of the two favored distributions. Set C represents three structures with the third-most optimal distribution found by a special classical GA (explained above). Sets D and E represent 12 configurations that had either of the two preferred vacancy distributions, but different Sm positions than were found in the GA. Set F contains four configurations with an alternative vacancy distribution, which has two vacancy pairs separated by 4.7 Å. The idea behind this was to test the relative energies of various structures with the two favored vacancy distributions, and compare their relative energies with structures of other, similar, vacancy distributions. Overall, the last two columns of Table 4 suggest there is reasonable agreement between the classical and DFT+U energies for 18.5% SDC structures. The lowest-energy structures had all vacancies separated by at least 6 Å and so are more likely to exist 709

dx.doi.org/10.1021/jp2074682 |J. Phys. Chem. C 2012, 116, 704–713

The Journal of Physical Chemistry C

ARTICLE

performed on 10.3% SDC slabs with 2  2 unit cell surfaces of 4, 6, and 8 OMO layers deep. In each GA, the population consisted of 500 structures and was evolved for a large number of generations. To determine convergence, we waited until the defect distribution among the best three configurations remained the same for at least 100 generations. This was deemed sufficiently converged based on previous benchmarks.2 For all (111) slab sizes, we found a Smvacancy distribution similar to that of the low-energy bulk structures. Specifically, the 6 Å separation between all vacancies is maintained at the surface, and most Sm lie NNN to the vacancies. If the optimal bulk structures are cut to form (111) slabs, their defects are spread evenly throughout the slabs. However, in the (111) slabs recovered from the 6 OMO and 8 OMO GAs, we observed defect segregation to the surface (see Figure 2). To quantify the level of segregation, an “effective dopant concentration” is calculated for the middle of the slab as well as the surface. The effective concentration in the slab middle was determined from the number of dopant atoms in the two central OMO layers of the slab. The effective concentration at the surface was determined from the number of dopant atoms in the top-most OMO layer on each side of the vacuum region. These correspond to the top- and bottom-most OMO layers in Figure 2. In the 4 OMO slab, no segregation was observed due to the relatively small depth of the slab. In the 6 and 8 OMO slabs, the effective dopant concentration is, on average, 0.8% in the middle of the slabs, and 18.5% at the surface. As the number of OMO layers was increased from 6 to 8, the degree of segregation across the slab did not change. Thus, larger slabs were not evaluated. Despite the high effective dopant concentration at the surface, oxygen vacancies are still separated from one another by at least 6 Å in the entire slab. Thus, it seems like the defects tend to segregate to the surface, but not so much as to create paired vacancies, which are known to be energetically unfavorable. This suggests that a segregation limit exists at the 10.3% SDC surface, which allows vacancies to be mobile by remaining 6 Å apart. The degree of segregation has not been measured by experiment. However, our results are in agreement with spectroscopic measurements of the 5.3% SDC surface, which showed some level of defect segregation.1 Experiments on GDC also showed segregation,23 which further illustrates the similarity of these two materials, alluded to in earlier works.2,15 We find that SDC differs from YSZ in that the NN dopantvacancy interaction is disfavored in SDC and favored in YSZ,19 although dopant segregation is found in both materials. Interestingly, the vacancies do not appear right at the surface of the material. While there are several close to the surface, they always exist in the bottom part of the top-most OMO layer. Similar subsurface vacancy aggregation was found by a Monte Carlo study of the (111) surface of YSZ, which is another common fuel cell electrolyte.19 It should be noted that the slabs found by the GA differ from the optimal bulk positions slightly. In the slabs cut from the optimal bulk positions, the dopants and vacancies spread throughout the cell, and there is no segregation. In contrast, the surface GA finds slabs with defect segregation, resulting in slightly more Smvacancy NNN interactions and slightly fewer NN interactions than in the slabs made from the optimal bulk positions. Thus, it could be that when SDC surfaces are prepared, the ions rearrange to their optimal slab positions, which, according to our GA, include defect segregation and slightly fewer NN Smvacancy interactions.

Figure 2. Lowest-energy 6 OMO (left) and 8 OMO (right) 10.3% SDC (111) slabs recovered from classical GAs.

at low temperatures as compared to structures that have a vacancy pair 4.7 Å apart. Yet there do exist some structures with the 4.7 Å vacancy pair that are as low as 0.17 eV above the lowestenergy structure. Because the energy ranges of structures with different defect distributions overlap, it is likely that 18.5% SDC contains structures with more than one defect distribution, at low temperatures. The third-best distribution favored at the classical level (found with the special GA) turns out to be lower in energy in DFT+U than expected. This result suggests that the correlation between the two methods is not as strong at the 18.5% concentration as it is at lower concentrations. It is interesting that structures with paired vacancies have energies as low as 0.05 eV/cell. At 873 K, the Boltzmann ratio of a configuration 0.05 eV/cell relative to the lowest-energy structure is 1:2, which suggests that it will be populated significantly at this temperature. Thus, the DFT+U energetics predict that paired vacancies are present in the material. Furthermore, these results also indicate that structures with vacancies 3.8 Å apart are favored over those with vacancy pairs at 4.7 Å (for example, see DFT+U energies of set C versus set B). This supports the work of Ou et al., who found that the Æ110æ pairs (3.8 Å) are more stable than the Æ111æ pairs (4.7 Å) in doped ceria.38,44 The above finding is distinct to this work and was not found in our prior study of 17.4% and 20.0% SDC.2 In our previous work, we were unable to determine the most favorable distance for vacancy pairs in highly concentrated SDC solely on the basis of classical energetics. This demonstrates the importance of using first-principles DFT+U calculations, to ensure an accurate representation of structures’ relative energies. The set F distribution, with two pairs of vacancies separated by 4.7 Å, was higher in energy (g0.66 eV/cell) than the classically favored defect distributions. This particular vacancy distribution was not found in any classical GA, and it has a correspondingly high energy in DFT+U. This suggests that the top three classically favored distributions may be the most stable in DFT+U. Thus, the correlation is still enough as not to compromise the validity of the classical GAs in searching for lowenergy structures.

4. SURFACE STRUCTURE OF SDC 4.1. 10.3% SDC Surface. The structure of the SDC (111) surface was explored by the classical GA. Searches were 710

dx.doi.org/10.1021/jp2074682 |J. Phys. Chem. C 2012, 116, 704–713

The Journal of Physical Chemistry C

ARTICLE

Table 5. Relative Energies and Structural Information of 10.3% 6 OMO Slabs structure retrieval

ΔEDFT+U/eV per cell

difference from optimal slab

ΔEclassical/eV per cell

optimal slab from GA

none

0

selected manually

more segregation

0.48

0 1.39

selected manually

more NN (fewer NNN)

0.66

3.98

from optimal bulk positions

no segregation

0.67

4.97

selected manually

more NN (fewer NNN) and more segregation

1.77

6.16

Figure 3. Lowest-energy 14.3% (left) and 18.5% (right) SDC (111) slabs recovered from classical GAs. The color scheme is the same as in Figure 2. Figure 4. Effective dopant concentration (% Sm2O3) as a function of depth in the 6 OMO surface slabs of 10.3%, 14.3%, and 18.5% SDC. The effective concentration was calculated by the average number of Sm atoms per layer in the three lowest-energy structures recovered from each GA. x = 1 represents the two surface layers, x = 2 the two subsurface layers, and x = 3 the two central layers of the 6 OMO slab. The lines of best fit show how the effective concentration changes from the surface (x = 1) to the middle of the slab (x = 3).

Even when no segregation is found, for the smaller 4 OMO slab, the validity of the GA can be tested by comparison with a bulk-cut surface. To show that the GA result from classical simulation is also favored in DFT+U, we optimized two 4 OMO surface slabs with the DFT+U method, one that was found by the GA, and the other that was cut from the optimal bulk positions. The slab from the surface GA had no Sm NN to any vacancy, but 28 NNN interactions, as compared to the other that had 4 NN and 22 NNN Smvacancy interactions. DFT+U calculations found the bulk surface slab to be higher in energy by 0.3 eV per cell, which indicates that the slab from the surface GA is heavily favored. Five 6 OMO slabs, each consisting of 279 atoms, were also optimized with the DFT+U method to confirm the GA result. Some were created manually and chosen to exhibit particular characteristics that differed from the GA-favored structures. A slab cut from the optimal bulk positions of 10.3% SDC was also tested. Table 5 summarizes the DFT+U and classical energies of these slabs. Overall, the results illustrate that the correlation between classical and DFT+U potentials holds even for surfaces. First, the slab recovered from the GA (which contained one undoped OMO layer at the center; see Figure 2) had the lowest energy in DFT+U. Second, the slab cut from the optimal bulk positions proved to be 0.67 eV higher in energy than that from the GA, most likely because it did not exhibit defect segregation. A slab with the same amount of segregation, but more Sm NN to the vacancies, was 0.66 eV higher in energy, in relative agreement with the classical potentials. The slab with more segregation than found by the GA (that is, with an effective surface concentration higher than 18.5%) was also high energy with both classical and DFT+U methods. Finally, a slab was made with two characteristics that differed from the optimal slabs: more NN Smvacancy interactions and more defect segregation, which had an additive effect on the structure energy, as both characteristics are unfavorable. In light of the remarkable correlation of classical and

DFT+U potentials at this concentration, we are confident that the classical GAs found the true limit to the segregation and the low-energy defect distributions that are present at the (111) surface. 4.2. 14.3 and 18.5% SDC Surfaces. It is interesting that, although the vacancies segregate to the surface along with the dopants in 10.3% SDC, the vacancies remain 6 Å apart. This suggests that the segregation in SDC, at the optimal doping level, continues until the point where paired vacancies would form. At higher doping levels, the number of defects must become so high that it forces vacancies to pair at the surface, as seen in experiment38 and in bulk SDC (section 3). To analyze if this hypothesis holds, we performed GA searches on 14.3% and 18.5% SDC (111) surfaces. We investigated only the 6 OMO slab, because it was large enough to observe segregation previously and still small enough to provide reasonable calculation times. Classical GAs were run on 14.3% and 18.5% (111) SDC surfaces, to find the low-energy defect distributions. The best structures from the surface GAs are shown in Figure 3. The effective dopant concentration as a function of the OMO layer in the 10.3%, 14.3%, and 18.5% SDC slabs is plotted in Figure 4. The average of the number of Sm atoms per layer in the three lowest-energy structures from the GA at each concentration was used to construct the points in the plot. Figure 4 allows for comparison of the level of segregation at various concentrations of SDC. For instance, the segregation in 14.3% SDC is just as evident as in the 10.3% SDC slab, as the slopes of the two lines are nearly identical. However, at 18.5%, less segregation is visible, due to the increased number of defects in the cell. 711

dx.doi.org/10.1021/jp2074682 |J. Phys. Chem. C 2012, 116, 704–713

The Journal of Physical Chemistry C

ARTICLE

At 14.3%, the effective dopant concentration is, on average (over the best three structures), 6.7% in the middle of the slab and 21.5% at the surface, which indicates significant defect segregation. Interestingly, at the 14.3% SDC surface, there exists a vacancy pair 3.8 Å apart, which was not seen in the bulk at the same concentration. Thus, in contrast to the 10.3% SDC surface, it appears that segregation produces vacancy pairs in 14.3% SDC, which could account for the decrease in ionic conductivity observed experimentally above 11.1%.5 At 18.5%, segregation is not as obvious as it is at lower concentrations, although still present. The middle of the slab has an effective dopant concentration of 16.4%, while the surface is 20.8%. The size of the slab could be limiting the observation of defect segregation here, and so 20.8% may be an underestimate of the true value. However, the purpose of these 14.3% and 18.5% SDC calculations was simply to examine if vacancy pairs formed, not to assess the number. Thus, larger slabs were not evaluated. In 18.5% SDC, vacancy pairs were found in the Æ110æ sites (i.e., the 3.8 Å separation), which is consistent with the findings of Ou et al., who claimed that vacancy pairs prefer Æ110æ sites over Æ111æ sites in doped ceria.38,44 We find no vacancy pairs exist in the Æ111æ sites in our slab calculations. From these results, we cannot determine the limit of segregation in 14.3% and 18.5% SDC materials, as larger slabs were not evaluated with the GA. However, the conclusions that (1) vacancy pairs exist at the surface, and (2) Æ110æ pairs are favored over Æ111æ pairs, are still consistent with the experiments of Ou et al.38,44 For these reasons, we feel the small slab size does not invalidate the results. The distribution of defects at the surface is similar to the bulk, in that the NNN Smvacancy interaction is dominant. This finding is somewhat counterintuitive, as one would expect that segregation of dopants and vacancies to the surface would cause an increase in NN interactions. However, the lack of NN dopant vacancy interactions in doped ceria has been suggested to play a key role in the increased mobility of the oxygen ions.37 That is, in materials where the NN site is heavily favored by the dopant ion, vacancies become trapped in these sites, thereby limiting oxygen conductivity.12,37,45 As SDC contains mostly NNN dopant vacancy interactions, the vacancies can move quickly, at the surface as well as in the bulk. By this rationale, our surface structure evaluation predicts that conductivity would not be diminished at the surface, relative to the bulk, in 10.3% SDC. However, in 14.3% and 18.5% SDC, more vacancy pairing is observed at the surface than in the bulk (as a result of segregation), which would limit conductivity, despite the favorable placement of the dopants.

5. CONCLUSION We have provided for the first time a DFT+U evaluation of the structure of concentrated SDC. This work confirms the results of classical simulations published earlier, which uncovered a relationship between the conductivity of SDC and oxygen vacancy separation.2 Our rigorous structure evaluation provides insight on the atomistic interactions that influence ionic conductivity and defect segregation. For instance, it was determined that Sm ions lie preferentially in the NNN site around vacancies, in 10.3%, 14.3%, and 18.5% SDC, in contrast to at lower concentrations.15 These findings indicate a lack of NN Smvacancy interactions, which is consistent with EXAFS measurements.37 It has been suggested that the lack of NN dopantvacancy interactions in certain doped ceria materials allows vacancies to remain mobile.37

In addition, we found a favorable 6 Å separation between oxygen vacancies in SDC, which may play a crucial role in the conductivity mechanism in this material. All vacancies in the bulk 10.3% SDC material were found to be separated by at least 6 Å, although at higher concentrations, vacancy pairs (separated by 3.8 Å) are expected to be present, most likely in nanosized domains of concentrated dopant. Because of the homogeneous distribution of dopants in our periodic model, vacancy pairing started to be observed in the bulk at 18.5%, and not at 14.3%. The above conclusions, which initially arose from classical simulations, are now confirmed by DFT+U calculations. Thus, for most concentrations of SDC (up to 18.5%), the classical potentials of Balducci35 and Senyshyn36 accurately reproduce the DFT+U-calculated energies for SDC configurations. A new result that emerged from solely DFT+U calculations was that Æ110æ vacancy pairs are favored over Æ111æ pairs, and this is consistent with experimental findings for doped ceria.38,44 In this work, we have also studied the surface structure of SDC at the classical and DFT+U level. The surface GA investigation found that defect segregation is present in 10.3%, 14.3%, and 18.5% SDC, which is consistent with experimental work on 5.3% SDC.1 Thus, we conclude that segregation is present at all doping levels. More importantly still is that our surface study sheds new light on the atomistic interactions at the SDC surface, which were not obtained by experiment. For instance, it was determined that segregation in SDC at the optimal doping level occurs until the point where paired vacancies would form. Above the optimal doping level, defect segregation in SDC causes vacancies to pair up, predominantly in the Æ110æ sites. The crucial decrease in vacancy separation from 10.3% to 14.3% SDC most likely represents the origin of the optimal doping level of 11.1% in SDC. In other words, above a doping level of 11.1%, the vacancies become so numerous that pairs are formed, which decreases conduction. At the surface, Sm dopants prefer to remain in NNN sites around the oxygen vacancies, which is the same as in the bulk. This result is important because it shows that, despite segregation, there are very few Sm NN to vacancies at the surface, which are believed to limit conductivity.37 Finally, it should be noted that the models used in this work do not account for certain phenomena that can occur in reality, such as grain boundaries, and inhomogeneity in the bulk material. The inhomogeneity was addressed in the surface study, as the slab was made larger to observe segregation. For the grain boundaries, it has been shown that conductivity in grain boundaries is markedly different from that in other parts of the material.46 Yet, this study was primarily concerned with outlining the reasons for the overall conductivity of the material, and understanding segregation in SDC; by neglecting the grain boundary effect, we are still able to reproduce available experimental results and explain the observed changes in conductivity as a function of dopant concentration.

’ ASSOCIATED CONTENT

bS

Supporting Information. GA search details. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected] (J.B.G.); [email protected] (T.K.W.). 712

dx.doi.org/10.1021/jp2074682 |J. Phys. Chem. C 2012, 116, 704–713

The Journal of Physical Chemistry C

ARTICLE

’ REFERENCES

(38) Ou, D. R.; Mori, T.; Ye, F.; Zou, J.; Auchterlonie, G.; Drennan, J. Phys. Rev. B 2008, 77, 024108. (39) Mori, T.; Drennan, J.; Lee, J.-H.; Li, J.-G.; Ikegami, T. Solid State Ionics 2002, 154155, 461. (40) Mori, T.; Drennan, J.; Wang, Y.; Auchterlonie, G.; Li, J.-G.; Yago, A. Sci. Technol. Adv. Mater. 2003, 4, 213. (41) Ou, D. R.; Mori, T.; Ye, F.; Takahashi, M.; Zou, J.; Drennan, J. Acta Mater. 2006, 54, 3737. (42) Ye, F.; Mori, T.; Ou, D. R.; Takahashi, M.; Zou, J.; Drennan, J. J. Electrochem. Soc. 2007, 154, B180. (43) Ye, F.; Mori, T.; Ou, D. R.; Zou, J.; Auchterlonie, G.; Drennan, J. Solid State Ionics 2008, 179, 827. (44) Li, Z. P.; Mori, T.; Ye, F.; Ou, D. R.; Zou, J.; Drennan, J. J. Chem. Phys. 2011, 134, 224708. (45) Nakayama, M.; Martin, M. Phys. Chem. Chem. Phys. 2009, 11, 3241. (46) Hong, S. J.; Mehta, K.; Virkar, A. V. J. Electrochem. Soc. 1998, 145, 638.

(1) Guo, M.; Lu, J. Q.; Wu, Y. N.; Wang, Y. J.; Luo, M. F. Langmuir 2011, 27, 3872. (2) Hooper, J.; Ismail, A.; Giorgi, J. B.; Woo, T. K. Phys. Chem. Chem. Phys. 2010, 12, 12969. (3) Inaba, H.; Tagawa, H. Solid State Ionics 1996, 83, 1. (4) Muthukkumaran, K.; Bokalawela, R.; Mathews, T.; Selladurai, S. J. Mater. Sci. 2007, 42, 7461. (5) Yahiro, H.; Eguchi, Y.; Eguichi, K.; Arai, H. J. Appl. Electrochem. 1988, 18, 527. (6) Andersson, D. A.; Simak, S. I.; Skorodumova, N. V.; Abrikosov, I. A.; Johansson, B. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 3518. (7) Fu, Y. P.; Wen, S. B.; Lu, C. H. J. Am. Ceram. Soc. 2008, 91, 127. (8) Mogensen, M.; Sammes, N. M.; Tompsett, G. A. Solid State Ionics 2000, 129, 63. (9) Dong, D. H.; Li, D.; Zhang, X. Y.; Chai, Z. L.; Wang, K.; Li, C. Z.; Zhao, D. Y.; Wang, H. T. J. Mater. Chem. 2010, 20, 1122. (10) Lin, B.; Chen, J.; Ling, Y.; Zhang, X.; Jiang, Y.; Zhao, L.; Liu, X.; Meng, G. J. Power Sources 2010, 195, 1624. (11) Zhang, X. G.; Robertson, M.; Yick, S.; Deces-Petit, C.; Styles, E.; Qu, W.; Xie, Y. S.; Hui, R.; Roller, J.; Kesler, O.; Maric, R.; Ghosh, D. J. Power Sources 2006, 160, 1211. (12) Wei, X.; Pan, W.; Cheng, L.; Li, B. Solid State Ionics 2009, 180, 13. (13) Dholabhai, P. P.; Adams, J. B.; Crozier, P.; Sharma, R. J. Chem. Phys. 2010, 132, 094104. (14) Dholabhai, P. P.; Adams, J. B.; Crozier, P.; Sharma, R. Phys. Chem. Chem. Phys. 2010, 12, 7904. (15) Ismail, A.; Hooper, J.; Giorgi, J. B.; Woo, T. K. Phys. Chem. Chem. Phys. 2011, 13, 6116. (16) Hooper, J.; Ismail, A.; Giorgi, J. B.; Woo, T. K. Phys. Rev. B 2010, 81, 224104. (17) Dalach, P.; Ellis, D. E.; van de Walle, A. Phys. Rev. B 2010, 82, 144117. (18) Lee, H. B.; Prinz, F. B.; Cai, W. Acta Mater. 2010, 58, 2197. (19) Mayernick, A. D.; Batzill, M.; van Duin, A. C. T.; Janik, M. J. Surf. Sci. 2010, 604, 1438. (20) Bernasik, A.; Kowalski, K.; Sadowski, A. J. Phys. Chem. Solids 2002, 63, 233. (21) Hughes, A. E.; Badwal, S. P. S. Solid State Ionics 1990, 401, 312. (22) Hughes, A. E.; Sexton, B. A. J. Mater. Sci. 1989, 24, 1057. (23) Scanlon, P. J.; Bink, R. A. M.; van Berkel, F. P. F.; Christie, G. M.; van Ijzendoorn, L. J.; Brongersma, H. H.; Van Welzenis, R. G. Solid State Ionics 1998, 112, 123. (24) Kresse, G.; Hafner, J. Phys. Rev. B 1993, 47, 558. (25) Kresse, G.; Furthmuller, J. Comput. Mater. Sci. 1996, 6, 15. (26) Blochl, P. E. Phys. Rev. B 1994, 50, 17953. (27) Perdew, J. P.; Burke, K.; Erzenhof, M. Phys. Rev. Lett. 1996, 77, 3865. (28) Pacchioni, G. J. Chem. Phys. 2008, 128, 182505. (29) Nolan, M.; Verdugo, V. S.; Metiu, H. Surf. Sci. 2008, 602, 2734. (30) Yang, Z.; Luo, G.; Lu, Z.; Hermansson, K. J. Chem. Phys. 2007, 127, 074704. (31) Yang, Z.; Luo, G.; Lu, Z.; Woo, T. K.; Hermansson, K. J. Phys.: Condens. Matter 2008, 20, 035210. (32) Gale, J. D. J. Chem. Soc., Faraday Trans. 1997, 93, 629. (33) Ewald, P. P. Ann. Phys. 1921, 64, 253. (34) Gale, J. D. Philos. Mag. B 1996, 73, 3. (35) Balducci, G.; Kaspar, J.; Fornasiero, P.; Graziani, M.; Islam, M. S. J. Phys. Chem. B 1998, 102, 557. (36) Senyshyn, A.; Oganov, A. R.; Vasylechko, L.; Ehrenberg, H.; Bismayer, U.; Berkowski, M.; Matkovskii, A. J. Phys.: Condens. Matter 2004, 16, 253. (37) Yoshida, H.; Deguchi, H.; Miura, K.; Horiuchi, M.; Inagaki, T. Solid State Ionics 2001, 140, 191. 713

dx.doi.org/10.1021/jp2074682 |J. Phys. Chem. C 2012, 116, 704–713