Distribution of Characteristic Times: A High-Resolution Spectrum Approach for Visualizing Chemical Relaxation and Resolving Kinetic Parameters of Ionic-Electronic Conducting Ceramic Oxides

Abstract

Surface exchange coefficient (k) and bulk diffusion coefficient (D) are important properties to evaluate the performance of mixed ionic-electronic conducting (MIEC) ceramic oxides for use in energy conversion devices, such as solid oxide fuel cells. The values of k and D are usually estimated by a non-linear curve fitting procedure based on electrical conductivity relaxation (ECR) measurement. However, the rate-limiting mechanism (or the availability of k and D) and the experimental imperfections (such as flush delay for gaseous composition change, τ_f) are not reflected explicitly in the time–domain ECR data, and the accuracy of k and D demands a careful sensitivity analysis of the fitting error. Here, the distribution of characteristic times (DCT) converted from time–domain ECR data is proposed to overcome the above challenges. It is demonstrated that, from the DCT spectrum, the rate-limiting mechanism and the effect of τ_f are easily recognized, and the values of k, D and τ_f can be determined conjunctly. A strong robustness of determination of k and D is verified using noise-containing ECR data. The DCT spectrum opens up a way towards visible and credible determination of kinetic parameters of MIEC ceramic oxides.

1. Introduction

The high-temperature energy conversion technologies, such as solid oxide fuel cells (SOFCs) and solid oxide electrolysis cells (SOECs), are gaining worldwide research interest. One reason is that these technologies are conceptually highly efficient and ecofriendly, however they also face challenges in delivering high performance and longevity using low-cost materials [1,2,3]. Another reason is that they are a very good platform for the fundamental research of solid state ionics and electrochemistry [4,5,6,7,8]. One representative example is the development of mixed ionic-electronic conducting (MIEC) materials for use in the electrodes. The performance of MIECs are usually evaluated by the rate of the oxygen reduction/evolution reactions on the surface, and the diffusivity of oxygen ions in bulk, represented respectively by the surface exchange coefficient (k, cm·s⁻¹) and bulk diffusion coefficient (D, cm²·s⁻¹) [9,10,11,12,13,14,15]. The measurements of k and D are usually based on ex-situ experiments, such as electrical conductivity relaxation (ECR) and isotope exchange depth profiling [12,13,15,16,17]. It has been shown that k and D determine the in-situ device-level performance, e.g., area specific resistance of electrodes [18]. In principle, the objective is to obtain high values of k and D, while considering the trade-off with longevity [2,19,20,21]. Researchers have developed atomic level descriptors that correlate to k and D for rational design of MIECs [3,22]. Therefore, the development and understanding of MIECs demands correct measurements of k and D.

Currently, the ECR method has been widely used to measure k and D, by fitting analytical solution to the relaxed conductivity data in response of a step change of atmosphere oxygen pressure. However, there is significant (orders of magnitude) discrepancy in the reported values of k and D for materials with the same nominal composition [23]. One reason is the incorporation of imperfections, such as the relaxation in initiating the driving force (e.g., the flush delay in change of oxygen pressure in atmosphere, τ_f), the nonlinearity if the sample is driven far away from equilibrium, the inhomogeneity of the sample and the instability at high temperatures [24]. Another reason for the unreliable chemical relaxation (CR) results derives from the trial-and-error approach. Reliable values of k and D can only be determined theoretically when the half thickness of the sample (L) close to a critical length (L_c = D/k), or in other words, the Biot number (Bi = Lk/D) being close to “1”. Otherwise, the relaxation is governed either by bulk diffusion (Bi >> 1), or by surface exchange (Bi << 1). Therefore, prior to measurements, the dimension of the sample should be designed properly based on the rational guess of k and D. This depends largely on experience. The Biot number (or the rate-limiting mechanism) is difficult to be identified from the ECR data, since the analytical solutions to the various rate-limiting mechanisms fit to the ECR data equally well. Efforts have been made to improve the reliability. For example, the flush delay can be shortened to a subsecond using small reactors [16]; The time–domain relaxation data can be converted into a frequency–domain electrochemical impedance spectroscopy (EIS), which is effective in capturing the flush delay and the rate-limiting step [25]; the confidence of the estimated k and D can be evaluated by a sensitivity analysis [26]. All the approaches are analogous to the EIS analysis methods based on equivalent circuit models. While there have been several emerging high-resolution methods for EIS analysis, such as the distribution of relaxation times (DRT) [27,28,29], there is a need to develop new methods independent of the trial-and-error approach for resolving of the CR kinetics and parameters.

