Revisiting the Temperature Dependent Ionic Conductivity of Yttria Stabilized Zirconia (YSZ)

The temperature dependent conductivity of yttria stabilized zirconia (YSZ) exhibits a bending in Arrhenius’ plots which is frequently discussedintermsoffreeandassociatedoxygenvacancies.However,theveryhighdopingconcentrationinYSZleadstosuchastrongdefectinteractionthattheconceptoffreevacanciesbecomeshighlyquestionable.Therefore,thetemperaturedependentconductivityofYSZisreconsidered.TheconductivityofYSZwithdifferentdopingconcentrationwasmeasuredinabroadtemperaturerange.ThedataareanalyzedintermsoftwodifferentbarrierheightsthathavetobepassedalonganaveragepathofanoxygenvacancyinYSZ(twobarriermodel).For8–10mol%yttria,thetwobarriersareintherangeof0.6eVand1.1–1.2eV,respectively.Theconductivityandthusthebarrierheightsalsodependonthecoolingrateafterahightemperaturepre-treatment.Thisindicatesthatdifferentfrozen-indistributionsofdopantsaffectthevacancymotionbydifferentenergylandscapes.Temporarilyexistingdefectconﬁgurations,possiblywithastrongeffectofrepulsiveoxygenvacancyinteraction,aresuggestedasthereasonofhighbarriers.Futuredynamicab-initiocalculationsmayrevealwhetherthismodiﬁedmodeloftheYSZconductivityismechanisticallymeaningful.©TheAuthor(s)2017.PublishedbyECS.ThisisanopenaccessarticledistributedunderthetermsoftheCreativeCommonsAttribution4.0License(CCBY,http://creativecommons.org/licenses/by/4.0/),whichpermitsunrestrictedreuseoftheworkinanymedium,providedtheoriginalworkisproperlycited.[DOI:10.1149/2.0641707jes]Allrightsreserved. stabilized zirconia (YSZ) is among the most important ion conducting solids and acts as a kind of model material representing fast oxide ion conductors. Owing to this model character of YSZ, but also due to its application in solid oxide fuel cells (SOFCs), solid oxide electrolysis cells (SOECs) and oxygen sensors, a vast amount of papers can be found dealing with its oxide ion conduction. The ionic conductivity is based on the motion of oxygen vacancies, introduced by Y 3 + ions replacing Zr 4 + . For concentrations above ca. 8 mol%, yttria doping also stabilizes the cubic structure down to room temper- ature. A detailed review of the science of YSZ and related materials

Yttria stabilized zirconia (YSZ) is among the most important ion conducting solids and acts as a kind of model material representing fast oxide ion conductors. Owing to this model character of YSZ, but also due to its application in solid oxide fuel cells (SOFCs), solid oxide electrolysis cells (SOECs) and oxygen sensors, a vast amount of papers can be found dealing with its oxide ion conduction. The ionic conductivity is based on the motion of oxygen vacancies, introduced by Y 3+ ions replacing Zr 4+ . For concentrations above ca. 8 mol%, yttria doping also stabilizes the cubic structure down to room temperature. A detailed review of the science of YSZ and related materials is far beyond the scope of this paper but a few important facts regarding the ionic conductivity of zirconia-based solid electrolytes can be briefly summarized as follows: 1-9 i) Doping concentrations above ca. 8 mol% Y 2 O 3 lead to a decrease of the conductivity, despite increasing oxygen vacancy concentration. ii) The conductivity not only depends on the vacancy concentration but also on the kind of dopant. For example, Sc-doped zirconia shows significantly higher conductivity than YSZ. iii) The temperature dependence of the conductivity cannot be described by a single activation energy (E act ) but shows higher E act values at lower temperatures.
Numerous theoretical studies were performed in order to understand these experimental observations and to get a deeper insight into defect thermodynamics and kinetics of doped zirconia and of the closely related ceria-based ion conductors, see e.g. Refs. 10-29. Those model studies employed different simulation approaches such as molecular dynamics (MD), density functional theory (DFT) and kinetic Monte Carlo. Particularly the effect of dopant concentration and dopant ion on the conductivity was often in the focus of research. The conductivity maximum with increasing dopant content was reproduced in many simulation studies on zirconia and ceria, e.g. in Refs. [10][11][12][13][14][15][16][17]19,[23][24][25][26]28 and is often discussed in terms of spatially varying migration barriers. For example, the vacancy migration barrier strongly depends on the two cations of the tetrahedral edge, which is crossed by an oxide ion during a jump to a nearest neighbor site of the anion sub-lattice: [13][14][15]20,23,26 Refs., 13,14,26 report barriers between 0.3 and 0.67 eV for a Zr-Zr edge and 0.85-1.29 eV for a Y-Zr edge. An increased number of such high barrier edges lowers the effective mobility for increasing doping concentrations.
In the course of such calculations, many data were collected on defect formation energies, defect interaction energies and defect migration barrier heights. Simulations on defect interaction often con-sidered defect pairs and found that energies are lower for an Y ion in the NNN (next nearest neighbor) cation shell of an oxygen vacancy, compared to Y in the NN (nearest neighbor) cation tetrahedron of the vacancy or interaction-free vacancies, even though calculated absolute values of defect interaction energies vary considerably. 11,13,18,19,22,24 This defect interaction also strongly depends on the dopant ion. 15,23,24 Defect energies were further calculated for more extended defects of several dopant ions and oxygen vacancies, [30][31][32] and such studies as well as experimental data 33 showed the relevance of vacancy-vacancy interaction. The importance of vacancy-vacancy interaction was also emphasized in Refs. 13,15,28,29,34,35 and ab initio metadynamic simulations revealed surprising effects due to complex (multi-)defectdefect interaction such as locally unstable sites for vacancies (several ten percent of the sites) and concerted vacancy jumps. 36 Interestingly, only few numerical simulation studies dealt with the temperature dependent activation energy of YSZ, e.g., 13,14,26,27,37,38 An activation energy lowering at high temperatures was reported in Refs. 27,37, even though calculated changes were smaller than the experimentally measured ones. In Ref. 37 this activation energy change was traced back to asymmetric barriers in YSZ due to different energies of final and initial states. Also a slight bending of the conductivity curve into the opposite direction was reported in some cases, 14,26,38 For doped ceria, an activation energy change similar to the experimentally observed one was found in kinetic Monte Carlo simulations using several types of defect-defect interactions and barrier heights. 28 The temperature dependence of the conductivity in YSZ (and doped ceria) seems to be a result of a complex energy landscape caused by multiple defect interactions including several neighboring shells and possibly larger defect clusters.
This absence of a clear picture of the temperature dependence in numerical simulations is in contrast to the fact that the different activation energies at high and low temperatures is often described as well understood in experimental studies. It is interpreted in terms of a defect interaction model with free (mobile) vacancies and immobile vacancies trapped in defect complexes (clusters, associates) such as (Y F791 shows a non-Arrhenius-type temperature dependence since the concentration of mobile (silver) vacancies (n) is temperature dependent. In Eq. 1, u is the mobility of the free vacancies, z denotes the (absolute) charge number of the vacancies and e 0 represents the elementary charge.
More specific, at low temperatures a large fraction of the vacancies does not contribute to the conductivity since they are trapped in immobile associates. With increasing temperature more silver vacancies become liberated and thus the ionic conductivity increases not only because of the thermally activated mobility u but also due to the temperature dependence of the free vacancy concentration n. At high temperatures virtually all silver vacancies are free and the activation energy of σT only reflects the energy barrier (E a ) of the (free) vacancy migration, provided the total vacancy concentration is fixed by doping. For silver halides very detailed conductivity measurements with varying doping levels allowed determination of defect chemical parameters such as defect formation enthalpy and entropy, association enthalpy and entropy, or defect migration energies of silver vacancies and interstitials. [44][45][46][47] This analysis partly includes further long range electrostatic interaction between defects, similar to the defect cluster model with Debye Hückel-type electrostatic interaction applied to doped NaCl. 48 In analogy with these studies, defect association of a vacancy with one or two Y ions is often assumed to immobilize a large number of oxygen vacancies in YSZ at low temperatures. At high temperatures, on the other hand, oxygen vacancies in YSZ are believed to be free. Hence, also for YSZ the measured activation energies were used to determine migration barriers of free vacancies and defect association enthalpies. Numerical studies frequently compared their simulated migration and interaction energies with these supposedly "measured" association energies.
However, even though such a defect association model is excellently suited to describe dilute systems such as slightly doped halides, it becomes highly questionable for YSZ or ceria with high doping concentrations. This problem was already mentioned in Refs. 12,49 and is discussed in more detail for ceria in Ref. 50 (in ceria, dilute solutions are accessible without changing the crystal structure). A doping of 8 mol% Y 2 O 3 , for example, corresponds to 16/108 Y ions on cation sites, i.e. 14.8% Y cations, compared to a few 10 ppm doping used in the defect interaction studies on halides. In case of a random distribution of these Y ions, about 47% of all cation tetrahedrons include an Y ion 10 and this number further increases when considering repulsive Y-Y interaction. 11,19 In a random case, about 92% of all oxygen sites have at least one Y ion in the NN or NNN cation shell and interaction of vacancies with Y in the NN as well as the NNN cation shell is very pronounced. Even interaction far beyond the second neighbor shell is important. 13,26 All this indicates that in highly doped oxides interaction-free oxygen vacancies hardly exist at all. Almost all vacancies exhibit significant interaction with Y dopants and neither a trapped nor a free vacancy is a well-defined state in YSZ. Over and above, the presumably very important vacancy-vacancy interaction is neglected in the simple defect association model. Therefore, it is more appropriate to consider vacancy motion in YSZ as a motion between sites with differently pronounced but always present defect interaction. Accordingly, also the meaning of the energies deduced from temperature dependent experimental conduc-tivities has to be reconsidered, see e.g. Ref. 26,28 . Neither does the  high temperature activation energy represent the barrier height of a  free vacancy jump nor does the difference to the low temperature activation energy yield the interaction enthalpy of a single vacancy with  one or two Y ions. Consequently, information on energies of isolated  vacancies or simple vacancy-dopant complexes cannot be deduced  from measurements on cubic YSZ. Molecular dynamics calculations, kinetic Monte Carlo simulations or ab-initio metadynamic modeling studies on YSZ do not distinguish mobile and immobile vacancies. Rather, they consider ion motion in YSZ as a dynamic process of all vacancies (or oxide ions) in a complex energy landscape. This landscape may even change with time if vacancy-vacancy interaction is included. In terms of Equation 1 this corresponds to an approach that considers all vacancies mobile (n = total concentration of vacancies n V,tot ). The defect mobility is then not a simple property of oxygen vacancies but varies in space and time; only a kind of effective mobility can be defined. However, this view on ion motion in YSZ is usually not applied by experimentalists when analyzing measured temperature dependent conductivity data, see above.
It is the scope of this paper to reconsider measured temperature dependent conductivities of YSZ by employing this different view on vacancy motion in terms of a modified model for data analysis. In this model, all vacancies are assumed to be mobile and they pass, on average, two types of barriers while diffusing in YSZ (two barrier model). The model is applied to a series of temperature dependent conductivity measurements (220-950 • C) on differently doped single and polycrystalline YSZ samples (8-11.5 mol% yttria). Conductivities were also analyzed after high temperature treatments (up to 1550 • C), followed by slow or fast cooling. This pre-treatment may affect the Y distribution in YSZ and thus the energy landscape for vacancy jumps, similar to aging of YSZ at temperatures around or above 1000 • C. [51][52][53][54][55] Future simulation studies with spatially and temporarily varying defect interaction and defect motion are encouraged to use the fit parameters of our two barrier analysis for the sake of comparison. Those could also clarify whether a dynamic bimodal barrier height distribution can indeed approximate the vacancy motion in YSZ, and whether our experimentally determined activation energies and prefactors have simple meanings.

Experimental
Single crystalline YSZ samples were purchased from different suppliers (CrT = CrysTec, Berlin, Germany, MaT = MaTeck, Jülich, Germany, UW = University Wafer, Boston, USA) with sizes of 5 × 5 or 10 × 10 mm 2 and a thickness of 0.5 mm. In case of CrT, single crystals from three different shipments (CrTA, CrTB, CrTC) were investigated. Additional samples (0.6 mm thickness) were cut from a single crystal grown at Max Planck Institute of Solid State Research, Stuttgart, Germany (MPS). Polycrystals (Tos) were prepared from nominally 8 mol% doped YSZ powder (Tosoh, Japan) and sintered at 1700 • C for 12 hours. Grain sizes (calculated by the line intersection method on a polished surface) were about 18 μm. The exact doping concentration was analyzed by X-ray fluorescence spectroscopy (Axios Advanced, PANalytical, Eindhoven, Netherlands) and the results are summarized in Table I. Lattice parameters were calculated  from X-ray diffraction measurements obtained on a X'Pert Pro Powder Diffractometer (PANalytical, Eindhoven, Netherlands). Values are also given in Table I and show a slight increase for increasing dopant concentration, in accordance with literature. 6,56,57 For conductivity measurements on single crystals, Pt paste was painted onto the two sides and annealed at 1000 • C for two hours in air. In case of YSZ polycrystals, the grain boundary and Pt electrode arcs showed a significant overlap in impedance spectra and therefore La 0.6 Sr 0.4 CoO 3-δ (LSC) thin film electrodes were employed. Those are known to exhibit very large chemical capacitances, 58 which leads to a much better separation of the grain boundary and electrode impedance contributions. LSC layers of 200 nm thickness were prepared by pulsed laser deposition (PLD) at 0.4 mbar O 2 background pressure and a substrate temperature of 650 • C. 59 Despite possible solid state reactions between LSC and YSZ, the electrodes turned out to be usable in the entire temperature range (up to ca. 950 • C) within the time frame employed in this study. Additional Pt paste on the LSC layer (annealed at 700 • C for 5 hours) ensured that the limited in-plane conductivity of the LSC layer did not cause a significant resistance.
Impedance spectra of the samples were measured in the temperature range of ca. 220-950 • C by a Solartron SI 1260 impedance analyzer in the frequency range from 0.1 Hz to 3 MHz. The exact temperature range used for impedance measurements varied slightly between the different samples. In all cases, a complete measurement cycle included impedance measurements with increasing as well as decreasing temperature. A spectrum was taken 10 minutes after reaching a temperature stability of ±0.1 • C. The temperature was measured by a type K thermocouple positioned very close to the sample. For eliminating impedance contributions of the wires, four point measurements were performed: Two wires were attached to each Pt sheet used for contacting the electrodes. The remaining ohmic resistance due to wiring was below 0.05 at all temperatures-the corresponding systematic error can thus be safely considered to be lower than 2%.

Impedance Spectra and Determination of Conductivities
Single crystalline YSZ.-At low and intermediate temperatures, impedance spectra of single crystalline YSZ samples exhibit a high frequency arc in the complex impedance plane, which reflects ion conduction in YSZ bulk, and a low frequency arc due to electrochemical oxygen exchange at the electrodes. Fig. 1 displays examples of such spectra. For higher temperatures, increasingly less of the high frequency arc is visible within the given frequency range (Figs. 1a-1d) while the contribution of the electrode impedance feature becomes more complete. At temperatures above ca. 650 • C, inductive contributions dominate at high frequencies and at low frequencies the electrode impedance almost reaches the real axis of the impedance plot (Fig. 1e). It is also clearly visible that the electrode impedance cannot be described by a single arc.
At first sight, determination of the ionic conductivity from such spectra seems straightforward. However, at least a few minor problems exist: i) The depressed arc of the YSZ bulk impedance in the spectrum shown in Fig. 1b (310 • C) can be quantified by a simple fit to a resistor in parallel with a constant phase element (CPE). Since the data points close to the |Z | minimum between bulk and electrode arc are affected by both electrode and bulk impedance, they have to be excluded from such a simple fit of the first arc. Therefore, some ambiguity remains which points to include in the fit of the bulk arc. ii) At lower temperatures a fit using one bulk R-CPE becomes rather non-ideal (Fig. 1a). Moreover, for all temperatures the constant phase element of the bulk arc has an exponential factor of about 0.85-0.90 and is thus rather far from purely capacitive ( = 1). An interpretation is beyond the scope of this paper and also in literature not much discussion exists on this dispersion in YSZ, see e.g. in Refs. 60-63. iii) In the high temperature regime without any bulk arc (Figs. 1d and 1e), one would need an accurate circuit model of the electrode impedance for exactly determining the serial YSZ bulk resistance.
This illustrates that a very precise determination of the ionic conductivity from the given data is not straightforward. However, keeping in mind the very strong temperature dependence of the conductivity and a certain inaccuracy of the temperature measurement, errors of a few percent are acceptable in the conductivity analysis and not relevant for any conclusion drawn in this paper. In order to avoid ambiguity in the choice of fit models (e.g. one or two R-CPE elements for describing the bulk impedance, varying number of R-CPE-elements for quantifying the electrode impedance) we have chosen the following pragmatic approach for extrapolating the YSZ bulk conductivity: i) As long as a part of the bulk arc is visible, the minimum value of the imaginary part between bulk and electrode impedance was identified and the corresponding real part of the impedance was used as bulk resistance. This is illustrated in Fig. 1f where it is also shown that the resulting resistance value is very close to that obtained from a fit using a single R-CPE element to quantify the high frequency part of the spectrum. ii) At high temperatures with a substantial inductive contribution, the intercept of the spectrum with the real axis was taken as the bulk resistance. iii) In the small temperature range without a minimum between bulk and electrode arc and without a real axis intercept (around 550-600 • C, see spectrum in Fig. 1d), the smallest real part of the impedance spectrum (i.e. that of the highest frequency) was taken as bulk resistance. From these bulk resistances (R) the YSZ bulk conductivity σ was calculated according to σ = L/(R · A) with L and A denoting the sample thickness and area.
Polycrystalline YSZ.-It is very often reported in literature (see e.g. Refs. 52,64-67) that at lower temperatures polycrystalline YSZ samples exhibit two arcs in the complex impedance plane, in addition to the electrode arc(s). The high frequency arc can be interpreted in terms of ion conduction in the grain interior (YSZ bulk), while the intermediate frequency arc originates from resistive grain boundaries. Two serial R-CPE elements are often used to quantify the corresponding resistances and a brick layer model [67][68][69] is then employed to determine bulk and grain boundary conductivities. Also in our measurements these two arcs are partly or completely visible below ca. 500 • C, but the bulk arc is always much larger than the grain boundary arc, see Fig. 2a. At low frequencies the impedance response of the LSC electrodes is visible. For higher temperatures, the onset of the bulk arc disappears (Fig. 2b) and above ca. 650 • C also the onset of the grain boundary arc can no longer be measured (Fig. 2c). This raises the question of how to analyze the ionic bulk conductivity of our polycrystalline YSZ in this case.
A separation into grain and grain boundary is possible at lower temperatures, even though some accuracy problems come again into play, similar to those discussed for single crystals: Two serial R-CPE elements do not exactly fit the transition between grain bulk and grain boundary arc and neglect any electrode contribution in the transition range from intermediate to low frequency (electrode) arc. Moreover, for higher temperatures the bulk arc completely disappears and quantities obtained by fitting with a (bulk) resistor in series with an R-CPE element for grain boundaries (R1+R2-CPE-fit, see Fig. 2d) depend somewhat on the frequency range chosen for the fit analysis. Above ca. 600 • C a meaningful separation into bulk and grain boundary contribution becomes impossible at all, since the remaining part of the grain boundary arc is too small to allow for a reasonable fit to a R1+R2-CPE model. Thus, only the total (effective) conductivity of polycrystalline YSZ can be obtained.
In this study, it is of essential importance to use one and the same analytical equation to describe the conductivity of YSZ for all temperatures from ca. 220 to 950 • C (see below). It would be inconsistent to consider the grain bulk conductivity of polycrystalline YSZ at lower temperatures and the total (effective) conductivity at higher temperatures. Hence, we have chosen the approach of considering only the total conductivity of polycrystalline YSZ, which is accessible for all temperatures. This total conductivity was determined as for single crystals (see above) except that here the |Z | minimum between grain boundary arc and the electrode impedance was used (instead of the minimum between bulk and electrode arc for single crystals). One may argue that the relevance of a comparison of single and polycrystals is then limited. However, in our case -also owing to the rather large grains obtained by 12 hours sintering at 1700 • C-grain boundary effects are small and thus also the difference between total and bulk conductivity is small, particularly compared to the huge conductivity variations caused by temperature. This is shown in Fig. 3, where grain bulk (σ grain ) and total conductivity (σ total ) of polycrystalline YSZ are compared up to 600 • C, i.e. in the temperature range for which a separation into bulk and grain boundary is indeed possible in terms of the brick layer model (R1+R2-CPE circuit). Absolute values and activation energies of bulk conductivity and total conductivity are so close that none of the effects discussed in this paper is compromised by our simplified approach. Fig. 4 Arrhenius plots (log(σT) vs. 1/T, T = temperature) are shown for all samples in the entire temperature range (Fig. 4a) and emphasizing the high and low temperature ranges (Figs. 4b and 4c). The often reported change of the activation energy with temperature is clearly visible also for our samples. Accordingly, a temperature independent acti-vation energy cannot describe the data. Moreover, the different doping concentrations lead to substantially different conductivities. Differences are less pronounced at high temperatures but amount to almost a factor of six at the lowest temperatures considered here (Fig. 4c), despite doping concentration differences of only a factor of 1.4. Owing to the strong temperature dependence of the conductivity, a precise analysis of the dopant concentration dependence of σ requires data at the very same temperature. Exactly matching sample temperatures were usually not achieved in our temperature cycles since set temperature of the furnace and sample temperature slightly differed. Therefore, the conductivities of all samples were parameterized in the entire temperature range by an analytical equation and the resulting fit parameters were then used to interpolate the conductivity to any temperature. Details on this parameterization are given below. In Fig. 5, interpolated conductivities are shown for 1000 • C, 600 • C and 400 • C.

