Optimized Process Parameters for a Reproducible Distribution of Relaxation Times Analysis of Electrochemical Systems

: The distribution of relaxation times (DRT) analysis offers a model-free approach for a detailed investigation of electrochemical impedance spectra. Typically, the calculation of the distribution function is an ill-posed problem requiring regularization methods which are strongly parameter-dependent. Before statements on measurement data can be made, a process parameter study is crucial for analyzing the impact of the individual parameters on the distribution function. The optimal regularization parameter is determined together with the number of discrete time constants. Furthermore, the regularization term is investigated with respect to its mathematical background. It is revealed that the algorithm and its handling of constraints and the optimization function signiﬁcantly determine the result of the DRT calculation. With optimized parameters, detailed information on the investigated system can be obtained. As an example of a complex impedance spectrum, a commercial Nickel–Manganese–Cobalt–Oxide (NMC) lithium-ion pouch cell is investigated. The DRT allows the investigation of the SOC dependency of the charge transfer reactions, solid electrolyte interphase (SEI) and the solid state diffusion of both anode and cathode. For the quantiﬁcation of the single polarization contributions, a peak analysis algorithm based on Gaussian distribution curves is presented and applied.

The EIS provides the frequency-dependent impedance of an electrochemical system by applying a sinusoidal current, measuring the voltage response (galvanostatic EIS) or vice versa (potentiostatic EIS), and applying Ohm's law. The impedance spectra are interpreted regarding their high-frequency inductive behavior [11] caused by geometry, wiring and the test setup, their low-frequency diffusion behavior [12] and mainly their interphase behavior. Numerous studies focus on battery impedance, e.g., correlating aging to changes in the impedance spectrum [13] or explaining the relaxation behavior in graphite anodes [14]. Analyzing the impedance without applying further techniques is complex, though. Effects overlap, processes are non-ideal, causing deviations between the expected and the measured spectra. Even the number of single processes cannot be obtained unambiguously. Therefore, the DRT became of interest in the 1990s [1] and was applied first for SOFCs [15,16]. Later on, the DRT was adapted for the characterization of PEMFC and batteries [2,7,17]. As the result of a mathematical transformation, the polarization contributions at predefined time constants of the investigated system are derived. Exemplary, in Figure 1a the impedance spectrum of two RC elements (resistor and capacitor in parallel) connected in series is shown: a deformed semi-circle with no distinct maximum of − Im(Z) results. This shape might be caused, e.g., by a non-ideal process due to large inhomogeneities within the electrode or by separate processes. Investigating the impedance spectrum, no further conclusions can be drawn. Applying the DRT algorithm as shown in Figure 1b, the two involved processes can be separated easily. Therefore, the contribution of single processes to the system's impedance can be calculated. The DRT offers a model-free, universal approach for the detailed analysis of electrochemical systems or any other data obtained from impedance measurement. This means that no a-priori assumptions, e.g., electrical, physical or chemical models or knowledge are necessary, which is an essential advantage when dealing with complex or hardly understood specimen such as fuel cell stacks or new materials with novel conducting mechanisms.
On the other hand, the DRT analysis strongly depends on the choice of its process parameters which are, in many publications, neither revealed nor discussed. Thus, the DRT algorithm containing those parameters is described in detail within this paper. Furthermore, an extensive parameter study is performed to evaluate the impact and correlation of those process parameters on the calculated distribution function. The study includes the optimization of the regularization parameter, the number of predefined discrete time constants as sampling points of the distribution function and the mathematical formulation of the regularization term. It is investigated whether real parts, imaginary parts of the complex impedance or both should be used for the calculation of the DRT. Finally, the impact of the optimization function and optimization algorithm is shown. As an instructive example of a complex impedance spectrum comprising inductive, resistive, capacitive and diffusive overlapping polarization contributions of multiple components within an electrochemical system, the algorithm is applied on impedance data from a commercial 3.3 Ah Nickel-Manganese-Cobalt-Oxide (NMC) lithium-ion pouch cell by Kokam and the results are discussed in detail. Furthermore, the features of the impedance spectra are assigned to the individual peaks observed in the distribution functions. For more insight into the single peaks within the DRT, a novel post-processing approach is presented, adapting modified Gaussian peaks to the distribution function resulting in a quantification of every single process within the spectrum.
By reporting all relevant information about the used algorithm, the results presented in this publication can be easily reproduced by other researchers. The information includes • the optimization algorithm, • the error function, • the type of data used, i.e., using real or imaginary parts of the complex impedance, • measurement parameters, e.g., current or frequency range, • regularization parameter, • number of time constants, • the pre-and post-processing routine, • the minimum and maximum time constants for the DRT.

