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

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.


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

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], where the Biot numbers are given by, The dimensionless parameters β k are the kth roots of the following equations, 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 ), 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], By treating τ f as a new characteristic time, Equation (5) can be rewritten as, where ∞ i=1 χ 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 << τ 1 , 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, subject to, +∞ −∞ χd log 10 τ = 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 10 τ versus χ. The DCT function χ(log 10 τ) 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 10 τ versus χ represents a new type of spectrum, showing the kinetic details that are hardly detectable by a time-domain representation of the ECR data.

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 −4 cm·s −1 , and 1 × 10 −5 cm 2 ·s −1 , respectively. L x serves as a variable parameter to simulate the effect of Biot number (Bi ≡ L x k/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 10 τ-χ 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.

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 (2Lx) is much smaller than the other dimensions. k and D are given as 1 × 10 −4 cm·s −1 , and 1×10 −5 cm 2 ·s −1 , respectively. Lx serves as a variable parameter to simulate the effect of Biot number (Bi ≡ Lxk/D). The analytical solution of DCT, represented by a series of CTs and the corresponding strengths (τi, Ai) was obtained directly from the analytical ECR solution. On the other hand, the DCT represented by the log10τ-χ 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 colimited) 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-Ai 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 −1 ), 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 −1 < 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 (SP1) 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, 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, In the co-limited regime, both the values of k and D are obtainable. In this case, both SP1 and τP1 are needed, following The first-order root of Equation (13), α1 approaches to 0, and to π/2 for the surface-exchangelimited and the bulk-diffusion-limited kinetics, respectively. Obviously, the accuracy of SP1 and τP1 revealed by the DCT spectrum is vital to the determination of k and/or D. Figure 3 shows that the values of SP1 and τP1, resolved from the DCT spectra in Figure 2a, fit excellently to the analytical values (A1 and τ1). It is shown that, for the surface-exchange-limited kinetics (e.g., Bi < 10 −1 ), 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 −1 < 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, 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, 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 The first-order root of Equation (13), α 1 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 1 and τ 1 ). Coatings 2020, 10, x FOR PEER REVIEW 6 of 11 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, Therefore, the determination of k and D demands the correct estimation of SP1, τ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 −4 cm·s −1 , D = 10 −5 cm 2 ·s −1 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 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, 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.
Coatings 2020, 10, x FOR PEER REVIEW 6 of 11 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, Therefore, the determination of k and D demands the correct estimation of SP1, τ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 −4 cm·s −1 , D = 10 −5 cm 2 ·s −1 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  The ECR data was generated using Equation (5) with τ f = 3 s, k = 10 −4 cm·s −1 , D = 10 −5 cm 2 ·s −1 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 Coatings 2020, 10, 1240 7 of 11 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.
Coatings 2020, 10, x FOR PEER REVIEW 7 of 11 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 SP1 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 SP1 and τP1 resolved from the reconstructed DCT agreed nicely to the analytical values. 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.
Coatings 2020, 10, x FOR PEER REVIEW 7 of 11 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 SP1 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 SP1 and τP1 resolved from the reconstructed DCT agreed nicely to the analytical values.    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/τ1 = 0 (Figure 7d), the ratelimiting mechanism can be captured by the DCT spectra. For τf/τ1 = 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 (τ2/τ1 = 1/9), so that the second-order peak was altered 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 /τ 1 = 0 (Figure 7d), the rate-limiting mechanism can be captured by the DCT spectra. For τ f /τ 1 = 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 (τ 2 /τ 1 = 1/9), so that the second-order peak was altered to be negative. Another reason is the strong interaction between τ f and τ 2 that cannot be properly Coatings 2020, 10, 1240 9 of 11 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 τ 1 (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 (τ 2 ≈ τ 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 /τ 1 = 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.

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.

Conflicts of Interest:
There is no conflict of interest.