Temperature and doping dependence of the conductivity.-In
The general trend is the same for all temperatures: The highest conductivity is found for 8 mol% Y 2 O 3 (polycrystalline sample) and the conductivity decreases with increasing doping level. Differences between maximal and minimal σ-values become less pronounced with increasing temperature, in accordance with literature: 6,8 At 400 • C the 8 mol% doped YSZ exhibits more than a factor of four higher conductivity than 11.5% doped YSZ while at 1000 • C the conductivity variation reduces to a factor of two. Moreover, when only considering 9-11 mol% yttria (single crystalline samples), the trend of decreasing conductivity with increasing doping level is almost absent at 1000 • C (dashed line in Fig. 5a).

Effect of thermal pre-treatment on the conductivity.-Different
Y distributions are expected to cause different energy landscapes for oxygen vacancy motion and might therefore affect the ionic conduc- Figure 3. Arrhenius plot of the effective (total) conductivity and the grain conductivity evaluated for a polycrystalline YSZ sample. Differences are marginal due to the very small grain boundary contribution. tivity of YSZ. 27 The exact Y ion distribution in a sample should depend on the temperature at which the cation sub-lattice becomes frozen-in during processing Hence, the conduction properties might depend on annealing conditions at high temperatures, cf. Ref. 6. This was investigated in the following measurement series: A YSZ single crystal (UW) was repeatedly heated to 1550 • C several times, annealed for 90 minutes and cooled to lower temperatures at different cooling rates. After each cooling experiment the conductivity of the sample was measured at temperatures much lower than the supposed freezing-in temperatures of the Y ion distribution; sample changes on the time scale of the conductivity measurements were therefore not observed.
The conductivity vs. T −1 curves were measured three times: after annealing and slow cooling from 1550 • C (0.1 • C/min), after an additional annealing step with subsequent fast cooling (10 • C/min), and after a third annealing followed by slow cooling. The results are plotted in Fig. 6 in the high and low temperature range. Obviously, the fast cooling procedure changed the sample: YSZ became better conductive at very low temperatures but less conductive at higher temperatures. This conductivity change was reversible and the conductivity originally obtained after slow cooling could be restored by another annealing and slow cooling. Especially this reversibility supports the assumption that Y ions become mobile at the annealing temperature and that their spatial distribution affects the ionic conductivity of YSZ. This spatial distribution is determined by the temperature at which the cation sub-lattice becomes frozen-in during cooling.