Deriving DRT from EIS Data
The main idea of the DRT is to express the complex-valued, frequency-dependent impedance of an electrochemical system as an infinite series of differential RC-elements in series with an ohmic resistor R 0 [7,15,17]: The integral is often normalized by introducing which leads to The normalized Equation (3) has two advantages compared to Equation (1): • The total polarization resistance of the system gets apparent by R pol . • Systems with impedances in different orders of magnitude can be compared easily.
The term 1+jωτ corresponds to the relative differential contribution of a single ohmic-capacitive element as τ = RC and g(τ) = R /R pol results in R 1+jωRC , which is the well-known equation for an RC-element.
The function g(τ) is computed numerically. Therefore, the integral in Equation (3) has to be discretized yielding the finite sum The range of time constants [τ 1 , τ n ] has to be chosen a priori based on the available frequency range of the impedance data. Usually, the interval of time constant matches the interval of measured frequencies. The number of time constants n has to be large compared to the expected numbers of processes within the investigated system. The determination of the optimal number is part of our process parameter analysis in Section 3.
Equation (4) describes the resistive-capacitive behavior of the impedance spectrum as a sum of RC elements connected in series with predefined time constants. Neither assumptions are made regarding the number of processes involved nor model concepts on electrical, physical or chemical processes are required as in case of the model-based analysis of electrochemical impedance spectroscopy. Hence, the DRT analysis is regarded as a model-free method for the characterization of electrochemical systems. Furthermore, it has to be stated that due to the ohmic-capacitive nature of Equation (4), the DRT is not valid for inductive systems. Purely capacitive contributions cannot be described neither. Mathematically, these restrictions regarding the examined complex impedance can be expressed as

Dealing with Artefacts at the Boundaries of the Measured Frequency Range
If the measured impedance spectrum is not converging towards the real axis at the boundaries of the measured frequency range as demanded by Equations (5) and (6), as is it the case e.g., for diffusion processes at low frequencies, the DRT analysis will yield artefacts, i.e., singularities or diverging, steep slopes at the boundaries of the interval of time constants corresponding to the extremes of the measured frequencies. To deal with this issue we propose to extend the predefined vector of time constants beyond the minimum or maximum time constant corresponding to the according maximum and minimum frequency. Mathematically this is legitimate since the choice of the predefined time constants does not have to match the measured frequencies compulsorily as can be seen in the definition of the optimization problem later. This approach does not extrapolate the measured spectrum beyond the measured frequency range. Not even new assumptions are made on the behavior of the device under test beyond the measured frequency range. Instead, the optimization problem maps a wider vector of time constants onto the exactly same impedance information in frequency domain and therewith helps to better extract all information contained in the measured frequency spectrum. A back calculation of the identified DRT into frequency domain would be the same as fitting an RC circuit to the measured spectrum. As for fitting an equivalent circuit model to a measured impedance spectrum, the model is only valid for the measured frequency range. Since real systems will never have a limited bandwidth and a measurement system will never be able to measure an unlimited frequency range, we show in the results that broadening the predefined vector of time constants is beneficial for the analysis of electrochemical systems, especially if polarization processes are not abated at the boundaries of the measured frequencies as it the case for the solid state diffusive behavior of lithium-ion batteries.

Calculation of g(τ)
For the equations above, g(τ) has to be calculated from an impedance spectrum. Measured or model-based, computer-generated spectra are limited in their number of frequency points. Typically, the frequency vector of the treated impedance data is logarithmically scaled, with 6 [14], 7 [18], 10 [7,19] or 16 [20] steps per decade, respectively. Often, this information is missing in publications [2,21,22]. However, this relatively small amount of m data points (i.e., the measured data including frequency and impedance values) leads to the issue that there are usually less data points than sampling points for the time constants of the DRT analysis (m < n).
Therefore, the resulting optimization problem with the matrix A and the vector b as introduced below is ill-posed and might be ill-conditioned. Therefore, it cannot be solved analytically. One way to overcome the ill-posed optimization problem in Equation (8) is to apply numerical regularization. Among the regularization techniques commonly used for calculation of DRTs of electrochemical systems, Tikhonov regularization has proven to be a powerful technique [23,24]. The optimization function is extended by a regularization term ||λx|| 2 where λ is the regularization parameter. Its influence on the resulting distribution will be discussed in Section 3. The parameter x represents the discrete values of the distribution function g(τ).
The parameter λx can be extended to M · λx, where M is the regularization matrix. In the simplest case, M is defined as I, representing the n × n identity matrix.
In addition, the constraint has to be taken into account, since negative contributions to the impedance are physically not feasible. The solution of the problem stated in Equation (9) can be obtained by multiple numeric algorithms. As the optimization problem is ill-posed with no unique solution and as commonly used least-squares algorithms are only able to find local optima, it is crucial to be aware of the used optimization algorithm. Otherwise, the results of the DRT can hardly be reproduced. Ciucci et al. [25] and Wan et al. [26] give insight into the optimization problem including a detailed discussion on the regularization matrix. Nevertheless, to the authors' knowledge there is no publication which gives insight into the used optimization algorithm.
Considering Equation (4), the matrix A consists of A i,k in i-th line and k-th column where represent normalized RC-elements with a time constant τ k evaluated at the angular frequency ω i . The vector b consists of the m measured or generated impedance values. This may happen in various ways:

