Nonrandom Point Defect Configurations and Driving Force Transitions

Nov 6, 2014 - The energetically favorable spatial configuration of M3+ ions and oxide-ion vacancies near a symmetrical grain boundary (GB) in cubic zi...
0 downloads 6 Views 4MB Size
Download PDF
Article pubs.acs.org/Langmuir

Nonrandom Point Defect Configurations and Driving Force Transitions for Grain Boundary Segregation in Trivalent Cation Doped ZrO2 T. Yokoi,*,† M. Yoshiya,†,‡ and H. Yasuda†,§ †

Department of Adaptive Machine Systems, Osaka University, 2-1 Yamadaoka, Suita, Osaka 565-0871, Japan Nanostructures Research Laboratory, Japan Fine Ceramics Center, 2-4-1 Mutsuno, Atsuta, Nagoya 456-8587, Japan § Department of Materials Science and Engineering, Kyoto University, Yoshidahonmachi Sakyo, Kyoto 606-8507, Japan ‡

S Supporting Information *

ABSTRACT: The energetically favorable spatial configuration of M3+ ions and oxide-ion vacancies near a symmetrical grain boundary (GB) in cubic zirconia is determined for various trivalent species M3+ (M = Al, Sc, Y, Gd, La), and the driving force for grain boundary segregation (GBS) quantitatively examined using atomistic Monte Carlo simulations in conjunction with static lattice calculations. For a high concentration of ∼10 mol %, it is found that point defects near a GB plane preferentially occupy specific sites to minimize total lattice energy, rather than being randomly distributed. Systematic analysis shows that energetically stable configurations of segregants vary depending on their ionic radii. Analysis of the driving force for GBS as a function of dopant concentration reveals that three important factors govern GBS. First, occupation of specific sites by point defects is necessary to minimize the total lattice energy; enrichment of point defects near the GB plane with random configuration does not decrease the total lattice energy significantly because of strong Coulombic interactions. Second, the factors governing GBS change with increasing dopant concentration. At dilute concentrations, relief of bond strain is the dominant factor, while at high concentrations Coulombic interactions, which depend strongly on the specific arrangement of defects, become another dominant factor. Third, the stabilization of matrix cations, Zr4+ ions, is the dominant factor to lower the driving force for GBS at all concentrations. In contrast, the stabilization of M3+ ions does not necessarily contribute to GBS of point defects at high concentrations. These findings suggest practical ways to control GBS to enhance materials’ properties or minimize detrimental effects.

1. INTRODUCTION Interfaces, such as surfaces, grain boundaries (GBs), and heterointerfaces, often trigger interfacial segregations of extrinsic and intrinsic point defects including solute elements, interstitials, antisite defects, and vacancies. Grain boundary segregation (GBS) is one of the important segregation phenomena, which frequently modifies the macroscopic properties of a material, even though the deviation from the bulk crystal structure and chemical composition, with the exception of space charge,1,2 may be confined to only a few nanometers on either side of the GB itself. GBS is also typically associated with detrimental effects, such as toughness in metals, with a typical example being embrittlement and lower fracture of sulfur GBS in iron and nickel.3−5 As many studies have shown, however, GBS and increased GB density can have beneficial effects on electrical,6,7 magnetic,8−10 and diffusional properties.11,12 To understand the effect of GBS and obtain insights into GBS mechanisms, continuum-matter-based theories, such as the Langmuir−McLean theory, have been proposed. However, the changes induced by GBS cannot © XXXX American Chemical Society

always be adequately described within continuum-matter-based theories alone: Recent studies have revealed that both detrimental and beneficial effects may result from the microscopic, or atomistic, changes near a GB, such as changes to local crystal structures,13−15 bond lengths,16 and electronic structures.17 When tailoring a material for a specific application, it is necessary to suppress detrimental effects and to acquire beneficial effects simultaneously, which requires an understanding of microscopic/atomistic structure and behavior of GBS, as well as to establish clearly the connection between underlying chemistry and physics and a material’s macroscopic properties. M-doped ZrO2 (MDZ), where M is an aliovalent cation, is commonly used as the solid-state electrolyte in solid oxide fuel cells,18−20 as well as a model material for studies of oxide-ion conductivity. GBs and GBS are known to have an important Received: August 21, 2014 Revised: November 6, 2014

A

dx.doi.org/10.1021/la503338x | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

Figure 1. (a) Schematic diagram of GB model. (b) Spatial configuration of Zr4+ and O2− ions along the [001] projection.

influence on the macroscopic O2− ionic conductivity21,22,24−26 in polycrystalline systems. In these materials, divalent or trivalent cations are typically substituted for Zr4+ ions and O2− vacancies created to maintain charge neutrality. In the case of MO- or M2O3-doped ZrO2, the defect reaction can be generally described using Kröger−Vink notation as × MO → M″Zr + V •• O + OO

columns. In addition, it is currently impossible to experimentally determine the precise position of O2− vacancies when O2− sublattices near GB planes are severely distorted and segregation of O2− vacancies brings about disorder in O2− columns. On the other hand, a number of studies, including those of Oyama et al.,30 Yoshiya and Oyama,31 and Lee et al.53 have examined specific atomic configurations of Y3+ and O2− vacancy in the vicinity of the Σ5 (310)/[001] symmetric tilt GB and other high-symmetry GBs31 in Y2O3-doped ZrO2.30 These studies also investigated the driving force for GBS based on results of atomistic simulations, and suggested the importance of O2− vacancy for Y segregation. Mao et al.35 and Marinopoulos36 employed ab initio quantum mechanical methods to calculate GBS energies at dilute dopant concentrations and suggested that the free volume at a GB is the dominant factor influencing GBS. Most of these computational studies have been concerned with dilute dopant concentrations, where the interactions between point defects are weak or nonevident. In realistic cases, however, interactions between many point defects cannot be ignored due to high dopant concentrations of approximately 8−10 mol %, even before defect enrichment triggered by GBS. To the best of our knowledge, few previous studies paid attention to the driving force for GBS at high defect concentrations, with most studies focusing only on dilute concentrations. Therefore, it is still unclear whether the knowledge based on dilute concentrations can be applied even to high defect concentrations. In order to understand GBS in a wide spectrum of materials, it is important to determine the energetically stable configurations of point defects near GBs at high concentrations and the driving force that triggers GBS phenomena, even if they are different from those for dilute concentrations. In order to understand GBS in a wide range of materials, it is important to determine the energetically stable configurations of point defects near GBs at high concentrations and the driving force behind GBS phenomena, even if they are different from those for dilute concentrations. In addition, most studies have looked at Y2O3-doped ZrO2 only. To obtain a more comprehensive understanding that can be applied to other systems, it is necessary to understand the relationship between point defect distributions and dopant species as well as GB structures.