Two Barrier Model for Data Analysis
A quantitative analysis of the temperature dependent conductivity as well as of its dependence on doping concentration and heattreatment requires a model for data parameterization (here σT vs T −1 ). In the introduction it was discussed why we consider the standard defect association model not as an appropriate approach for analyzing the experimental data of highly doped ion conductors. In the following, we employ an alternative analytical model in which all vacancies  the barriers being successfully passed along an average vacancy path. For example, much more high barriers may be present than suggested by their factor β i since oxygen vacancies may avoid most of those barriers by detours.
A successful jump across a certain type of barrier i takes place at a jump frequency ν i which depends on temperature T according to The pre-factor ν 0 i includes the attempt frequency of the corresponding jump, but also a migration entropy term and correlation effects, k = Boltzmann's constant. Barriers of the same height but very different ) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 128.131.44.51 Downloaded on 2020-01-08 to IP pre-factors might be considered separately. For an average pathway of a vacancy, an effective jump frequency υ eff can be defined by the total number (N) of successful jumps of a vacancy per time t: . [3] Symbol N i denotes the number of successful jumps across barriers i and τ i is the residence time, i.e. the time needed until a successful jump takes place (τ i = ν −1 i ). With we may write Eq. 3 as An effective vacancy diffusion coefficient D eff can be defined from this effective jump frequency by Symbol a 0 denotes the average jump distance ( = 1/2 of the lattice constant of YSZ). From Eqs. 5, 2 and Nernst -Einstein's relation with effective mobility u eff u eff = ze 0 D eff kT = ze 0 a 2 0 ν eff kT [7] we obtain the (effective) conductivity: σ = z · e 0 · u eff · n V,tot = z 2 e 2 0 a 2 0 n V,tot kT Symbol n V,tot denotes the total vacancy concentration. Please note that correlation effects may be present, despite absence of an explicit correlation factor in these equations; rather they are included in the ν 0 i pre-factors. When defining γ i = z 2 e 2 0 a 2 0 n V,tot ν 0 i kβ i [9] Eq. 8 becomes [10] This equation may be used to fit experimental conductivity data with 2m independent fit parameters (E a,i and γ i ) provided E a,i and γ i are temperature independent; m denotes the number of different barrier types.
In this terminology the conductivity of a virtual material with only barriers of height E a,i (i.e. β i = 1) is given by [11] and the conductivity of YSZ can therefore also be expressed as [12] Hence, our material behaves like a sample with a series connection of regions with different conductivities σ i . The fraction β i determines how much each "region" contributes to the total conductivity. However, it should be emphasized that the vacancy motion is not onedimensional and thus the β i values do not have the simple meaning of a thickness fraction as one might conclude from a 1D interpretation of Eq. 12.
Still, it has to be decided how many barrier types should be used for analyzing an experimental data set. The analysis of our YSZ data sets showed that two different barriers and thus four parameters are already sufficient for an excellent fit to all measured conductivity data, see Fig. 7a. Eq. 10 then simplifies to [13] For E a,1 < E a,2 and γ 2 >> γ 1 , Eq. 13 leads to σT-T −1 plots with a temperature dependent activation energy: E a,1 determines the high temperature range and E a,2 is found at low temperatures. Please note that for similar γ i values the entire temperature range would be characterized by the larger E a,2 . When considering the ratio γ 2 /γ 1 one should keep in mind that it not only depends on the barrier fraction β 1 /β 2 but also on the jump frequency pre-factors ν 0 i . According to Eq. 9 it reads Still γ 2 >> γ 1 most probably correlates with β 1 >> β 2 since ν 0 i pre-factors hardly differ by several orders of magnitude and thus much more low barriers are successfully passed on an average path of a vacancy in YSZ with γ 2 >> γ 1 . In such a case, the two limiting activation energies for high and low temperature can be easily understood: At high temperatures the much fewer high barriers are irrelevant due to the high thermal energy of the jumping ions and the ) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 128.131.44.51 Downloaded on 2020-01-08 to IP activation energy of σT is simply the lower barrier height E a,1 . At lower temperatures it becomes increasingly difficult to pass the high barriers (strong increase of the e −E a,2 /kT factor), and even though only a small number has to be passed, the high barriers determine the low temperature activation energy, which is thus E a,2 .
If both barriers had the same pre-factor ν 0 i factor one could use the fit parameters γ i to determine the ratio of barrier fractions from γ 1 /γ 2 = β 2 /β 1 . Then γ 2 /γ 1 reflects the ratio of high and low barriers within an average sequence of successful vacancy jumps in YSZ. However, β 2 /β 1 does not provide much quantitative information on the density of each barrier in the three dimensional energy landscape. High barriers can often be avoided by a detour via low barriers and thus β 2 /β 1 ratios might strongly underestimate the true number of high barriers in the 3D system. A discussion on possible atomistic meanings of these two types of serial barriers is given below.
Please note that, despite applicability of Eq. 13, a significant amount of barriers with an even lower activation energy than E a,1 might be present in an average vacancy path without affecting the temperature dependence. Owing to their larger e −Ea,i/kT term in Eq. 10 they hardly hinder the measurable ion transport as long as their β i fraction is not much larger than β 1 . Moreover, it should be emphasized that applicability of Eq. 13 to an experimental data set does not validate the existence of such a two barrier energy landscape.

