Ab initio kinetic Monte Carlo model of ionic conduction in bulk yttria-stabilized zirconia

The kMC method is a variant of the Monte Carlo method and can simulate kinetic processes because time is also updated during the simulation [32]. YSZ maintains the cubic phase at intermediate temperature when the doping concentration is over 6 mol% [59]. The unitcell of cubic YSZ is shown in figure 1. Each unitcell contains four sites to be occupied by cations (Zr⁴⁺ or Y³⁺), and eight sites to be occupied by oxygen anions or vacancies. The cations occupy the cubic corners and the face centers, and the oxygen anions and vacancies occupy the center of the eight sub-cubes inside the unit cell. Therefore, the anion sites form a simple cubic sublattice. The simulation supercell is subjected to periodic boundary conditions (PBCs) in all three directions and consists of N_x × N_y × N_z unitcells. For simplicity, all lengths are expressed in units of 1/4 of the lattice constant a₀. In this way, the coordinates of each cation and oxygen site can be specified by three integers. The coordinates of cations are even numbers and the coordinates for oxygen sites are odd numbers. The lattice constant is chosen to be 5.14 Å, which corresponds to the experimental value for 8 mol% YSZ [33].

Zoom In Zoom Out Reset image size

Since the vacancy mechanism determines the ionic diffusion in YSZ, we simulate the oxygen vacancy diffusion instead of oxygen ion diffusion. Each oxygen vacancy can possibly jump to one of their six adjacent sites (±2, ±2, ±2) if the adjacent site is occupied by an oxygen anion. The jump unitcell is also shown in figure 1. For every vacancy, the probability rate j = ν exp(−E_b/k_BT) is calculated for every jump direction. Each possible jump is called an event. E_b is the energy barrier of event and is calculated by the energy difference between the current state and the transition state during the event. ν is the jump trial frequency¹. Vacancy diffusion is simulated by selecting one event from all possible events and updating the vacancy positions and the physical time accordingly at each interaction [32, 34].

Since selecting an event at each kMC step depends on the probability rate of each event, accurate calculation of the energy barrier of each event is essential for an accurate kMC simulation. Even though DFT calculation is usually accepted as a much more accurate method compared with empirical potential models, it is very challenging to calculate the energy barrier from DFT because of its heavy computational cost. Thus we need to develop an energy barrier model that is efficient but maintains an accuracy similar to DFT calculations.

Krishnamurthy et al [21] developed an energy barrier model, in which the energy barrier is determined by the chemical species of the two cations on the transition plane that the diffusing vacancy passes through (ions 3 and 4 in the jump unit cell in figure 1) and the energy barrier is increased as more Y ions occupy these two sites. This model predicts the optimal doping concentration (8–10 mol% at high temperature, T > 1000 K) consistent with experimental observations. We will call this model the non-interacting model because the energy barriers for each vacancy jump and its reverse jump are the same, meaning that the system energy remains the same before and after each jump. However, this is not true in ionic solids which have strong interactions between ions. Due to the ionic interactions, the energy is different before and after vacancy jump, unless the ionic configurations before and after vacancy jump are equivalent by symmetry. The ionic interaction makes the energy barrier of vacancy jump different in the forward and backward directions. Hence, to build a more physical energy barrier model, the ionic interactions must be considered.

The KRA energy barrier model has been proposed for this purpose [27, 28]. In the KRA model, energy barrier is divided into a kinetic part and a thermodynamic part. The thermodynamic part is defined as the energy difference of system before and after the vacancy jump, ΔE, and the kinetic part describes the energy barrier when ΔE is zero, $E_{\rm b}^0$. The kinetic part and the thermodynamic part are assumed to be independent each other in the KRA model. Thus the energy barrier E_b can be written as a function of $E_{\rm b}^0$ and ΔE.