•
Using complex values: In this case, a solver is needed which is able to handle complex values as well as constraints.

•
Using the real parts only: Ignoring imaginary parts is legitimate due to the Kramers-Kronig relationship [27,28], as long as its conditions are fulfilled. It can be assumed that a properly measured spectrum complies with these. The linear Kramers-Kronig test is one option to prove the Kramers-Kronig validity of the examined impedance data.

•
Using the imaginary parts only: equivalent to Equation (13) substituting real parts by imaginary parts.
which provides more, in case of a noise-free measurement redundant, information. If not stated otherwise, this method is used within this publication.
With both real and imaginary parts used, a number of time constants n > 2 × m leads to an under-determined system. With n < 2 × m, the problem is overdetermined. Both cases lead to an ill-posed problem as there is no unique solution. If n = 2 × m, the system is well-posed. For all three cases, Tikhonov regularization is applicable.
A non-negative least squares (NNLS) solver is used as proposed by Lawson and Hanson ([29], p. 161). As this algorithm can only deal with matrix-vector-systems, Equation (9) is re-written according to Kaipio [30] and Wu [31]: The solver returns h(τ k ) which turns into g(τ k ) after normalization.

Pre-Processing of Measurement Data
As stated above, the DRT can only deal with resistive-capacitive systems. The frequencyindependent R 0 has to be subtracted from the impedance before running the DRT routine. In addition, electrochemical systems are never only resistive-capacitive. There are inductive contributions due to the measurement setup or the specimen itself [11] as well as diffusive behavior e.g., when investigating batteries. Therefore, it is critical that such influences are removed from the impedance data beforehand.
There are two proper ways to deal with inductive and capacitive branches: the first is to fit an appropriate equivalent circuit (EC) to the measurement data using the model-and-reduce approach. This EC may contain a resistor, an inductor, a capacitor, RL-elements (resistor and inductivity in parallel), RC-elements, ZARC-elements (Z ZARC = R 1+(jωRC) ϕ ), constant phase elements (CPE) (Z CPE = 1 1+(jω) ϕ C ) and Warburg-elements (planar and spherical diffusion) all connected in series. The contribution of non-resistive-capacitive EC elements is subtracted from the measured spectrum to fulfil the RC-constraint of the DRT (Equations (5)- (7)) This pre-processing step is in opposition to the advantage of not needing any a priori knowledge for the DRT at first glance as the EC fit needs to be precise. Only non-RC elements have to be identified and characterized precisely to be subtracted properly while the exact characterization of electrode and interface processes is not essential within this step. Imprecise EC fits lead to peaks without any physical meaning, offering space for misinterpretations. For an appropriate fit, Equations (5) and (6) should be fulfilled after subtracting non-inductive-resistive components.
Choosing an EC requires experience by the scientist. A by-inspection analysis of the impedance spectrum while keeping in mind the spectra of the specific elements is beneficial. The EC is then fitted using a trust-region method. The two-dimensional error function to be minimized is given by with Z mod representing the impedance function of the EC whose parameters are to be optimized. The second option, called cut-and-shift approach, does not need any a priori assumptions: first, all data points with Im(Z) > 0 are cut off the spectrum. As a second step, the ohmic offset is removed by subtracting R 0 = min{Re(Z)} from the spectrum.
The entire workflow is summarized in Figure 2. The validity of the measured spectrum is proven by Kramers-Kronig test: in accordance to literature, a spectrum is said to be valid as long as the residual between the measured and the Kramers-Kronig reconstructed impedance is below 1% for every data point. Within the pre-processing step non-resistive-capacitive components are removed. Finally, the DRT is calculated from the reduced spectrum and g(τ k ) is obtained.
Tikhonov regularization is not the only possible approach for calculating the DRT. Fourier transformation is used as well [15,17]. A comparison between regularization and Fourier transform is not part of this paper, though a study on benefits and disadvantages of either method might be interesting in the future.

Post-Processing of Result
In general, once g(τ k ) has been calculated, the distribution function can simply be interpreted in a phenomenological manner: the number of processes involved can be identified as well as changes e.g., during ageing or relaxation or differences between used materials, morphologies or cell designs. For a proper quantitative analysis, post-processing is required.
We quantify the polarization of each peak by fitting the numerical distribution to a sum of modified, non-normalized Gaussian probability density functions with adjustable height H, standard deviation σ and skew s using a logarithmic time constants scale. The introduction of the skew allows asymmetric distributions: Thus, the area under each specific peak corresponds to the impedance contribution of the underlying process. Furthermore, the standard deviation corresponds to the ideality or uniformity of a process, e.g., a narrow peak indicates a process with a small band width, e.g., a uniform particle size. Wide peaks point to distributed processes, e.g., a wide spread particle size distribution or a geometric influence of a large-scale specimen. The skew indicates the distribution of the non-ideal process, e.g., whether there are more small-sized or large particles within the electrode. Yet, an overall quantification or correlation to the system's attributes is not easily achievable, but a quantitative comparison between various DRTs can be achieved.
The entire algorithm, as described within this chapter, is available for download at http://www. ec-idea.uni-bayreuth.de.