Quantitative Conductivity Data Analysis
Analysis of the dopant dependence.-Eq. 13 was used to parameterize all measured σT-T −1 curves with four fit parameters (E a,1 , E a,2 , γ 1 , γ 2 ). Fig. 7a shows that an excellent quantitative agreement of fit curves and experimental data is found. All fit parameters are summarized in Table II and these parameters were also used to calculate the conductivity values for 400, 600 and 1000 • C mentioned in the analysis of the doping dependence (Fig. 5). Fig. 8a displays the dopant concentration dependence of the two activation energies. High energy barriers are around 1.1-1.25 eV and the low barrier height varies between 0.57 and 0.83 eV. In both cases, the values increase for increasing dopant concentration and this effect is slightly more pronounced for the lower barrier which determines the high temperature range. The much stronger conductivity variation found for lower temperatures when increasing the doping concentration (see Fig. 5) is largely due to the fact that a change in activation energy has a much stronger effect at lower temperatures. The significant barrier height difference ( E a = E a,2 -E a,1 ) of ca. 0.42-0.6 eV reflects the substantial bending of the curves in the Arrhenius plots. In Ref. 70 the same analysis based on Eq. 13 was performed for 10 mol% doped YSZ, and the fit results were very similar (E a,1 = 0.682 eV, E a,2 = 1.168 eV, γ 1 = 1.1 · 10 5 K S/cm, γ 2 = 2.2 · 10 7 K S/cm); however, an interpretation in terms of our two barrier model was not made.
The two γ i factors as well as their ratios are plotted versus doping concentration in Fig. 9. The γ factors of the lower barriers (γ 1 ) are more than two orders of magnitude smaller than γ 2. Supposed the meaning of all factors in γ i is that suggested in our two barrier model (see Eq. 9), ν 0 1 values only depend on lattice vibration frequency and the correlation factor; therefore we assume that all ν 0 i values are within one order of magnitude. Hence, the large γ 2 /γ 1 ratios is mainly due to a much larger number of low barriers passed in an average sequence of successful vacancy jumps (β 1 >> β 2 ), see Eq. 14. We may thus assume β 1 ∼ = 1, and ν 0 1 values then resulting from Eq. 9 are shown in Table II for measured lattice constants and n V,tot being half of the measured Y concentration. Pre-factors ν 0 1 around 10 13 s −1 are found, which is within the order of magnitude of typical lattice vibration frequencies in solids.
The γ 1 factors show a strong increase for doping concentrations above ca. 10 mol% yttria (Fig. 9a). For β 1 ∼ = 1 one might therefore conclude that ν 0 1 strongly increases for higher doping, possibly due to correlation effects (cf . Table II). However, one should keep in mind that possibly also barriers with activation energies lower than E a,1 are present, see above. Hence, β 1 <1 cannot be excluded and also an increasing fraction of additional (invisible) barriers with E a,i < E a,1  (i.e. a decrease of β 1 ) could explain the trend. The factor γ 2 shows an almost linear increase with increasing doping concentration (see Fig. 9b). However, since the unknown barrier fraction β 2 as well as the pre-factor ν 0 2 enter γ 2 this is not further interpreted.
Analysis of the effect of thermal pretreatment.-Also the heattreatment dependent conductivities were analyzed by the two barrier model (Eq. 13) and the results are summarized in Table III. After the first heat-treatment with slow cooling, absolute conductivity values for all temperatures and hence also fit parameters were only slightly changed compared to the as-received sample; only γ 2 showed an increase by about a factor of two. Annealing and fast cooling, however, changed all fitting parameters: γ 1 increased by a factor of two while γ 2 decreased to ca. 1/3 of its original value. Accordingly, the ratio γ 2 /γ 1 decreased by a factor of 6. The change of the lower activation energy (E a,1 ) is also very pronounced, with values after fast cooling being lowered by 0.15 eV. The high barrier E a,2 is only slightly changed by ca. 0.05 eV. Restoration of the original conductivity curve by another slow cooling step was confirmed by the very similar fit parameters of the two data sets obtained after slow cooling.
Such changes are in accordance with the assumption that different cooling rates lead to different temperatures at which the cation sublattice becomes frozen-in and thus also to different energy landscapes. Most probably the cation distribution for a fast cooling rate is more random-like and the changed parameters (e.g. lower E a,1 value) are a consequence of this fact. A more detailed interpretation, e.g. of changes in γ i is again difficult, since those may be caused either by β i and/or ν 0 i contributions. Some temperature gradient or spatially varying cooling rates during the single crystal growth might also be the reason why nominally identical single crystals (CrTA, CrTB, CrTC) not only exhibit slightly different doping levels, but partly also deviate from trends shown in Figs. 8, particularly the sample with the lowest doping level (CrTA). However, here a smaller activation energy E a,1 comes along with a smaller value of γ 1 in contrast to the effect of the fast cooling rate employed in our study. Similar variations might exist also for the other single crystals used here, but further systematic investigations have not been performed yet.
Comparison with the associate model.-For comparison, all measured σ · T vs. T −1 curves were also analyzed by the standard defect association model. The most frequently discussed defect complex (cluster) is that of a singly charged associate between Y and an oxygen vacancy, which forms according to with mass action law K = n Y−V n Y,free · n V,free . [16] Symbols n Y,free , n V,free , n Y−V , and n Y,tot denote concentrations of free Y dopants, free vacancies, associates and the total Y concentration, respectively. Mass balance and charge neutrality requires n Y,tot = n Y,free + n Y−V and n Y,free = 2n V,free + n Y−V , respectively. Hence, results from Eq. 16 with the solution n V,free = − 1 2 n Y,tot 2 [18] (This treatment is in accordance with Ref. 41 except that the last term in Eq. 18 has a positive sign here.) In accordance with Eq. 1, Nernst-Einstein's equation u = ze 0 D/kT and D = D 0 e −Ea/kT [19] the ionic conductivity can thus be expressed by with mass action constant K = K 0 · e − Has / kT [21] ( H as = standard association enthalpy, E a = migration barrier height, D 0 = high temperature limit of the diffusion coefficient D). For given doping concentration, four parameters (E a , H as , K 0 , D 0 ) can be determined by fitting Eq. 20 to the measured σT curves. The fit quality is again excellent (see Fig. 7b) and from this point of view both models are equally applicable. The resulting fit parameters are given in Table IV (Table IV). In particular, the observed variations of D 0 and E a contradict the basic model assumption of isolated free vacancies.
In this analysis the total activation energy for low temperatures is simply the sum of E a and | H as |. This sum only slightly increases with the dopant concentration and not surprisingly coincides excellently with E a,2 of the two barrier model. Also the experiments with different cooling rates can be quantified by the associate model and results are shown in Table IV. The supposed vacancy migration energy E a drops from 0.77 eV to 0.60 eV for fast cooling; also D 0 decreases by about a factor of five. This is again in contradiction with the assumption of free vacancies in YSZ.
Migration energies and association enthalpies reported in literature partly differ significantly from our values, despite considering the same association model. Association enthalpies in literature are often smaller, for example 0.26 eV (determined for 9.5 mol% yttria measured up to 800 • C), 41 0.36 eV (12 mol% yttria, up to 1000 • C), 71 0.33 eV 7 or 0.15 eV 72 for 8 mol% yttria up to 1000 • C, 0.13-0.19 eV (8.4−12 mol% yttria up to 500 • C), 73 0.32 eV for 9.4 mol% yttria up to 600 • C. 74 Reported migration energies, on the other hand, are often higher than our fit values, e.g. 0.84 eV for 9 or 9.4 or 9.5 mol% yttria 52,74,41 0.92 eV 72 or 0.79 eV 75 or 0.74 eV 7 for 8 mol% yttria, 0.93 eV for 8.7 mol% yttria, 76 0.91 eV for 12 mol% yttria. 71 This discrepancy to our values may be simply caused by the way how activation energies are usually determined. Mostly, the conductivity data are not parameterized in the entire temperature range by a single equation (e.g. Eq. 20) but are simply split into two parts: High  temperature data are approximated by a linear line, the slope of which is interpreted in terms of E a , and the same is done for low temperature data, leading to | H as | + E a . This approach is exemplified in Fig. 10 for one YSZ single crystal (MPS). First, this analysis suffers from the ambiguity with respect to the temperature range chosen for extrapolation. In Table V we summarize effective activation energies of sample MPS deduced from linear fits to the σT values of restricted temperature ranges. The low temperature activation energy is not very sensitive and values only change from 1.14 to 1.12 eV when increasing the upper limit of the temperature window from 441 to 570 • C while keeping the low temperature limit at 204 • C. (The exact value from the fit to Eq. 20 is 1.15 eV.) Deviations are somewhat larger when neglecting the low T conductivity values, cf. data in Ref. 77 For high temperatures, however, the activation energy strongly depends on the chosen temperature range. Values between 0.69 eV and 0.89 eV were deduced from such an analysis (Table V). However, the correct value from the true fitting analysis using Eq. 20 is still lower than all approximate values, namely 0.54 eV. This shows that the transition between the two activation energy regimes is very broad and that only temperatures far above 1000 • C reflect the true high temperature activation energy without any influence of the bending of the curve. Such high temperature data are often not available and not surprisingly the reported supposed activation energies E a are often too large. Accordingly, the supposed association enthalpies (difference between low T activation energy and E a ) are too small. Indeed an accurate analysis in a broad temperature range lead to values that are very similar to that in Table IV, e.g. | H as | = 0.57 eV and E a = 0.73 eV for 12 mol % yttria. 40 In Ref. 40 it is also shown that the transition  from the low temperature to the high temperature regime ranges from about 800 K to 1600 K.