A simple linear function has been proposed to correct the energy barrier by a half of the energy difference after vacancy jump ΔE [27, 28], i.e. $E_{\rm b} = E_{\rm b}^0+{\Delta E}/2$. This function means that the saddle point remains the same (along the reaction coordinate) regardless of ΔE. This is a good approximation only when ΔE is much smaller than $E_{\rm b}^0$. Since in our DFT calculations for YSZ, ΔE is on the same order of magnitude as $E_{\rm b}^0$, we designed a model applicable in the extended range of ΔE.

The energy profile usually can be approximated by the summation of several sinusoidal functions. When ΔE is zero, the energy profile during vacancy jump can be approximated by $E(x)=E_{\rm b}^0\cos^2(\pi x)$, where x is the reaction coordinate; −0.5 indicates the state before vacancy jump and +0.5 indicates the state after vacancy jump. The energy barrier for the forward direction is the energy difference between x = 0 (saddle) and x = −0.5 (local minimum). The energy barrier for the backward direction is the energy difference between x = 0.5 and x = 0. If the energy changes after vacancy jump, we can assume that the energy profile is linearly shifted, i.e. $E(x)=E_{\rm b}^0\cos^2(\pi x)+\Delta E x$.

This changes the positions of the saddle point and local minimum as illustrated in figure 2(a). The energy barrier for the forward direction is obtained by finding a new local minimum and a new saddle point in the energy profile. If ΔE is too large, the energy barrier predicted by this approach is less than zero. This means that there is no energy minimum before the vacancy jump. To avoid the negative energy barrier in the kMC simulations, the energy barriers less than 0.001 eV are all set to 0.001 eV. Once the energy barrier is obtained for the forward direction $(E_{\rm b}^{\rm FW})$, the energy barrier for the backward direction $(E_{\rm b}^{\rm BW})$ is set to $E_{\rm b}^{\rm FW}-\Delta E$ to satisfy detailed balance. Using this numerical method, the energy barrier as a function of ΔE can be obtained as shown in figure 2(b).

Zoom In Zoom Out Reset image size

Now that we have obtained function $f(E_{\rm b}^0,\Delta E)$, the task of the energy barrier calculation has been converted to the task of computing $E_{\rm b}^0$ and ΔE between metastable states. For $E_{\rm b}^0$, the values in the non-interacting model [21] are used in this paper. However, ΔE still needs to be calculated at every kMC step.

DFT calculations can be used to compute the energy of different ionic configurations by solving for the electronic ground state of the system. DFT calculations are considered to be highly accurate and widely applicable, but are also computationally very expensive. For a system containing N electrons, the computational time of DFT scales as O(N³).

Even though the recent progress in the combination of high-performance computing and DFT methods enables us to calculate the energy for larger systems (up to a few hundred ions), the computation is still very expensive. For example, one DFT-based structural relaxation calculation for about 300 ions in this paper takes about 3 days on 160 CPUs. Thus, it is impossible to do one such DFT calculation at every kMC step and we need an efficient method to predict the energy of arbitrary ionic configurations based on a limited set of DFT data. The CEM enables us to reduce the computational load dramatically by avoiding the necessity of performing one DFT calculation for every ionic configuration.

In the CEM, the configurational properties of system are expanded with the orthonormal basis functions which are defined by the ionic configurations on the lattice [35]. Once we fit the expansion coefficients to a limited number of DFT calculations, these coefficients can be used to predict the energy of new ionic configurations. The CEM provides the basis functions for the expansion.

Consider an ionic configuration of YSZ which contains N lattice sites. In YSZ, four kinds species can occupy the lattice sites—they are host cation (Zr), dopant cation (Y), anion (O) and anion vacancy (VA) in YSZ. For simplicity we can assume that cations and anions exist only on their own sublattices: Zr and Y exist on the cation sublattice while O and VA exist on the anion sublattice. The ion type in each sublattice is labeled by a spin σ_i: +1 for Zr, −1 for Y, +1 for O and −1 for VA in YSZ. Then the ionic configuration of a YSZ single crystal can be represented by a vector $\vec{\sigma}=(\sigma_1,\sigma_2,\ldots,\sigma_N)$, where N is the number of lattice sites.

