Extended Distribution of Relaxation Time Analysis for Electrochemical Impedance Spectroscopy

Abstract

The distribution of relaxation time (DRT) is increasingly investigated as a novel analytical method for electrochemical impedance spectroscopy. However, this method has not yet been generalized, as it cannot be applied to a spectrum influenced by an isolated capacitator and a resistor connected in parallel with an inductor. To generalize the DRT analysis, we propose a novel principal relation and actual calculation methods using an iterative elastic net regularization algorithm. The elastic net regularization was implemented using the Levenberg–Marquardt method. The proposed algorithm aids in analyzing the DRT for spectra that cannot be solved using conventional methods. We compared the DRT results in different regularization methods obtained from the proposed program with those obtained from other open-source programs. Additionally, the realistic problems of the DRT analysis are discussed considering the calculated results and their theoretical basis. Different DRT curves can be obtained depending on the particular program and regularization methods, even when the spectrum is the same. Moreover, overconfidence in the electrochemical DRT method can be avoided with a clear understanding of DRT basics.

1. Introduction

The distribution of relaxation time (DRT) analysis of the impedance (Z) spectrum is increasingly researched as a new method for electrochemical impedance spectroscopy analysis.^1–31 This can be attributed to the result generated by the DRT, which is objective in comparison with that obtained from a conventional impedance analysis using an equivalent circuit model,³² as a circuit model need not be set for fitting. Over the last decade, several types of DRT algorithms have been proposed,^1–26 which exhibit certain constraints, such as the imaginary part of Z (Z_imag) should approach 0 at high- and low-frequency limits, for applicable spectra. These constraints indicate that the DRT analysis cannot be used for the spectra collected from electrochemical devices, namely, batteries, capacitors, and corrosive systems. This is because the imaginary part of these spectra diverges with decreasing frequency.³³

Recently reported DRT programs have solved the problem of Z_imag divergence at the high-frequency limit. These methods involve the addition of an isolated inductor,^2–4,6,19 wherein Eq. 1 is used as the principal equation for DRT analysis.   

\begin{equation} Z (j\omega) = L_{0}j\omega + R_{\infty} + \int_{-\infty}^{\infty}\frac{g(\ln \tau)}{1 + \tau j\omega}\mathrm{d}\ln \tau, \end{equation} (1)

where j, ω, and Z(jω) denote the imaginary unit, angular frequency, and impedance at ω, respectively; and L₀, R_∞, τ, and g(ln τ) indicate the inductance, resistance at high-frequency limit, relaxation time, and DRT in logarithmic timescale, respectively. As the unit scale for τ is arbitrary, log τ can be used rather than ln τ in Eq. 1. However, we used the ln τ scale in this study. Despite using this equation, the DRT cannot be applied to a spectrum with diverging Z_imag at the low-frequency limit, which contains an isolated capacitor. Further, a spectrum including inductive loop owing to a resistor connected in parallel with an inductor is also nonapplicable. Therefore, the principal equation and calculation algorithm of DRT must be extended to develop a generalized method.

Although several DRT programs with different algorithms have recently been developed as open-source software,^11,19,34,35 DRT results have not been compared considering identical spectrum data. This is important because DRT curves may differ based on the program used.

In this study, we discuss the mathematical aspects and calculation methods to generalize the DRT analysis. Particularly, we discuss the solution methods of the residual problems; namely, the application of DRT when it is influenced by an isolated capacitor and inductive loop. We propose using pseudo-negative DRT considering the basic property of regularized complex nonlinear least-squares (CNLS). Additionally, we analyze the dependencies of the DRT results by employing regularization methods and programs. We then address the existing realistic problems of the DRT method.

2. Theory and Algorithm

2.1 Supporting method for Z_imag diverging spectrum at the low-frequency limit

The problem of Z_imag divergence at the low-frequency limit can be avoided by adding a capacitor in Eq. 1. In this case, the modified equation is represented as   

\begin{equation} Z (j\omega) = L_{0}j\omega + \frac{1}{C_{0}j\omega} + R_{\infty} + \int_{-\infty}^{\infty}\frac{g(\ln \tau)}{1 + \tau j\omega} \mathrm{d}\ln\tau, \end{equation} (2)