Interpretation of the Two Barrier Model
It was shown above that the two barrier model is able to quantitatively describe the measured temperature dependences of the ionic conductivity. We also introduced the model of an average vacancy path with different barrier heights that includes the two barrier model as a limit. However, so far these considerations do not give any specific atomistic interpretation of the different energy barriers. Dynamic numerical simulations with all types of long range defect interactions and a detailed analysis of sequences of successful vacancy jumps could validate or disprove our hypothesis of two types of important barriers. Such simulations are far beyond the scope of our paper but some first remarks on possible meanings of the barriers are given in the following. Two extreme cases leading to spatially varying barrier heights are sketched in Fig. 11: Higher barriers can either be caused by energetically unfavorable transition states (a) or by low levels of ground states (b).
We first consider the case that a certain fraction of vacancy jumps takes place across energetically unfavorable transition states. Strongly different energies of transition states were already calculated for jumps across Y-Zr and Zr-Zr tetrahedral edges, e.g. an activation energy difference of 0.5-0.7 eV was found between the energetically lower Zr-Zr tetrahedral edge and a Y-Zr edge. 13,14,26 This energy further varies when allowing presence of additional Y ions in the two relevant tetrahedrons. 26 Hence, different (parallel) paths are possible for vacancies moving in the complex energy landscape of YSZ: paths with some jumps across high barriers and paths bypassing all high barriers, possibly at the expense of more jumps. This situation is similar to that of two parallel reaction paths in electrochemical electrode reactions and in Ref. 78 it is shown that such parallel paths always lead to a temperature dependence with higher activation energy at F801 higher temperatures. Indeed kinetic Monte Carlo simulations with two different barrier heights and high yttria concentrations revealed a slight increase of the activation energy for increasing temperature. 14,26 However, this is in contradiction to our experimental observations and this parallel path interpretation of our two barrier model is therefore not reasonable.
Only if it is unavoidable to pass both barrier heights in series, the higher activation energy is found at lower temperatures, cf. Ref. 78. Above the upper percolation limit of the low barriers such a situation certainly results and oxygen vacancies can then never avoid all high barriers. However, this upper percolation limit of paths only consisting of energetically favorable Zr-Zr edges is by far not reached for the dopant concentrations used in our study: 10 For a random distribution of Y about 22 mol % Y 2 O 3 would be required to block essentially all paths consisting of low Zr-Zr edges only. Therefore, a model explaining the high activation energy simply by unavoidable jumps across Y-Zr edges is also not possible.
Accordingly, additional effects have to come into play, which might enforce jumps across high transition states much below the nominal percolation limit of (low) Zr-Zr edges. Oxygen vacancies repel each other 15,28,29,34 and all jumps leading to a decreasing distance to the next vacancies may thus have higher activation energies, irrespective of the cation edge. Moreover, in Ref. 36 it was shown by ab-initio metadynamic simulations that vacancy-vacancy interaction causes surprising phenomena such as unavailable sites for oxygen vacancy jumps or collective oxygen vacancy jumps. The exact characteristics of unavailable oxygen sites or of defect configurations leading to collective jumps remained unknown since both are an outcome of complicated multi-defect interactions. Cooperative motion of oxygen vacancies was also reported in Ref. 11. Moreover, it was shown that long range defect interaction may modify the local barrier height by about 0.4 eV even for identical nearest neighbors. 29 Hence, multi-defect interactions may easily increase barrier heights by about 0.5 eV, i.e. by the value suggested in our study.
These complex effects of multi-defect interaction on vacancy migration depend on the location of all vacancies and are thus highly dynamic in nature. Temporarily, situations may arise in which a vacancy faces high barriers in all possible jump directions even though it is not surrounded by high energy Zr-Y (or Y-Y) edges only. We may call this a "dynamic high barrier". A sophisticated interplay of defects might therefore cause additional high barriers that cannot be avoided by vacancies while moving in YSZ. Accordingly, already for the dopant concentrations considered here we may have a situation with all possible paths of vacancies necessarily including at least some high barriers. Even if only one percent or less of all vacancy jumps cannot avoid such "dynamic high barriers", this can lead to a high activation energy at low temperatures, cf. γ 2 /γ 1 in Fig. 9c with ν 0 1 ≈ ν 0 2 . The second extreme case with serial high activation energies corresponds to a situation with vacancies becoming trapped in deep ground states (Fig. 11b). Formally, this is similar to the standard association model presented above. However, identifying this trapped vacancy with a (Y x complex in a static, dilute solution of trap sites in YSZ is not appropriate due to the very high Y concentration. Moreover, we also have vacancy-vacancy interaction and thus the energy landscape is highly dynamic. Still we may face situations in which locally a kind of "dynamic trap" with low energy emerges for specific local defect configurations and only disappears after the motion of one or several vacancies. Such "dynamic traps" may resemble the large defect clusters considered in several modeling studies [30][31][32] and most probably involve numerous Y / Zr and V •• O defects. Also an energetically advantageous alignment of oxygen vacancies along the <201> or <111> direction 21,79,35 or curved chains 31 was reported. Our two barrier model could thus also approximate a situation in which from time to time vacancies find themselves in "dynamic traps" and can only leave those traps by overcoming a high barrier. In contrast to the defect associate of Eq. 15 such "dynamic traps" only emerge temporarily for certain local defect configurations. Formally, this can again be treated by a thermodynamic approach, despite the different meaning of the traps compared to Eq. 15. We may write this trapping as a reaction according to [22] with mass action law K = K 0 · e − Htrap / kT = n V,trapped n V,untrapped . [23] Symbols n V,trapped and n V,untrapped denote the concentrations of trapped and untrapped vacancies, respectively and H trap is the enthalpy difference between an "average vacancy" and a vacancy in a dynamic trap. The average vacancies are those exposed to a kind of averaged defect interaction (outside the traps) and should not be confused with free vacancies in YSZ; the latter are by definition not exposed to any defect interaction. The "dynamic traps" do not enter the mass action law, since their concentration is neither fixed nor does it limit the number of trapped vacancies. It is shown in the appendix that also this situation leads to the conductivity equation of the two barrier model (Eq. 13).
The assumption of such "dynamic traps" would also solve the conflict between a pre-treatment dependent association enthalpy (cf . Table IV) and the conventional associate model. These are the two extreme cases that may lead to a two barrier model but possibly both a lowering of the ground state and an increased energy of the transition state play together and cause the unavoidable serial high barriers. Moreover, a "bimodal" distribution of barrier heights may exist, rather than two very specific and welldefined energy barriers, cf. Fig. 11c. Hence, our two specific activation energies may reflect averaged values of two groups of barrier heights and do not exclude existence of other barrier heights. In any case, within the framework of our model the majority of the successful jumps has to take place across the lower barriers with an activation energy of approximately E a,1 . These majority jumps take place between states of strong but still similar defect interactions, i.e. between similar initial and final states and across similar barriers. Such similar energy levels may be caused by the large number of defects in several neighboring shells that affect local energies. Thus, effects of the exact locations of all defects may be largely levelled out and a kind of average jump barrier results. Only from time to time local situations are such that a much higher barrier has to be passed.
However, we have to emphasize that the excellent agreement between all measurement data and the two barrier model is far from being a proof of its validity. Further analysis of the temperature dependence of the YSZ conductivity by ab-initio calculations may reveal the true distribution of barrier heights in an average oxygen vacancy path. Provided such simulations reproduce the activation energy change in the log(σT) vs. T −1 plot, one could test whether the two E a,i values of a two barrier analysis find their counterpart in atomistically meaningful barrier heights.

Conclusions
Owing to the very high defect concentration in yttria stabilized zirconia (YSZ), essentially all oxygen vacancies are exposed to pronounced interaction with dopants and other oxygen vacancies. Therefore simple defect associate models with a temperature dependent fraction of (interaction-)free, mobile vacancies are considered as inappropriate to interpret the temperature dependent activation energy of the ionic conductivity. An alternative model is suggested, based on the assumption that an average path of an oxygen vacancy in YSZ includes (at least) two types of barriers with significantly different barrier heights. The corresponding analytical equation excellently describes the temperature dependent conductivities measured for YSZ of different dopant content.
The lower barriers determine the activation energy at high temperatures and we suggest that this barrier type reflects the average defect interaction oxygen vacancies in YSZ are exposed to.

F802
Journal of The Electrochemical Society, 164 (7) F790-F803 (2017) Occasionally, however, defect configurations may arise which lead to high barriers, either due to strong attractive energies in local and temporarily existing defect clusters or due to defect configurations with pronounced repulsive oxygen vacancy interaction. Even if only a few jumps within an average path of an oxygen vacancy in YSZ have to take place across such high barriers, a much higher activation energy of the conductivity results at low temperatures. However, additional dynamic ab-initio calculations are required to test this hypothesis of two barrier types in YSZ, and those simulations have to include long range vacancy-vacancy interaction.
Supposed this model of two serial barriers is meaningful, we can conclude from our data analysis that low barriers are in the 0.6 eV range for 8-10 mol% yttria and high barriers are of the order of 1.0-1.1 eV. Both barrier heights slightly increase with increasing doping content. Conductivities and thus also barrier heights in YSZ also depend on the thermal prehistory: The exact conductivity curve can be reproducibly changed by modifying the cooling rate after high temperature annealing at 1550 • C. This is most probably due to different distributions of Y dopants in YSZ, resulting from different temperatures at which the cation sub-lattice becomes frozen-in.
The conductivity measurements showed the expected dopant dependence of the YSZ conductivity, with the lowest conductivities found for the highest doping concentration (11.5 mol%). However, while at low temperatures the conductivity decreased by up to a factor of six when changing from 8 to 11.5 mol %, the effect was much smaller at 1000 • C. Finally, it was shown that a simple analysis of conductivity curves in terms of two specific temperature regimes with different activation energies may lead to severe misinterpretation. The transition range from one to the other activation energy is very broad and thus particularly the high temperature activation energy can be easily overestimated in such a simplified analysis.