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

**Fuyao Yan <sup>1</sup> , Yiheng Wang <sup>1</sup> , Ying Yang <sup>1</sup> , Lei Zhu <sup>1</sup> , Hui Hu <sup>1</sup> , Zhuofu Tang <sup>1</sup> , Yanxiang Zhang 1,\*, Mufu Yan 1,\*, Changrong Xia 2,\* and Yueming Xu <sup>3</sup>**


Received: 13 November 2020; Accepted: 10 December 2020; Published: 17 December 2020

**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, τ*<sup>f</sup>* ) 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 τ*<sup>f</sup>* are easily recognized, and the values of *k*, *D* and τ*<sup>f</sup>* 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.

**Keywords:** distribution of characteristic times; electrical conductivity relaxation; surface exchange coefficient; bulk diffusion coefficient

## **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–3]. Another reason is that they are a very good platform for the fundamental research of solid state ionics and electrochemistry [4–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 −1 ) and bulk diffusion coefficient (*D*, cm<sup>2</sup> ·s −1 ) [9–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–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–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, τ*<sup>f</sup>* ), 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*<sup>c</sup> = *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–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 2*L<sup>x</sup>* yields the solution [13],

$$\sigma(t) = 1 - \sum\_{k=1}^{\infty} \frac{2Bt\_x^2 \exp\left(-\beta\_k^2 Dt / L\_x ^2\right)}{\beta\_k^2 \left(\beta\_k^2 + Bt\_x^2 + Bi\_x\right)}\tag{1}$$

where the Biot numbers are given by,

$$\text{Bi}\_{\text{x}} = \text{L}\_{\text{x}} \text{k} / \text{D} \tag{2}$$

The dimensionless parameters β*<sup>k</sup>* are the *k*th roots of the following equations,