Regularization Parameter and Number of Time Constants
For reliable and meaningful results, the regularization parameter λ and the number of time constants n have to be chosen adequately. In Equation (9), the condition λ > 0 forces x towards small ||x|| as a direct consequence of the minimization problem. Other publications have noted that an increase of λ results in a smoother shape of the DRT. This might lead to the risk that narrow peaks merge or small peaks disappear while small λ may cause oscillations or peaks with no physical meaning [7,32,33]. It is intuitive that small values of ||x|| lead to a smoother shape of g k and vice versa. Small ||x|| are achieved by a vector whose elements differ little, as outliers deteriorate the norm in the second order.
The expectations from literature are perfectly matched by the results shown in Figure 3. For small λ, oscillations occur for low time constants and around the second peak. If λ is chosen too large, the two peaks begin to merge and cannot be separated anymore. Furthermore, the amplitude decreases while the width of the peaks increases. The total area below the graphs does not change significantly. Several publications were dealing with the issue of optimizing the regularization parameter for their purposes. For this study, we investigated three approaches made by Saccoccio et al. [34] and Hansen [35]: using the so-called discrepancy method, the term is minimized by varying λ. It calculates the DRT using real parts of the impedance values within A and b returning x Re and using imaginary parts only resulting in x Im . The difference is meant to be a measure for the quality of the DRT as a function of λ [34]. Due to Kramers-Kronig relation, the difference should be zero. Thus, the regularization parameter is expected to fit best when the above expression is closest to zero. In this case, noise impact is minimized but there is no over-regularization.
The second criterion, cross-validation [34], is based on a similar idea: Equation (19) takes x from imaginary parts and A, b containing real parts and vice versa. Using Kramers-Kronig relation again, both terms should be zero in theory. As for the discrepancy method, noise and discretization effects lead to deviations from the ideal result. In this case, λ is chosen appropriately if the DRT calculated from the imaginary part of the measurement data matches the real part of the measurement data and vice versa. The aim is similar to that of the discrepancy method: Kramers-Kronig being fulfilled best without causing over-regularization.
An L-shaped graph emerges with a distinct edge. The regularization parameter leading to the data point at the edge is used as the optimum. It can be seen as a compromise between the precision of the result (abscissa) and the norm (ordinate) which is forced to be small by the regularization: λ smaller than that value lead to large x by only marginally increasing the precision; large λ reduce the precision at a slight decrease in ||x||.
As can be seen from Table 1, the optimum value for λ differs strongly with the used method. Taking Figure 3 into account, oscillations occur at λ = 0.01. Thus <1 × 10 −5 and 6 × 10 −4 cannot be well-chosen. The parameter λ = 0.13 seems suitable, though. Therefore, for the RC-ZARC element used for this analysis, the discrepancy method is the only useful method. Modifying the ZARC to a second RC, keeping R 2 and C 2 constant, leads to a suggested regularization parameter of 2 and thus to over-regularization. There are two possible explanations why none of the proposed algorithms works satisfactorily: first, this could be due to another algorithm used as explained in Section 2 having a different sensitivity towards λ. Second, the magnitude of the investigated system's impedance could play a role as λ · x scales linearly with x. The latter could be ruled out by scaling R and C of our example system where the results coincide.  Since none of the above criteria lead to physically meaningful values for λ, the accuracy based on the sum of squared errors sse = ∑ i (y i −ŷ i ) 2 of the reconstructed spectrum is used. Considering the sse as a quality criterion for the choice of λ, a co-dependency from the number of considered time constants n must be assumed. This is reasonable because a higher amount of predefined time constants should eventually lead to a better match of the original and the reconstructed spectrum. Figure 4a,b reveal that the spectrum obtained by the DRT correlates disproportionately less to the original impedance spectrum with growing λ. Starting at around 0.2, the slope increases significantly. For n ≥ 2 m, there is no significant reduction of the error below λ = 0.1 while oscillations increase (cf. Figure 3). Thus, 0.1 ≤ λ ≤ 0.2 should be chosen. This range was verified by a parameter variation of the simulated test circuit and by the application of the proposed DRT algorithm to measurement data. These values proved adequate in most of our use cases but possibly have to be adapted when applying the algorithm to different impedance data. Hereinafter, λ is set to 0.1. Even if smaller or larger values for λ might lead to reasonable distributions, the later interpretation of the polarization processes would be less conclusive.
Analyzing the number of time constants n, Figure 4a,b show that the sse decreases significantly until n = 3 m and does not change much for larger numbers of time constants within the preset range for λ. Figure 4c displays g k for n = m up to n = 5 m. Up to n = 2 m, the graph is edgy, especially in the time constants range of the huge, narrow peak (approximately 1 × 10 −3 s) caused by the large ∆τ between τ k and τ k+1 . Furthermore, the peaks differ in height although g k is normalized. The reason for this is that the normalization is dependent on the number of time constants the plot consists of. So, more time constants would lead to a greater ∑ k g k if the height was equal. It is therefore expected that the height scales linearly with n. This is proven in Figure 4d where each graph, consisting of c · n time constants, is multiplied by c. The resulting curves overlap. Though, having only few time constants, the peaks lack height, but are wider as there are not enough time constants to precisely form the sharp peak. In addition, the characteristic time constant for the oscillation at low time constants shifts slightly to higher values with increasing n. This has no physical meaning, but is caused by the optimization algorithm and the regularization.
Summarizing, the number of time constants should be three times the number of frequency sampling points of the impedance data. The DRT does not change significantly with a further increase of sampling points, but the equation system enlarges and therefore the computation time rises from milliseconds to few seconds on a standard PC. As stated when analyzing the regularization parameter, there is no n where the DRT could be called wrong, but the cost-benefit ratio can be suboptimal.
In Figure 5, the parameter dependency of the condition number cond(A Reg ) is investigated. The non-regularized optimization problem is ill-conditioned as the condition number is in the range of 10 19 . Increasing λ, the condition number is reduced significantly, resulting in a linear relationship on a double logarithmic scale. The smaller the number of time constants within a constant frequency range, the lower the condition number. A dependency on the absolute values of τ min and τ max was observed as well as a correlation to the width of the interval [τ min τ max ]. Thus, Tikhonov regularization has also a substantial impact on the noise sensitivity of the DRT which is mainly determined by the condition number.