The cluster functions are defined as² φ_α = ∏_i∈ασ_i, where α indicates a cluster composed of multiple lattice sites including empty cluster, N point clusters,_NC₂ pair clusters,_NC₃ three-body clusters, ..., and_NC_NN-body cluster where $_NC_m=\frac{N!}{N!(N-m)!}$. These cluster functions satisfy the orthonormality condition, $\langle \phi_{\alpha}\phi_{\beta}\rangle=\frac{1}{2^N}\sum_{\vec{\sigma}}\phi_{\alpha}\phi_{\beta}= \delta_{\alpha\beta}$ where the sum is over all possible choices of $\vec{\sigma}$. Because the cluster functions satisfy the condition of an orthonormal basis, the energy of a system can be expressed as a unique expansion of the cluster functions φ_α.

$$\begin{equation} E(\vec{\sigma}) = \sum_{\alpha}V_{\alpha}\phi_{\alpha}(\vec{\sigma}), \end{equation} \tag{ 1 }$$

where the coefficient V_α associated with φ_α is called the effective cluster interactions (ECIs) of cluster α, and the sum is over all 2^N cluster functions. We can obtain φ and E from pre-computed DFT datasets. ECI then can be calculated from the least-square solution of equation (1). Once ECIs are obtained, the energy for any given $\vec{\sigma}$ can be predicted by equation (1).

To prepare pre-computed DFT datasets, we performed DFT calculations for 140 different ionic configurations at the same doping concentration at 9 mol%, using the Vienna Ab-initio Simulation Package (VASP) [37–40]. GGA-PW91 scheme [41, 42] is used for exchange-correlation function and pseudopotentials with projector augmented wavefunctions (PAWs) [43, 44] are used. All ions are allowed to move freely during the conjugate gradient relaxation.

The size of supercell used in this work is 3[1 0 0] × 3[0 1 0] × 3[0 0 1], which contains 108 cations (Cat) and 216 anions (An). This size is much larger than the supercell in previous DFT calculations of YSZ [21, 23, 45], and it can cover the 9th nearest neighbor (9NN) for Cat-An, 9NN Cat-Cat and 18NN Cat-An. Because the size of supercell is large, gamma point is enough for the k-point sampling. The charge smearing factor SIGMA = 0.05 is used. Each configuration usually takes 3 days for full relaxation (volume and shape relaxation) on 160 CPUs. This high computational cost in DFT calculations is the main reason that we want to develop a cluster expansion model to predict the energies of arbitrary ionic configurations. The mean value of the energy in these 140 configurations is −2974.15 eV and the standard deviation is 1.06 eV.

To solve equation (1) by the least-square method, the number of dataset must be greater than or equal to the number of coefficient. There are 2¹⁰⁸ × 2²¹⁶ = 2³²⁴ ≃ 3.417 × 10⁹⁷ clusters in the 3[1 0 0] × 3[0 1 0] × 3[0 0 1] YSZ supercells. Since the number of datasets, 140, is much smaller than the number of cluster functions, we need to truncate the clusters. The clusters in the supercell can be classified by their translational and rotational symmetries. PBCs limit the distance between two ions in a cluster within the half length of the supercell in all directions. Thus, a cut-off distance of 1.5a₀ is applied to all three dimensions. We also limit the number of ions in the cluster to be less than 4.

There are two types of point clusters, one is the cation point cluster and the other is the anion point cluster. For the anion–anion pair clusters, the independent clusters (not related by symmetry) are specified not only by the distance between anions but also by the number of cations between the anions pair. The three-body clusters can be represented as a triangle, consisting of three line segments. Because each segment can be associated with a pair cluster, three-body clusters can be labeled by the set of three pair clusters. All pair clusters and one example three-body cluster are shown in the auxiliary material.

As a result of the above truncation schemes, the total number of independent ECI is 192. For robust fitting, the number of pre-computed datasets needs to be several times larger than the number of clusters [46]. Therefore we need to truncate the clusters further.