(1)

and × M 2O3 → 2M′Zr + V •• O + 3OO

(2)

for divalent and trivalent dopants, respectively. The vacancies formed allow O2− ions to migrate through the crystal at elevated temperatures. In the case of Y2O3, the maximum conductivity is achieved for a dopant concentration of 8−10 mol %.22,23 Recent studies have suggested that point defect distributions at the atomic level have a strong effect on ionic conductivity. For example, Pornprasertsuk et al.54 revealed that the specific spatial configuration of Y3+ defects surrounding a given O2− migration path affects its energy barrier. In addition, polycrystals of MDZ, rather than its single crystals, are generally used for practical applications so that careful attention must be paid to the interactions not only between point defects but also between GBs, point defects, and any segregations that may occur. According to experimental reports, the ionic conductivity of GBs in polycrystals is several orders of magnitude lower than the bulk ionic conductivity even in nominally pure specimens.26−29 These results imply that the effects of GBs and GBS on ionic conductivity are intrinsic. However, influence of GBS on GB ionic conductivity remains unclear. In order to reveal their effects, first of all, it is necessary to understand atomic-level GB structure and point defects distribution near GBs. With the development of aberration-corrected transmission electron microscopy (TEM) and scanning TEM (STEM) with atomic resolution, spatial configurations of cations near GBs have been determined for an increasing number of materials.37−39 Dickey et al.,40 Shibata et al.,41−43 Lei et al.,44 and recently An et al.52 have successfully observed directly the detailed arrangements of cations in Y2O3-doped ZrO2 using TEM/STEM techniques. However, it is still a considerable challenge to experimentally distinguish between different elements with similar atomic numbers, for example, Y (Z = 39) and Zr (Z = 40) within the same or neighboring crystal B

dx.doi.org/10.1021/la503338x | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

2.2. GBS Energy and Distribution of Point Defects. In this study, we examined the mechanisms of GBS at both dilute and high dopant concentrations to identify factors that drive GBS. First, we calculated GBS energies at dilute concentrations, ΔEGBS dilute, as a function of distance between dopant ions and the Σ5 GB plane to determine the dependence of ΔEGBS dilute on dopant site occupancy. We assume that ′ and one V•• a defect pair consisting of one MZr O vacancy, written in • Kröger−Vink notation as (M′Zr − V•• O ) , is introduced into the bicrystal model with uniform background charges to maintain charge neutrality. The influence of the uniform background charges is canceled out in •• the calculation of ΔEGBS dilute. VO is assumed to reside on the next-nearestneighbor sites relative to M′Zr, in accordance with experimental observations45,46 and previous atomistic analyses,47,48 although again GBS its influence is canceled out in the calculation of ΔEGBS dilute. Here ΔEdilute is defined as follows:

In this study, we performed atomistic simulations for a number of different trivalent dopant species to determine the energetically favorable configurations of point defects after segregation, and to clarify the driving force and the mechanisms behind such a phenomenon. The trivalent species considered were Al, Sc, Y, Gd, and La. In order to determine the energetically favorable configurations of defects in these systems, atomistic Monte Carlo simulations involving site exchange were carried out in conjunction with static lattice calculations using empirical potential models. Then, based on optimized configurations of Y3+ and O2− vacancy, the driving force for GBS was analyzed, focusing on how specific spatial configurations affect the total lattice energy of a GB-containing system.

GBS ΔEdilute = E(r ) − E(∞)

2. COMPUTATIONAL PROCEDURE

• where E(r) is the total lattice energy when (M′Zr − V•• O ) is positioned at distance r from the GB plane, and E(∞) is the reference energy, • which is assumed to be equal to the lattice energy when (M′Zr − V•• O) is located at the midpoint between two GBs in the supercell. When •• • ΔEGBS dilute is negative, segregation of (M′Zr − VO ) to the GB is energetically favorable. Second, to examine the case of high dopant concentrations typically used in real materials, we prepared three sets of GBs with different in dopant distributions as well as local concentration near the GBs: These we refer to as (1) the high concentration segregation (HCSG) model; (2) the high concentration random (HCRD) model; and (3) the overall random (OR) model. In the HCSG model, M3+ ions and O2− vacancies occupy specific sites to minimize the total lattice energy of a system. Such specific spatial configurations of M3+ and O2− vacancy are obtained by MC simulation, as explained in the next section. We defined the GB region to be the region within which M3+ ions and O2− vacancies segregate in the HCSG model. In this case, the grain interior remains undoped.30,58 The thickness of the GB region is determined from the magnitude of enrichment of the segregants in the HCSG model. In the HCRD model, these point defects are randomly distributed within the GB region, whereas enrichment of both M3+ and O2− vacancy is assumed to occur to the same magnitude within the GB region. In the OR model, it is assumed that GBS, enrichment or depletion of M3+ and O2− vacancy, does not occur at all; these segregants are randomly distributed throughout the supercell. Random distributions of M3+ and O2− vacancy result in different lattice energies in both the OR and HCRD models. To obtain a representative value, we prepared 30 different configurations and averaged their lattice energies. At high dopant concentrations, the driving force for GBS, ΔEGBS high , is defined as the difference in the lattice energy of the HCSG model relative to that of the OR model. At high concentration, ΔEGBS high is given by