where C₀ denotes the capacitance. Although this equation is theoretically accurate, it is unsuitable for regularized CNLS, which is used for actual DRT calculations. This is because non-determinable parameters converge to zero in the case of regularized nonlinear least-squares,³⁶ and the capacitor term (1/C₀jω) diverges when C₀ is non-determinable. To avoid divergence, the capacitor term is replaced by a reciprocal capacitor ($C'_{0}$), as follows:   

\begin{equation} Z (j\omega) = L_{0}j\omega + \frac{C'_{0}}{j\omega} + R_{\infty} + \int_{-\infty}^{\infty}\frac{g(\ln \tau)}{1 + \tau j\omega} \mathrm{d}\ln \tau. \end{equation} (3)

After the optimization of the Z spectrum data using Eq. 3, the value of C₀ is calculated as 1/$C'_{0}$. When $C'_{0}$ is determined to be zero, the output is presented using a specific number of a computer, namely, “Not a number.”

2.2 Supporting method for the inductive loop

If the concept of conventional DRT relationship is applied to an inductive loop (Appendix A), Eq. 3 can be extended as follows:   

\begin{align} Z (j\omega) &= L_{0}j\omega + \frac{C'_{0}}{j\omega} + R_{\infty} + \int_{-\infty}^{\infty}\frac{g(\ln \tau)}{1 + \tau j\omega} \mathrm{d}\ln \tau \notag\\ &\quad+ \int_{-\infty}^{\infty}\frac{\gamma (\ln \tau')\tau' j\omega}{1 + \tau' j\omega} \mathrm{d}\ln\tau', \end{align} (4)

where τ′ and γ(ln τ′) denote the relaxation time and DRT in the logarithmic timescale of the inductive loop, respectively. However, this equation is invalid for optimization owing to numerous parameters even use of regularized CNLS. Therefore, a different method is required to support the inductive loop.

Considering the symmetry of the spectra between the two equivalent circuits illustrated in Figs. 1a and 1b, we determined that the inductive loop spectrum that is calculated using Fig. 1b can be alternatively reproduced by combining the resistors and capacitors (Fig. 1c) when a negative resistor and capacitor are added. These results are confirmed by the example spectra in Figs. 1d and 1e.

Figure 1.

(a) to (c) Equivalent circuits used for calculating the impedances illustrated in (d) to (f). The used element values are (a) R₁₁ = 10 Ω and C₁₁ = 100 µF, (b) R₂₁ = 10 Ω and L₂₁ = 10 mH, (c) R₃₀ = 10 Ω, R₃₁ = −10 Ω, and C₃₁ = −100 µF. The spectra obtained from the circuits in (b) and (c) using the aforementioned values overlap completely.

The impedance spectrum in Fig. 1b can be calculated as   

\begin{equation} Z_{R_{21}\|L_{21}} = \frac{R_{21}(L_{21}/R_{21})j\omega}{1 + (L_{21}/R_{21})j\omega} = \frac{R_{21}\tau j\omega}{1 + \tau j\omega}, \end{equation} (5)

and the spectrum in Fig. 1c can be obtained as   

\begin{align} Z_{R_{30} + R_{31}\|C_{31}} &= \frac{( R_{30} + R_{31}) + R_{30}R_{31}C_{31}j\omega}{1 + R_{31}C_{31}j\omega} \notag\\ &= \frac{(R_{30} + R_{31}) + R_{30}\tau j\omega}{1 + \tau j\omega}. \end{align} (6)

When R₃₀ + R₃₁ equals 0 and both τ and R₃₀ are positive values, Eqs. 5 and 6 are considered equivalent. This condition is established when R₃₀ > 0, R₃₁ = −R₃₀ < 0, and C₃₁ < 0.

Equation 6 indicates that the inductive loop spectrum can be treated by the DRT method if negative g(ln τ) and shift of R_∞ caused by the negative g(ln τ) are accepted. Based on the relation R₃₁ = −R₃₀, the value of R₂₁ in Fig. 1b can be calculated using the negative resistance generated by the DRT. This method is valid due to that g(ln τ) and γ(ln τ′) are not independent parameters but depended ones when negative g(ln τ) is used. The dependencies are principally given by Eqs. 5 and 6.

2.3 Calculation of DRT in a logarithmic timescale

We used CNLS to optimize a parameter vector (p) based on the iterative calculation and minimized the sum of least-squares (S).³³   

\begin{align} S &= \sum\nolimits_{m=0}^{m_{\max}}\biggl[ \frac{\{\mathit{Re}(Z[m]) - \mathit{Re}(Z^{\text{mdl}}[ j\omega_{m},\mathbf{p}])\}^{2}}{w_{m}^{\text{Zr}}} \notag\\ &\quad+ \frac{\{\mathit{Imag}(Z[m]) - \mathit{Imag}(Z^{\text{mdl}}[j\omega_{m},\mathbf{p}])\}^{2}}{w_{m}^{\text{Zi}}}\biggr], \end{align} (7)

where m denotes the data line number; Re() and Imag() indicate functions to output the real and imaginary parts of the input object; $w_{m}^{\text{Zr}}$ and $w_{m}^{\text{Zi}}$ represent the data weights of the real and imaginary parts, respectively; and Z^mdl[jω_m, p] denotes the impedance calculated by the model.³³

To obtain the optimized g(ln τ), we require an integral term in Eq. 3 for the numerical calculation after appropriately discretizing ln τ within the finite range. Therefore, it is necessary to introduce new parameters, such as the maximum and minimum values of ln τ, namely, ln τ_max and ln τ_min; divided ln τ under equally spaced ln τ_n between ln τ_max and ln τ_min, where n denotes the index number of divided ln τ; and equally divided space, Δ ln τ. Additionally, the following integral equation is assumed to generate a true answer (A + jB).   

\begin{equation} \int_{\log \tau_{n} - \frac{\Delta \ln \tau}{2}}^{\log \tau_{n} + \frac{\Delta \ln \tau}{2}} g (\ln \tau) \frac{1}{1 + \tau j\omega}\mathrm{d}\ln \tau = A + jB, \end{equation} (8)

wherein two uniquely determinable values (g^Zr(ln τ_n) and g^Zi(ln τ_n)) exist, which can be obtained as   

\begin{equation} g^{Z\text{r}}(\ln \tau_{n}) = \frac{A}{\Delta \ln\tau}\ \text{and}\ g^{Z\text{i}}(\ln \tau_{n}) = \frac{B}{\Delta \ln \tau}. \end{equation} (9)

Here, the DRT in logarithmic timescale can be redefined as   

\begin{equation} g(n) = \frac{g^{Z\text{r}}(\ln \tau_{n}) + g^{Z\text{i}}(\ln \tau_{n})}{2}. \end{equation} (10)

Subsequently, g(n) can be optimized because these values are determined uniquely. This treatment is equivalent such that g(ln τ), as a continuous function of ln τ, is replaced by a statistical average in each Δ ln τ in the range of (g(n)). Although this interpretation cannot fundamentally eliminate the inconsistency between discrete and continuous τ, both cases of τ can be included for numerical calculation as a tentative method. When g(n) is introduced, Z^mdl[jω_m, p] can be expressed as   

\begin{align} Z[j\omega_{m},\mathbf{p}] &= L_{0}j\omega_{m} + \frac{C'_{0}}{j\omega_{m}} + R_{\infty} \notag\\ &\quad+ \sum\nolimits_{n=0}^{n=n_{\max}}g(n)\frac{1}{1 + \tau_{n}j\omega_{m}}\Delta \ln \tau. \end{align} (11)

Herein, n need not be equal to the value of m because ln τ and g(ln τ) are treated as parameters in Eq. 3 (Appendix B).

Solving Eqs. 3 and/or 11 using CNLS is nearly impossible owing to an ill-posed problem caused by the instability of the optimized answer.^1–31,36 Several studies used the Tikhonov regularization method to avoid this problem.^{1–22,27,30,31} In this study, we implemented an advanced regularization method, namely, the elastic net regularization method,^{11,20,21,36,37} which is a parameter optimization algorithm that adds a penalty function to S. In the case of elastic net regularization, the penalty function (pf) is presented using the regularization parameter (λ) and partition coefficient (α) as follows:^21,38   

\begin{equation} pf = \lambda \{\alpha\|\textbf{p}\| + (1 - \alpha)\|\textbf{p}\|^{2}\}, \end{equation} (12)

where ||p|| denotes the norm of the p vector; ||p|| and ||p||² are obtained using p elements (p₀, p₁, …, p_n) as   

\begin{align} \|\mathbf{p}\| &= |\text{p}_{0}| + |\text{p}_{1}| + \cdots + |\text{p}_{n}|\ \text{and}\ \|\mathbf{p}\|^{2} \notag\\ &= |\text{p}_{0}|^{2} + |\text{p}_{1}|^{2} + \cdots + |\text{p}_{n}|^{2}. \end{align} (13)

where |p_n| indicates the absolute value of p_n. Therefore, the target function for optimizing the parameters using the elastic net regularization (S_enr) can be presented as   

\begin{equation} S_{\text{enr}} = S + pf. \end{equation} (14)

By the regularized CNLS, not only g(n) but also R_∞, C₀, and L₀ are optimized by searching minimum value of S_enr.

Elastic net regularization is a hybrid method that combines the Tikhonov and Lasso regularization methods. When α is set to 0, the analysis uses the Tikhonov regularization method, and the same program executes the Lasso regularization when α is set to 1.

2.4 Differences among DRT programs

We were unable to review the detailed differences among existing DRT programs as actual differences can be confirmed only through code analysis and their analyzing ability. Therefore, we compared the codes of open-source DRT programs with that of the proposed extended-DRT program. To the best of our knowledge, the two DRT program codes, namely, DRTtools^19,34 and EIS-DRT,^11,35 are open-source programs. In this study, we used the Python version of DRTtools (pyDRTtools) for comparison.³⁴ The primary differences between our program and the others are the support of the capacitor term and inductive loop spectrum.

Considering the principal relation indicated in Eq. 14, DRTtools employed different penalty functions and various interpolation methods for calculating the integral term in Eq. 1. The penalty function (pf_DRTTL) was implemented as a selectable function as   

\begin{equation} pf_{\text{DRTTL}} = \left\| \frac{\mathrm{d}g(\ln \tau)}{\mathrm{d}\ln \tau} \right\|\ \text{or}\ \left\| \frac{\mathrm{d}g(\ln \tau)}{\mathrm{d}\ln \tau} \right\|^{2}. \end{equation} (15)

Another difference is the method used to calculate the integral term. In DRTtools, intermediate values between g(ln τ_n) and g(ln τ_n+1) are interpolated using a function, such as linear or Gaussian functions,¹¹ which should be selected by the user before the calculation. Although this method is highly accurate, certain situations, such as the scenario where g(ln τ_n) is constructed by single and/or multiple impulse signals, are implicitly excluded. Although this scenario appears excludable, it can be executed based on the first principal equation as explained based on Eq. A1 in Appendices A and B.³⁹ In the case of our interpretation (Eqs. 8 to 10, the situation can be included if deconvolution within the Δ ln τ space is excluded.

According to Boukamp,⁴⁰ the equivalent calculation of our proposed method is possible using a Gaussian function with a full width at half maximum value of 0.15. This is referred to as the mollifier of the Gaussian function.³⁹ However, the DRT calculation diverged under this setting when DRTtools was used. This can be attributed to the peculiar penalty function, as explained in Eq. 15. Based on this result, the implicit assumption that DRT is a continuous function is included in the DRTtools results.

Although EIS-DRT implements the same penalty function as our program, the DRT is calculated using the imaginary part (Z_imag) of the impedance data, whereas the real part (Z_real) is used only for the optimization of R_∞. This method is employed in other programs that do not provide open-source codes.^{12,15,26,28,29} As Z data are constructed using Z_real and Z_imag, the DRT analysis should ideally use both Z_real and Z_imag under equal weighting. However, the analysis in EIS-DRT weighs toward Z_imag data.

3. Experimental

3.1 Implementation of the extended-DRT

The proposed extended-DRT function was coded on Igor Pro (WaveMetrics, Inc., USA) version 7 or higher. The DRT was implemented as a novel function of our integrated analysis program for electrochemical impedance.^41,42 The built-in Levenberg–Marquardt algorithm of Igor Pro was used for the optimization of DRT, L₀, C₀, and R_∞. This algorithm ensures calculation stability and high optimization ability for nonlinear functions. Furthermore, the optimized parameter set is searched through iterative calculation.

To perform flexible analysis, the values of ln τ_max, ln τ_min, Δ ln τ, α, and λ were used as the settable parameters on the graphic user interface.^41,42 Additionally, two conditions were set as user selectable. First, the constraints of the parameter values were set, wherein they can be only positive or no constraints should be applied. Second, either both parameters L₀ and C₀ or none should be added. In terms of $w_{m}^{\text{Zr}}$ and $w_{m}^{\text{Zi}}$ in Eq. 14, the parameter values were fixed at unity to stabilize the calculation.

3.2 Comparison of results

3.2.1 Simulated spectrum

To confirm the validity of the spectrum, including the influences of the capacitor and inductive loop, the simulated impedance spectrum was analyzed using our program. Figure 2 depicts the equivalent circuit and element values. The simulated impedance spectrum was generated between 1 MHz and 10 mHz, and the frequency values were spaced equally on a logarithmic scale with 10 points per order of frequency.

Figure 2.

Equivalent circuit used for calculating the simulated spectrum. The element values used for the calculation are also indicated.

Furthermore, this spectrum was analyzed using three programs, namely, the proposed program employing Tikhonov and Lasso regularization; DRTtools combined with Tikhonov regularization, piecewise linear discretization, and first-order derivative penalty function;^19,34 and EIS-DRT with Tikhonov regularization.^11,35 In the case of EIS-DRT, the programs were coded based on each target spectrum analyzed by the authors.³⁵ This implies that other users must modify the source code when they intend to analyze their own data. Owing to the difficulty in modifying the EIS-DRT source code, the calculation was performed using one file for MATLAB code, namely, (Example_Experiment_LiB.m),³⁵ on GNU Octave version 6.3.0.⁴³ The source code file was originally used for the spectrum collected by a lithium battery.^11,19 The details of the data are explained in a previous study.²¹ The simulated spectrum was analyzed based on two modifications, namely, a change in loading the data file and calculating a basis matrix (A) using the “unifrnd” function. Although basis matrix A was imported from two files by default, the matrix size relied on the impedance data size. Therefore, the size of matrix A was calculated considering the simulated impedance data.

3.2.2 Measured spectrum

In studies on DRTtools^11,34 and EIS-DRT,^19,35 the same spectrum collected by a commercial lithium-ion battery was analyzed.^11,19 In a previous study,²¹ the impedance spectrum was collected using a commercial lithium-ion battery (LiCoO₂-Ansmann 18650) at a 25 % state of change under 300 K. These data are bundled with DRTtools,³⁴ which were analyzed using our program. The data analyzed in terms of this spectrum were compared with the results obtained from DRTtools and EIS-DRT.

4. Results and Discussion

4.1 Simulated spectrum affected by the capacitor and inductive loop

Figure 3a illustrates the fitting of the simulated spectrum, wherein the proposed extended-DRT program reproduced the simulated spectrum using both the Tikhonov and Lasso regularization methods. Additionally, the figure indicates that DRTtools and EIS-DRT failed to reproduce the target spectrum because they do not support a spectrum influenced by an isolated capacitor and inductive loop.

Figure 3.

Results obtained from the proposed extended-DRT that are fitted using Tikhonov and Lasso regularization methods. (a) Impedance plot and (b) DRT curves. The impedance data used for the analysis were calculated using the circuit illustrated in Fig. 2. The fitted results in the impedance plot in (a) overlap completely. The arrow marks in (b) represent the satellite peaks that appeared during Tikhonov regularization by applying the proposed extended-DRT. The results calculated using DRTtools and EIS-DRT are plotted for comparison. These calculations used λ = 10⁻⁴, and the range between log τ_min and log τ_max was identical to the frequency range.

Despite adequate concurrence with the fitted impedance spectrum irrespective of regularization methods used in the proposed extended-DRT program (Fig. 3a), the DRT curves exhibit large differences (Fig. 3b). This can be attributed to the influence of the penalty function in Eqs. 14 and 15. In comparison with the true answer, the DRT peaks were broader in the case of Tikhonov regularization. Additionally, several small satellite peaks were observed around the large peaks in the case of Tikhonov regularization in the proposed program, as depicted by the arrow marks in Fig. 3b. However, these satellite peaks do not appear clearly in the case of Lasso regularization. Almost no satellite peaks appear in the results of DRTtools and EIS-DRT owing to the constraint, positive DRT values. Therefore, the appearance of the satellite peaks in the extended-DRT program was not caused by program bugs.

Table 1 lists the additional optimized parameters, namely, R_∞, L₀, and C₀. The R_∞ values output by the extended-DRT program differ from the true answer in Fig. 3. This difference can be attributed to the influence of the negative resistance evaluation that reproduces an inductive loop resistance, as explained in Eqs. 5 and 6. The R₃ value in Fig. 2 can be calculated using the absolute DRT area of the negative peak in Fig. 3b. Additionally, the corrected R_∞ was calculated as R_∞ − R₃, which was evaluated considering the area. The deviation between the true R₃ value (5 Ω) and evaluated value is larger in the case of Tikhonov regularization (Table 1), which is caused by the overlap of R₄ on the DRT and satellite peaks.

Table 1. List of additional optimized parameters, calculated inductive loop resistance, and corrected R_∞ in Fig. 3.

Programs	Proposed program		DRTtools	EIS-DRT
Regularization method	Tikhonov	Lasso	Tikhonov	Tikhonov
R∞/Ω	10.1	10	11.6	17.9
L0/µH	1.0	1.0	0.7	0.5
C5/F	1.0	1.0	NA	NA
R3 (Area of negative peak)/Ω	6.1	5.2	NA	NA
R∞ − R3/Ω	4	4.8	NA	NA

In the case of the simulated spectrum, true DRTs are combinations of impulse signals (Fig. 3b) when the resistor in the inductive loop represents the negative value. The DRT peaks in the DRT analyzed using our program with Tikhonov regularization are clearly broad. Conversely, the DRT in the case of Lasso regularization is close to the answer. Based on these results, we concluded that the sparse ability of Lasso regularization is stronger than that of Tikhonov regularization.

4.2 Measured spectrum

In the case of the lithium battery spectrum,²¹ the impedance spectrum can be accurately reproduced using the proposed extended-DRT program under Tikhonov and Lasso regularizations (Fig. 4a). In the case of this spectrum, expanding the ln τ_max range was preferred over data frequency range expansion. This is because this spectrum was terminated in the middle of the depressed arc around the lower limit frequency (5 mHz).

Figure 4.

Results obtained from the proposed extended-DRT that are fitted using Tikhonov and Lasso regularization methods based on the measured spectrum of lithium battery.^20,32 (a) Impedance plot and (b) DRT curves. The fitted results in the impedance plot in (a) obtained from our extended-DRT program overlap completely. The results calculated using DRTtools and EIS-DRT are plotted for comparison. Based on the enlarged plot of Z_real = 0.112–0.120 in (a), certain deviations of results from the data are observed in the case of DRTtools and EIS-DRT. These calculations considered λ = 10⁻⁴, and the τ_max was set to 1/(f_min/10), where f_min denotes the minimum frequency of the measured impedance spectrum (5 mHz).

This spectrum can be reproduced either by adding a capacitor at low-frequency termination or without adding the capacitor. This is because extremely large capacitance values are optimized by the extended-DRT, as listed in Table 2. In the case of the frequency range of the data spectrum, a substantially large capacitor has almost no influence on the spectrum. Therefore, capacitor termination was considered unnecessary in the case of these data.

Table 2. List of additional optimized parameters in Fig. 4.

Programs	Proposed program		DRTtools	EIS-DRT
Regularization method	Tikhonov	Lasso	Tikhonov	Tikhonov
R∞/Ω	0.11	0.11	0.11	0.11
L0/µH	0.75	0.70	0.56	0.63
C5/F	7975	4686	NA	NA

In comparison with the results obtained from DRTtools and EIS-DRT, the degree of fit to the data spectrum in the proposed extended-DRT was better in the high-frequency range, as illustrated in the enlarged plot in Fig. 4a. Additionally, a large deviation was observed in the case of EIS-DRT in the low-frequency range (Fig. 4a). In the case of the DRT curves (Fig. 4b), differences were observed around log τ ≈ 1 to 2 and −3 to −2.5. Based on the detailed degree of fit to the data spectrum (Fig. 4a), our program appeared better than the other programs.

4.3 Existing problems of the DRT analysis

We determined that the DRT curves exhibit large differences based on the program used, as depicted in Figs. 3 and 4. Additionally, the DRT curves change depending on the regularization method. Based on these results, a simple comparison of the DRT curves is difficult among those analyzed using different programs and regularization methods. This aspect is essential when certain tendencies are analyzed by machine learning using several DRT curves. However, this is not a particular disadvantage of the DRT method because a similar tendency exists in conventional equivalent circuit analysis.

As the DRT curve indicates the secondary information, it cannot be considered as the primary data. This is because the DRT curve is calculated using a model, as indicated in Eqs. 2 or 3, although the users do not set the model. Furthermore, the regularization method is selected based on the users’ objective, such as whether they require a smooth curve or highly sparse curve. This implies that the DRT result is influenced by the user’s objective, which invalidates the belief that the DRT result does not involve an arbitrary choice made by users.

If users intend to exhibit a change in the impedance spectra, plotting Z_imag versus log f is preferable. This is because the change in spectra should be reflected in the Z_imag − log f plot if changes occur in the DRT curves. The reliability of the change in Z_imag − log f is higher than that observed in the DRT curve because the Z_imag − log f data serve as the primary data. Considering that satellite peaks appear in the DRT spectra, the changes in DRT curves must be analyzed to understand their influence.

However, almost no satellite peaks were observed (Fig. 3) in the case of Lasso regularization because the DRT values converged to a small value in the case of Lasso regularization owing to its highly sparse ability compared with the Tikhonov regularization. During the impedance analysis, continuous and smooth DRT curves are not always the correct answers, as explained in Section 4.1. Conversely, the DRT curves become smooth in the case of Tikhonov regularization. This indicates that the deconvolution resolution of Tikhonov regularization is lower than that of Lasso regularization. However, it is difficult to objectively determine the preferred result between Tikhonov and Lasso regularizations. This is because both DRT results can reproduce the impedance spectrum data, as depicted in Figs. 3a and 4a, owing to the different rules of optimized state introduced using the penalty function in the case of regularized least-squares.^36–38 Furthermore, different treatments on probability and statistics are included depending on the penalty function.⁴⁴

DRT curves exhibit large dependencies in terms of the regularization methods and programs used. Therefore, the DRT curves can be used only as a basis for constructing equivalent circuit models. A detailed discussion requires the relative change of the fitted parameters based on the equivalent circuit analysis. This is because the analysis program dependencies of the equivalent circuit model are weaker than those of the DRT curves if the same equivalent circuit model is used. Therefore, overconfidence in the electrochemical DRT method can be avoided with a clear understanding of the basics.

The DRT method can be summarized as follows. The DRT analysis is performed by optimizing the model based on Eqs. 3 or 4, which include numerous parameters. The optimized answer is unstable owing to the presence of numerous parameters. Although regularized CNLS is employed to avoid the instability of the answer, the answer becomes ambiguous. Instead, the answer should be used as an important hint for the construction of an equivalent circuit model.⁵

5. Conclusions

In this study, we expanded the DRT analysis by introducing a capacitor term and negative DRT for analyzing the impedance spectrum. Additionally, its influence on the capacitor and inductive loop was investigated. We considered the basic property of regularized nonlinear least-squares and simple alternative expression of resistors connected in parallel to inductor elements using resistors and capacitors with negative values. Furthermore, we successfully implemented an iterative elastic net regularized CNLS using the Levenberg–Marquardt algorithm.

We compared the DRT results obtained from the proposed extended-DRT program with those generated by other open-source programs, namely, DRTtools and EIS-DRT. The results verify that our program reproduces both the simulated and measured spectra adequately. Moreover, we determined that the DRT results rely on the programs and regularization methods used. Therefore, we invalidated the belief that DRT results do not involve arbitrary choices made by users. Based on this analysis, we concluded that it is better to use the DRT results as a basis to set an equivalent circuit model. As DRT results are not primary data and are secondarily analyzed results that serve as secondary information, it is better to understand the DRT basics to avoid overconfidence in the DRT results.

Meanwhile, two problems have been identified in our program. One is the high calculation load caused by the many iterations of the Levenberg–Marquardt method. The second is the use of expensive commercial software (Igor Pro). To use our method widely, it is necessary to develop another open-source program.

Acknowledgments

This work was supported in part by JSPS KAKENHI Grant Numbers 17H01317 and 19H02804, and JST, CREST Grant Number JPMJCR1996, Japan.

CRediT Authorship Contribution Statement

Kiyoshi Kobayashi: Conceptualization (Equal)

Tohru S. Suzuki: Conceptualization (Equal)

Conflict of Interest

The authors declare no conflict of interest in the manuscript.

Funding

Japan Society for the Promotion of Science: 17H01317

Japan Society for the Promotion of Science: 19H02804

Japan Science and Technology Agency: JPMJCR1996

Footnotes

T. S. Suzuki: Equal Contribution

K. Kobayashi: ECSJ Active Member

References

• 1)  T. Paul, P. W. Chi, P. M. Wu, and M. K. Wu, Sci. Rep., 11, 12624 (2021).
• 2)  E. Quattrocchi, T. H. Wan, A. Belotti, D. Kim, S. Pepe, S. V. Kalinin, M. Ahmadi, and F. Ciucci, Electrochim. Acta, 392, 139010 (2021).
• 3)  J. Liu, T. H. Wan, and F. Ciucci, Electrochim. Acta, 357, 136864 (2020).
• 4)  J. Liu and F. Ciucci, Electrochim. Acta, 331, 135316 (2020).
• 5)  S. Dierickx, A. Weber, and E. Ivers-Tiffée, Electrochim. Acta, 355, 136764 (2020).
• 6)  F. Ciucci, J. Electrochem. Soc., 167, 126503 (2020).
• 7)  S. V. Pereverzev, S. G. Solodky, V. B. Vasylyk, and M. Žic, Comput. Methods Appl. Math., 20, 517 (2020).
• 8)  M. Žic, S. Pereverzyev, Jr., V. Subotić, and S. Pereverzyev, GEM Intnl. Geomath., 11, 2 (2020).
• 9)  N. Schlüter, S. Ernst, and U. Schröder, ChemElectroChem, 6, 6027 (2019).
• 10)  E. Quattrocchi, T. H. Wan, A. Curcio, S. Pepe, M. B. Effat, and F. Ciucci, Electrochim. Acta, 324, 134853 (2019).
• 11)  X. Li, M. Ahmadi, L. Collins, and S. V. Kalinin, Electrochim. Acta, 313, 570 (2019).
• 12)  K. Kobayashi and T. S. Suzuki, J. Phys. Soc. Jpn., 87, 094002 (2018).
• 13)  A. Mertens and J. Granwehr, J. Energy Storage, 13, 401 (2017).
• 14)  E. Ivers-Tiffée and A. Weber, J. Ceram. Soc. Jpn., 125, 193 (2017).
• 15)  A. L. Gavrilyuk, D. A. Osinkin, and D. I. Bronin, Russ. J. Electrochem., 53, 575 (2017).
• 16)  M. B. Effat and F. Ciucci, Electrochim. Acta, 247, 1117 (2017).
• 17)  Y. Zhang, Y. Chen, M. Li, M. Yan, M. Ni, and C. Xia, J. Power Sources, 308, 1 (2016).
• 18)  Y. Zhang, Y. Chen, M. Yan, and F. Chen, J. Power Sources, 283, 464 (2015).
• 19)  T. H. Wan, M. Saccoccio, C. Chen, and F. Ciucci, Electrochim. Acta, 184, 483 (2015).
• 20)  F. Ciucci and C. Chen, Electrochim. Acta, 167, 439 (2015).
• 21)  M. Saccoccio, T. H. Wan, C. Chen, and F. Ciucci, Electrochim. Acta, 147, 470 (2014).
• 22)  P. Büschel, T. Günther, and O. Kanoun, 2014 IEEE Intern. Instr. Meas. Tech. Conf. Proceedings (IEEE Instr. Meas. Soc., Uruguai), pp. 901–904 (2014).
• 23)  D. Klotz, J. P. Schmidt, A. Kromp, A. Weber, and E. Ivers-Tiffée, ECS Trans., 41, 25 (2012).
• 24)  S. Hershkovitz, S. Tomer, S. Baltianski, and Y. Tsur, ECS Trans., 33, 67 (2011).
• 25)  A. B. Tesler, D. R. Lewin, S. Baltianski, and Y. Tsur, J. Electroceram., 24, 245 (2010).
• 26)  A. Leonide, V. Sonn, A. Weber, and E. Ivers-Tiffée, J. Electrochem. Soc., 155, B36 (2008).
• 27)  F. Dion and A. Lasia, J. Electroanal. Chem., 475, 28 (1999).
• 28)  T. Hörlin, Solid State Ionics, 107, 241 (1998).
• 29)  T. Hörlin, Solid State Ionics, 67, 85 (1993).
• 30)  J. Weese, Comput. Phys. Commun., 69, 99 (1992).
• 31)  C. Elster, J. Honerkamp, and J. Weese, Rheol. Acta, 31, 161 (1992).
• 32)  H. Schichlein, A. C. Müller, M. Voigts, A. Krügel, and E. Ivers-Tiffée, J. Appl. Electrochem., 32, 875 (2002).
• 33)  E. Barsoukov and J. R. Macdonald, Impedance spectroscopy, Wiley, New Jersey, Chaps. 1, 3, and 4 (2018).
• 34)  T. H. Wang, DRTtools, https://sites.google.com/site/drttools/ referred on 20 October (2021).
• 35)  X. Li, EIS-DRT, https://github.com/nonmin/EIS-DRT referred on 20 October (2021).
• 36)  Y. Murakami, Gyakumondai no Kangaekata, Kodansha, Tokyo, Chap. 7 (2014). [in Japanese]
• 37)  H. Zou and T. Hastie, J. R. Stat. Soc. B, 67, 301 (2005).
• 38)  T. Tada, Atarashii Jinkochino no Kyokasho, Shoeisha, Tokyo, Chap. 4 (2016). [in Japanese]
• 39)  T. Mori, Gendai no Kotenkaiseki, Chikuma Gakugei Bunko, Tokyo, Chaps. 8, 11, 13, and 21 (2006). [in Japanese]
• 40)  B. A. Boukamp and A. Rolle, Solid State Ionics, 314, 103 (2018).
• 41)  K. Kobayashi, Y. Sakka, and T. S. Suzuki, J. Ceram. Soc. Jpn., 124, 943 (2016).
• 42)  K. Kobayashi and T. S. Suzuki, Kagaku Kogyo, 72, 556 (2021). [in Japanese]
• 43)  J. W. Eaton, GNU Octave, Scientific Programming Language, https://www.gnu.org/software/octave/index referred on 20 October (2021).
• 44)  S. B. McGrayne, The Theory That Would Not Die: How Bayes’ Rule Cracked the Enigma Code, Hunted Down Russian Submarines, and Emerged Triumphant from Two Centuries Controversy, Yale Univ. Press, Chaps. 1 and 3: (Japanese translation) Itan no Toukeigaku Bayes, Soshiha, Tokyo, Chaps. 1 and 3 (2011). [in Japanese]
• 45)  E. R. v. Schweidler, Ann. Phys., 329, 711 (1907).
• 46)  C. J. F. Böttcher and P. Bordewijk, Theory of Electric Polarization. Vol. II, Dielectrics in Time-Dependent Fields, 2nd completely revised ed., Elsevier, Amsterdam, Chap. 8 (1996).
• 47)  K. W. Wagner, Ann. Phys., 345, 817 (1913).
• 48)  D. D. Macdonald, Electrochim. Acta, 51, 1376 (2006).

