Application of finite Gaussian process distribution of relaxation times on SOFC electrodes

Electrochemical impedance spectroscopy (EIS) is a powerful tool in characterisation of processes in electro-chemical systems, allowing us to elucidate the resistance and characteristic frequency of physical properties such as reaction and transport rates. The essence of EIS is the relationship between current and potential at a given frequency. However, it is often the case that we do not understand the electrochemical system well enough to fit a meaningful physical model to EIS data. The distribution of relaxation times (DRT) calculation assumes an infinite series of relaxation processes distributed over a characteristic timescale. The DRT calculation may identify the number of processes occurring, as well as their respective resistivity and characteristic timescale, and may resolve processes which have relatively similar timescales. Using a nonparametric tool known as Gaussian process (GP) regression, we showcase a method of finding a unique solution to the ill-posed DRT problem by optimising kernel hyperparameters as opposed to ad-hoc regularisation. In this work, we use finite GP regression under inequality constraints (fGP) to analysed EIS data generated by a (Ni/CGO | CGO | YSZ | Reference Cathode) solid-oxide fuel cell in a gas mixture of 0.5 bar H 2 /0.5 bar H 2 O and at a temperature of 600 ◦ C. By varying the current density, we can characterise the current-voltage relationship of the electrode and shed light on the reaction mechanism governing charge transfer at the solid-gas interface. Our findings also show that even at relatively high current densities ( ± 600 mA cm (cid:0) 2 ) the electrode process is limited by charge transfer.


Introduction
Amongst many looming environmental disasters, the concentration of CO 2 in our atmosphere is widely recognised to be driving the most catastrophic climate devastation.The need for a sustainable, reliable, and affordable energy economy has never been more urgent.The solid oxide fuel cell (SOFC) is a highly efficient chemical-to-electrical and electrical-to-chemical energy conversion technology compatible with both existing fuel (i.e.natural gas) and future fuel (i.e.renewably sourced hydrogen) infrastructures [1][2][3][4].Electrochemical devices, such as solid-oxide fuel cells (SOFCs), allow for reversible chemical to electrical energy conversion, with efficiency surpassing that of the combustion engine [5].The Faradaic reactions at the fuel (i.e.H 2(g) , or CO (g) ) and air electrodes of an SOFC can be given simply as [6,7] Fuel electrode : and Air electrode : respectively, where O 2− and e − represent oxide ions in the electrolyte and electrons in the electrode, respectively.At the electrodes, electrons are not directly transferred between the molecules, but are instead passed through a circuit whereby the chemical energy is used to drive external work.The open circuit voltage (V 0 ) of the SOFC at operational temperatures is approximately 1V.As current is drawn, resistances act by changing the cell voltage (V), this change in voltage is described as the overpotential (η) [8][9][10][11][12].Gadolinium doped ceria (CGO) displays considerably faster oxygen transport kinetics and high electronic conductivity under reducing conditions.To enhance the overall electronic conductivity of the E-mail address: nw7g140@gmail.com(N.J. Williams).electrode and to increase the density of reactive triple phase boundaries (TPB), it is beneficial to mix a metallic phase with the MIEC to create a cermet electrode i.e. nickel/gadolinium doped ceria (Ni/CGO) [1,[13][14][15][16].The operation of the Ni/CGO electrode under electrical bias has been studied for decades, however a unifying model for hydrogen electro-oxidation or water electrolysis has not yet been agreed upon [14,[17][18][19][20][21][22][23][24][25][26][27].

Contents lists available at ScienceDirect
A convenient way to characterise the reaction mechanisms at electrodes is to derive its current-voltage relationship using electrochemical impedance spectroscopy (EIS) [28][29][30][31][32][33].However, the complex relationship between charge transfer and transport means that to model the porous electrode, we often use the transmission-line-type (TL) model which encompasses many assumptions regarding the homogeneity of the electrode [13,34].In such a case, it is likely that the constraints derived from the equivalent circuit fitting will not be representative of the electrode.Therefore, it is often the case that we do not understand the electrochemical system well enough to fit a meaningful physical model to EIS data.The distribution of relaxation times (DRT) calculation assumes an infinite series of relaxation processes distributed over a characteristic timescale [35].The DRT calculation may identify the number of processes occurring, as well as their respective resistivity and characteristic timescale.A key utility of DRT is that it may resolve processes which have relatively similar timescales, and would otherwise appear to be overlapping after EIS analysis [7,36].
The optimisation problem in general is ill-posed, meaning there are an infinite number of solutions which satisfy the calculated DRT [37].Consequently, regularisation is frequently implemented to make the solution physically meaningful [35].Ciucci et al. have investigated a variety of machine learning tools to deconvolute the DRT and to predict the EIS values at unmeasured frequencies [38,36,[39][40][41][42].One such method is known as Gaussian process (GP) regression, which is a computationally cheap, nonparametric Bayesian framework for modelling stochastic processes [43,44].However, GP regression of the DRT has two drawbacks, i) only the imaginary part of the impedance can be predicted, and ii) positivity of the DRT is not obeyed [41].Moreover, they provide less realistic uncertainties when physical systems satisfy inequality constraints such as linear inequalities, gradient inequalities or convex inequalities [45,46].To overcome these issues, we follow the works of Maradesa et al. by implementing Gaussian processes with inequality constraints, otherwise known as finite Gaussian process (fGP) regression [42].

Distribution of Relaxation Times
EIS describes polarisation behaviour as a complex transfer function [47,48,7].The DRT problem looks at taking a Fourier transform of the linear Fredholm integral of the first kind, to analyse at the same data in the time domain (Derivation in SI) [49] where γ(lnτ) represents the non-negative distribution function, and ∫ (1 + iωτ) − 1 dlnτ represents the functional-space kernel [35].γ(lnτ) is approximated as a sum of radial basis functions (RBF), ψ k (lnτ), centred around the timescale τ k and can be given in the weight-space view [50] γ(lnτ) = where x k are unknown scalars to be computed numerically, and ψ k (lnτ) is often a radial basis function kernel [51].Note that Eq. 4 is the Fredholm integral of the second kind if R ∞ is allowed to vary.By combining Eq. 3 and 4 the DRT is computed numerically, as such, the DRT integral must be discretised yielding the finite sum, and separated into real and imaginary parts and can be given in the functional-space view [36,39,50] where 1⩽n⩽N is the total number of experimental measurements.The optimisation problem in general is ill-posed, meaning there are an infinite number of solutions which satisfy Eq. 5 [37].Consequently, regularisation is required to make the solution physically meaningful [35].The most popular method is Tikhonov regularisation, where the optimisation problem forces the distribution to smooth [52][53][54][55]48,36].
x can be solved with respect to x by minimising the loss function, S(x), by means of regularisation where Z exp is the experimental impedance data, and Ω ′′ are weighted factors (kept at unity for this study), λ is the regularisation parameter and ||L m x|| 2 2 is the regularisation term with L m a suitable differentiation matrix [40,39,36,37,41].