2.1. GB Model. A Σ5 (310)/[001] symmetric tilt GB model was chosen to study segregation in cubic zirconia because a bicrystal containing such a GB can be experimentally fabricated.40,52 It is known that this high symmetry GB has a slightly higher energy31 than many others that have been studied, which should promote GBS of point defects. Experimental measurements have confirmed that GBS of point defects occurs in Σ5 (310)/[001] bicrystals of Y2O3-doped ZrO2.40,52 In order to reduce computational time, we utilized the property of point symmetry of the computational cell as shown in Figure 1a. The elements within the region 0.0 ≤ x ≤ 0.5 and 0.0 ≤ y ≤ 0.5 are inversed in point P (x,y) = (0.5,0.5) in fractional coordinates. Figure 1b depicts the GB model which contains two equivalent GBs upon application of three-dimensional periodic boundary conditions. The undoped GB model consists of 236 Zr4+ ions and 512 O2− ions. The GB planes are initially separated by 48 Å so that the interaction between two equivalent GBs in the supercell is negligible (Figure1b). In the undoped GB model, we modified the initial configuration of Zr4+ ions and O2− near the GB plane from the one based on a simple CSL construction so that ions located too close to each other were removed to avoid unrealistically high ionic repulsions.30 In order to calculate the total energy of a system, we used the static lattice calculations implemented in the GULP code57 to relax the structure until zero-net forces act on each atom. An empirical potential model was used, as ab initio calculations based on density functional theory, for example, are too computationally demanding at this stage. At least 106−107 trial moves were performed using the Monte Carlo (MC) method to obtain a converged minimum-energy structure at high dopant concentrations. We note that previous studies have successfully reproduced the experimentally observed configurations of cations in Y2O3-doped ZrO2.30 A Buckingham-type potential was used with parameters that successfully reproduce the crystal structure of Y 2 O3 -stabilized ZrO 2 , as well as those of the various M 2 O 3 sesquioxides.47 In this model, the total lattice energy is represented as the summation of short-range interaction energy Eshort and longrange Coulombic interaction energy ECoulomb,

E

total

=E =

short

1 2

+E

⎧ ⎪

i

j>i

GBS ΔE high = EHCSG − E OR HCSG

⎫ ⎛ −r ⎞ C ⎪ ij ⎟ − ij ⎬ + 1 ⎟ r6 ⎪ ⎝ ρij ⎠ ⎭ 2

⎩

∑∑ i

j>i

(5)

OR

and E are the lattice energies of the HCSG and OR where E models, respectively. It should be noted that the number of M3+ and O2− vacancy species is the same in these three GB models at a given concentration so that we can compare the lattice energies of these GB models directly. If ΔEGBS high is negative, the specific spatial configuration of M3+ ions and O2− vacancies in the HCSG model is energetically more favorable than the OR model. The magnitude of ΔEGBS high indicates to what extent M3+ ions and O2− vacancies form a specific configuration in the HCSG model compared with the random distribution in the OR model. The driving force for GBS in the HCRD model was also evaluated in the same manner to understand the influence of M3+ and O2− vacancy enrichment in the vicinity of the GB. The influence of specific configurations can be also quantified by comparing ΔEGBS high of the HCSG model with that of the HCRD model. To obtain greater insights into GBS phenomena, ΔEGBS high is divided into two components according to eq 3, namely a short-range interactions, ΔEshort, and Coulombic interactions, ΔECoulomb. Furthermore, ΔEGBS high can also be separated into components from each of

Coulomb

∑ ∑ ⎨⎪Aij exp⎜⎜

(4)

qiqj rij (3)

where rij is the distance between ions i and j; Aij, ρij, and Cij are empirically derived potential parameters for a particular combination of ion species i and j; and qi is the formal charge of ion i. For shortrange interactions, the parameters reported by Minervini et al. were used.47 Although the potential parameters for Zr4+ and O2− cannot reproduce monoclinic or tetragonal phases of ZrO2, which are the thermodynamically stable phases when either dopant content or temperature is low, in this study, we were concerned only with dopant concentrations near the GB sufficient to stabilize the cubic fluorite structure. The cutoff radius for short-range interaction is 20 Å. C

dx.doi.org/10.1021/la503338x | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

the constituent elements, ΔEZr for Zr4+, ΔEO for O2−, and ΔEM for M3+; such an analysis allows the contribution of these elements to the GBS overall ΔEGBS high to be assessed. Thus, ΔEhigh can be expressed as GBS ΔE high = ΔE short + ΔE Coulomb

= ΔEZr + ΔEO + ΔEM

(6)

A similar separation is also applied to ΔEGBS dilute for dilute dopant concentrations. 2.3. Monte Carlo Simulations. In order to obtain minimal-energy configurations of M3+ ions and O2− vacancies, we used the Metropolis MC method, which is a well-established optimization algorithm based on random sampling.49 Only a short description of the procedure as applied to this study will be given here, as more comprehensive explanations are available elsewhere.30 First, M3+ ions are randomly substituted for Zr4+ ions and O2− vacancies are introduced into an initial supercell while maintaining charge neutrality. Then, M3+ ions and O2− vacancies are exchanged with Zr4+ ions and O2− ions, respectively, within MC cutoff radius of 6 Å to generate a new configuration. After calculating the lattice energies of the initial and new configurations by the static lattice method, a decision is made whether the new configuration is accepted or rejected according to the Metropolis algorithm. The acceptance probability, P, is

⎧1 if ΔE ≤ 0 p=⎨ ⎩ exp( −ΔE/kBT ) if ΔE > 0 ⎪ ⎪

(7)

where kB is Boltzmann’s constant, T is absolute temperature, and ΔE is the difference in lattice energies between the current and new configurations. This procedure is repeated for a given segregant content until the lattice energy converges to a minimum value, consistent with the magnitude of thermal fluctuation indicated by eq 7. Once the minimum energy and configuration of the segregants are obtained, another defect complex is introduced to each GB region and the procedure repeated until segregants begin to move into the grain interior, which is assessed by calculating the GBS energy for dilute dopant concentrations as described above. This procedure corresponds to obtaining an equilibrium concentration of point defects in the grand canonical ensemble.

3. RESULTS AND DISCUSSION 3.1. Dilute Concentrations. In this section, we discuss the energetics and mechanism of GBS at dilute defect concentrations for five species of dopants. The model used in this section is close to the ideal case of infinite dilution, but its simplicity enables us to identify and understand the factors influencing the segregation mechanism more easily. First, according to Figure 2a and f, which shows ΔEGBS dilute as a function • GBS of distance between (YZr ′ − V•• O ) and the GB plane, ΔEdilute • tends to decrease as (Y′Zr − V•• ) approaches the GB plane, O although it is not a simple function of distance from the GB • plane. The overall tendency of ΔEGBS ′ − V•• dilute for (YZr O) , especially in the vicinity of the GB plane, is in good agreement with the results of previous studies,30,31 even though the supercell in this study is much larger, with more degrees of freedom, than previously reported models. In order to identify • GBS the dominant factors for GBS of (YZr ′ − V•• O ) , ΔEdilute is divided into ΔEshort and ΔECoulomb components as shown in Figure 2 a. This shows that ΔEshort is the main contributor to the decrease in ΔEGBS dilute, driving GBS of the defect pair. In other words, the relief of the bond strain is responsible for the decrease in ΔEshort, since ΔEshort is strongly correlated with the magnitude of bond strain. On the other hand, comparison of the components of each element, ΔEZr, ΔEO, and ΔEY shown in Figure 2f, suggests that stabilization of the O2− sublattice does not contribute