Among the 192 clusters, we select a small subset of clusters by minimizing the cross validation score (cvs) [47]. As the quality of the selected clusters depends on the size of subset, k, we need to decide the optimal value for k itself. When the number of coefficient is too large and approaches the number of pre-computed datasets, the clusters' ECI can be over-fitted meaning that the fitting error is very small but the fitted model will make huge errors in making new predictions. (See the online supplementary material for the demonstration of the dependence of the predictive power on k (stacks.iop.org/MSMSE/20/065006/mmedia).) Thus it is necessary to check whether the fitted coefficients have any predictive power. We can measure the predictive power from the error in the prediction of the energy of new ionic configurations which are not included in the fitting. Among the 140 datasets, 100 datasets are used for fitting of ECI and the remaining 40 datasets are used for checking the predictive power of the fitted ECI.

A reasonable approach to assess whether or not over-fitting occurs is to compare the standard deviation in the error (E_DFT − E_CEM) for the two datasets used in fitting and in testing real predictions, respectively. For a fixed k, Monte Carlo algorithm is used to search a set of k clusters providing the lowest cvs. This process is repeated by changing k systemically. Standard deviations of the error in fitting and in real prediction for different k are obtained as shown in figure 3. From this result, we conclude that k = 9 is optimal, because the summation of two standard deviations is minimized when k = 9. It also implies that in this case the optimized ratio between the size of cluster set and the number of datasets is close to 0.1(= 9/100). When the selected nine clusters are used to predict the energy of 40 datasets which are not included in the fitting, the result is given in figure 4. The root-mean-square (rms) error in real prediction is 0.0049 eV per cation, which is consistent or even smaller compared with previous works [45].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One of the nine clusters is a cation point cluster, representing the cohesive energy of the YSZ crystal. The ionic positions of remaining eight clusters are given in table 1. It is interesting that all the selected clusters (except the point cluster) are three-body clusters. In empirical potential models, the ionic interaction has been usually considered to be dominated by pair interactions and three-body terms are often considered as higher order corrections. The previous works using the CEM also reported that pair clusters are selected more frequently than three-body clusters and have stronger ECIs [45].

The difference from previous works may originate from the fact that the number of selected clusters k in previous works is substantially larger than what we use (k = 9) in this paper. The other possible reason is the purpose of the fitted CEM model. Previous works aim to the energy and structure of the ground state, while the aim in this paper is to predict the energies of all the metastable states encountered in kMC simulations of vacancy diffusion. Pair clusters may be more important for predicting ground-state energies and three-body clusters may be more important for metastable states.

The ECIs obtained for the selected clusters, except the cation point cluster, are given in table 1. The ECI corresponding to the cation point cluster is not relevant for kMC simulations of ionic conduction because it contributes the same constant to all metastable states. Similarly, the ECI corresponding to CCC cluster also cancels out in the energy difference, because cations do not move in our kMC simulations. ECIs of the selected clusters are not strongly correlated with the size of the clusters.

Using the fitted ECIs, the binding energy between the Y ion and oxygen vacancy can be obtained. The energy of the system with one Y ion and one vacancy (E_1V1Y) is calculated by changing the position of vacancy with respect to the Y ion (d_YV). This energy is then compared with the energies of a crystal containing one Y ion (E_1Y), of a crystal containing one oxygen vacancy (E_1V) and of a pure zirconia crystal (E₀₀).

$$\begin{equation} E_{\rm bind}(x) = E_{\rm 1V1Y}(x) - E_{\rm 1Y} - E_{\rm 1V} + E_{\rm 00}. \end{equation} \tag{ 2 }$$

The result is plotted in figure 5. We can see that the vacancy prefers to be located at the second nearest neighbor site of the Y ion. This behavior of binding energy is consistent with the previous works where the binding energy was calculated by fitting to several DFT calculations [22, 23, 48] and is also consistent with experiments [49–51].

Zoom In Zoom Out Reset image size