By the existing ECR method, the analytical solutions to the ECR data are an infinite series of exponential terms represented by the characteristic times (CTs) and the corresponding pre-exponential factor (or strength) of each CT. Therein, the k and D are imbedded in the CTs and strengths in a non-linear, non-analytical manner. Therefore, the estimation of k and D using non-linear curve fitting procedure face challenges, although it is technically feasible. Recently, our team developed a new method namely distribution of characteristic times (DCT) to analyze the chemical relaxation in porous MIECs [30]. In that work, the time–domain ECR data is converted to a spectrum where the strength is plotted as a continuous function of logarithmic CT. For porous MIECs, the DCT spectrum was interpreted by a DCT model, and the values of k were estimated using the characteristics of the surface exchange-dominated peak in the DCT spectra. In principle, the DCT method is applicable to dense bulk MIEC bars/sheets/disks that are usually used in the literature. The DCT method provides a different manner to present the ECR data. As compared to the time–domain representation of ECR data, the DCT spectrum can visualize the CTs and their strengths without any preknowledge of the system. From the DCT spectrum, the intrinsic information of the ECR process can be visualized. In this work, the feasibility of estimation of k and D using the DCT spectrum reconstructed from the ECR data of dense bulk MIEC is demonstrated. It will be shown that the DCT method has a clear advantage over the time–domain analysis method.

2. The Theory of the Distribution of Characteristic Times

To deduce the DCT theory for dense bulk MIEC materials, it is necessary to review the analytical solution to the ECR data. For example, a dense MIEC sheet with a thickness 2L_x yields the solution [13],

$$\sigma \left(t\right)=1-{\displaystyle \sum _{k=1}^{\infty }\frac{2B{i}_x^2\exp \left(-{\beta }_k^{2}Dt/L_x{}^2\right)}{{\beta }_k^2\left({\beta }_k^2+B{i}_x^2+B{i}_x\right)}}$$

(1)

where the Biot numbers are given by,

$$B{i}_x=L_{x}k/D$$

(2)

The dimensionless parameters β_k are the kth roots of the following equations,

$${\beta }_k\tan {\beta }_k=B{i}_x$$

(3)

The infinite series solution is obtained from the eigenfunctions for space and time by separation of the variables method. By inspection of Equation (1), we can see that this solution can be rewritten as a series of exponential functions of CTs (τ_i) and strengths (λ_i),

$$\sigma \left(t\right)=1-{\displaystyle \sum _{i=1}^{\infty }{\lambda }_i\exp \left(-\frac{t}{{\tau }_i}\right)}$$

(4)

where ${\displaystyle \sum _{i=1}^{\infty }{\lambda }_i}=1$. It is noted that, if there is a relaxation in the change of atmosphere, say a flush delay of τ_f, Equation (4) is altered to [31],

$$\sigma \left(t\right)=1-{\displaystyle \sum _{i=1}^{\infty }{\lambda }_i\frac{{\tau }_i}{{\tau }_i-{\tau }_f}\exp \left(-\frac{t}{{\tau }_i}\right)}-\left(1-{\displaystyle \sum _{i=1}^{\infty }{\lambda }_i\frac{{\tau }_i}{{\tau }_i-{\tau }_f}}\right)\exp \left(-\frac{t}{{\tau }_f}\right)$$

(5)

By treating τ_f as a new characteristic time, Equation (5) can be rewritten as,

$$\sigma \left(t\right)=1-{\displaystyle \sum _{i=1}^{\infty }{\chi }_i\exp \left(-t/{\tau }_i\right)}$$

(6)

where ${\displaystyle \sum _{i=1}^{\infty }{\chi }_i}=1$. Obviously, if the change of atmosphere is instantaneous (τ_f = 0), we have χ_i > 0, else if there is a flush delay subject to τ_f << τ₁, we have χ_i ≤ 0 for τ_i ≤ τ_f, and χ_i > 0 for τ_i > τ_f. It is natural to see that the right-hand-side of Equation (6) can be generally represented by an integral,

$$\sigma \left(t\right)=1-{\displaystyle {\int }_{-\infty }^{+\infty }\left[\chi \exp \left(-t/\tau \right)\right]d{\log }_{10}\tau }$$

(7)

subject to,

$${\displaystyle {\int }_{-\infty }^{+\infty }\chi d{\log }_{10}\tau }=1$$