• Figure 2. GBS energy as a function of distance of (MZr ′ − V•• O ) from the GB plane in a Σ5 (310)/[001] ZrO2 bicrystal.

significantly to GBS; ΔEO exhibits positive values when (Y′Zr • − V•• O ) is close to the GB plane. The relatively unstable nature of the O2− sublattice relative to the cation sublattice may be related to the high diffusivity of O2− in this well-known oxideion conductor. In contrast, ΔEZr and ΔEY decrease as (Y′Zr − • V•• O ) approaches the GB plane. Therefore, the stabilization of 4+ Zr ions and Y3+ ions, rather than O2− ions, is the dominant • factor for GBS of (Y′Zr − V•• O ) . This may be attributed to the 4+ strong bonds by Zr ions in the bulk crystal. It is likely that disorder on the Zr4+ sublattice carries with a much higher energy penalty compared to the O2− sublattice. If ΔEO were negative, probably the activation enthalpy for O2− migration would be too high for O2− diffusion. This simple model, where ′ and one V•• only one YZr O are introduced, should help us understand the more complex phenomena in systems with high dopant concentrations, as will be discussed in subsection 3.3. • The same tendency as (Y′Zr − V•• O ) can be seen for the other dopant species, as shown in Figure 2b−e, although the absolute value of ΔEGBS dilute depends strongly on the dopant species, most likely because of the different size mismatches between Zr4+ D

dx.doi.org/10.1021/la503338x | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

and M3+ ions. Figure 2f−h shows that ΔEZr also plays an important role in decreasing ΔEGBS dilute for trivalent cations other • than Y; overall, ΔEZr is negative and decreases as (MZr ′ − V•• O) approaches the GB plane. This feature of the stabilization of Zr4+ sublattice appears to be common to the various trivalent cations at dilute concentrations. Needless to say, the nature of the bonding depends on the main phase elements, so that different behavior can be expected for systems other than zirconia. Nevertheless, the stabilization of the host elements in response to the distribution of segregant ions and vacancies plays an important role in minimizing ΔEGBS dilute in dilute systems. It should be noted that the geometric constraint of (M′Zr − • V•• O ) , which was described in section 2, is imposed in order to quantify the magnitude of ΔEGBS dilute relative to grain interior. The A site in Figure 2, which exhibits the lowest ΔEGBS dilute, is not necessarily the segregation site for high defect concentrations. In fact, there are some sites that exhibit similar values of ΔEGBS dilute to that at the A site. Thus, it is not surprising that, when the geometric constraint is removed, other sites near a GB become energetically favorable segregation sites depending on dopant concentrations, as will be discussed in subsection 3.2. 3.2. GBS at High Dopant Concentrations. In this subsection, the spatial configuration of M3+ ions and O2− vacancies at high concentrations is discussed. Figure 3 shows

alternately along the [001] direction. Experimentally, it is difficult to identify precisely the positions of Y3+ in Y2O3-doped ZrO2. However, the Y3+ configuration calculated here has been indirectly verified by STEM observation of segregation at GBS.30,40 Thus, this configuration appears to be not only favorable, but also exhibited in real zirconia materials. Identifying the precise distribution of oxide-ion vacancies within the material is not straightforward, since the O2− sublattice near GBs is highly distorted. Instead, we examined the overall one-dimensional distribution of O2− using Gaussian broadening with a full width at half-maximum of 1.0 Å, as shown in Figure 4. The valleys in the O2− profile indicate

Figure 4. Distribution of ions near GB in the case of (a) the HCSG model and (b) the HCRD model. Blue, red, and green lines correspond to Zr4+, O2−, and Y3+ ions, respectively. Dashed line shows the position of the GB plane.

enrichment of O2− vacancy. The O2− vacancy distribution in the HCSG model (Figure 4a) is also nonuniform while that in the HCRD model (Figure 4b) is uniform. In addition, while the distributions of point defects in the HCRD model are almost identical, those in the HCSG model are different from each other. It appears that O2− vacancies also occupy different sorts of specific sites near the GB plane to accommodate the different sizes of trivalent cations. O2− vacancies might thus also form an ordered arrangement, as has been intensively investigated in single-crystal models.55,56 Unfortunately, it is very difficult to identify positions of O2− ions or columns of ions near the GB. Figure 5 shows the spatial configurations for other trivalent cations examined in this study, Al3+, Sc3+, Gd3+, and La3+. Three of the dopants, Y3+, Gd3+, and La3+, are larger than Zr4+, while Al3+ and Sc3+ are smaller. The arrangement of the dopants can be seen to vary with the type of trivalent cations. For example, the larger cations (Figure 5c and d) occupy the same A sites, while the smaller cations (Figure 5a and b) do not occupy these sites. In contrast, the smaller cations occupy the B sites while the larger cations do not. In addition to the A and B site, some columns are fully occupied by M3+ while some are not in the [001] direction. The overall point defects distributions are also not completely random as shown in Figure 6: The dopant distributions are not uniform like the well potential but are not uniform over the GB region, either. In addition, the O2− vacancy distributions seem to be spatially correlated with those of the individual dopants in Gd- and La-doped ZrO2, as shown in Figure 6c and d. In Al- and Sc-doped ZrO2, the O2− vacancy distributions look like they are random, but Figure 5a and b shows that the columns on or in the vicinity of the GB plane decrease the number of O2−. Therefore, O2− vacancies occupy the energetically favorable sites even in these cases.

Figure 3. Distribution of Y3+ ions near a Σ5 GB. Blue, red, and green balls represent Zr4+, O2−, and Y3+ ions, respectively. The dashed line indicates the GB plane.