We recall that E_b is predicted as a function of the calculated ΔE and $E_{\rm b}^0$, using the function illustrated in figure 2(b). In this paper, $E_{\rm b}^0$ is taken as the same values in the non-interacting model [21], where $E_{\rm b}^0$ is assumed to depend on the two cations on the transition plane: $E_{\rm b}^0$ is 0.58 eV for Zr–Zr pair, 1.29 eV for Y–Zr pair and 1.86 eV for Y–Y pair, respectively. We examined three sampled cases to compare the E_b prediction by our model with that by a direct DFT calculation. For a direct DFT calculation, the constrained minimization method and the 3[1 0 0] × 3[0 1 0] × 3[0 0 1] supercell were used. The result, given in table 2, shows that the discrepancy between our model and a direct DFT calculation is within ∼0.12 eV.

Table 2. Comparison between the energy barrier prediction by a direct DFT calculation $(E_{\rm b}^{\rm DFT})$ method and our model (E_b), in three sampled cases.

	ΔE	$E_{\rm b}^0$	$E_{\rm b}^{\rm DFT}$	Eb
1	0.32	0.58	0.65	0.67
2	−0.52	1.29	1.10	0.98
3	−0.44	1.86	1.40	1.45

In the first hypothetical distribution, a fraction x of Y ions segregate forming a spherical cluster and the remaining Y ions are randomly distributed in the supercell. We calculate the impedance by changing the fraction of segregated Y ions. The initial positions of vacancies are randomly distributed and data collection starts after a transient period of 6 × 10⁵ kMC jumps to reach equilibrium. 80 million kMC jumps at 1800 K are performed for each distribution. The results are shown in figure 10(a). We observe that segregation of Y ions into a sphere only increases the electrical impedance (i.e. decreases the ionic conductivity) of YSZ. To gain insight on the microscopic mechanisms, we examine the dependence of various parameters of the equivalent circuit on dopant concentration.

Zoom In Zoom Out Reset image size

When the Nyquist plots are fitted to a two-arc model, R₃ is always much smaller than R₂ (by two orders of magnitude). This supports that the Nyquist plot has a single arc shape. In the following, we still use the equivalent circuit with two arcs to fit the impedance spectrum. R₁ corresponds to the left end of the arc and R₂ + R₃ corresponds to the width of the arc. The equivalent circuit parameters as functions of the fraction of segregated Y ions are plotted in figure 10(b). We can see that both R₁ and R₂ + R₃ increase when more Y ions are segregated.

The hypothetical Y-segregation was intended to enhance the ionic conductivity by providing Y-free regions, which has low energy barrier for the diffusing vacancy. KMC simulations based on the non-interacting model have predicted that such a structure would enhance ionic conductivity [52]. However, now the interacting model predicts that Y-segregation forming a spherical cluster is detrimental to ionic conductivity. The origin of this discrepancy between the predictions by the interacting model and by the non-interacting model can be explained in terms of the vacancy distribution. Recall that the vacancy concentration and the effective energy barrier together determine the conductivity [59]. In the case of the non-interacting model, all metastable states have identical energy, which leads to a uniform vacancy distribution despite the existence of the interaction between the vacancy and cations. However, in the interacting model, the vacancy distribution will no longer be uniform if the Y ions are not randomly distributed.

We obtained the vacancy distribution using the simulated annealing method and observed that the energetically favorable positions of the vacancies are at the first nearest neighbor of outer Y ions of the spherical Y cluster. To diffuse over long distances and contribute to the ionic conductivity, vacancies must detach from the Y clusters. This requires overcoming a binding energy of about 0.12 eV (figure 5), which reduces the ionic conductivity. Therefore the reduction in ionic conductivity is caused by the increase in spatial variation of the interaction energy. This finding explains why the non-interacting model is quite successful when the Y ion distribution is completely random, but experiences problems as soon as this is not the case.

