Computation of distribution of relaxation times by Tikhonov regularization for Li ion batteries: usage of L-curve method

In this paper, the distribution of relaxation times (DRTs) functions are calculated numerically in Matlab for synthetic impedance data from single parallel \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$RC$$\end{document}RC circuit and two parallel \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$RC$$\end{document}RC circuits connected in series, experimental impedance data from supercapacitors and α-LiFeO2 anode based Li ion batteries. The quality of the impedance data is checked with the Kramers–Krönig (KK) relations. The DRTs are calculated within the KK compatible regime for all the systems using Tikhonov regularization (TR) method. Here we use a fast and simple L-curve method to estimate the TR parameter (λ) for regularization of the Fredholm integral equations of first kind in impedance. Estimation of the regularization parameters are performed effectively from the offset of the global corner of the L-curve rather than simply using the global corner. The physical significances of DRT peaks are also discussed by calculating the effective resistances and capacitances coupled with peak fitting program. For instance, two peaks in the DRTs justify the electrical double layer capacitance and ion diffusion phenomena for supercapacitors in low to intermediate frequencies respectively. Moreover, the surface film effect, Li/electrolyte and electrode/electrolyte charge transfer related processes are identified for α-LiFeO2 anode based Li-ion batteries. This estimation of the offset of the global corner extends the L-curve approach coupled with the Tikhonov regularization in the field of electrochemistry and can also be applied in similar process detection methods.

Electrochemical impedance spectroscopy (EIS) is useful to separate contributions of different electrochemical phenomena related to polarization losses, basic ion transport and kinetic parameters of various electrochemical cells, and batteries [1][2][3][4][5][6] . EIS is often analyzed with complex non-linear least squares fit, which is computed by equivalent circuit modelling. However, equivalent circuits are not easy to find and it is generally difficult to have proper circuit for complex electrochemical phenomena.
It is of utmost interest of researchers to find alternative approach to interpret the EIS data. Distribution of relaxation times (DRT) provide a rich analysis of EIS by transforming the data from frequency ( f ) domain to time (τ ) domain. DRTs exhibit peaks or single peak on a log (τ ) axis (or f axis) indicating some electrochemical processes. The primary advantage of DRTs is its representation being model free and do not need much information about the system. Up to now several methods have been adopted to invert impedance data to a DRT such as Tikhonov Regularization (TR) [7][8][9][10][11][12][13][14][15][16] , Maximum Entropy (ME) 17 , Fourier transform (FT) [18][19][20] . Concerning TR, ridge and Lasso methods need an adjustment of the regularization parameters to get reasonable DRTs and different values lead to different number of peaks hindering true electrochemical phenomena 21 . The DRTs of yttria and calcium stabilized cubic zirconia have been computed by choosing regularization parameter based on root mean square error of impedance 22 as well as by least square method coupled with Tikhonov regularization 15 . Combination of Tikhonov regularization and multi-(RQ) CNLS-fits technique is used for solid oxide fuel cells and steam electrolysers to determine the analytic form of DRTs 18 . DRTs based on ME method has been computed for NASICON materials 17 . Both ionic and electronic contribution are determined from DRTs of solid oxide fuel cell cathodes by Fourier transformation of the EIS data 23 . Least Absolute Shrinkage and Selection Operator 21 based regularization and multi-(RQ) CNLS-fits 24,25 are used for commercial Li-ion batteries for EIS data synergistically to compute the DRTs. On the other hand, another approach has been made by the group of Tsur 26-28 by constructing a parametric form of DRTs. The advantage of this method is, it does not need any smoothing parameter or window function. DRTs based on this approach have been computed for oxide ion conducting electrolytes 29 of solid oxide fuel cells, ceramic actuators 30 as well as Fe and Ni based catalysts 31 .
The goal of this work is to provide a simple and fast method to calculate DRTs using Matlab application with valid regularization parameter which is determined by L-curve. The offset of the global corner of the L-curve method leads to better estimate results for all the studied examples. To demonstrate this new method, we have worked on the synthetic impedance data from single unit of RC circuit as well as series connection of RC circuits. DRTs are also calculated for graphite-based supercapacitors and LiFeO 2 anode based Li-ion batteries. A Matlab based curve fitting program is also used to estimate the R and C 's from the DRTs quantitatively.