the energetically favorable configuration of Y3+ ions and O2− vacancies, where 48 Zr4+ ions are substituted by Y3+ ions, and 24 O2− vacancies are introduced. The number of Y3+ ion corresponds to a concentration of ∼20 atom % in the GB region, and the local concentration near the GB plane is enriched to ∼30 atom %. The dopant content remains almost constant above this concentration, and addition of even more dopant simply results in the high dopant concentration in the grain interior. The deviation of dopant and vacancy concentrations from their bulk values can be easily understood in terms of chemical potential. It should be noted that the local concentration of dopants and vacancies exceeds the solubility limit of the cubic phase of zirconia. Local deviations of the thermodynamic states relative to the bulk are found for other systems including stabilization of intergranular glassy layers which are unstable in the bulk phase.32−34 It is found that there are the energetically favorable and unfavorable cation sites for segregation of Y3+. For instance, in Figure 3, no Y3+ occupies the A sites while the B sites are fully occupied by Y3+, and the C sites are partially occupied by Y3+. This trend suggests that the spatial configurations of Y3+ are nonrandom. In addition, an ordered arrangement of Y3+ also forms within ca. 6 Å either side of the GB plane. There are many columns where Zr4+ and Y3+ occupy cation sites E

dx.doi.org/10.1021/la503338x | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

Figure 5. GB structures projected onto a plane perpendicular to the GB plane for M2O3-doped ZrO2, where M = Al, Sc, Gd, and La.

minimum ΔEGBS dilute in Figure 2 due to the geometric constraint. In order to compare directly the results at dilute and high concentrations, we have carried out a set of extra calculations without the geometric constraint at dilute concentrations. As a result, it is found that for larger M3+ a dopant occupies the A site, while the B site is most favorable for smaller M3+. This result is in agreement with the results in Figure 5. The detailed discussion is described in the Supporting Information. 3.3. Driving Force for GBS at High Concentrations. In this subsection, we discuss the driving force for GBS at high dopant concentrations, based on calculations of ΔEGBS high , with the aim of determining the mechanism responsible for GBS. In particular, our attention is given to the dominant factors causing GBS that vary with increasing dopant concentration. First, we focus on the Y2O3-doped ZrO2 system and then discuss the other systems to obtain a more general understanding. Figure 7 shows ΔEGBS high as a function of the number of defect complexes, that is, two Y3+ ions and one O2− vacancy, at various concentrations, from 3.5 to 11.3 mol %, which was defined as the fraction of the number of Y2O3 to that of ZrO2 in the GB region. It is noted that the concentration in the OR model is uniform throughout the simulation cell, while those in the

Figure 6. Gaussian distribution as a function of position from the GB plane in the case of (a) Al, (b) Sc, (c) Gd, and (d) La doped ZrO2.

Except for Y2O3-doped ZrO2, the dopant configuration near the GB plane is yet to be verified by experiment and thus further study is needed. These results clearly show that the distribution of oxide-ion vacancies, in addition to dopant cations, differs depending on the type of dopant species. The profile for the larger cations is clearly correlated with that of the anion vacancies, while the smaller cations show the greatly different spatial profile from O2− vacancy. This difference in vacancy distributions between trivalent cations larger or smaller than Zr4+ undoubtedly affects O2− diffusivity in the GB plane region. It should be emphasized that dopant and vacancy profiles differ from one another, meaning that, at the atomic level, local charge neutrality will be violated to some extent. Although the distributions for the Σ5 (310)/[001] ZrO2 GB are different from that suggested by the space charge theories applied for zirconia systems by Guo et al.,50,51 a nonuniform ionic charge distribution is indeed created when electronic segregation is neglected in the classical model used in this study. As discussed in subsection 3.1, the segregation site at high concentrations is not necessarily the same as the site having

Figure 7. Driving force as a function of the number of defect complexes in the case of (a) the HCSG model and (b) the HCRD model. Red and purple lines correspond to the total value, ΔEshort, and ΔECoulomb, respectively. F

dx.doi.org/10.1021/la503338x | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

dopant. For cations smaller than Zr4+ (Figure 9a and b), ΔEshort remains the dominant factor at all concentrations, while for larger cations (Figure 9c and d), ΔECoulomb decreases with increasing dopant concentration, becoming lower than ΔEshort. The difference in the dominant component of ΔEGBS high between smaller and larger cations is a result of their different preferred spatial configurations near the GB, as shown in Figures 5 and 6. For the smaller dopants, defects segregate within a narrower region than do the larger cations, as seen in Figures 5a,b and 6a,b. On the other hand, the larger cations spread out over a wider region with a certain degree of order, as shown in Figures 5c,d and 6c,d. The energetically stable configurations of dopants and vacancies are likely determined by a balance or competition between the magnitudes of Eshort and ECoulomb at high dopant concentrations. On the other hand, Figure 9e−h shows that ΔEZr is entirely negative for all the cases examined irrespective of the dopant species, while the contribution of ΔEO and ΔEM to ΔEGBS high depends on the type of dopant. Therefore, we conclude that the stabilization of Zr4+ sublattice also makes a dominant contribution to ΔEGBS high regardless of the type of trivalent dopant.

HCSG and the HCRD models are not uniform, with Y3+ ions and O2− vacancies concentrated in the vicinity of the GB plane. Figure 7a indicates that ΔEGBS high for the HCSG model is significantly lower than that of the OR model. In contrast, Figure 7b shows that ΔEGBS high of the HCRD model is nearly zero for all concentrations. These results clearly indicate that GBS produces a specific spatial configuration of Y3+ ions and O2− vacancies in the GB region. Substitution on specific sites near the GB plane minimizes ΔEGBS high , while a random distribution of 3+ 2− defects does not produce low ΔEGBS high , even if Y ions and O vacancies are concentrated near the GB plane. short Examining components of ΔEGBS is the dominant high , ΔE factor at low concentrations resulting in minimum values of ΔEGBS high in both the HCSG and HCRD models, as shown in Figure 7a and b, respectively. This indicates that segregation of defects to the GB results in decrease of ΔEshort . In addition, the magnitudes of ΔEshort in the HCSG and HCRD models are similar to each other. Thus, increased segregation of defects leads to a decrease in Eshort whether they are randomly distributed or not. However, the difference in ΔEGBS high between the HCSG and HCRD models stems from the magnitude of ECoulomb. ΔECoulomb in the HCSG model is lower by ∼1 eV than that in the HCRD model, and ΔECoulomb is determined by the distribution of Y3+ ions and oxide-ion vacancies. The GBS contribution of ΔECoulomb to ΔEhigh is enhanced by the 3+ occupation of specific sites by Y ions and oxide-ion vacancies. Consequently, ΔEGBS high for the HCSG model is the lowest in magnitude of all three models, and the calculated configuration is the one closest to that observed experimentally. ΔEGBS high can also be divided into the components of each element, ΔEZr for Zr4+, ΔEO for O2−, and ΔEY for Y3+. Figure 8a shows that ΔEZr is strongly negative, while ΔEO and ΔEY are