Let us return to the change in the equivalent circuit parameters with the fraction of segregated Y ions in spherical distribution. As more Y ions are segregated, vacancy concentration becomes higher near the Y-sphere, which leads to more vacancy–vacancy and vacancy–yttrium interactions. We propose that the increased ionic interactions is the cause of increase in R₂ + R₃ with increasing Y segregation (figure 10(b)), similar to the increase in R₂ with increasing doping concentration (figure 6(b)). We propose that $R_1^{-1}$ is proportional to the vacancy concentration in Y-free region (i.e. mobile vacancy concentration).

Following the reasoning given above, we might expect a different result when Y ions are segregated to a layer, because intuitively it provides the diffusing channels without requiring oxygen vacancies to detach from isolated Y clusters. In these channels, vacancies jumping along a direction parallel to the Y layer experience a reduced spatial variation of the energy. We tested configurations where Y ions are segregated to a layer parallel to the x-axis (conduction direction) in 3[1 0 0] × 3[0 1 0] × 3[0 0 1] supercell. To minimize the reduction in spatial variation of the energy, all cation sites in the segregation layer are completely occupied by Y ions, which corresponds to 9 mol% of doping concentration. 800 million kMC steps are performed at T = 1800 K. The conductivity is predicted to be 0.21 S cm⁻¹, while it is 0.51 S cm⁻¹ and 0.37 S cm⁻¹ for the random distribution and the spherical cluster, respectively. Unfortunately, the interacting model predicts that the ionic conductivity when Y ions segregate to a layer is still lower than when the Y ions are randomly distributed.

Using the simulated annealing method, the optimal vacancy distribution is obtained when Y ions are segregated to a layer. We observe that vacancies are located in the first anion layer next to the Y-segregation layer. The vacancy distribution at finite temperatures will not be very different from the minimum-energy vacancy distribution. Hence, when a vacancy jumps in the x direction, it mostly experiences a Y–Zr pair $(E_{\rm b}^0=1.29\,{\rm eV})$ at its saddle configuration, which results in high E_b. This region can be called slow diffusing channels while the region where a jumping vacancy experiences Zr–Zr pair $(E_{\rm b}^0=0.58\,{\rm eV})$ can be called fast diffusing channels.

When all the cation sites in one layer are occupied by Y, we have seen that excessive Y-segregation in this layer attracted vacancies to the first anion layer and this lowers ionic conductivity. The number of vacancies in the first anion layer could be reduced by reducing the number of Y ions in the segregation layer. Another factor to determine the conductivity of layer structure is the distance between adjacent Y layers. As more Y-segregation layers are accommodated in the supercell, the vacancy concentration increases and ionic conductivity may be enhanced. However, too short inter-layer distance may cause more vacancies to get too close to the Y layers and to reduce ionic conductivity.

Various layered structures are examined to verify the hypothesis that two factors determine the ionic conductivity in a layer structure: the number of Y ions in the segregation layer and the period of the segregation layer. Layer structures in a 3[1 0 0] × 3[0 1 0] × (p/a₀)[0 0 1] supercell are prepared by changing the number of Y ions in the segregation layer and by changing the inter-layer distance p. In all cases Y ions are randomly distributed in the segregation layer. Since Y ions are located in the segregation layer only, doping concentration changes according to the number of Y ions in the segregation layer. 800 million kMC jumps are performed at 1800 K and the results are shown in figure 11.

Zoom In Zoom Out Reset image size

We observe that the ionic conductivity is maximized at 9 mol%, when p = 1.5a₀. At this doping concentration, half of cation sites in the segregation layer are occupied by Y ions. The effective energy barrier of the entire structure can be obtained from the temperature dependence of ionic conductivity. The optimal distribution, the Y-segregation layer structure at 9 mol% with p = 1.5a₀, is examined. The effective energy barrier of this structure is obtained as 0.74 eV for T > 900 K and 0.85 eV for T < 900 K, which is very close to the effective energy barrier in random Y distribution case. This similar effective energy barrier results in very small enhancement of the ionic conductivity. The conductivity is enhanced by only 20% at T = 1800 K and 17% at T = 500 K, compared with random Y distribution at 8 mol% of doping concentration.