Setup of A and b
As mentioned in Section 2, in our case the matrix A as well as b consist of real and imaginary parts following the idea that using as much data as possible reduces noise or discretization effects and leads to more reliable results. To prove that assumption, the DRT of various RC-ZARC elements were calculated using real parts only, imaginary parts only and both real and imaginary parts. Furthermore, the time constants were varied from sub-µs to several seconds to cover a wide band of typical time constants in electrochemical systems. The sse was calculated for all parameters and the three matrices. As can be seen in Table 2, the error is minimal for all parameters when using both real and imaginary parts. This is supposed to be caused by the Kramers-Kronig relation not fulfilled perfectly due to numeric deviations for this mathematical equivalent circuit model. As the sse returns an absolute error, the absolute sse value varies strongly with respect to the total polarization resistance of the system. The capacities do not influence the sse significantly. Thus, by using real and imaginary parts, containing redundant information when neglecting numerical issues and measurement noise, the result can be improved by a factor of 1.2 to 1.6. Though, it has to be considered that an enlarged A is accompanied by an increase of the size of the equation system and thus by an increase in computation time.

Impact of Optimization Function and Solving Algorithm
The impact of the algorithm was investigated by comparing our approach ("Algorithm 1") with another one using an interior-point method as suggested by Byrd et al. [37] ("Algorithm 2"). The latter is based on the direct evaluation of Equation (9) which is possible due to the structure of the algorithm. For this investigation, only the real parts were taken into account for A as using both real and imaginary parts is not possible within Algorithm 2. Again, the identical RC-ZARC element was used setting λ = 0.1 and n = 3 m.
Four major differences can be obtained from Figure 6a. First, the height of the RC-peak increases while its width decreases. Second, the minimum between the peaks does not reach zero for Algorithm 2. Third, the time constant of the ZARC element shifts slightly to higher time constants coinciding with an increase in height and decrease in width. Fourth, the small peak at low time constants is not visible using Algorithm 2. Furthermore, Algorithm 2 causes g k to rise when reaching the highest time constants and |g k | is higher than for algorithm 1 when investigating the time constants where no specific process is visible. Concerning the sse, Algorithm 1 (4.57 × 10 −9 ) outmatches Algorithm 2 (9.78 × 10 −9 ) by a factor of approximately 2. This is in accordance with the analysis of the reconstructed impedance spectra (not displayed) underlining the hypothesis made in Section 2 that the DRT is strongly dependent on the exact algorithm and optimization function used and that spectra can neither be compared nor reproduced when the underlying mathematics are not revealed. Nevertheless, both algorithms yield a polarization resistance of 12.1 mΩ. This absolute value can therefore be used as a cross-link between different approaches on how to obtain the DRT. Three differences between the algorithms might cause the different DRT: • Convergence towards different local minima. This can be resolved using proper initial values and analyzing the reconstructed impedance spectrum. • Different optimization functions might lead to differing results, e.g., by using relative or absolute errors. • Treatment of non-negativity constraint. Various algorithms handle constraints differently which can influence the result.
The latter two cannot be resolved nor separated as they coincide. Thus, it has to be taken into account that, although both algorithms converge, the results may differ for numerical reasons. This must not be seen as a flaw of the DRT but more as an issue that has to be considered when comparing DRT of different sources.