4. CONCLUSIONS We have determined the energetically stable spatial configurations of trivalent dopants and O2− vacancies and GBS mechanisms at various dopant concentrations by means of atomistic simulations. First, we evaluated ΔEGBS dilute and identified the dominant factors for GBS at dilute dopant concentrations. Energetically stable configurations of the same defect species were also identified at high dopant concentrations. With the obtained stable configuration of M3+ ions and O2− vacancies for each dopant species, we have evaluated and examined ΔEGBS high at high concentration. Based on these results and findings, the following conclusions can be drawn: (1) At dilute dopant concentration, ΔEGBS dilute deceases as a • (MZr ′ − V•• O ) becomes closer to the GB plane. This result is in good agreement with previous studies. The dominant factor for GBS is short-range interactions between cations and anions, which is related to bond strain. (2) Stabilization of Zr4+ and M3+ on energetically favorable sites is the dominant factor for GBS at dilute concentration. The results of the five trivalent cations shows that the most stable sites of M3+ is the same regardless of dopant species. We supposed that this trend of a segregated site and mechanism behind it can be applied to other dopant species having different valency. (3) At high dopant concentrations, M3+ ions occupy the specific sites near the GB plane to minimize total lattice energy. According to the Monte Carlo simulation, O2− vacancy distributions are also nonuniform, rather than random distributions. Comparison of the results among the five trivalent cations revealed that the energetically stable configuration of M3+ ions and O2− vacancies strongly depends on dopant species. (4) Enrichment of GB regions with M3+ ion and O2− vacancy distributed randomly does not result in significant decrease in ΔEGBS high due to strong Coulombic interactions. Well-defined spatial distribution of M3+ ions and O2− vacancies that decreases the Coulombic interaction is essential to minimize GBS ΔEhigh . The reduction of Coulombic interactions and stabilization of host cation species on specific sites are the key factors leading to GBS of both dopant and vacancy species.

Figure 8. Driving force as a function of the number of defect complexes in the case of (a) the HCSG model and (b) the HCRD model.

positive at all concentrations in the HCSG model. These results indicate that the stabilization of matrix cations, that is, Zr4+, plays a dominant role in GBS, while the stabilization of O2− and Y3+ ions has a slightly negative impact on ΔEGBS high . These same factors influencing GBS of point defects hold for the other dopant species as well. Comparing results from the HCSG model with those from the HCRD model reveals that ΔEZr is lower for the HCSG model than the HCRD model, and thus the difference in ΔEGBS high between the HCSG and HCRD model stems from the difference in magnitude of ΔEZr, that is, the degree of stabilization of the host cations. The same analysis for Y3+ was carried out for the other dopants. Figure 9a−d shows that the dominant component controlling GBS gradually changes with the ionic radius of the G

dx.doi.org/10.1021/la503338x | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

Figure 9. Driving force as a function of the number of defect complexes for M2O3-doped ZrO2, where M = Al, Sc, Gd, and La. These figures show the results of the HCSG model in all cases.

■

Notes

ASSOCIATED CONTENT

The authors declare no competing financial interest.

S Supporting Information *

■

Discussion of extra calculations without the geometric constraint at dilute concentrations to compare directly the results at dilute and high concentrations; figure showing site occupation rate as a function of the number of dopant. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

ACKNOWLEDGMENTS This work was supported by Grant in Aid for Scientific Research for Exploration of Nanostructure-Property Relationships for Materials Innovation. Fruitful discussion with Profs. Y. Ikuhara, N. Shibata, Dr. T. Nakagawa, and Mr. B. Feng at the University of Tokyo are gratefully acknowledged. We would like to acknowledge valuable discussion with Dr. C. A. J. Fisher.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: +81 6 6879 7476. Fax: +81 6 6879 7476. E-mail: tatsuya. [email protected].

REFERENCES

(1) Ikeda, J. A. S.; Chiang, Y.-M. Space Charge Segregation at Grain Boundaries in Titanium Dioxide: I, Relationship between Lattice

H

dx.doi.org/10.1021/la503338x | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