The reason for this low enhancement can be explained by the destruction of fast diffusing channels. We removed some Y ions in the segregation layer to reduce vacancy attraction to the first anion layer, but random removal of Y ions increases the spatial variation of the potential energy in the ionic conduction direction. The energy fluctuation in the interaction energy increases the energy barrier of the jumping vacancies.

If Y ions are distributed in a line in the x-direction, a vacancy jumping in the x-direction will not experience any energy variation as well as the number of Y ions in the segregation layer is reduced. All cation sites in the segregation line need to be completely occupied by the Y ion for the energy variation to vanish. The energy of a hypothetical system, which has one Y-segregation line in the x-direction and one oxygen vacancy, is calculated by changing the vacancy position in (y, z) space. The vacancy concentration in (y, z) space can be approximated by the Boltzmann factor using this energy. The result at T = 500 K is plotted in figure 12(a). We can clearly see that the vacancy concentration is high in fast diffusing channels, where a vacancy jumping in the x-direction experiences a Zr–Zr pair only and ΔE by vacancy–cation interactions is zero. As more Y-segregation lines are accommodated in the supercell, more vacancies are created and the ionic conductivity should increase.

Zoom In Zoom Out Reset image size

However, some minimum spacing between Y-segregation lines is required. Otherwise, we get to the same planar segregation in the limit of inter-rod spacing going to zero. In that case, as shown in figure 12(b), the vacancy concentration is the highest at the first anion layer, in which vacancy experiences Zr–Y pairs (high $E_{\rm b}^0)$ when jumping in the x-direction. On the other hand, if the spacing is too wide, the ionic conductivity will also be low due to insufficient oxygen vacancy concentration. Thus there must exist an optimal spacing between Y-segregation lines. From figure 12(a), we can infer that the promising candidates of optimal spacing are 1.5a₀ and 2a₀.

To search for the optimal microstructure, several 2D superlattices are tested. We examined three cases of spacing between Y-segregation lines—(1.5a₀,1.5a₀), (1.5a₀,2a₀) and (2a₀,2a₀) for the y- and z-directions. After 800 million kMC steps, the conductivity at each spacing is calculated as 69.71, 58.72 and 35.27 (S cm⁻¹) at T = 1800 K and 0.0029, 0.0059 and 0.0012 (S cm⁻¹) and at T = 500 K. This means that the optimal spacing is (1.5a₀,1.5a₀) at T = 1800 K and (1.5a₀,2a₀) at T = 500 K. We will call these structures A and B, respectively. The vacancy distribution is obtained by statistically collecting the positions of all vacancies at all kMC steps. The results are shown in figure 13. We can see that in both structures the vacancy concentration is high at fast diffusing channels.

Zoom In Zoom Out Reset image size

Structures A and B correspond to 6 mol% and 4.5 mol% of doping concentration, respectively. Lower doping concentration in structure B leads to less vacancy–vacancy interactions, which helps us to decrease the effective energy barrier in structure B. The temperature dependence of conductivity for each structure is given in figure 14. Structure A has higher conductivity at T > 900 K, while structure B has higher conductivity T < 900 K. From this result, the effective energy barriers can be extracted and are listed in the table 4. We can see that the effective energy barrier of structure B is always less than that of structure A as expected. Both structures A and B show enhanced ionic conductivity compared with the structure with random Y distribution at 8 mol% of doping concentration as given in table 5.

Zoom In Zoom Out Reset image size

It may seem difficult to manufacture these two structures using today's synthesis techniques. However, The field of nano-fabrication is advancing very fast. Many bottom-up or top-down techniques are continuously added to our repertoire to fabricate more and more complex nano-structures. Moreover, the same method described here can be applied to predict the optimal structure of other solid electrolyte material such as gadolinium-doped ceria (GDC). If the optimal structure in other electrolyte materials were such that dopant ions segregate into planes, it would then be quite straightforward to synthesize by existing technologies such as molecular beam epitaxy (MBE) and atomic layer deposition (ALD) [60–62].