Regularization Matrix
Usually, the regularization parameter is multiplied by an identity matrix. Nevertheless, e.g., Weese [24], uses a second derivative operator depicted in matrix form. Figure 6b displays the impact of the second derivative by Weese as well as a first derivative matrix following the same scheme. Considering the first derivative, the RC-peak decreases in height and increases in width while the ZARC peak behaves vice versa. Using the second derivative, both peaks decrease in height and increase in width showing a behavior similar to an increase of the regularization parameter. The distribution function between these peaks not becoming zero allows the same conclusion. In opposition, oscillations at high time constants increase significantly when using derivative matrices as would be expected by low regularization parameters.
Applying regularization using the first derivative can be interpreted as a regularization towards low slopes: at time constants where g k is steep, the regularization has a large effect compared to areas where g k is flat. This is exactly what is observed for the shape and position of the two main peaks in Figure 6b. The interpretation for the second derivative is conducted analogue to the first derivative considering its mathematical nature: The regularization term forces the curvature to be small. Hence, the peaks should widen and the tip should become less sharp. Again, this can be witnessed well in Figure 6b. The τ 1 peak shows a smoother shape at the local maximum leading to a lower overall height. The peaks becoming wider and therefore less easily separable coincides with oscillations leading to physically non-meaningful peaks. Using the derivative matrices for the regularization is not advantageous in our example, but might be appropriate when handling different data, e.g., when dealing with shallow slopes. Figure 6c underlines the hypothesis above, showing the regularization term λ · M · x: the regularization terms are large for large ||x|| using the identity, for large slopes using the first derivative and for large curvatures when using the second derivative. Thus, each of the matrices focuses on regularizing specific sections of the distribution function to a larger extent than others leading to the deviations visible in Figure 6b.