$$
\beta\_k \tan \beta\_k = Bi\_\chi \tag{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 (τ*<sup>i</sup>* ) and strengths (λ*<sup>i</sup>* ),

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

where <sup>P</sup><sup>∞</sup> *i*=1 λ*<sup>i</sup>* = 1. It is noted that, if there is a relaxation in the change of atmosphere, say a flush delay of τ*<sup>f</sup>* , Equation (4) is altered to [31],

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

By treating τ*<sup>f</sup>* as a new characteristic time, Equation (5) can be rewritten as,

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

where <sup>P</sup><sup>∞</sup> *i*=1 χ*<sup>i</sup>* = 1. Obviously, if the change of atmosphere is instantaneous (τ*<sup>f</sup>* = 0), we have χ*<sup>i</sup>* > 0, else if there is a flush delay subject to τ*<sup>f</sup>* << τ1, we have χ*<sup>i</sup>* ≤ 0 for τ*<sup>i</sup>* ≤ τ*<sup>f</sup>* , and χ*<sup>i</sup>* > 0 for τ*<sup>i</sup>* > τ*<sup>f</sup>* . It is natural to see that the right-hand-side of Equation (6) can be generally represented by an integral,

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

subject to,

$$\int\_{-\infty}^{+\infty} \chi d \log\_{10} \tau = 1 \tag{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 log10τ versus χ. The DCT function χ(log10τ) 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 log10τ 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 (2*Lx*) is much smaller than the other dimensions. *k* and *D* are given as 1 <sup>×</sup> <sup>10</sup>−<sup>4</sup> cm·s −1 , and 1 <sup>×</sup> <sup>10</sup>−<sup>5</sup> cm<sup>2</sup> ·s −1 , respectively. *L<sup>x</sup>* 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 (τ*<sup>i</sup>* , *A<sup>i</sup>* ) 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. *Coatings* **2020**, *10*, x FOR PEER REVIEW 4 of 11 **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 (2*Lx*) is much smaller than the other dimensions. *k* and *D* are given as 1 × 10−4 cm·s−1, and 1×10−5 cm2·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

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<sup>i</sup>* 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. 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.

**Figure 1.** Comparison of the time domain representation of the electrical conductivity relaxation (ECR) data and the reconstructed distribution of characteristic times (DCT) spectra. (**a**,**c**,**e**) The time– domain representation of ECR data calculated from analytical solution of the mixed ionic-electronic conducting (MIEC) sheet with *k* = 10<sup>−</sup>4 cms−1, *D* = 10<sup>−</sup>5 cm2·s−1 and a Biot number of 0.01 (**a**), 1 (**c**) and 100 (**e**). (**b**,**d**,**f**) The DCT spectra calculated from (**a**,**c**,**e**) and the corresponding analytical solutions of DCT. **Figure 1.** Comparison of the time domain representation of the electrical conductivity relaxation (ECR) data and the reconstructed distribution of characteristic times (DCT) spectra. (**a**,**c**,**e**) The time–domain representation of ECR data calculated from analytical solution of the mixed ionic-electronic conducting (MIEC) sheet with *k* = 10−<sup>4</sup> cm·s −1 , *D* = 10−<sup>5</sup> cm<sup>2</sup> ·s <sup>−</sup><sup>1</sup> and a Biot number of 0.01 (**a**), 1 (**c**) and 100 (**e**). (**b**,**d**,**f**) The DCT spectra calculated from (**a**,**c**,**e**) and the corresponding analytical solutions of DCT.

(*A*1 and *τ*1).

*Coatings* **2020**, *10*, x FOR PEER REVIEW 5 of 11

**Figure 2.** The reconstructed DCT spectra (**a**) and the analytical solution of DCT (**b**) of MIEC sheets at various Biot numbers, with *k* = 10<sup>−</sup>4 cm·s−1 and *D* = 10<sup>−</sup>5 cm2·s−1. The strength of characteristic time (CT) scales with the color (**a**,**b**) and the size of scatter (**b**). **Figure 2.** The reconstructed DCT spectra (**a**) and the analytical solution of DCT (**b**) of MIEC sheets at various Biot numbers, with *k* = 10−<sup>4</sup> cm·s <sup>−</sup><sup>1</sup> and *D* = 10−<sup>5</sup> cm<sup>2</sup> ·s −1 . The strength of characteristic time (CT) scales with the color (**a**,**b**) and the size of scatter (**b**).

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. It is shown that, for the surface-exchange-limited kinetics (e.g., *Bi* < 10−<sup>1</sup> ), 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−<sup>1</sup> < *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, 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\_{\rm x} / \tau\_{\rm P1} \tag{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, 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\_{\rm Pl} \right) \tag{10}$$

are needed, following 2 *Lx* 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\_{\rm Pl} = \frac{L\_{\rm x}^{\prime 2}}{D\alpha\_1^{\prime 2}} \tag{11}$$

$$S\_{\rm P1} = \frac{2Bi^2}{\alpha\_1^{-2}(\alpha\_1^{-2} + Bi^2 + Bi)}\tag{12}$$

$$Bi \equiv \frac{L\_{\chi}k}{D} = \alpha\_1 \tan \alpha\_1 \tag{13}$$

$$\dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots \quad \dots$$

1 1 tan *<sup>D</sup>* ≡ =α α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 *S*P1 and *τ*P1 revealed by the DCT spectrum is vital to the determination of *k* and/or *D*. Figure 3 shows that the The first-order root of Equation (13), α<sup>1</sup> 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*<sup>1</sup> and τ1).

values of *S*P1 and *τ*P1, resolved from the DCT spectra in Figure 2a, fit excellently to the analytical values

*Coatings* **2020**, *10*, x FOR PEER REVIEW 6 of 11

**Figure 3.** (**a**) The major peak strength of the reconstructed DCT (the integral area of the peak P1 in Figure 2a, *S*P1) and the analytical strength of *τ*1 (*A*1), and (**b**) the CTs of the major peak in the reconstructed DCT (*τ*P1) and the analytical values of *τ*1 of MIEC sheets with *k* = 10<sup>−</sup>4 cms−1 and *D* = 10<sup>−</sup><sup>5</sup> cm2·s−1, as a function of the Biot number. **Figure 3.** (**a**) The major peak strength of the reconstructed DCT (the integral area of the peak P1 in Figure 2a, *S*P1) and the analytical strength of τ<sup>1</sup> (*A*<sup>1</sup> ), and (**b**) the CTs of the major peak in the reconstructed DCT (τP1) and the analytical values of τ<sup>1</sup> of MIEC sheets with *k* = 10−<sup>4</sup> cm·s <sup>−</sup><sup>1</sup> and *D* = 10−<sup>5</sup> cm<sup>2</sup> ·s −1 , as a function of the Biot number. **Figure 3.** (**a**) The major peak strength of the reconstructed DCT (the integral area of the peak P1 in Figure 2a, *S*P1) and the analytical strength of *τ*1 (*A*1), and (**b**) the CTs of the major peak in the reconstructed DCT (*τ*P1) and the analytical values of *τ*1 of MIEC sheets with *k* = 10<sup>−</sup>4 cms−1 and *D* = 10<sup>−</sup><sup>5</sup> cm2·s−1, as a function of the Biot number.

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, 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 (τ*<sup>f</sup>* ), denoting the relaxation time of the gaseous composition change. By considering the effect of τ*<sup>f</sup>* , Equation (12) is revised to be, 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,

$$\mathcal{S}\_{\rm P1} \frac{\tau\_{\rm P1} - \tau\_f}{\tau\_{\rm P1}} = \frac{2Bi^2}{\alpha\_1^2 (\alpha\_1^2 + Bi^2 + Bi)}\tag{14}$$

2

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. Therefore, the determination of *k* and *D* demands the correct estimation of *S*P1, τP1 and τ*<sup>f</sup>* . It is hard to detect the effect of τ*<sup>f</sup>* on the ECR data from the time–domain representation. However, as shown by Equation (5), the strengths of the CTs below τ*<sup>f</sup>* are theoretically negative values. Therefore, the DCT spectrum is expected to be a visible manner identifying the flush delay. Figure <sup>4</sup> demonstratesan example that how a rational DCT spectrum and a correct value of <sup>τ</sup>*<sup>f</sup>* can be determined from aflush-delay-imbedded ECR data. 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.

**Figure 4.** Reconstruction of DCT from the analytical ECR solution of an MIEC sheet with a flush delay *τf* = 3 s, *k* = 10<sup>−</sup>4 cm·s−1, *D* = 10<sup>−</sup>5 cm2·s−1 and a Biot number of 0.01. (**a**) The standard deviation (std) between DCT fitting and the original ECR data as a function of *τf* for use in DCT reconstruction. (**b**) The comparison between the original ECR data and the DCT fitting with different values of *τf*. (**c**) The reconstructed DCT spectrum (the left axis), and the analytical DCT (the left axis). **Figure 4.** Reconstruction of DCT from the analytical ECR solution of an MIEC sheet with a flush delay *τf* = 3 s, *k* = 10<sup>−</sup>4 cm·s−1, *D* = 10<sup>−</sup>5 cm2·s−1 and a Biot number of 0.01. (**a**) The standard deviation (std) between DCT fitting and the original ECR data as a function of *τf* for use in DCT reconstruction. (**b**) The comparison between the original ECR data and the DCT fitting with different values of *τf*. (**c**) The reconstructed DCT spectrum (the left axis), and the analytical DCT (the left axis). **Figure 4.** Reconstruction of DCT from the analytical ECR solution of an MIEC sheet with a flush delay τ*<sup>f</sup>* = 3 s, *k* = 10−<sup>4</sup> cm·s −1 , *D* = 10−<sup>5</sup> cm<sup>2</sup> ·s <sup>−</sup><sup>1</sup> and a Biot number of 0.01. (**a**) The standard deviation (std) between DCT fitting and the original ECR data as a function of τ*<sup>f</sup>* for use in DCT reconstruction. (**b**) The comparison between the original ECR data and the DCT fitting with different values of τ*<sup>f</sup>* . (**c**) The reconstructed DCT spectrum (the left axis), and the analytical DCT (the left axis).

The ECR data was generated using Equation (5) with *τf* = 3 s, *k* = 10−4 cm·s−1, *D* = 10−5 cm2·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 The ECR data was generated using Equation (5) with *τf* = 3 s, *k* = 10−4 cm·s−1, *D* = 10−5 cm2·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 The ECR data was generated using Equation (5) with τ*<sup>f</sup>* = 3 s, *k* = 10−<sup>4</sup> cm·s −1 , *D* = 10−<sup>5</sup> cm<sup>2</sup> ·s −1 and *Bi* = 0.01. For reconstruction of DCT, a guess of τ*<sup>f</sup>* must be given, since DCT is constricted to be negative/positive for CTs below/above τ*<sup>f</sup>* . To determine the value of τ*<sup>f</sup>* , a parametric sweep study on τ*<sup>f</sup>*

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 τ*<sup>f</sup>* shows an L-shaped curve. The corner point of the L-shaped curve can be used as a good estimator of τ*<sup>f</sup>* , 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 τ*<sup>f</sup>* 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 τ*<sup>f</sup>* with different values are shown in Figure 5, showing that the values of τ*<sup>f</sup>* 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. *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 *<sup>τ</sup>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 *<sup>τ</sup>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 *<sup>τ</sup>f* with different values are shown in Figure 5, showing that the values of *<sup>τ</sup>f* can be estimated correctly by the L-shaped curve of standard deviation.

**Figure 5.** Determination of *τf* for the various flush-delay-containing ECR data by the evolution of standard deviation (std) with a sweep of *τf*. The red dots represent the coordinates of the calculated *τf* and the theoretical *τf* determined as the point where the std does not decrease obviously with *τf*. The 45° line is for eye-guide. **Figure 5.** Determination of τ*<sup>f</sup>* for the various flush-delay-containing ECR data by the evolution of standard deviation (std) with a sweep of τ*<sup>f</sup>* . The red dots represent the coordinates of the calculated τ*<sup>f</sup>* and the theoretical τ*<sup>f</sup>* determined as the point where the std does not decrease obviously with τ*<sup>f</sup>* . The 45◦ line is for eye-guide. **Figure 5.** Determination of *τf* for the various flush-delay-containing ECR data by the evolution of standard deviation (std) with a sweep of *τf*. The red dots represent the coordinates of the calculated *τf* and the theoretical *τf* determined as the point where the std does not decrease obviously with *τf*. The 45° line is for eye-guide.

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. 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 τ*<sup>f</sup>* , 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. 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 *<sup>τ</sup>f*, showing that a correct value of 3 s can be estimated. To check if the correct values of *<sup>S</sup>*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.

datasets with the same theoretical solution imbedded with 1% standard deviation level random Gaussian noise), (**b**) the major peak strength of the reconstructed DCT (*S*P1) and the analytical strength of *τ*1 (*A*1) and (**c**) the CTs of the major peak in the reconstructed DCT (*τ*P1) and the analytical values of *τ*1 for noise-containing ECR data of MIEC sheet with a flush delay *τf* = 3 s, *k* = 10<sup>−</sup>4 cm·s−1 and *D* = 10<sup>−</sup><sup>5</sup> **Figure 6.** Determination of the flush delay ((**a**) the three lines represent the results of three ECR datasets with the same theoretical solution imbedded with 1% standard deviation level random Gaussian noise), (**b**) the major peak strength of the reconstructed DCT (*S*P1) and the analytical strength of *τ*1 (*A*1) and (**c**) the CTs of the major peak in the reconstructed DCT (*τ*P1) and the analytical values of *τ*1 for noise-containing ECR data of MIEC sheet with a flush delay *τf* = 3 s, *k* = 10<sup>−</sup>4 cm·s−1 and *D* = 10<sup>−</sup><sup>5</sup> **Figure 6.** Determination of the flush delay ((**a**) the three lines represent the results of three ECR datasets with the same theoretical solution imbedded with 1% standard deviation level random Gaussian noise), (**b**) the major peak strength of the reconstructed DCT (*S*P1) and the analytical strength of τ<sup>1</sup> (*A*<sup>1</sup> ) and (**c**) the CTs of the major peak in the reconstructed DCT (τP1) and the analytical values of τ1 for noise-containing ECR data of MIEC sheet with a flush delay τ*<sup>f</sup>* = 3 s, *k* = 10−<sup>4</sup> cm·s <sup>−</sup><sup>1</sup> and *D* = 10−<sup>5</sup> cm<sup>2</sup> ·s −1 . A random Gaussian noise with a 1% standard deviation level is added in the ECR data at each Biot number.