(8)

Equation (7) suggests two manners to represent the ECR data. The left-hand-side of Equation (7) corresponds to a time–domain representation, such as a plot of t versus σ. By this plot, it is easy to present the experimental ECR data, however, it is hard to see the rate-limiting step of the ECR experiment that is vital to the determination of k and D. In contrast, the right-hand-side of Equation (7) suggests a DCT representation, such as a plot of log₁₀τ versus χ. The DCT function χ(log₁₀τ) can be reconstructed from the experimental σ(t) data, by converting Equations (7) and (8) to a standard form of bound constrained quadratic programming problem based on a Tikhonov regularization [30]. A plot of log₁₀τ versus χ represents a new type of spectrum, showing the kinetic details that are hardly detectable by a time–domain representation of the ECR data.

3. Results and Discussion

A systematic study for the feasibility of resolving the values of k and D from ECR data was verified based on artificial ECR experiments. The ECR data was generated using the analytical solution of the MIEC sheet, for which the thickness (2L_x) is much smaller than the other dimensions. k and D are given as 1 × 10⁻⁴ cm·s⁻¹, and 1 × 10⁻⁵ cm²·s⁻¹, respectively. L_x serves as a variable parameter to simulate the effect of Biot number (Bi ≡ L_xk/D). The analytical solution of DCT, represented by a series of CTs and the corresponding strengths (τ_i, A_i) was obtained directly from the analytical ECR solution. On the other hand, the DCT represented by the log₁₀τ–χ relationship was reconstructed using the synthetic ECR data. The following contents first show the feasibility and merits of DCT spectroscopy in visualizing the rate-limiting step and the flush delay. Then, the robustness against noise and the accuracy in determining the values of k and D were demonstrated.

Figure 1 shows the comparison between the various representations of the synthetic ECR data with different Biot numbers. From the time–domain representation of the ECR data (Figure 1a,c,e) under Biot numbers of 0.01 (surface-exchange limited), 1 (surface-exchange & bulk-diffusion co-limited) and 100 (bulk-diffusion limited), it is hard to distinguish the rate-limiting mechanisms. From the reconstructed DCT spectra (the black curves in Figure 1b,d,f), it is shown clearly that the DCT shows a single peak for Bi = 0.01. For Bi = 1 and 100, the second-order peak appears in the DCT, with the strength increasing with Bi. It is noted that, for Bi = 100, the higher-order peaks appeared, which cannot be reconstructed properly by the DCT. These incorrect peaks may be derived from the Tikhonov regularization, since their strengths are negligible as compared to the first-order peaks. The similar problem is also encountered in reconstructing the DRT from the EIS spectrum by Tikhonov regularization [32]. Although the third and higher order peaks may be miscalculated, the first two major peaks with significant strength play a dominant role in visualizing the rate-limiting step. It is shown by the τ_i–A_i plot that the first- and second-order peaks of the reconstructed DCT agree well with the theoretical solutions (the orange stems in Figure 1b,d,f). Therefore, the reconstructed DCT spectrum can be used to identify the rate-limiting mechanism. To provide a whole picture that how the major peak(s) of DCT evolves with rate-limiting regimes, a contour plot of the reconstructed DCT in the logarithmic Bi–τ axis system is presented in Figure 2a.

It is shown that, for the surface-exchange-limited kinetics (e.g., Bi < 10⁻¹), the DCT presents only the first-order peak (P1). For the bulk-diffusion-limited kinetics (e.g., Bi > 10), the second-order peak (P2) appeared. In this condition, the ratio of the CTs of P1 and P2 was τ_P1/τ_P2 = 9. For the colimited kinetics (e.g., 10⁻¹ < Bi < 10), the ratio was subject to τ_P1/τ_P2 > 9. These features were verified by comparing with a scatter plot of the analytical DCT shown in Figure 2b. Therefore, the rate-limiting mechanism, which was hardly revealed from the time–domain representation of the ECR data, can be visualized explicitly in the DCT spectrum.

Then, we turned to the issue of determining the values of k and D from the DCT spectrum. Mathematically, k and D can be determined using two quantitative characteristics of the DCT, for example the CT (τ_P1) and the strength (S_P1) of P1. In the surface-exchange limited regime, the DCT shows one peak (P1). In this case, only the value of k can be determined, given by,

$$k=L_x/{\tau }_{\mathrm{P}1}$$

(9)