Analysis of Measurement Data
The algorithm was applied on measured impedance spectra to validate the method and to quantify the predominant polarization effects in a real electrochemical system. As specimen, a commercial Li-ion battery, Kokam SLPB526495 with a nominal capacity of 3.3 Ah in pouch cell format with a high energy NMC cathode was investigated. The impedance spectrum was measured using a Zahner ZENNIUM pro electrochemical workstation applying galvanostatic EIS at an amplitude of 120 mA within a frequency range of 20 mHz to 5 kHz with 21 steps per decade. The SOC (state of charge) was varied between 20% and 100% in 20% steps. For the variation, the cell was charged to 100% SOC as recommended by the manufacturer. After 3 h of relaxation, the impedance was measured. Then, the cell was repeatedly discharged at 0.5C for ∆SOC = 20%. Again, after 3 h of relaxation, the impedance was measured. The five resulting spectra are displayed in Figure 7a. The upper frequency limit was set so that there was no inductive behavior within the spectrum, while interface processes and the diffusion branch were visible. The ohmic resistance increased monotonously from 5.9 mΩ to 6.1 mΩ with decreasing SOC. The high-frequency feature of the impedance spectrum is typically associated with the process at the current collector/active material interface process [4]. Though, this assignment has to be studied carefully as there might be inductive contributions within that frequency range. The mid-frequency feature was strongly dependent on the SOC. It can be assumed that the underlying process was related to the charge transfer reaction [4]. The diffusion, which occurs at low frequencies, did not seem to show an obvious SOC-dependency within the Nyquist plot. A quantitative comparison is hard due to the offset in Re(Z) and possibly overlapping electrode processes between the different graphs. In Figure 7b, the distribution functions of the spectra in Figure 7a are shown. The DRT have been obtained by applying the cut-and-shift approach introduced in Section 2. In addition, the range of time constants has been extended by three decades towards low frequencies compared to the measured range also introduced within Section 2. The branch itself was fitted to the time constants for the measured frequencies, while semi-circles which were caused by contributions of lower frequencies in h k are necessary to bypass the resistive-capacitive restriction. This procedure was proven to be suitable by comparing the reconstructed spectrum to the original data in Figure 7c as the graphs do not deviate strongly at 20% SOC. The reconstructions of the other SOCs have been of equal quality (not shown here). Only for high frequencies, the graphs do not entirely match due to the onset of superimposed inductive contributions. This influence must be neglected for the DRT anyway. For the data point at 20 mHz there was a deviation due to the non-RC restriction which cannot be compensated perfectly, but in a sufficient manner. Figure 7d shows the absolute residual between measured and reconstructed impedance within the measured frequency range where Z > 0. Due to the cut-and-shift approach, the data point for the highest displayed frequency is close to zero. Thus, even a small difference would lead to a huge relative error which, in this case, is of no specific relevance. The good quality of the fit is proven by deviations below 0.2 mΩ for the entire spectrum, except the very small time constants. As stated above, this was caused by a beginning inductive influence, which can neither be avoided nor subtracted precisely. For the highest time constants, the difference begins to grow as the limitations regarding Equation (5) begin to impact the reconstructed spectrum.
Analyzing Figure 7b in more detail, the absolute distribution function h is used instead of g, enabling both a quantitative and qualitative analysis as the distribution functions are all within the same magnitude. The value h 1 , representing the contribution of the lowest time constant, is not equal to zero. This is caused by a combination of imperfect subtraction of the purely ohmic resistance and the presence of some minor inductive contributions at high frequencies. The absolute value of h 1 is strongly dependent on the number of measurement points, the regularization parameter and the number of time constants used within the DRT. This peak can be ignored for further investigations as it is not related to an electrode process. At high time constants, the contribution of the solid-state diffusion is dominant. Preliminary studies with constant phase and Warburg elements have shown that the three peaks closest to the main peak are also diffusion-related and necessary for a precise reconstruction. This is in accordance to a recent publication by Boukamp [38] who also derives a formula for the characteristic time constants of the peaks obtained in relation to the exact time constant calculated from finite length Warburg element formula. The results of our preliminary study are in very good accordance to the formula given as [38] where the main diffusion peak is referred to as τ 1 and the smaller ones as τ 2 , τ 3 etc. with decreasing time constant, respectively. With synthetic Warburg data, τ 0 calculated from τ 1 to τ 3 is equal. Increasing the regularization parameter, the calculated τ 0 differ increasingly, especially for high indices of τ, and shift towards lower time constants. Investigating the battery as a whole, the match of the theoretical and practical peak positions τ 0 differs up to a factor of eight between the three small and one large diffusion-related peaks. Furthermore, at 60% there is a small additional peak visible. This non-ideal behavior may be caused by the distributed solid-phase diffusion throughout the porous electrode material and by the simultaneous presence of anode and cathode diffusion process. These result in non-ideal Warburg behavior in combination with the regularization-based effect as mentioned above. Further, the polarization resistances of the four peaks differ from their theoretical values for the same reason. In contrast to Figure 7a, a significant change in diffusive behavior becomes visible when investigating the DRT. For intermediate SOCs, the diffusion shows a smaller peak in the DRT than for high or low SOCs within the measured frequency range. As previously seen in Figure 7a, the slowest interface process is strongly dependent on the electrodes' SOC and therefore on the lithium surface concentration of the particles. Starting with the fully charged cell, the polarization and the characteristic time constant decreases strongly and is almost constant for 80% and 60%. When discharging further, the polarization rises by more than a factor of two and the time constant shifts towards higher values. At around τ ≈ 1 × 10 −3 s there are two more processes visible. Due to their characteristic time constants, it can be assumed that these peaks represent electrode or interface processes. The faster process shows the highest polarization resistance at 20% and 100%, while the value is lower and close to constant for intermediate SOCs. The fastest process shows a constant shape and height. For further investigations, a quantitative analysis of the three electrode process-related peaks is inevitable.
The peak analysis shown in Figure 8 allows a quantitative analysis of each separate distribution in Figure 7b. As mentioned above, peak 1 does not show a noteworthy SOC-dependency and has a polarization resistance of 1.8 mΩ. It can be assumed that this peak represents the particle-particle or particle-current collector interface [4].   Peaks 2 and 4 show a similar SOC dependency: For high and low SOCs, the polarization resistance is large, while it is low for intermediate SOCs. The higher the resistance, the higher is furthermore the time constant. The width of the peaks does only change due to the change in height. As there are two charge transfer reactions expected within the battery, it can be assumed that peaks 2 and 4 represent those. Peak 3 behaves differently: it significantly changes its width and the characteristic time constant at 60% is lower than could be expected when compared to the other SOCs. This peak could hence be assigned to the SEI, while the exact reason for this behavior remains unclear.
Considering these findings, peak 4 can be assigned to the anodic charge transfer reaction based on the following indicators:

1.
The time constant of the anodic charge transfer is typically higher than that of the SEI [39,40] 2.
The resistance increases for high and low values reaching a minimum at medium SOCs. Figure 9 underlines this observation. This is in accordance with the Butler-Volmer kinetics, where the exchange current density is small at high and low SOCs [14,41]. To differentiate the process from the cathodic reaction, Figure 9 shows that the resistance is highest at low SOCs, whereas the increase is only shallow at 80% and 100%. It has to be considered that in commercial cells, the anode is over-dimensioned so that lithium plating is avoided by restricting high lithium concentrations in the anode. This leads to an disproportionate distribution of the anode's degree of lithiation. Thus, its SOC is not well-proportionate but shifted towards smaller concentrations. This leads to the conclusion that peak 3 is caused by the anodic charge transfer reaction. The qualitative investigation of the diffusive process in Figure 9 is of interest as well. The absolute value is calculated by adding the resistances of peaks 5-8. Summing up the values of charge transfer and diffusion obtained e.g., for 20% SOC, approx. 37 mΩ are achieved. However, the maximum value of |Z| is below 30 mΩ in Figure 7a. As the assumption of convergence towards Im(Z(jω min )) = 0 for frequencies below the measured range is not necessarily fulfilled, e.g., when a limited reservoir leads to capacitive behavior for very low frequencies, the polarization contribution obtained for the diffusive branch must neither be interpreted as an absolute value for an overall diffusion resistance nor be used to draw conclusions of processes outside the measured frequency range. Instead, the values can be taken into account for a comparison of different spectra to visualize changes in diffusive behavior, if the frequency range and the regularization parameter is equal for the investigated spectra.
The qualitative behavior of the diffusive values can be interpreted in two ways: in a descriptive manner, the diffusion is limited when the concentrations are either close to the maximum while there are still more ions to be intercalated and transported towards the particles' center or when the concentration is close to zero, but there is a demand of ions from the inner particle to be deintercalated at the surface. Considering electrochemistry, diffusion leads to a high polarization when the differential capacity is small, i.e., a small change of concentration leads to a large change in the open circuit voltage when the surface concentration differs from the average concentration. This is the case at high SOCs, but even stronger at low SOCs for the investigated NMC cell. Both explanations lead to the same result: As the resistance increases for high and for low SOCs in a similar way, it can be assumed that anodic and cathodic solid state diffusion overlap within the diffusion peak showing similar time constants without the possibility of separation.
Analyzing the DRT, we are able to gain significantly more information on the dominant polarization mechanisms. The SOC dependency of solid state diffusion could be assessed quantitatively and the dominating electrode processes could be assigned to peak 2-4.

Conclusions
The DRT offers a model-free characterization for electrochemical systems. Its information value is strongly dependent on a proper choice of process parameters which are needed to solve the ill-posed problem. Choosing the regularization parameter λ too small leads to oscillations and a large number of peaks not connected to physical processes. A large regularization parameter leads to a merger of peaks so that independent processes cannot be separated properly. In our case, λ should be chosen between 0.1 and 0.2. This value is strongly dependent on the exact algorithm used as there are publications using λ differing by magnitudes. The measurement data and its signal-to-noise-ratio also influences the optimum value. The number of sampling points of time constants for the DRT analysis should be large enough for a smooth result. Exceeding a reasonable amount increases the computational effort without significant improvements to the result, though. The threefold number of measured frequency points has shown to be appropriate. A smaller number increases the difference between measured and reconstructed impedance spectrum notably, while a higher number does not reduce the error significantly. Choosing the highest time constant three decades above the minimum measured frequency, processes where lim τ→0 (Im(Z)) = 0 can be analyzed.
The regularization parameter can be used as a linear factor on the distribution function itself or on the first or second derivative. Dependent on that, either the absolute height, the slope or the curvature is regularized leading to contrasting effects on g k . The implementation of those approaches verifies this hypothesis. We recommend to use the non-derived regularization matrix for most cases. Furthermore, the optimization algorithm used, the optimization function fed into the algorithm and the way the algorithm deals with the non-negativity constraint affect the outcome. This constraint is physically inevitable as there must not be negative contributions to the overall polarization. With the parameters set as above, processes with similar polarization resistances can be separated if their time constants differ by a factor of two to four. The exact factor cannot be determined in advance but is dependent on the shape of the individual peaks. If the polarization resistance of peaks differ strongly, the gap between the characteristic time constants needs to be higher to be able to observe the smaller peak. The algorithm used is available at http://www.ec-idea.uni-bayreuth.de within a software package based on the Matlab runtime environment.
The introduced peak analysis algorithm facilitates the separation of the peaks and enables the calculation of the individual polarization resistances. One large and several smaller peaks at high time constants can be assigned to solid state diffusion. A hypothesis was established to correlate three specific peaks to the anodic and cathodic charge transfer and the SEI due to their time constants and SOC behavior. For both identified processes, diffusion and charge transfer, a non-monotonous behavior emerges with a maximum polarization resistance for low SOCs, decreasing for medium and again increasing resistances for high SOCs. Investigating the impedance spectrum only, e.g., by fitting equivalent circuits, there is no evidence on the number of involved processes. The DRT offers the benefit of a detailed insight into the electrochemical mechanisms taking place in the specimen and their characteristic time constants. Acknowledgments: This publication was funded by the German Research Foundation (DFG) and the University of Bayreuth in the funding programme Open Access Publishing.

Conflicts of Interest:
The authors declare no conflict of interest.