Appendices

Appendix A. Theoretical Basis for Conventional DRT Equation and Expansion in Terms of the Inductive Loop Spectrum

The original interpretation of the dielectric spectrum^45,46 was translated into an impedance spectrum. When the influence of the inductor is negligible, the non-ideal impedance spectrum can be reproduced from the mathematical equivalent of the dielectric relationship using Eq. A1.^12,45,46   

\begin{equation} Z(j\omega) = R_{\infty} + \sum\nolimits_{k}\frac{R_{k}}{1 + \tau_{k}j\omega} \end{equation} (A1)

where R_k and τ_k denote the k-th resistance and relaxation time, respectively. In this case, τ_k is obtained using the k-th capacitance (C_k) as   

\begin{equation} \tau_{k} = R_{k}C_{k}. \end{equation} (A2)

Figure 1c illustrates a part of the equivalent circuit of Eq. A1. However, Eq. A1 cannot be transformed into another equation analytically, the following expression was introduced to be “sufficiently general to permit an adequate description.”^46,47   

\begin{equation} Z(j\omega) = R_{\infty} + \int_{0}^{\infty}\frac{R(\tau)}{1 + \tau j\omega} \mathrm{d}\tau. \end{equation} (A3)

where R(τ) denotes the relaxation time on a linear timescale. In the case of impedance measurements, the frequency range is several orders wide. Therefore, the use of a logarithmic timescale is convenient, and Eq. A3 can be transformed using the relaxation time in logarithmic timescale [R(ln τ)] as   