In the bulk-diffusion limited regime, the DCT shows two or more peaks subject to τ_P1/τ_P2 = 9. In this case, only the value of D can be determined, given by,

$$D=4L_x{}^2/\left({\pi }^2{\tau }_{\mathrm{P}1}\right)$$

(10)

In the co-limited regime, both the values of k and D are obtainable. In this case, both S_P1 and τ_P1 are needed, following

$${\tau }_{\mathrm{P}1}=\frac{L_x{}^2}{D{\alpha }_1{}^2}$$

(11)

$$S_{\mathrm{P}1}=\frac{2B{i}^2}{{\alpha }_1{}^2\left({\alpha }_1{}^2+B{i}^2+Bi\right)}$$

(12)

$$Bi\equiv \frac{L_{x}k}{D}={\alpha }_1\tan {\alpha }_1$$

(13)

The first-order root of Equation (13), α₁ approaches to 0, and to π/2 for the surface-exchange-limited and the bulk-diffusion-limited kinetics, respectively. Obviously, the accuracy of S_P1 and τ_P1 revealed by the DCT spectrum is vital to the determination of k and/or D. Figure 3 shows that the values of S_P1 and τ_P1, resolved from the DCT spectra in Figure 2a, fit excellently to the analytical values (A₁ and τ₁).

The above discussion shows that it is theoretically feasible to determine the values of k and D from the DCT spectrum. However, in practice, non-ideal factors are usually imbedded in the ECR data that can influence the quality of DCT spectrum, therefore the determination of k and D. One major factor is the flush delay (τ_f), denoting the relaxation time of the gaseous composition change. By considering the effect of τ_f, Equation (12) is revised to be,

$$S_{\mathrm{P}1}\frac{{\tau }_{\mathrm{P}1}-{\tau }_f}{{\tau }_{\mathrm{P}1}}=\frac{2B{i}^2}{{\alpha }_1{}^2\left({\alpha }_1{}^2+B{i}^2+Bi\right)}$$

(14)

Therefore, the determination of k and D demands the correct estimation of S_P1, τ_P1 and τ_f. It is hard to detect the effect of τ_f on the ECR data from the time–domain representation. However, as shown by Equation (5), the strengths of the CTs below τ_f are theoretically negative values. Therefore, the DCT spectrum is expected to be a visible manner identifying the flush delay. Figure 4 demonstrates an example that how a rational DCT spectrum and a correct value of τ_f can be determined from a flush-delay-imbedded ECR data.

The ECR data was generated using Equation (5) with τ_f = 3 s, k = 10⁻⁴ cm·s⁻¹, D = 10⁻⁵ cm²·s⁻¹ and Bi = 0.01. For reconstruction of DCT, a guess of τ_f must be given, since DCT is constricted to be negative/positive for CTs below/above τ_f. To determine the value of τ_f, a parametric sweep study on τ_f is performed to obtain the evolution of standard deviation between the original ECR and the ECR determined by the reconstructed DCT, as shown by Figure 4a. It is shown that the evolution of standard deviation with τ_f shows an L-shaped curve. The corner point of the L-shaped curve can be used as a good estimator of τ_f, which was estimated to be 3.07 s for this case (very close to the theoretical value of 3 s). Figure 4b shows the comparison of the original ECR and the DCT-fitted ECR, showing that the ECR could be perfectly fitted by the DCT when a correct value of τ_f is given. The corresponding DCT spectrum, as shown by Figure 4c, shows a negative peak at the CT of 3 s. To check the effectiveness of the procedure, a case study of τ_f with different values are shown in Figure 5, showing that the values of τ_f can be estimated correctly by the L-shaped curve of standard deviation.

Another factor that influences the quality of ECR data is the noisy perturbations. Therefore, the robustness of DCT reconstruction using noise-containing ECR data must be verified for practical usage. Based on the flush-delay-imbedded ECR, the effect of noise on the DCT reconstruction was studied by adding a random Gaussian noise with a standard deviation of 1% in the ECR data in Figure 4. Figure 6a shows the L-shaped curves of standard deviation for three-times repeats of parametric sweep on τ_f, showing that a correct value of 3 s can be estimated. To check if the correct values of S_P1 and τ_P1 can be resolved from the DCT, the noise-containing ECR data with different Biot numbers were simulated. Figure 6b,c shows that the values of S_P1 and τ_P1 resolved from the reconstructed DCT agreed nicely to the analytical values.