Method of simulation
The data inversion in calculating DRTs needs to solve a Fredholm integral equation of first kind supposing the impedance spectrum as Z * (ω) = Z ′ (ω) − jZ ′′ (ω) , where Z * (ω) is the complex impedance and Z ′ and Z ′′ are the real and imaginary parts of complex impedance, respectively. Equation (1) as per Ref. 32 presents the relation between impedance data and its frequency distribution, where R ∞ is the series resistance, R p is the total polarization resistance, G is the DRT and τ = RC is the relaxation time,R is the effective resistance and C is the effective capacitance, and ω is the angular frequency. Considering ω = 0 the normalization condition can be obtained as Decomposing both real and imaginary parts of the impedance, Eqs. (2) and (3) can be obtained and where F(τ ) = R p · G log (τ ) is the DRT spectrum. Finding F(τ ) using the above-mentioned equations is known as ill-posed inverse problem as in principle, many solutions are possible. The usage of TR is widely accepted for such problems. As solution of the right-hand sides for Eq. (2) and (3), construction of linear systems of equations is performed by discrete mesh { τ n , n = 1, . . . , N }, where N is the total number of f points. Presently, we have used τ n τ n , = 1 ω n and ω n = τ n+1 − τ n . Furthermore, we can modify Eqs. where h τ = dlog(τ ) . In terms of matrix notation For physically acceptable solution a non-negative constraint x n ≥ 0∀n is applied. The b matrix is constructed by both real and imaginary parts of given impedance data such as min xǫR n �Ax − b� 2 + 2 �x� 2 , www.nature.com/scientificreports/ L-curve is a popular graphical method for determining the suitable regularization parameter. L-curve 33 shows the negotiation between the two quantities as described in Eq. (6) and is used to find the regularization parameter. Secondly, this technique is fast and straightforward. Thirdly, computing the distribution of relaxation times is an ill-posed inverse problem and it is of utmost necessity to use some tool to reduce this deficiency. Presently, we use L-curve to find the global corner (with maximum curvature) to the residual norm. Since this does not need the information of the solution norm, provides the necessary information of the residual norm. However, offsets of the global corner are more efficient in estimating the regularization parameters as discussed in later section. The package for determining the DRTs is developed in Matlab environment and for single value decomposition Matlab based script 33 is used. The overall code takes one to five seconds to calculate a DRT spectrum on a 1.99 GHz PC depending on the number of frequency points.
The quantification and isolation of each peak from the DRTs are done by peak fitting program assuming Gaussian distribution. A peak fitting that justifies the central position with the characteristic time constant τ c ; the height, an indication of the polarization; the standard deviation, which indicates the width of the peak. Presently, the peak fitting is performed using a Matlab based script peakfit 34 . Time consumes in peak fitting also as several iterations are necessary. Overall, the polarization contribution or the R at a characteristic time constant is equal to the normalized area under each peak multiplied by R p . Then the effective capacitance ( C eff ) is computed by τ c R relation. Thus, we are able not only in identifying different processes but also R and C eff 's corresponding to these processes from the peaks.

Results and discussion
In this context, impedance spectra of different electrochemical systems are analysed using the TR technique. Before processing the impedance data the quality of data has been checked by conventional KK relations using Lin-KK tool 35 .
Single RC circuit. A simulated impedance arc along with the computed arc having R = 50 Ω and C = 10 × 10 −6 F are shown in Fig. 1a between 1 MHz to 1 mHz. It has been observed that the KK relations are valid in the entire frequency region (Fig. S1). The computed DRT spectrum as plotted as a function of frequency ( f = 1 2πτ ) has a single broad peak, implying single electrochemical process (Fig. 1b). Considering both the global corner (4.44 × 10 −5 ) and the offset of the global corner (9 × 10 -3 ), smoother variation for the latter in the DRT especially at high frequencies is achieved (Fig. 1b). It has been observed that the global corner is not providing proper choice of a regularization parameter (Fig. S2). Similar failure of the global corner in estimation of the regularization parameter is also discussed for ill-posed inverse problem of electrical resistance tomography 36 . Nevertheless, the R and C eff are obtained as 50.997 Ω and 1.70 × 10 -6 F respectively from the area calculation of the peak fit considering the offset of the global corner. www.nature.com/scientificreports/ Double RC circuits. The simulated impedance arc for two parallel R and C 's as connected in series having resistances and capacitances as 1 kΩ, 100 Ω and 10 × 10 -6 F, 1 × 10 -6 F respectively is plotted in Fig. 2a between 1 and 1 MHz. As obtained from Fig. S3, the residual errors using KK relations is null and the entire frequency range can be used for DRT calculations. Figure 2b shows the DRT spectrum corresponding to the computed impedance plot. Although global corner is obtained using the L-curve method (Fig. S4), still the DRT becomes noisy at high frequencies restricting precise estimation (Fig. S5). Further peak profile analysis leads to R and C 's 998.82 Ω, 78.29 Ω and 3.55 × 10 -6 F, 20.86 × 10 -9 F respectively with = 0.04 (offset of the global corner). However, we have obtained an onset of a weak pear around the marked region which is vanished when we simulate the DRT for more closed semicircle at the high frequencies (viz. 50 MHz).