\begin{equation} Z(j\omega) = R_{\infty} + \int_{-\infty}^{\infty}\frac{R(\ln \tau)}{1 + \tau j\omega} \mathrm{d}\ln \tau. \end{equation} (A4)

If a non-ideal shape spectrum exists in an inductive loop, its spectrum can be reproduced under the same analogy as indicated in Eq. A1, as follows:   

\begin{equation} Z(j\omega) = R_{\infty} + \sum\nolimits_{l}\frac{R_{l}\tau'_{l}j\omega}{1 + \tau'_{l} j\omega} \end{equation} (A5)

where R_l and $\tau '_{l}$ denote the l-th resistor and relaxation time, respectively, and $\tau '_{l}$ is calculated using the l-th inductor (L_l) as   

\begin{equation} \tau'_{l} = L_{l}/R_{l}. \end{equation} (A6)

A part of the corresponding equivalent circuit is depicted in Fig. 1b. The same analogy can be applied to transform Eq. A6, as follows:   

\begin{equation} Z(j\omega) = R_{\infty} + \int_{-\infty}^{\infty}\frac{\gamma (\ln \tau')\tau' j\omega}{1 + \tau' j\omega} \mathrm{d}\ln \tau', \end{equation} (A7)

where γ(log τ′) denotes the relaxation time in the logarithmic timescale for the inductive loop.

Appendix B. Problems of Continuous Relaxation Time and DRT

Mathematically, Eqs. A3, A4, and A7 are the Fredholm integral equations of the first kind,³⁶ and those include generalized functions such as the Dirac delta function. The transformations from Eqs. A1 and A5 to A3, A4, and A7 are established when τ and τ′ are the sets of discretized impulse signals.³⁹ If the number of elements in the sets is sufficiently large, continuous DRT becomes a valid model. Contrary, continuous DRT is derived analytically from Eq. A4 with the implicit assumption that τ and τ′ are continuous integer numbers, which can cause a slight excess expansion for the impedance spectrum without accurate validation of the assumption.

Based on impedance theory,⁴⁸ the free variable of Z is jω, and τ_k and R_k in Eq. A1 and τ, R(τ), and R(ln τ) in Eqs. A3 and A4 are the parameters. This property has already been established indirectly³² because the authors³² explained that the relaxation time and DRT can be calculated analytically and not measured directly. Although τ and ln τ appear to be treated as free variables in Eqs. A3 and A4, the parametric property is maintained in these equations owing to the fixed interval of the integral. It is clear that τ and/or ln τ are eliminated after the definite integral from 0 to ∞ in the case of Eq. A3, and from −∞ to ∞ for Eq. A4. This indicates that the parametric property is maintained owing to mapping by integral calculus from a function to another function under a fixed interval.³⁹ The free variable of one function must differ from that of another function based on the mapping.

The integrand in Eq. A4 can be expressed as I and presented using its free variable and parameters.   

\begin{equation} I(\ln \tau; R(\ln \tau), j\omega) = \frac{R (\ln \tau)}{1+\exp (\ln \tau)j\omega}. \end{equation} (B1)

In function I, the free variable is implicitly set by ln τ, and R(ln τ) and jω are treated as the parameters. However, the integral term in Eq. A4 is presented using Z_{ln τ} similar to I, and this function is expressed as   

\begin{equation} Z_{\ln \tau} (j\omega; \ln \tau, R (\ln \tau)) = \int_{-\infty}^{\infty} \frac{R(\ln \tau)}{1 + \exp (\ln \tau)j\omega} \mathrm{d}\ln \tau. \end{equation} (B2)

Herein, the free variable is identical to that of Z, which is jω, and the parameters are ln τ and R(ln τ). This indicates that the free variable is exchanged between ln τ and jω by integral calculus within the interval −∞ to ∞.

Considering that ln τ and R(ln τ) are the parameters of a model, the belief that continuous ln τ and R(ln τ) are natural would be weakened. Therefore, the validity of continuous DRT in impedance spectroscopy remains ambiguous.