Defect Chemistry and Space Charge Potential. J. Am. Ceram. Soc. 1993, 76, 2437. (2) Ikeda, J. A. S.; Chiang, Y.-M.; Garratt-Reed, A. J.; Sande, J. B. V. Space Charge Segregation at Grain Boundaries in Titanium Dioxide: II, Model Experiments. J. Am. Ceram. Soc. 1993, 76, 2447. (3) Mulford, R. A. Grain Boundary Segregation in Ni and Binary Ni Alloys Doped with Sulfur. Metall. Trans. A 1983, 14A, 865. (4) Liu, C. T.; White, C. L.; Horton, J.A. Effect of Boron on GrainBoundaries in Ni3Al. Acta. Metall. 1985, 33, 213. (5) Rice, J.R.; Wang, J.-S. Embrittlement of Interfaces by Solute Segregation. Mater. Sci. Eng., A 1989, A107, 23. (6) Yan, Y.; Jiang, C.-S.; Noufi, R.; Wei, S.-H.; Moutinho, H. R.; AlJassim, M. M. Electrically Benign Behavior of Grain Boundaries in Polycrystalline CuInSe2 Films. Phys. Rev. Lett. 2007, 99, 235504. (7) Tsen, A. W.; Brown, L.; Levendorf, M.P.; Ghahari, F.; Huang, P. Y.; Havener, R. W.; Ruiz-Vargas, C. S.; Muller, D. A.; Kim, P.; Park, J. Tailoring Electrical Transport Across Grain Boundaries in Polycrystalline Graphene. Science 2012, 336, 1143. (8) Gupta, A.; Gong, G. Q.; Xiao, G.; Duncombe, P. R.; Lecoeur, P.; Trouilloud, P.; Wang, Y. Y.; Dravid, VP; Sun, J. Z. Grain-Boundary Effects on the Magnetoresistance Properties of Perovskite Manganite Films. Phys. Rev. B 1996, 54, 629. (9) Lin, M. N.; Hsu, H. S.; Lai, J. Y.; Guo, M. C.; Lin, C. Y.; Li, G. Y.; Chen, F. Y.; Huang, J. J.; Chen, S. F.; Liu, C. P.; Huang, J. C. A. Enhanced Ferromagnetism in Grain Boundary of Co-doped ZnO Films: A Magnetic Force Microscopy Study. Appl. Phys. Lett. 2011, 98, 212509. (10) Hsu, H. S.; Huang, J. C. A.; Chen, S. F.; Liu, C. P. Role of Grain Boundary and Grain Defects on Ferromagnetism in Co:ZnO Films. Appl. Phys. Lett. 2007, 90, 102506. (11) Chokshi, A. H. Diffusion, Diffusion Creep and Grain Growth Characteristics of Nanocrystalline and Fine-Grained Monoclinic, Tetragonal and Cubic Zirconia. Scr. Mater. 2003, 48, 791. (12) Swaroop, S.; Kilo, M.; Argirusis, C.; Borchardt, G.; Chokshi, A. H. Lattice and Grain Boundary Diffusion of Cations in 3YTZ Analyzed Using SIMS. Acta Mater. 2005, 53, 4975. (13) Rittner, J. D.; Seidman, D. N. Solute-Atom Segregation to ⟨110⟩ Symmetric Tilt Grain Boundaries. Acta. Mater. 1997, 45, 3191. (14) Lejcĕk, P.; Hofmann, S.; Paidar, V. Solute Segregation and Classification of [100] Tilt Grain Boundaries in α-Iron: Consequences for Grain Boundary Engineering. Acta Mater. 2003, 51, 3951. (15) Millett, P. C.; Selvam, R. P.; Bansal, S.; Saxena, A. Atomistic Simulation of Grain Boundary EnergeticsEffects of Dopants. Acta Mater. 2005, 53, 3671. (16) Braithwaite, J. S.; Rez, P. Grain Boundary Impurities in Iron. Acta Mater. 2005, 53, 2715. (17) Yamaguchi, M.; Shiga, M.; Kaburaki, H. Grain Boundary Decohesion by Impurity Segregation in a Nickel−Sulfur System. Science 2005, 307, 393. (18) Hui, S. R.; Roller, J.; Yick, S.; Zhang, X.; Decës-Petit, C.; Xie, Y.; Maric, R.; Ghosh, D. A Brief Review of the Ionic Conductivity Enhancement for Selected Oxide Electrolytes. J. Power Sources 2007, 172, 493. (19) Mogensen, M.; Lybye, D.; Bonanos, N.; Hendriksen, P. V.; Poulsen, F. W. Factors Controlling the Oxide Ion Conductivity of Fluorite and Perovskite Structured Oxides. Solid State Ionics 2004, 174, 279. (20) Steele, B. C. H.; Heinzel, A. Materials for Fuel-Cell Technologies. Nature 2001, 414, 345. (21) Verkerk, M. J.; Middelhuis, B. J.; Burggraaf, A. J. Effect of Grain Boundaries on the Conductivity of High-Purity ZrO2-Y2O3 Ceramics. Solid State Ionics 1982, 6, 159. (22) Badwal, S. P. S. Zirconia-Based Solid Electrolytes: Microstructure, Stability and Ionic Conductivity. Solid State Ionics 1992, 52, 23. (23) Arachi, Y.; Sakai, H.; Yamamoto, O.; Takeda, Y.; Imanishi, N. Electrical Conductivity of the ZrO2-Ln2O3 (Ln = Lanthanides) System. Solid State Ionics 1999, 121, 133.