Bayesian Inference
In Bayesian statistics, the probability represents the degree of belief in an event, rather than the relative frequency with which the event occurs [37,[56][57][58].Before data are introduced to the model we express our belief about the coefficients through a prior probability distribution, known as the prior, π(γ|ξ).The prior beliefs are then updated based on the data obtained to yield the posterior distribution, π(γ|Z exp , ξ), where we have used the notation given by Ciucci et al. such that ξ = logf = − logτ, and where the full data-set is ξ = [ξ 1 , …, ξ K ] T [59].Updating the probability of a hypothesis as more evidence becomes available using Bayes' Theorem is the essence of Bayesian inference.In the context of the DRT problem, Bayes' Theorem states [37] π(γ|Z where the likelihood, π(Z exp |γ, ξ), quantifies how well the model fits the data and the evidence (or marginal likelihood), π(Z exp |ξ), is a normalizing constant that is independent of the coefficients x.It is possible to interpret optimisation problems from the Bayesian perspective as an expression of a prior belief about the inverse function.In order to derive the ordinary least squares loss function we assume that the measured impedance is the sum of the model impedance and the error term where ε ∼ N (0, σ 2 n I) is a normal distribution.As such, the magnitude of the regularization may be taken to correspond to the "expected error" of the fit [50].

Finite Gaussian Processes Distribution of Relaxation Times
Gaussian processes (GPs) are a nonparametric Bayesian framework for modelling stochastic processes [43,44].To maintain positivity of the DRT solution, we follow the works of Maradesa et al. by implementing Gaussian processes with inequality constraints [42].Note that in this derivation we neglect contributions from inductive effects where Z ′′ > 0. To handle the conditional distribution incorporating both interpolation conditions and inequality constraints, Maatouk et al. proposed using a finite-dimensional approximation of Gaussian processes N.J. Williams et al. [45,42] where Ψ = (Ψ 0 , …, Ψ K ) T is a vector of piecewise (or hat) basis functions [60] Using the model proposed by Maatouk et al. means that the conditional GP is reduced to simulating the Gaussian vector γ(ξ) within the the space of bounded functions described by a convex set [45].The coefficients x k are drawn from an independent, identically distributed Gaussian probability distributions, giving the Gaussian prior, π(γ|ξ) ∼ N (0,K), where K is defined in the SI and depends on σ 2 f , the variance of weights, x k .This leads to a distribution of DRT estimators, known as a Gaussian process (GP) where m(ξ) is the mean and k(ξ, ξ ′ ) is the kernel (full derivation in SI) [61].The smoothest kernel is the Radial Basis Function kernel (RBF), also known as a square exponential kernel [62,63,61] where σ f determines the average distance between your function and the mean, and l determines the length of the "wiggles" of the function, both are hyperparameters [57].If we now take a vector of function values such that [γ(ξ 1 ), γ(ξ 2 ), …, γ(ξ K ))] =: γ(ξ) we arrive at normal distribution function where K is shorthand for the Gram matrix k(ξ, ξ) of dimension K × K with ijth entry k(ξ j , ξ j ) (otherwise known as the kernel matrix) [62,63,61,43].The DRT is an inversion of the real data, and therefore it is not possible to condition the GP to the "true" DRT, as in the nature of the ill-posed problem.GPs are closed under linear transformations, therefore Z DRT is also a GP [64] Z where are the real and imaginary operators, respectively (defined in SI) [42].Eq. 8 can be modified to account for random errors for both real and imaginary impedance [42] where x ∼ N (0, Γ) and Γ = diag(σ 2 R , K) and K is the kernel Gram matrix.
We now have all hyperparamaters defined Θ = [σ n , σ R , σ f , l] T .A key utility of using GPs to model stochastic processes is the ability to make predictions on an unconditioned location, or a test set of K discrete points, using the joint distribution (covariance terms are derived in SI) [42] π Regression to determine the unmeasured output variables requires a conditional distribution.We can now condition x on Z exp with positivity constaintans [0, ∞] we get a truncated multi-normal distribution [65,66,56,57] π where the mean and covariance for the test function are [61,67] and Empirical Bayesian approaches such as maximising the marginal likelihood allow us to use continuous optimisation [56].The marginal likelihood π(Z exp |ξ,Θ), or evidence, implicitly integrates over all possible function values at unconditioned locations, and is calculated as the integral over the likelihood and prior [67].This integration can be computationally expensive and therefore avoided in other Bayesian methods such as maximum likelihood estimation and maximum a posteriori estimation [44].However, the advantage of using GPs is that the marginal likelihood can be integrated relatively easily [68,57,45].The marginal likelihood allows us to compute models by balancing between the fit with the data and the complexity of the model [56].The negative marginal log likelihood (NMLL) is illustrated in Fig.
The hyperparameters are then learned by minimising the NMLL via Hamiltonian Monte Carlo (HMC) sampling [70].HMC is used rather than the standard Gibbs sampler [71] since it has shown to be more efficient in higher dimensions, whilst also reducing mixing time for the Markov chain Monte Carlo (MCMC) algorithm.[58,65].Thus, making hyperparameter selection more robust.With the optimised hyperparameters it is then possible to calculate the mean and covariance of the inversion.

Results and Discussion
By conditioning the fGP model onto experimental data we observe an excellent fit to both the real and imaginary data, illustrated in Fig. 2a  and b, respectively.Experimental EIS data was collected from a (Ni/ CGO-CGO-YSZ-Reference Cathode) solid-oxide fuel cell in a gas mixture of 0.5 bar H 2 :0.5 bar H 2 O and at a temperature of 600 • C. By varying the current density, the impedance of the low-frequency process was modulated.We note that due to commercial sensitivity the ticks labels for the impedance have been censored.Since the fGP model allows us to use both parts of the complex transfer function, we can also illustrate the inversion results in the form of a Nyquist plot in Fig. 2c.Here we see that the convoluted high-frequency process is symmetric and appears consistent with the experimental data.The DRT in Fig. 2d shows four distinct peaks, combined with a 95% credibility interval.The convoluted peak in 2d has a slight asymmetry which is unexpected given the RC character of the convoluted transfer function in Fig. 2c.The DRT peak centred at 10 − 5 s represents the medium-frequency process in Fig. 2c.Characterising the processes responsible for these DRT features is outside the scope of this study.We believe that the relatively small peak escorting the large low-frequency process is caused by charge transport in the electrode, and is analogous to the asymmetry of the lowfrequency arc displayed in the Nyquist plot.Finally, we have assigned the large low-frequency DRT peak centered around 10 − 2 Hz to the general charge transfer reaction at the electrode-gas interface.Importantly, the non-negativity constraint has been satisfied making the DRT physically meaningful.
When a net current is drawn/applied an overpotential is generated.In this work the fuel cell has been driven to current densities between +600 mA (fuel cell mode) and -600 mA (electrolysis mode), as illustrated in Fig. 2e and f, respectively.Additionally, the low frequency process of the calculated DRT at negative and positive current densities is displayed in Fig. 2g and h, respectively.For both positive and negative current densities, all other DRT peaks are invariant as a function of overpotential.However, the large low-frequency process appears to be significantly dependent on the magnitude of the overpotential.We assign this process to the collection of electrode processes generally described by the transmission-line-type model [13].We note in the transmission-line-type model, the lowest frequency process becomes more pronounced from the distribution of escorting peaks when the charge transfer resistance is significantly higher than the transport resistances [13].This is observed in the Nyquist plot, where for increasing negative current densities, the charge transfer resistance is increased (Fig. 2e), where the transmission-line-type model appears to resemble an RC (or RQ) element.By integrating the area under the DRT peak using a Gaussian curve, we determined the resistance of the process.Using Ohm's law, V = IR, we calculated the overpotential as a function of current density, displayed as a Tafel plot in Fig. 2i.Tafel plots are key to understanding the physics governing charge transfer reactions for electrochemical systems [11,8].At high current densities the Ohmic resistance of the electrolyte phase often becomes limiting.However, the electrolyte layers in these commercial SOFC cells are very thin and therefore we do not observe significant transport limitations.
The characteristic relaxation time is calculated as the resistance multiplied by the capacitance, τ = R⋅C.For the low frequency electrode process, the chemical capacitance of the bulk phase is dependant on the relationship between defect concentration and the oxygen partial pressure (derived in SI), where the gradient ∂logC chem /∂logp O2 = − 0.25 has  been reported for CGO electrodes using electrochemical and thermogravimetric methods [72,13,73,74].The capacitance as a function of overpotential is displayed in Fig. 2j, where the relationship between the "effective oxygen partial pressure" and overpotential is p O2,eff = p O2 exp(4eη/k B T) [72].The black circles in Fig. 2j represent the total capacitance extracted directly from the DRT calculation, which is the sum of the chemical capacitance (C chem ) of the CGO bulk and interfacial capacitance (C int ) of the double layer at the CGO-gas interface, C total = C chem + C int .The interfacial capacitance is a manifestation of the high defect concentration stored in the first few atomic layers from the CGO surface [75][76][77].Following the derivation of the interfacial capacitance given in the SI, it is possible to correct the total capacitance to yield the chemical capacitance of the bulk of the CGO phase shown, illustrated by the red circles in Fig. 2j.The apparent gradient ∂logC chem /∂logp O2 = − 0.25 was determined to be equal to the expected theoretical value [75,72].This demonstrates that over the full range of current densities, the process responsible for the low frequency DRT peak is consistent, as mass transport limitation would constrain the timescale and would not give the characteristic ∂logC chem /∂logp O2 = − 0.25 relationship.
Overall, the application of the finite Gasussian process algorithm to find the solution to the DRT problem for SOFC test cell data appears to be satisfactory.The signature for bulk chemical capacitance (∂logC chem / ∂logp O2 = − 0.25) was realised, and the truncated data set was able to undergo the mathematical inversion by learning the distribution function, something which is not possible using maximum a posteriori estimation.It must be noted that the low-frequency arc displayed in Fig. 2c is tilted, and it is therefore not ideal to use a series of parallel RC elements to carry out the inversion.A more appropriate model would include Warburg or Gerischer elements to better capture the physical processes as shown by Fu et al., where the low-frequency arc corresponding to oxygen electrocatalysis at SOFC cathodes was conclusively shown to be a Gerischer type element [78].Effendy et al. has demonstrated the use of different circuits with random parameter distributions in a meta-optimisation study known as gEIS [35].Future development of DRT problem solving should therefore endeavour to unify the robust advantages of Gaussian process regression with more physically meaningful models as introduced by Effendy et al..
The fGP-DRT model developed by Maradesa et al. is an improvement on the maximum a posteriori estimation (Tikhonov regularization) method which is frequently used for carrying out mathematical transformation of electrochemcial data [42] The key advantages of using the fGP-DRT formalism is the ability to condition the model to truncated data-sets, automatic and unbiased optimisation, and the absence of a prior which penalizes sharp gradients in the DRT, meaning that peaks are not broadened by regularisation and instead display the inherent distribution of time-sales related to the process of origin.These features novel to the fGP-DRT model mean that mathematical transformations carried out on electrochemical systems (batteries, fuel cells) with processes of similar timescales can be separated and characterised, without prior speculation of the system, with a higher level of precision than equivalent circuit modelling or DRT using Tikhonov regularization.

Conclusion
Gaussian process regression under inequality constraints was used to calculate the DRT of Ni/CGO electrodes as part of a full SOFC stack.This method of mathematical inversion allowed us to avoid using ad-hoc regularisation methods by optimising the kernel hyperparameters through maximising the evidence.By learning the constrained DRT function we successfully convoluted the truncated experimental EIS data and separated four electrochemical processes according to their timescales.The largest of the resistances we assigned to the low-frequency charge transfer process at the electrode-gas interface.When a net current was drawn/applied, the low-frequency DRT peak changed in timescales and resistance, while all other processes remained invariant.
The change in resistance was reported in the Tafel plot which displayed asymmetry between positive and negative overpotentials, and can be utilised to characterise the physics of charge transfer at the electrode-gas interface.The total capacitance was reported as a function of effective oxygen partial pressure.The bulk chemical capacitance was determined by removing the interfacial capacitance from the total capacitance, yielding a relationship ∂logC chem /∂logp O2 = − 0.25 which was found to agree with previous experimental reports.

Author Contributions
Nicholas J. Williams performed conceptualization, methodology, data curation, writing -original draft.Conor Osborne performed methodology, writing -original draft.Ieuan D. Seymour performed validation, writing-reviewing and editing.Martin Z. Bazant and Stephen J. Skinner performed reviewing and editing and supervision.

Conflict of Interest
The authors declare no conflict of interest.

Supporting Information
Figures S1 and S2.Supplementary notes: derivations of Fredholm integral of the first kind, covariance matrix, and the nonequilibrium effects to nonstoichiometry.

Fig. 1 . 2 R
Fig. 1.Visualisation of the Gaussian process distribution of relaxation times optimisation using EIS data (N data points).a, Marginal log likelihood with the covariance matrix Γ is a K + 1 × K + 1 matrix (plus one accounts for the additional hyperparameter σ 2 R ). b, hyperparameter optimisation procedure via Hamlitonian Monte Carlo.c, solution to the GP mean.

Fig. 2 .
Fig. 2. Experimental EIS data and fGP model for the (a) Real and (b) imaginary impedance as a function of measurement frequency, and (c) Nyquist plot where the red dots and solid black line represent the experimental observation and the model fitting, respectively.d, DRT calculated from the posterior, where the solid black line and the grey filled region represent the mean of the posterior and the 95% confidence interval, respectively.e-f, illustrate the EIS at negative and positive current densities, respectively, where the arrow represents the increasing absolute current density going from top to bottom.g-h, illustrate the DRT of the low frequency process for negative and positive current densities, respectively.i, Tafel plot illustrating the absolute current density as a function of overpotential for the low frequency process.j, capacitance of the Ni/CGO electrode in a full cell with varying polarisation to control the p O2,eff where black and red circles represent the observed and corrected capacitance, respectively.All measurements on the Ni/CGO electrode in full cell configuration are 0.5 bar H 2 , 0.5 bar H 2 O, 600 • C.

Table 1
Symbols used