A Kernel for Calculating PEM Fuel Cell Distribution of Relaxation Times

Impedance of all oxygen transport processes in PEM fuel cell has negative real part in some frequency domain. A kernel for calculation of distribution of relaxation times (DRT) of a PEM fuel cell is suggested. The kernel is designed for capturing impedance with negative real part and it stems from the equation for impedance of oxygen transport through the gas-diffusion transport layer (doi:10.1149/2.0911509jes). Using recent analytical solution for the cell impedance, it is shown that DRT calculated with the novel K 2 kernel correctly captures the GDL transport peak, whereas the classic DRT based on the RC-circuit (Debye) kernel misses this peak. Using K 2 kernel, analysis of DRT spectra of a real PEMFC is performed. The leftmost on the frequency scale DRT peak represents oxygen transport in the channel, and the rightmost peak is due to proton transport in the cathode catalyst layer. The second, third, and fourth peaks exhibit oxygen transport in the GDL, faradaic reactions on the cathode side, and oxygen transport in the catalyst layer, respectively.


INTRODUCTION
Electrochemical impedance spectroscopy provides unique opportunity for testing and characterization of PEM fuel cells without interruption of current production mode (Lasia, 2014). A classic approach to interpretation of EIS data is construction of equivalent electric circuit having impedance spectrum close to the measured one. However, a more attractive option provides the distribution of relaxation times (DRT) technique. Note that DRT is sometimes called distribution function of relaxation times.
In the context of PEM fuel cell studies, the idea of DRT can be explained as follows. To a first approximation, PEM fuel cell impedance Z can be modeled as impedance of a parallel RC circuit where ω is the angular frequency of applied AC signal. This approximation corresponds to the cell with ideal (fast) transport of reactants in all transport media (Kulikovsky, 2017). In that case, R describes Tafel resistivity of the oxygen reduction reaction (ORR), and C represents the superficial double-layer capacitance of the electrode. In general, to calculate transport contributions to cell impedance, one has to develop a transport model for the ORR reactants. However, if information on transport resistivities and characteristic frequencies suffices, DRT provides a simpler option. Denoting RC τ, multiplying the right side of Eq. 1 by nonnegative function c(τ) and integrating over τ, we get where R pol is the total polarization resistivity of the cell, and the function c is the DRT of impedance Z. Mathematically, Eq. 2 means expansion of Z(ω) into infinite sum of RC impedances, with the resistivity of each elementary RC circuit being R pol c dτ. The function under integral in Eq. 2 is usually called "Debye model" (Barsoukov and Macdonald, 2018). Eq. 2 can also be considered as integral transform of c(τ), which justifies the term "RC kernel" for Eq. 3. Quite evidently, DRT of a single parallel RC circuit is Dirac delta function c δ(τ − τ * ) positioned at τ * RC. This example illustrates the main feature of DRT: it converts any RC-like impedance into a single, more or less smeared on the τ scale δ-like peak.
All transport processes in a fuel cell eventually are linked to the double-layer capacitance in the catalyst layer; thus, it is usually assumed that impedance of every process is not far from impedance of a parallel RC circuit. That means that the DRT of a PEMFC is expected to consist of several delta-like peaks. As the regular frequency f 1/(2πτ), it is convenient to plot c(f) instead of c(τ). Position of each c(f) peak marks a characteristic frequency of the respective transport process and gives the contribution of process resistivity to the total cell polarization resistivity R pol . Here, τ n and τ n+1 are the peak boundaries on the τ scale. Fuel cell impedance is usually measured on equidistant in logscale frequency mesh {f n , n 1, . . . , N}, with ln (f n+1 ) − ln (f n ) being independent of n. From numerical perspective, it is beneficial to deal with the function G(τ) satisfying to where the term R ∞ is added to describe pure ohmic (highfrequency) fuel cell resistivity. Clearly, c G/τ, and Eq. 4 in terms of G takes the form where the frequencies f n 1/(2πτ n+1 ) and f n+1 1/(2πτ n ) mark the peak boundaries. Note that Eqs 4, 6 return a single peak resistivity if the DRT peaks are not overlapped; otherwise, R n would represent the sum of overlapped peak resistivities. The problem of overlapped peak resistivity is alleviated if the DRT in analytical form is found, as in, for example, the ISGP method (Avioz Cohen et al., 2021).
Setting in Eq. 5 ω 0 and taking into account that Z (0) − R ∞ R pol , we see that G obeys to normalization condition In the following, Eq. 5 will be discussed, as G is usually used instead of c in practical calculations.
DRT technique, Eq. 2, was invented by Fuoss and Kirkwood (1941) in the context of polymer materials impedance and brought to the fuel cell community seemingly by Schichlein et al. (2002). Since 2002, a lot of works from the group of Ivers-Tiffée have been devoted to deciphering of solid oxide fuel cell spectra by means of DRT [see a review (Ivers-Tiffee and Weber, 2017)]. Analysis of PEMFC impedance spectra using DRT is a relatively new field (Heinzmann et al., 2018;Avioz Cohen et al., 2021;Reshetenko and Kulikovsky, 2021;Wang et al., 2021). Heinzmann et al. (2018) measured impedance spectra of a small (1 cm 2 ) laboratory PEMFC and studied DRT peak behavior depending on cell temperature, relative humidity (RH), oxygen concentration, and current density. They obtained DRT with up to five peaks; the leftmost on the frequency scale peak P1 was attributed to oxygen diffusion in the gas-diffusion and cathode catalyst layers (CCLs). Avioz Cohen et al. (2021) performed impedance measurements of a 5-cm 2 cell varying temperature, RH, and current density. The calculated DRT consisted of four peaks, attributed (in ascending frequencies) to (1) oxygen transport in the GDL/CCL, (2) ORR, (3) proton transport in the CCL, and (4) proton transport in membrane. Note that Heinzmann et al. and Cohen et al. used different codes for DRT calculation. Wang et al. (2021) measured impedance spectra of application-relevant 25-cm 2 PEMFC and obtained a three-peak DRT; the lowest frequency peak was attributed to oxygen diffusion processes in the cell. In our recent work (Reshetenko and Kulikovsky, 2021), DRT spectra of a low-Pt PEMFC have been reported; we attributed the low-frequency peak to oxygen transport in the GDL and, possibly, in the channel.
In PEMFCs, the supplied oxygen (air) is transported through the four quite different media: channel, GDL, open pores of the CCL, and finally through Nafion film covering Pt/C agglomerates. One therefore could expect four corresponding peaks in the DRT spectra. However, in Heinzmann et al. (2018), Avioz Cohen et al. (2021), and Wang et al. (2021), a single oxygen transport peak has been reported. There are two options to explain this result: either some of the oxygen transport peaks overlap with each other (or with the ORR peak) and DRT is not able to separate them, or the codes used were unable to resolve all the transport processes. It is important to note that the code for DRT calculation of Wan et al. (2015) used in Heinzmann et al. (2018) and Wang et al. (2021), the ISGP code (Hershkovitz et al., 2011) used in Avioz Cohen et al. (2021, and our code using Tikhonov regularization (Tikhonov, 1995) in combination with NNLS solver (Kulikovsky, 2020a;Kulikovsky, 2021a) are based on the RC kernel, Eq. 5.
The real part of RC-circuit impedance, Eq. 1, is positive, and the imaginary part is negative. This imposes limits on functions Z, which could be represented by Eq. 5. For example, impedance of an inductive loop cannot be expanded in infinite series of RC impedances, as the imaginary part of inductive impedance is positive. Quite similarly, the impedance Z having negative real part in some frequency domain also cannot be represented by Eq. 5.
Below, we show that the DRT calculated with RC kernel could completely miss some of the transport peaks in PEMFC spectra. An alternative K 2 kernel better capturing oxygen transport processes in the cell cathode is suggested. The kernel is illustrated by calculation of DRT of the recent analytical PEM fuel cell impedance spectrum (Kulikovsky, 2021b). Finally, we show that the new kernel well separates the channel, GDL, and ORR peaks in the DRT spectra of a standard Pt/C-based PEM fuel cell.

MODEL: K KERNEL
Below, the following dimensionless variables will be used Here, t * is the characteristic time of double-layer charging c is the oxygen concentration, c in h is the reference oxygen concentration, C dl is the volumetric double-layer capacitance, J is the mean current density in the cell, η is the ORR overpotential, positive by convention, b is the ORR Tafel slope, i * is the volumetric ORR exchange current density, l t is the CCL thickness, D b is the oxygen diffusion coefficient in the GDL, and l b is the GDL thickness.
Analytical GDL impedanceZ gdl has been derived in Kulikovsky and Shamardina (2015). For the cell current densities well below the limiting current density due to oxygen transport in the GDL,Z gdl has the form Eq. 10 has been obtained from the general expression for the cathode side impedance in the limit of infinite air flow stoichiometry (Kulikovsky and Shamardina, 2015). Eq. 10 is a Warburg finite-length impedance divided by the factor (1 + iω/J). The factor describes the effect of double-layer charging by the cell current densityJ transported in the form of oxygen flux through the GDL to the attached CCL (Kulikovsky and Shamardina, 2015). In his classic work (Warburg, 1899), Warburg used static polarization curve to derive the boundary condition for calculation of transport impedance of a semiinfinite electrode [Eq. 4 of Ref. (Warburg, 1899)]. Later, Warburg model was extended for the case of finite-length transport layer, using the same boundary condition on the electrode side [see Lasia (2014), page 104]. However, account of capacitive term in the electrode charge conservation equation changes the boundary condition for the oxygen transport equation (Kulikovsky and Shamardina, 2015), leading to the  Table 1. FIGURE 2 | Nyquist spectra of the channel impedance, Appendix Eq. A1 and of impedance of oxygen transport in the CCL (Kulikovsky, 2017 (10) is shown in Figure 1A; frequency dependence of Re(Z gdl ) is depicted in Figure 1B.
As can be seen, between 20 and 200 Hz, the real part of Z gdl is essentially negative, and at higher frequencies Re(Z gdl ) tends to zero. As the real part of Warburg finite-length impedance is positive, the negative real part of impedance (10) is due to the term (1 + iω/J) in the denominator. Impedance of oxygen transport in channel, Appendix Eq. A1, and in the CCL (Kulikovsky, 2017) also exhibits negative real part ( Figure 2). Last but not least, in low-Pt cells, an important role plays oxygen transport through a thin Nafion film covering Pt/C agglomerates in the CCL (Greszler et al., 2012;Weber and Kusoglu, 2014;Kongkanand and Mathias, 2016). The spectrum of this transport layer is quite similar in shape to the spectrum in Figure 1 (Kulikovsky, 2021a). Thus, all the oxygen transport processes in a PEMFC cannot be described by the standard RC kernel, and another kernel suitable for description of impedance elements with the negative real part is needed.
In a standard PEMFC, unless the cell current density is small, the DRT peaks of oxygen transport in the GDL, CCL, and channel are expected to locate at frequencies below the frequency f ct of faradaic (charge transfer) processes. Thus, to capture the oxygen transport peaks, correction for the negative real part of impedance is needed in the range of frequencies f < f ct . The kernel K 2 suggested in his work consists, thus, of two parts: where f * ≃ f ct is the threshold frequency. Selection of optimal f * is discussed below. A function similar to impedance of a transport layer (TL) (10) forms the low-frequency part of K 2 ; the real part of this function is negative at ωτ > 1.81052. The high-frequency (RC) part of K 2 is the standard RC-circuit kernel. Switching between TL and RC kernels is necessary, as the TL part itself does not describe well RC-circuit impedance (see below). The real part of GDL impedance becomes negative at ≃20 Hz ( Figure 1B),  while the imaginary part at this frequency is quite large ( Figure 1A). This shape of impedance around and below 20 Hz is thus far from the RC-circuit impedance and the TL kernel is needed to "recognize" this shape. At higher frequencies, transport impedances tend to zero, and the Debye kernel works better. The idea behind Eq. 12 is thus to expand the lowfrequency components of cell impedance using the TL kernel and the high-frequency components using the standard RC kernel.
It is convenient to combine Eq. 12 into one function where α is a step function of the frequency f H(x) is the Heaviside step function, and ϵ 10 -10 is a small parameter to avoid zero division error. Parameter α therefore changes from 1 to 0 at the threshold frequency f f * . With α → 0, the Warburg factor in Eq. 13 tends to unity: and hence the α function serves as a switch between the two kernels in Eq. 12.
With Eq. 13, 5 takes the form Setting in Eq. 16 ω → 0, it is easy to show that G still obeys to normalization conditions, Eq. 7.    Figure 3 shows the imaginary part of the RC-circuit impedance for R 1, C 0.01 and the DRT calculated using Eq. 16, and α 1 for all frequencies (TL kernel). As can be seen, the main peak is positioned correctly at the frequency 1/(2πRC); however, the DRT spectrum exhibits three phantom peaks located to the right of the main peak. It is worth noting that quite similarly, the RC kernel generates several phantom peaks in the DRT of Warburg finite-length impedance (Kulikovsky, 2020a), that is, no single kernel is able to correctly represent DRT of all impedance components.

Synthetic Impedance Tests
In a standard PEMFC, the characteristic frequencies of channel and GDL impedance are approximately 1 and 10 Hz, respectively (Kulikovsky, 2021b). Thus, the typical value of the threshold frequency f * in Eq. 14 should be about 10 Hz; however, the exact value can always be selected simply looking at the calculated DRT spectrum. In standard PEMFCs operated at oxygen stoichiometry of 2, the faradaic DRT peak is located to the right of the GDL transport peak on the log-frequency scale (see below).
Analytical impedanceZ tot of the PEMFC cathode side has been obtained in Kulikovsky (2021b), assuming fast proton and oxygen transport in the CCL. Equation forZ tot includes three components: impedanceZ chan due to oxygen transport in channel, impedanceZ gdl of oxygen transport in the GDL, and faradaic (charge-transfer) impedanceZ ct The respective formulas are listed in the Appendix; these solutions allow us to check how well DRT from Eq. 16 captures the channel, GDL, and faradaic components in the total impedanceZ tot spectrum (Figure 4).
The spectrum ofZ tot has been calculated in the frequency range of 10 -2 to 10 4 Hz with 22 points per decade. Parameters for impedance calculation are listed in Table 1. Figure 4D depicts the imaginary part of the impedance components calculated using equations in the Appendix. Imaginary part ofZ tot has then been used for calculation of DRT using our recent algorithm based on Tikhonov's regularization and nonnegative least-squares (NNLS) solver (Kulikovsky, 2020a;Kulikovsky, 2021a). The NNLS method greatly outperforms projected gradient iterations suggested in Kulikovsky (2020a). In all the cases, the L-curve method (Hansen, 1992) gave the regularization parameter λ T ≃ 10 -3 . Variation of λ T in the range of plus-minus order of magnitude did not change the number and position of peaks. At 10 times larger λ T , the peaks get wider and started overlapping. Figure 4A shows the DRT spectrum ofZ tot , Appendix Eq. A6, calculated using the RC kernel. As can be seen, the standard kernel returns only two peaks corresponding to the channel (left peak) and faradaic (right peak) impedance. The GDL peak, which is clearly seen in Figure 4D, is completely missing. Note poor quality of reconstructed imaginary part (red open circles) between 0.1 and 20 Hz. This is a result of poor description of the cell impedance by Eq. 5 in the frequency domain where the real part of GDL impedance is negative. Figure 4B displays the DRT calculated with the "pure" TL kernel, Eq. 12. The GDL peak is resolved, and the quality of reconstructed imaginary part is much better; however, phantom high-frequency peaks to the right of the faradaic peak are clearly seen (cf. Figure 3). Figure 4C shows the DRT spectrum calculated using the K 2 kernel with the threshold frequency of 10 Hz; the GDL peak is well resolved, and the phantom peaks vanish. It is worth mentioning that the real part of the total cell impedance is always positive due to positive contributions of the faradaic and proton transport impedances. However, the Debye kernel fails to recognize the transport impedances having negative real part, as the whole shape of this impedance strongly deviates from the RC circuit one. The result in Figure 4C shows that in order to be recognized, the DRT peak corresponding to the process with negative real part must be fully "covered" by the TL kernel. The respective DRT peak frequency corresponds to the peak value of the imaginary part of impedance (cf. Figures 4C,D). Table 2 shows the resistivities, Eq. 6, corresponding to individual peaks in Figure 4. As can be seen, K 2 kernel provides good estimate of the channel and faradaic resistivities; however, the GDL resistivity R gdl is underestimated by 30%. Nonetheless, as the contribution of R gdl is small, the 30% accuracy could be tolerated.

Real PEMFC Spectra
A crucial check for the new kernel is calculation of DRT of a real PEM fuel cell. Impedance spectra of a standard Pt/C-based PEMFC have been measured in the frequency range of 0.1 to approximately 10 3 Hz with 11 points per decade. The cell geometrical parameters and operating conditions are listed in Table 3; note that the air flow stoichiometry was 2 in this set of measurements. The impedance points in the frequency range above ≃ 10 3 Hz have been discarded due to effect of cable inductance. More details on experimental setup and measuring procedures can be found in Reshetenko and Kulikovsky (2019). Figure 5 shows DRT spectra calculated with the real part of measured impedance using the RC kernel. Figure 6 shows the respective peak frequencies and resistivities. The DRT spectra in Figures 5A-C exhibit four peaks, whereas in Figure 5D, the most high-frequency peak disappears. This peak represents proton transport in the CCL, and at high cell currents, it shifts to frequencies that have been discarded. The characteristic frequency f 4 of proton transport in the CCL is given by Kulikovsky (2020b) With the typical values of σ p ≃ 0.01 S cm −1 and C dl ≃ 20 F cm −3 [Ref. (Reshetenko and Kulikovsky, 2019)], we get f p ≃ 700 Hz, which by the order of magnitude agrees with the proton peak position in Figures 5A-C. The growth of f 4 with the cell current and the respective decay of the peak resistivity R 4 ( Figure 6C) is due to growing amount of liquid water improving the CCL proton conductivity.
The leftmost peak in Figures 5A-D represents impedance due to oxygen transport in the cathode channel. The characteristic frequency f 1 of this peak linearly increases with the cell current density ( Figure 6A), which is a signature of channel impedance (Kulikovsky, 2021b).
The highest, second peak in the DRT spectra ( Figures  5A-D) represents the contributions of ORR and oxygen transport in the GDL (see below). In the absence of strong oxygen and proton transport limitations, the ORR resistivity is given by which follows the Tafel law. Qualitatively, the shape of second peak resistivity R 2 follows the trend of Eq. 18 due to dominating contribution of ORR resistivity to this peak ( Figure 6B). Note that the separate GDL peak is not resolved by the RC kernel. The third CCL peak in Figures 5A-D is most probably due to oxygen transport in the CCL pores. For the estimate. we take the Warburg finite-length formula for the transport layer where D is the oxygen diffusivity in the transport layer of the thickness l. Setting f W 35 Hz ( Figure 5A) and l l t , for the oxygen diffusion coefficient in the CCL, we get D ox ≃ 1.2 · 10 -4 cm 2 s −1 , which agrees with measurements (Reshetenko and Kulikovsky, 2019). With the growth of cell current, the peak frequency f 3 rapidly shifts to 200 Hz ( Figure 6C), corresponding to D ox ≃ 7 · 10 -4 cm 2 s −1 . This result correlates with the growth of D ox resulting from fitting of physics-based impedance model to the PEMFC cell spectra (Reshetenko and Kulikovsky, 2016). The mechanism of this growth yet is unclear. Figure 7 shows the DRT of the same impedance spectra calculated with the K 2 kernel. The properties of K 2 kernel are immediately seen: setting of the threshold frequency f * in Eq. 14 just to the left of the "ORR + GDL" peak in Figures 5A-D splits this peak into two well-resolved peaks ( Figures 7A-D). The left peak of this doublet corresponds to the GDL impedance and the right peak to the ORR impedance. With the K 2 kernel, the proton transport peak is seen only at the smallest cell current density ( Figure 7A), whereas at higher currents the peak vanishes, indicating its shift to the frequencies greater than 1 kHz ( Figures 7B-D).
Splitting the "CCL + GDL" peak into GDL and ORR peaks is confirmed by the behavior of peak resistivities in Figures 8B,C, respectively. The ORR peak resistivity R 3 follows the trend of Eq. 18 (solid line in Figure 8C) with the ORR Tafel slope b 30 mV, which is a typical value for Pt/C cells (Neyerlin et al., 2006). The GDL resistivity R 2 decreases in the range of cell currents 100 to 400 mA cm −2 and remains nearly constant at higher currents ( Figure 8B). Using again the Warburg formula Eq. 19, with l l b 230 · 10 -4 cm and the frequency between 10 and 30 Hz ( Figure 8B), for the GDL oxygen diffusivity, we get quite reasonable values of D b ≃ 0.013-0.033 cm 2 s −1 . The increase of D b in the range of 100 to 400 mA cm −2 is probably due to growing air flow velocity in the channel at the constant stoichiometry, which facilitates liquid droplets removal from the GDL.
Comparing the GDL peaks in Figures 7A,B one can see that the area under this peak increases with the cell current density. However, the area under the peak gives the fraction of the respective process resistivity in the total polarization resistivity of the cell R pol . The latter value decreases with the cell current, leading to decrease in the absolute value of R GDL ( Figure 8B).
K 2 kernel returns twice lower resistivity and about twice higher frequency of the CCL peak (cf. f 3 , R 3 in Figure 6C and f 4 , R 4 in Figure 8C). This shift leads to twice higher estimate of the CCL oxygen diffusivity, which is still acceptable (see above). Overall, confirmation of the CCL peak nature requires measurements at variable oxygen concentration and RH. An important issue is accuracy of peak resistivities calculated using the K 2 kernel. Impedance spectra fitting by physics-based model seems to be the only way to get an independent estimate of the resistivities for comparison.
Over the past years, large efforts have been directed toward development of universal code capable to calculate DRT based on RC kernel, not using any a priori information on the system [see a nice review of Effendy, Song and Bazant (Effendy et al., 2020)]. However, in PEMFC studies, it would be wasteful to ignore analytical results, showing that the RC kernel alone is not well suited for DRT description of the spectra.

CONCLUSION
Impedance of all oxygen transport processes in a PEM fuel cell exhibits negative real part in some frequency range. This makes it difficult to accurately calculate the respective DRT peaks using the standard RC kernel 1/(1 + iωτ). A novel kernel K 2 , Eq. 13, is suggested. K 2 combines the low-frequency transport layer kernel having a domain with negative real part and the standard RC kernel for description of faradaic and high-frequency processes in the cell. Calculation of DRT for analytical PEMFC impedance shows that K 2 kernel captures the peak due to oxygen transport in the gas-diffusion layer, whereas the RC kernel can miss this peak. Comparison of Pt/C PEMFC DRT calculated using RC and K 2 kernel shows that the K 2 kernel resolves the GDL oxygen transport peak, which otherwise is merged to the ORR peak when using the standard RC kernel. Overall, the K 2 spectra of a standard Pt/C PEMFC operating at the air flow stoichiometry λ 2 consist of five peaks. In the frequency ascending order, these peaks are due to (1) oxygen transport in channel, (2) oxygen transport in the GDL, (3) faradaic reactions, (4) oxygen transport in the CCL, and (5) proton transport in the CCL. If the CCL proton conductivity is high, the peak (6) shifts to the frequencies well above 1 kHz, and it may not be resolved due to inductance of measuring system.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.