Supercapacitors. Graphite based supercapacitors have been studied electrochemically at different bias
voltages with frequency range of 10 mHz to 10 kHz with an AC excitation voltage of 10 mV after 4000 cycles. To obtain a shift in the impedance plot different bias voltages are applied. It is well known that the typical Nyquist plot of impedance for supercapacitors consist of single semicircles at high to medium high frequency range but a spike like nature at low frequencies (Fig. S6). As proposed by De Levie 37 electrical double layer capacitance (EDLC) forms towards the porous electrode at lowest frequencies and ion diffusion occurs at intermediate frequencies 38 . For the sake of completeness, the complex impedance mode has been transformed to complex capacitance mode via the relation 39,40 C(ω) = 1 jωZ(ω) , where Z(ω) is the complex impedance and C(ω) is the complex capacitance. Here C ′ (ω) and C ′′ (ω) are the real and imaginary parts of the complex capacitance respectively (Figs. 3a,b). This transformed data show KK compatibility are within 2% as shown by the residual error plots (Figs. S7, S8). Interestingly, noise free DRTs are obtained using the regularization parameters which are obtained from the offset of the global corner of L-curve (see Figs. S9, S10). Although the impedance arcs show similar behaviour but the DRT peaks shift towards high frequency side indicating effect of charge relaxations due to increasing bias voltages (Figs. 3c,d). Presently, the total area of each peak is multiplied by the normalization factor to obtain C eff . The DRTs are comprised of two broad peaks around 10 Hz and 1 Hz with capacitances 0.06 F, 0.57 F, and 0.08 F, 0.57 F, for 1 V and 2 V respectively. The peaks around 10 Hz and 1 Hz can be associated with diffusion of electrolyte ions towards and inside the pores of the electrodes and EDLC respectively 40 . The diffusion process consumes maximum capacitance due to local polarizations 40 .
Li-ion battery. Lithium ion batteries (LIB) show complex nature in their impedance spectra due to different electrochemical processes. In the half cell configuration of LIBs, impedance response consists of contributions from anode with its solid electrolyte interphase, the cathode and electrolyte having polarization processes like ohmic losses, charge transfer processes, as well as diffusion processes 41   www.nature.com/scientificreports/ ured at room temperature with a minimum frequency of 10 mHz up to maximum frequency of 100 kHz with 71 measured frequencies in total. The AC excitation voltage was set as 10 mV, and the impedance data have been obtained after 1st, 5th and 10th cycles of cyclic voltammetry measurements. The total impedance increases with increasing cycling as observed in Figs. 4a,b. Since we have observed a spike like extension towards ω → 0 for Z ′′ (ω) vs f , (not shown) restricting the use of DRTs at low frequencies for LIBs. On the other hand, at present, the KK relations are valid up to 1 Hz for the impedance data, restricting the computation of the DRTs up to 1 Hz from high frequencies (Figs. S11, S12, residue errors are within 1%). This also explains why the fitting deviates significantly at low frequencies. As observed from Figs. S13 and S14 the global corner does not produce satisfactory DRT. Henceforth, we have used the offset of the global corner as a regularization parameter (λ = 0.0659) (Figs. S13, S15-S17). The calculated DRTs consist of three peaks as shown in Fig. 4c, Supplemental Fig. S18, and Fig. 4d after 1st, 5th and 10th cycles of cyclic voltammetry measurements respectively. The ohmic resistances are subtracted before the DRT calculations, thus the impedance spectra consist of polarization behaviour only. Furthermore, peak fitting analysis suggests that three peaks are located around 251 Hz, 2 Hz and 0.20 Hz named as P1, P2 and P3 respectively (Fig. 4c,d, Supplemental Fig. S19). With increasing cycles, the effective resistances for peaks P1 and P2 decrease while for P3 increases (see Table 1). Concerning the positions of the peaks and resistance values peaks P1, P2 and P3 correspond to surface film resistance as originated due to Li ion migration through the electrode surfaces, and α-LiFeO 2 /electrolyte interfacial charge transfer resistances respectively 42 . It is noted, DRT analysis has been investigated on LiFePO 4 /lithium cells 9 . The results show solid state diffusion in the cathode, and charge transfer resistance between cathode and electrolyte in between 10 3 and 1 Hz at different temperatures and state of charge 9 . The spike like extension at low frequency regime is avoided by pre-processing using equivalent circuit model in Ref. 9 . In our study, we have neglected this regime indirectly by choosing the KK compatible regime only. In this way, we have considered only the polarization behaviour and identified anode related phenomena as shown in Fig. 4c,d. As our interest lies within the KK compatible regime and for Li-ion batteries the electrochemical phenomena occur within the semicircular regime only. Considering these facts peak P1 and P2 are of utmost interest and since peak P3 is in the vicinity of the incompatible regime lower confidence level is obtained. As observed from Fig. S20, during the first discharge curve a sharp peak at 0.56 V attributes to the formation of solid electrolyte interface 43 as well as reduction of Fe 2+ and Fe 3+ to Fe 0 Ref. 44 . Secondly, during charging, a broad peak around 1.69 V corresponds to oxidation of Fe 0 to Fe 3+ Ref. 45 . To understand the effect of the lithiation and delithiation peaks on the ion transport additional EIS measurements were performed at 0.1 °C rate at different potentials and the results are shown in Fig. 5a,b. Distinguished features in the semicircular arcs as well as in the tails at low frequencies are noted for both charging and discharging potentials. Additionally, the DRTs www.nature.com/scientificreports/ are shown in Fig. 5c,d considering the offset of the global corner of the L-curve method (see Fig. S21). To better understand our observation, we have also determined the DRT of Li||Li symmetric cell as shown in Fig. S22. Comparing Fig. S22 with Fig. 5c,d, we can conclude that the peak S1 (between 10 3 and 10 Hz) is related to Li metal to electrolyte interface and continues till 2.5 V of charging as observed in Fig. 5d. Interestingly, during discharging, peak S1 decreases towards low frequency, whereas it shifts towards high frequency during charging exhibiting the effect of lithiation and delithiation respectively, which are also time bound processes. It can be concluded that these potentials activate the lithium related process but the direction of the potential sweep (OCV (1.16 V) → 0.01 V → 3 V) corresponds to the shifting of the peak position. An increasing trend in their effective resistances exhibits the aging effect as displayed in Fig. 5e. In addition, the peaks S2 and S3 occurring between 10 and 0.1 Hz correspond to Li ion motion for the α-LiFeO 2 /electrolyte interfacial charge transfer process, and the solid state Li ion diffusion respectively. Similar to the previous description, peak S3 is in the KK incompatible regime so that we do not carry out further analysis. The lithiation at 0.46 V (Fig. S20) has strong effect on peak S2 both in position as well as in peak shape. For example, it is observed that initially at OCV and 1 V of discharging the S2 is separable from S3 and not visible beyond 2 V during charging of the battery (Fig. 5 (c) and (d)). Since S2 evolves from the electrode/electrolyte interface as well as direction dependent of the potential sweep (OCV → 0.01 V → 3 V), it is expected to be related to the reduction and oxidation of Fe respectively. The effective resistance related to peak S2 also shows increasing behaviour as shown in Fig. 5f. It is noted that  www.nature.com/scientificreports/ the reduction of Fe 3+ and Fe 2+ to Fe 0 is related to the initial capacity drop in ferrite based LIBs 44 . Therefore, our results provide clues to better understand not only the different electrochemical processes, but also the different transitions of Fe in the α-LiFeO 2 electrode.

Conclusions
In this paper, a Matlab based fast application is developed to transform the EIS data into non-parametric DRT arrangement. The code has been developed with the help of Tikhonov regularization with special constraint for eliminating non-negative parts in DRTs. The global corner of L-curve is failed to compute meaningful DRTs whereas the offset of the global corner is found to be useful. Nevertheless, we are successful in identifying different electrochemical processes from the DRTs by peak fitting methods quantitatively. For instance, electrical double layer capacitance and ion diffusion phenomena are identified for supercapacitors in low to intermediate frequencies respectively. Moreover, we have used the method to identify the three processes contributing to the half-cell LIBs. These processes are related to the charge transfer reactions and surface reaction processes. More