(24) De Souza, R. A.; Pietrowski, M. J.; Anselmi-Tamburini, U.; Kim, S.; Munir, Z. A.; Martin, M. Oxygen Diffusion in Nanocrystalline Yttria-Stabilized Zirconia: The Effect of Grain Boundaries. Phys. Chem. Chem. Phys. 2008, 10, 2067. (25) Chen, X. J.; Khor, X. J.; Chan, S. H.; Yu, L. G. Influence of Microstructure on the Ionic Conductivity of Yttria-Stabilized Zirconia Electrolyte. Mater. Sci. Eng., A 2002, 335, 246. (26) Aoki, M.; Chiang, Y.-M.; Kosacki, I.; Lee, I. J.; Tuller, H.; Liu, Y. Solute Segregation and Grain-Boundary Impedance in High-Purity Stabilized Zirconia. J. Am. Ceram. Soc. 1996, 79, 1169. (27) Ramamoorthy, R.; Sundararaman, D.; Ramasamy, S. Ionic Conductivity Studies of Ultrafine-Grained Yttria Stabilized Zirconia Polymorphs. Solid State Ionics 1999, 123, 271. (28) Nakagawa, T.; Sakaguchi, I.; Shibata, N.; Matsunaga, K.; Yamamoto, T.; Haneda, H.; Ikuhara, Y. Oxygen Diffusion Blocking of Single Grain Boundary in Yttria-Doped Zirconia Bicrystals. J. Mater. Sci. 2005, 40, 3185. (29) Park, J. S.; Kim, Y. B.; An, J.; Prinz, F. B. Oxygen Diffusion Across the Grain Boundary in Bicrystal Yttria Stabilized Zirconia. Solid State Commun. 2012, 152, 2169. (30) Oyama, T.; Yoshiya, M.; Matsubara, H.; Matsunaga, K. Numerical Analysis of Solute Segregation at Σ5 (310)/[001] Symmetric Tilt Grain Boundaries in Y2O3-Doped ZrO2. Phys. Rev. B 2005, 71, 224105. (31) Yoshiya, M.; Oyama, T. Impurity and Vacancy Segregation at Symmetric Tilt Grain Boundaries in Y2O3-Doped ZrO2. J. Mater. Sci. 2011, 46, 4176. (32) Yoshiya, M.; Tanaka, I.; Adachi, H. Energetical Role of Modeled Intergranular Glassy Film in Si3N4SiO2 Ceramics. Acta. Mater. 2000, 48, 4641. (33) Yoshiya, M.; Tatsumi, K.; Tanaka, I.; Adachi, H. Theoretical Study on the Chemistry of Intergranular Glassy Film in Si3N4SiO2 Ceramics. J. Am. Ceram. Soc. 2002, 85, 109. (34) Yoshiya, M.; Tanaka, I.; Adachi, H. Atomic-Level Modeling and Computation of Intergranular Glassy Film in High-Purity Si3N4 Ceramics. J. Eur. Ceram. Soc. 2010, 101, 57. (35) Mao, Z.; Sinnott, S. B.; Dickey, E. C. Ab Initio Calculations of Pristine and Doped Zirconia Σ5 (310)/[001] Tilt Grain Boundaries. J. Am. Ceram. Soc. 2002, 85, 1594. (36) Marinopoulos, A. G. First Principles Study of Segregation to the Σ5(310) Grain Boundary of Cubic Zirconia. J. Phys.: Condens. Matter 2011, 23, 085005. (37) Lee, H. S.; Mizoguchi, T.; Mitsui, J.; Kang, J. L.; Yamamoto, T.; Ikuhara, Y. Characterization and Atomic Modeling of an Asymmetric Grain Boundary. Phys. Rev. B 2011, 83, 104110. (38) Sato, Y.; Mizoguchi, T.; Shibata, N.; Yamamoto, T.; Hirayama, T.; Ikuhara, Y. Atomic-Scale Segregation Behavior of Pr at a ZnO [0001] Σ49 Tilt Grain Boundary. Phys. Rev. B 2009, 80, 094114. (39) Sato, Y.; Mizoguchi, T.; Oba, F.; Ikuhara, Y.; Yamamoto, T. Arrangement of Multiple Structural Units in a [0001] Σ49 Tilt Grain Boundary in ZnO. Phys. Rev. B 2005, 72, 064109. (40) Dickey, E. C.; Fan, X.; Pennycook, S. J. Structure and Chemistry of Yttria-Stabilized Cubic-Zirconia Symmetric Tilt Grain Boundaries. J. Am. Ceram. Soc. 2001, 84, 1361. (41) Shibata, N.; Oba, F.; Yamamoto, T.; Sakuma, T.; Ikuhara, Y. Grain-Boundary Faceting at a Σ = 3, [110]/{112} Grain Boundary in a Cubic Zirconia Bicrystal. Philos. Mag. 2003, 83, 2221. (42) Shibata, N.; Oba, F.; Yamamoto, T.; Ikuhara, Y. Structure, Energy and Solute Segregation Behaviour of [110] Symmetric Tilt Grain Boundaries in Yttria-Stabilized Cubic Zirconia. Philos. Mag. 2004, 84, 2381. (43) Shibata, N.; Oba, F.; Yamamoto, T.; Ikuhara, Y.; Sakuma, T. Atomic Structure and Solute Segregation of a Σ = 3, [110]/{111} Grain Boundary in an Yttria-Stabilized Cubic Zirconia Bicrystal. Philos. Mag. Lett. 2002, 82, 393. (44) Lei, Y.; Ito, Y.; Browning, N. D.; Mazanec, T. J. Segregation Effects at Grain Boundaries in Fluorite-Structured Ceramics. J. Am. Ceram. Soc. 2002, 85, 2359. I

dx.doi.org/10.1021/la503338x | Langmuir XXXX, XXX, XXX−XXX

Langmuir

Article

(45) Li, P.; Chen, I.-W.; Penner-Hahn, J. E. X-ray-Absorption Studies of Zirconia Polymorphs. I. Characteristic Local Structures. Phys. Rev. B 1993, 48, 14. (46) Li, P.; Chen, I.-W.; Penner-Hahn, J. E. Effect of Dopants on Zirconia StabilizationAn X-ray Absorption Study: I, Trivalent Dopants. J. Am. Ceram. Soc. 1994, 77, 118. (47) Minervini, L.; Grimes, R. W.; Sickafus, K. E. Disorder in Pyrochlore Oxides. J. Am. Ceram. Soc. 2000, 83, 1873. (48) Zacate, M. O.; Minervini, L.; Bradfield, D. J.; Grimes, R. W.; Sickafus, K. E. Defect Cluster Formation in M2O3-Doped Cubic ZrO2. Solid State Ionics 2000, 128, 243. (49) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087. (50) Guo, X. Physical Origin of the Intrinsic Grain-Boundary Resistivity of Stabilized-Zirconia: Role of the Space-Charge Layers. Solid State Ionics 1995, 81, 235. (51) Guo, X.; Sigle, W.; Fleig, J.; Maier, J. Role of Space Charge in the Grain Boundary Blocking Effect in Doped Zirconia. Solid State Ionics 2002, 154-155, 555. (52) An, J.; Park, J. S.; Koh, A. L.; Lee, H. B.; Jung, H. J.; Schoonman, J.; Sinclair, R.; Gür, T. M.; Prinz, F. B. Atomic Scale Verification of Oxide-Ion Vacancy Distribution near a Single Grain Boundary in YSZ. Sci. Rep. 2013, 3, 2680. (53) Lee, B. H.; Prinz, F. B.; Cai, W. Atomistic Simulations of Grain Boundary Segregation in Nanocrystalline Yttria-Stabilized Zirconia and Gadolinia-Doped Ceria Solid Oxide Electrolytes. Acta. Mater. 2013, 61, 3872. (54) Pornprasertsuk, R.; Ramanarayanan, P.; Musgrave, C. B.; Prinz, F. B. Predicting Ionic Conductivity of Solid Oxide Fuel Cell Electrolyte from First Principles. J. Appl. Phys. 2003, 94, 174. (55) Bogicevic, A.; Wolverton, C.; Crosbie, G. M.; Stechel, E. B. Defect Ordering in Aliovalently Doped Cubic Zirconia from First Principles. Phys. Rev. B 2001, 64, 014106. (56) Dalach, P.; Ellis, D. E.; van de Walle, A. First-Principles Thermodynamic Modeling of Atomic Ordering in Yttria-Stabilized Zirconia. Phys. Rev. B 2010, 82, 144117. (57) Gale, J. D. GULP: A Computer Program for the SymmetryAdapted Simulation of Solids. J. Chem. Soc. Faraday Trans. 1997, 93, 629. (58) Due to the small gradient of chemical potential with respect to the dopant content in bulk phases, having an undoped grain inteiror does not have a noticeable impact on energetics of GBS as well as point defect configurations according to a previous study.30

J

dx.doi.org/10.1021/la503338x | Langmuir XXXX, XXX, XXX−XXX