*Coatings* **2020**, *10*, 1240 cm2·s−1. A random Gaussian noise with a 1% standard deviation level is added in the ECR data at each Biot number.

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 (τ*<sup>f</sup>* /τ<sup>1</sup> = 0, 0.2, 0.4) and Biot numbers (*Bi* = 10−<sup>2</sup> ~ 10<sup>2</sup> ), as shown in Figure 7a–c. 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*/*τ*1 = 0, 0.2, 0.4) and Biot numbers (*Bi* = 10−2 ~ 102), as shown in Figure 7a–c.

*Coatings* **2020**, *10*, x FOR PEER REVIEW 8 of 11

**Figure 7.** Determination of *k* and *D* from reconstruction of DCT. (**a**–**c**) The synthesized ECR data for the MIEC sheets with *k* = 10<sup>−</sup>4 cms−1, *D* = 10<sup>−</sup>5 cm2·s−1, std = 1%, Bi = 10−2 ~ 102 and *τf*/*τ*1 = 0 (a), *τf*/*τ*1 = 0.2 (**b**) and *τf*/*τ*1 = 0.4 (c). (**d**–**f**) The reconstructed DCT from the synthesized noise-containing ECR datasets with *τf*/*τ*1 = 0 (**d**), *τf*/*τ*1 = 0.2 (**e**) and *τf*/*τ*1 = 0.4 (**f**). (**g**–**i**) The evolution of the major peak strength of the reconstructed DCT (*S*P1) and the analytical strength of *τ*1 (*A*1) (**g**), the calculated Biot number (**h**) and the calculated *D* and *k* (**i**) as a function of Biot number from the reconstructed DCT spectra with different levels of *τf*. **Figure 7.** Determination of *k* and *D* from reconstruction of DCT. (**a**–**c**) The synthesized ECR data for the MIEC sheets with *k* = 10−<sup>4</sup> cm·s −1 , *D* = 10−<sup>5</sup> cm<sup>2</sup> ·s −1 , std = 1%, Bi = 10−<sup>2</sup> ~ 10<sup>2</sup> and τ*<sup>f</sup>* /τ<sup>1</sup> = 0 (a), τ*f* /τ<sup>1</sup> = 0.2 (**b**) and τ*<sup>f</sup>* /τ<sup>1</sup> = 0.4 (c). (**d**–**f**) The reconstructed DCT from the synthesized noise-containing ECR datasets with τ*<sup>f</sup>* /τ<sup>1</sup> = 0 (**d**), τ*<sup>f</sup>* /τ<sup>1</sup> = 0.2 (**e**) and τ*<sup>f</sup>* /τ<sup>1</sup> = 0.4 (**f**). (**g**–**i**) The evolution of the major peak strength of the reconstructed DCT (*S*P1) and the analytical strength of τ<sup>1</sup> (*A*<sup>1</sup> ) (**g**), the calculated Biot number (**h**) and the calculated *D* and *k* (**i**) as a function of Biot number from the reconstructed DCT spectra with different levels of τ*<sup>f</sup>* .

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 τ*<sup>f</sup>* /τ<sup>1</sup> = 0 (Figure 7d), the rate-limiting mechanism can be captured by the DCT spectra. For τ*<sup>f</sup>* /τ<sup>1</sup> = 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 τ*<sup>f</sup>* appeared. This resulted from the high level of flush delay that was higher than the CT of the second-order peak (τ2/τ<sup>1</sup> = 1/9), so that the second-order peak was altered to be negative. Another reason is the strong interaction between τ*<sup>f</sup>* and τ<sup>2</sup> 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 τ<sup>1</sup> (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 (τ<sup>2</sup> ≈ τ*<sup>f</sup>* ), 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 τ*<sup>f</sup>* 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 τ*<sup>f</sup>* /τP1. From Equation (14), we can see that *S*P1 is proportional to τ*<sup>f</sup>* /τ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 τ*<sup>f</sup>* , τ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 τ*<sup>f</sup>* /τ<sup>1</sup> = 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 τ*<sup>f</sup>* /τ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 τ*<sup>f</sup>* could be determined explicitly. The conjunct determination of *k*, *D* and τ*<sup>f</sup>* 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.

**Author Contributions:** Data curation, F.Y.; Formal analysis, F.Y.; Funding acquisition, Y.Z., M.Y., and Y.X.; Investigation, F.Y., Y.W., Y.Y., L.Z., H.H., Z.T., and C.X.; Methodology, Y.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key Research and Development Program of China (2018YFB2001901) and National Natural Science Foundation of China (21673062, 52,073,072 and 51972298), and the Key Research and Development Program of Sichuan Province (2020YFSY0026).

**Conflicts of Interest:** There is no conflict of interest.

## **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