Finally, we demonstrated the capability of determining the values of k and D using the DCT spectrum reconstructed from noise-containing ECR data with different levels of flush delay (τ_f/τ₁ = 0, 0.2, 0.4) and Biot numbers (Bi = 10⁻² ~ 10²), as shown in Figure 7a–c.

It is shown that the time–domain ECR curves present a similar shape, from which it is hard to identify the flush delay. The reconstructed DCT spectra are shown in Figure 7d–f, showing that the flush delay can be explicitly represented by the negative peak. For τ_f/τ₁ = 0 (Figure 7d), the rate-limiting mechanism can be captured by the DCT spectra. For τ_f/τ₁ = 0.2 and 0.4 (Figure 7e,f), it is noted that, for bulk-diffusion-limited kinetics, the second-order peak was not reconstructed properly and pseudo small peaks around τ_f appeared. This resulted from the high level of flush delay that was higher than the CT of the second-order peak (τ₂/τ₁ = 1/9), so that the second-order peak was altered to be negative. Another reason is the strong interaction between τ_f and τ₂ that cannot be properly decoupled by the present Tikhonov regularization method. The same challenge is also encountered in the reconstruction of DRT from impedance spectroscopy. That is, the small polarization resistance process is hard to be decoupled by the DRT from the high resistance process with a similar relaxation time [32]. Pseudo small DRT peaks are hard to be avoided by the present methods, such as Tikhonov regularization, and Fourier transform [24]. In analogous to DRT, the reconstruction of DCT requires high quality ECR data. In practice, the flush delay is typically several seconds, which is usually much smaller than τ₁ (e.g., thousands of seconds for dense bulk MIECs [11,17,33,34]). If the flush delay is significant so that the second-order peak interplays strongly with the flush delay peak (τ₂ ≈ τ_f), the DCT spectrum can still serve as a tool for identifying the flush delay. For determination of k and D, the first-order peak and the flush delay played a crucial role. Figure 7d–f shows that τ_P1 and τ_f could be correctly estimated. Figure 7g shows that the value of S_P1 could be properly estimated for the various flush delays and Biot numbers. It is shown that the value of S_P1 increased with the increase of τ_f/τ_P1. From Equation (14), we can see that S_P1 is proportional to τ_f/τ_P1 at a fixed Biot number. Equation (14) also indicates that, combining with Equation (13), the value of the Biot number can be estimated using the values of τ_f, τ_P1 and S_P1. Figure 7h shows the calculated Biot numbers as a function of the theoretical Biot number for the different levels of flush delay. It is shown that the Biot number could be determined properly above a theoretical value of “1”. The calculated Biot number was overestimated if the value was below “1”. This overestimation resulted from the low estimation of S_P1 at small Biot numbers, as shown by the case τ_f/τ₁ = 0.4 in Figure 7g. According to the estimated Biot number, the values of k and D could be estimated using Equations (11) and (13). It is shown in Figure 7i that the value of k could be properly estimated under Bi < 10, while the value of D could be properly estimated under Bi > 1. It is noted that, for the case τ_f/τ_P1 = 0, the value of D is obtainable under Bi > 0.1, for which the Biot number can be properly estimated (Figure 7h). Therefore, the DCT spectrum demonstrated the same capability with the time–domain analysis method in determining k and D when the flush delay is negligible [26]. For the case that flush delay was significant, the DCT spectrum was capable of determining the value of k and D when the estimated Biot number was below “10” and above “1”, respectively. More importantly, the value of τ_f could be determined explicitly. The conjunct determination of k, D and τ_f was not demonstrated by the existing time–domain analysis methods. The DCT spectrum provides a visual, high-resolution and reliable manner for revealing chemical relaxation processes and resolving the kinetic parameters.

4. Conclusions

In summary, the DCT spectrum converted from ECR data was demonstrated to be a visual, high-resolution and robust tool for revealing the chemical relaxation kinetics in the dense bulk MIEC materials. As compared to the time–domain representation of ECR data, the DCT spectrum was capable of visualizing the rate-limiting mechanism and the flush delay. The values of surface exchange and bulk diffusion coefficients and flush delay could be determined conjunctly using the characteristics of the first-order peak and the flush delay peak of the DCT spectrum. The robustness against noise of DCT reconstruction was verified to be sufficient for practical usage. The DCT method was generally applicable to other types of chemical relaxation measurements for revealing rate-limiting steps and determining kinetic parameters, such as carburizing and nitriding of steels.