Inferring solid-state diffusivity in lithium-ion battery active materials: improving upon the classical GITT method
A. Emir Gumrukcuoglu James Burridge Kieran O'Regan Jamie M. Foster




School of Mathematics and Physics, University of Portsmouth, Lion Terrace, Portsmouth PO1 3HF, UK

School of Mathematics and Physics, University of Portsmouth, Lion Terrace, Portsmouth PO1 3HF, UK

About:Energy, Labs Camden - Atrium The Stables Market, Chalk Farm Road, London, England, NW1 8A

School of Mathematics and Physics, University of Portsmouth, Lion Terrace, Portsmouth PO1 3HF, UK
The Faraday Institution, Quad One, Becquerel Avenue, Harwell Campus, Didcot, OX11
0RA, UK


The Galvanostatic Intermittent Titration Technique (GITT) is a ubiquitous method for determining the solid-state diffusivity in lithium-ion battery materials. However, it is notoriously time-consuming and relies upon assumptions whose applicability is questionable. We propose a novel methodology that allows inference of the diffusivity for a  more general class of data that is simpler and faster to harvest. We infer the diffusivity (as a function of stoichiometry) by minimising the residual sum of squares between data and solutions to a spherically-symmetric diffusion model in a single representative active material particle.

Using data harvested from the NMC cathode of a commercial LG M50 cell we first demonstrate that our method is able to reproduce the diffusivities inferred by the GITT, which requires ten days of galvanostatic intermittent titration data. We then demonstrate that our method reliably reconstructs diffusivity using significantly less data. Despite arising from quick-to-measure data, our method more accurately infers diffusivities. This work is a contribution towards developing faster and more reliable techniques in parameter inference for lithium-ion batteries.



Inferring solid-state diffusivity in lithium-ion battery active materials: improving upon the classical GITT method
    Jamie M. Foster
    April 29, 2024
===================================================================================================================





§ INTRODUCTION


The Doyle-Fuller-Newman (DFN) model <cit.> is widely acknowledged as the gold standard in physics-based modelling for lithium-ion batteries and offers a framework for understanding and predicting device behaviour at the level of individual and electrode pairs. Provided that the parameterisation has been performed accurately, the DFN model has been shown to reliably predict the electrochemical response of real batteries in realistic operating conditions; with observed and predicted voltages often matching to within a single percentage <cit.>. So ubiquitous is the DFN model that a wealth of literature exists on various model simplifications and extensions, many of which are discussed in a recent review by Brosa Planella et. al <cit.>. One prevalent class of simplified DFN models are the so-called single particle models (SPM) <cit.>. These are generally accurate for low and moderate C-rate operating (<1C, or thereabouts), in which all the particles in each electrode behave similarly and so the model can be reduced to that of lithium transport within one “representative” particle in each electrode. The distribution of lithium ions within these representative particles are assumed to obey a spherically-symmetric diffusion process.

As alluded to above, the predictive power of the DFN, and related physics-based models, is predicated upon knowledge of parameters that cannot be directly measured but must instead be inferred through physical principles and indirect experimental observations. Therefore, the process of parameterisation – the accurate inference of these parameters – is crucial for developing effective battery models, which in turn are key to advancing device improvement and control <cit.>.

 One of the most influential parameters in physics-based models (especially SPMs) is the solid state diffusivity of lithium-ions within the active materials, D(c), which varies (often rather strongly) with concentration c, and must be inferred. The most widespread inference approach is to collect data with the Galvanostatic Intermittent Titration Technique (GITT), where the electrode is forced with a small constant current for a short time, then allowed to relax until the lithium ions are fully diffused, with the forcing–relaxation cycle repeated until the state-of-charge range of interest is covered. This technique allows each forcing-relaxation cycle to be analysed separately as an independent diffusion problem with constant diffusivity and homogeneous initial concentration. Using the analytic solution of a semi-infinite slab (or, the Sand equation <cit.>), one can determine a diffusion constant corresponding to the small range of concentrations in the cycle <cit.>.

Despite its widespread application, the GITT approach presents several limitations, with one of the most prominent being the practical challenge of data collection.

For typical particle sizes R∼𝒪(1)μ m and diffusion constants D∼𝒪(10^-15) m^2/s, the constant current should be applied for 𝒪(10) seconds to accurately determine D(c). As a consequence, in order to adequately reconstruct D(c) across the relevant range of concentrations, it is necessary to apply low currents several hundred times. Each pulse is followed by a relaxation period, typically lasting around an hour, to restore uniform concentration, a prerequisite for applying the Sand equation to the next pulse. Thus the total data collection time can span weeks.
A second type of issue is that of the questionable embedded assumptions. The GITT inference is carried out using a solution that pertains for semi-infinite slabs, yet the models in which the parameters will be used contain diffusion equations in a sphere. Thus there is inconsistency in between inference and prediction.
The semi-infinite slab approximation is appropriate for pulse durations t_ pulse much shorter than the diffusion time scale t_d = R^2/D. The accuracy decreases as the ratio t_ pulse/t_d increases, with a 5% deviation from the spherical solution at t_ pulse/t_d ≈3× 10^-2 <cit.>. For t_ pulse>0.04 t_d, the Sand equation becomes inadequate for determining D(c) (a detailed discussion is presented in Appendix <ref>).


In addition to the GITT, other techniques to determine the diffusivity exist. Cyclic voltammetry allows estimation of the average (constant) diffusivity <cit.>;
electrochemical impedance spectroscopy <cit.>;

the potentiostatic intermittent titration technique <cit.>;

the intermittent current interruption technique <cit.>.


In this paper, we propose determining the diffusivity by inference from spherical diffusion model (ISDM). In the spirit of SPMs, we solve the fully nonlinear spherical diffusion model in a single representative particle for a known surface lithium flux (directly proportional to the known current applied to/drawn from the cell) and determine the concentration at the particle surface. Via the equilibrium overpotential, U_ eq(c), we calculate the cell voltage, V(t), and by comparing this quantity to the observed cell voltage, we infer D(c).
In contrast with the GITT, our approach offers
the following benefits: (i) it can be applied to data with any current profile as long as the concentration spans the range of state of charge of interest, thereby allowing for very significantly reduced data collection times; (ii) the model being used for inference is compatible with those in which the parameters will subsequently be used to make predictions, thus there are no embedded assumptions and it is therefore valid for all diffusion timescales; (iii) as we will show the parameters inferred by ISDM may offer improved predictive capability for many use cases of interest.

To validate our approach, we use half-cell data harvested from a commercial LG M50 NMC811 cathode from Ref.<cit.>
and show that by using data covering a significantly shorter time, we are able to match, and arguably improve upon, the parameter inference produced by the GITT.
The remainder of this paper is organised as follows. In Section <ref>, we outline the mathematical method underlying ISDM. We validate our approach in Section <ref> using both synthetic and experimental data for GITT current profiles. In Section <ref> we demonstrate the generality of our approach by inferring diffusivity from experimental data generated by applying a constant current. We then compare ISDM to the GITT in recovering a known diffusivity in a full DFN model simulation data in Section <ref>. We conclude with Section <ref> where we summarise and discuss our results.







§ THE ISDM APPROACH


In this section we present the system of equations (reminiscent of a SPM) that lie at the heart of the ISDM. We then discuss the details of the inference approach before moving onto validation in Section <ref>.

For the present study, our aim is to write down the simplest physics-based model that contains only the parameters that one infers with the GITT, thereby ensuring an even playing field in terms of what we are expecting to learn from ISDM vs. GITT. If confronted with data where the influence of other parameters are nonnegligible, then our ISDM methodology remains robust once one the inference model is extended to include those parameters concurrently with diffusivity.




 §.§ A single particle model

We consider a representative particle with radius R, where the molar Li ion concentration c(t,r) on the sphere follows the radially symmetric diffusion equation


    ∂ c/∂ t = 1/r^2 ∂/∂ r(r^2 D(c) ∂ c/∂ r) ,


subject to the initial condition


    c(0,r) = c_0 ,


and boundary conditions


    .∂ c/∂ r|_r=0 = 0 ,
        .-D(c) ∂ c/∂ r|_r=R = j ,


where j(t) is the molar surface flux, describing the rate of Li ions transferred across the particle's surface. Under the usual SPM assumption that the dynamics of each particle are very similar, we can (to a good approximation) apportion the current applied to/drawn from the cell, I(t), equally between particles, such that


    j = ±I/F n 4 π R^2 ,


where n is the number of active particles in the electrode, and F=e N_A=9.6485×10^4 C/mol is Faraday's constant. The sign above is negative for the cathode and positive for the anode, so that positive values of I(t) correspond to a lithiating cathode (and delithiating anode) and hence battery discharging.

In most experiments, one applies I(t) and measures the voltage V(t). Neglecting the reaction overpotential and internal resistance, the dominant contribution to the voltage comes from the equilibrium potential, that is,


    V = U_ eq(c_ surf) ,


where c_ surf(t) = c(t,R)
is the surface concentration of Li ions.



 §.§ Inference of diffusivity and the parameterisation scheme


Let 𝒟 = {t_i, I_i, V_i}_i=1^N  be an experimental data and let us assume that the particle number n, radius R and U_ eq(c) are known from independent measurements.  In order to infer concentration-dependent diffusivity, we parameterise it as D(c ; θ), where θ is the parameter vector.

The spatial derivatives in Eq. (<ref>) are treated using the conservative control-volume method <cit.> thereby reducing the PDE to a system of coupled ODEs for the concentrations at given collocation points, r_i, as functions of time only. We proceed to numerically solve the ODE system for given parameters θ. Using the solution on the particle surface c(t,R;θ), we define the loss function as the mean of the residual sum of squares:


    ℒ(θ) = 1/N∑_i=1^N (V_i - U_ eq(c(t_i,R;θ)))^2 .


The parameters θ can be estimated by minimising the loss:


    θ̂ = θ argmin ℒ(θ) ,


which is equivalent to a maximum likelihood estimate under the assumption of Gaussian errors.

The main challenge in this approach is the parameterisation of D(c). To model diffusivity, there exists an infinite number of choices of functional forms. Parameterisations with few number of parameters generally do not reflect the complexity of the GITT estimation. On the other hand, increasing the number of parameters leads to a high dimensional parameter space which makes it impractical to scan for a global minimum of the loss function.

In this proof-of-principle study, we model D(c) as a piecewise linear function, where the 𝒩 dimensional parameter vector θ is defined as the array of diffusivity values at predetermined knots. To estimate these parameters, we first bin the time series data into 𝒩 partitions. Each partition α has the corresponding data:


    𝒟_α = {t_α i, I_α i, V_α i}_i=1^N_α ,


where N_α is the number of data points in partition α∈{1,…,𝒩}.

We evolve the diffusion equation for each partition sequentially, assuming a series of Fickian diffusion processes with constant diffusivity θ_α in each segment, and solve the corresponding optimisation problem. Defining the local loss in partition α as
 
 
    ℒ_α(θ_α) = 1/N_α∑_i=1^N_α[V_α i - U_ eq(c(t_α i,R ; θ_α))]^2 ,


we estimate the constant diffusivity θ_α as


    θ̂_α = θ_α argmin ℒ_α(θ_α) .


The value of the corresponding knot ĉ_α is defined by evolving the diffusion equation for the partition using the estimated diffusivity θ̂_α and averaging it over time and volume:


    ĉ_α = 3/T_α R^3∫_0^R  ∫_t_α1^t_α1+T_α r^2 c(t, r ; θ̂_α)  dt dr ,


where T_α is the duration of the evolution that falls into partition α. At the end of this procedure, we obtain 𝒩 pairs of (ĉ_α, θ̂_α) values which can be seen as sample points from the full diffusivity function D(c) that is most likely to produce the experimental data in each partition.

It should be noted that this approach is not a comprehensive optimisation strategy in the global sense. For scenarios in which the concentration gradients within the particle are minimal, e.g. a galvanostatic intermittent titration (GIT) experiment with slow constant charge input and frequent relaxation periods, the local optimisation procedure yields highly precise outcomes. On the other hand, for cases with high overlap in c values across partitions, it is necessary to transition to a more global optimisation scheme. In this study, this transition is accomplished by our second optimisation phase, where D(c;θ) is initialised as a piecewise linear function with knots (ĉ_α, θ̂_α) obtained in the first phase; the values of θ̂_α are then improved towards the minimum of the global (or quasi-global) loss function using a numerical gradient algorithm. The updated parameter values are then used to construct the final inferred function D̂(c) = D(c;θ̂). The pseudocode for our two-stage optimisation approach is presented in Appendix <ref>.

For the success of the first (local) optimisation stage, the selection of the number and locations of the partitions is a crucial step. The optimum choice heavily relies on the form of input current and the density of data, such that each partition contains sufficient information about the diffusion of Li ions. For instance, a typical GIT experiment consists mostly of stabilisation data, where the potential does not noticeably change; these data points contain minimal information once the relaxation is well under way. For the specific example of GIT data, a natural choice is to select the partitions such that they each coincide with at least one pulse and relaxation period to democratically distribute the diffusion information across the partitions.

On the other hand, for an experiment where there is no obvious choice for partitioning, e.g. a battery sourced by a continuous constant current (see Section <ref>), one option is to bin the data such that each partition contains an equal amount of useful data points, which we define as points where the change in voltage Δ V is above a certain pre-set threshold value Δ V_ min. [If the current direction remains the same throughout the experiment, an alternative would be to define the partitions such that the transferred electric charge is uniformly distributed across the partitions.]


The parameterisation and optimisation schemes adopted in this study are only a subset of possible approaches for ISDM. For materials with well-characterised diffusion properties, adopting a physics-informed parameterisation
will closely align the model with the underlying physical principles, while reducing the number of inferred parameters θ significantly.
Additionally, the adjoint sensitivity method, which facilitates an effective use of gradient descent variants, can be adopted as an alternative optimisation strategy. For the present study, where a major focus is a GIT experiment with very large number of data points, a global parameterisation with adjoint methods proves to be computationally expensive and can become unstable due to the complexity of the diffusion equation. Instead, the two-stage approach adopted here allows us to analyse individual segments of data separately with parameters localised in each segment, providing a more practical application of ISDM in this initial study.





§ VALIDATION


In this Section, we show that ISDM produces a diffusivity that is consistent with the one obtained with the GITT. For the purposes of demonstration we will use a commercially relevant cathode comprised of NMC811. The measurements have already been reported in Ref. <cit.>. In the interests of brevity we refer the interested reader to <cit.> with a brief summary of the data collection method given below.

The data contains four cycles: the first one (cycle 0) is a charging cycle with constant current for about 10 hours, followed by (cycle 1) a discharge cycle of the same duration, but opposite current. Then there is a GIT charge cycle (cycle 2) consisting of 249 periods of constant pulse (150 seconds) and relaxation (1 hour). The last part of the data is a GIT discharge cycle (cycle 3) with the same parameters for 250 periods, with an opposite-sign current. For this experiment, a 15  mm diameter and 75.6  μ m thickness coin cell is used, and the constant current for all cycles is I=0.78  mA (or C/10).

The average particle radius is R=5.22 μ m.
The delithiation starts at minimum state of charge with fully diffused concentration c = 0.9084 c_ max. Maximum concentration is reported as c_ max = 51765 mol/m^3 <cit.>.

In this section, we establish the diffusivity derived from the conventional GITT procedure, denoted as D_ GITT, as the baseline diffusivity for comparison with the ISDM results. We calculate D_ GITT by applying the Sand equation to each forcing-relaxation period within the GIT charging data.



 §.§ Self-Validation


We start by demonstrating the efficacy of the ISDM method in accurately reproducing diffusivity.
We generate synthetic data by solving the spherical diffusion equation for a particle using D_ GITT(c) and U_ eq(c), both obtained from the conventional GITT approach. The diffusion equation is forced by the same current time series in the GIT charging data, with 249 periods of pulse and relaxation.

We note that this self-validation step focusses solely on the performance of the ISDM method and is not influenced by the validity of the semi-infinite slab approximation or the suitability of the spherical diffusion model for this context; within the scope of this synthetic data analysis, the D_ GITT(c) function is treated as the “true” diffusivity.

We divide the synthetic data into 𝒩=50 partitions, such that each partition contains 5 whole GIT periods, with the last one covering 4.
For the first stage of the optimisation, we create a uniform grid with 10 points for log_10(10^15/ m^2/s D) within the range [-2,2]. We evaluate the loss at these points to identify the minimum. If a minimum is located, the range encompassed by the grid point and its two adjacent points is subdivided into 10 finer points, and this process is repeated until desired accuracy is reached. If no minimum is found within this range, the number of grid points is expanded by a factor of 10, and the range is extended by 1 unit on both sides.
As a result we obtain 50 pairs of (ĉ_α, θ̂_α).

The second stage of the optimisation scheme is intended to change the inferred diffusivity towards the minimum of the global loss. To this goal, we adopt a variant of cyclic coordinate descent algorithm. We construct the stepwise linear D(c;θ) function with parameters θ initialised at the estimates θ̂ from the previous step. For the GIT data, there is minimal overlap in concentration across the partitions; the slow forcing and relaxation periods prevent the growth of c gradients. We therefore employ a quasi-global loss function for each partition, defined as:


    ℒ_α^(β)(θ_α) = 1/∑_i=α-β^α+β N_i ∑_i=α-β^α+β ∑_j=1^N_i[V_i,j - U_ eq(c(t_i,j, R; θ_α))]^2 ,


where β is the number of neighbours (left and right) for each partition α. For inference from synthetic data we use β=2, that is, for each θ_α parameter, we include 5 partitions to calculate the quasi-global loss (<ref>). We then estimate the gradient of ℒ_α^(β)(θ_α) with respect to θ_α  numerically and update the values of θ̂_α, sequentially going over each partition. We repeat this procedure until the global loss ℒ(θ) converges and we obtain our final D̂(c) = D(c;θ̂) from ISDM.


The results of ISDM procedure is shown in Fig.<ref> for

𝒩=50

partitions, along with the “true” value D_ GITT(c) which was used to generate the synthetic data. In the plot, we also show D̃_ GITT(c) which is obtained by applying the traditional GITT to this synthetic data for comparison.



For the synthetic data, the ISDM result D̂(c) is in good agreement with the original D_ GITT(c) used to generate the data.

To quantify the accuracy of the predictions from our model, we first note that an overwhelming portion of the data is determined by the correct choice of the equilibrium potential U_ eq(c). Therefore we expect the predictions to be very close to the data. Instead of directly quantifying how much of the data is explained, we instead quantify how much of the data in excess of a null model can be predicted by the model. To this end, we choose a null model corresponding to D(c)→∞, which predicts instantaneous diffusion, i.e.


    V^ null(t) = U_ eq(c_ av(t)) ,


where c_ av(t) is the average concentration across the particle at time t, and U_ eq is determined from the relaxation values in the GITT forcing-relaxation data.
Defining the voltage data in excess of the null model as


    Δ V_i = V_i - V^ null(t_i) ,


the predicted voltage deviation from the null model for a given D(c) is


    Δ V^ pred_i (D)= U_ eq(c_ pred(t_i, R; D)) - V^ null(t_i) ,


where c_ pred(t,R;D) is the concentration predicted by the specific D(c) model used.

We define the coefficient of determination R^2 as the proportion of variability in the dataset beyond that explained by the null model, accounted for by a spherical diffusion model with diffusivity D(c). Specifically,


    R^2(D) = 1 - ∑_i=1^N (Δ V_i - Δ V^ pred_i(D))^2/∑_i=1^N(Δ V_i - 1/N ∑_j=1^NΔ V_j)^2 .


Note that R^2 ∈(-∞, 1], where 1 indicates a perfect fit, values close to 0 suggest that the model does not improve upon the null model, and negative values reveal that the model adds more variability than the null model.

The coefficient of determination for the ISDM estimate D̂(c) is calculated as follows:


    R^2(D̂) = 0.997 ,


that is, ISDM yields a model that can account for almost all the variability in the data. For comparison we apply the traditional GITT, using the Sand equation on the synthetic data to infer the diffusivity D̃_ GITT(c). This model produces less accurate predictions compared to ISDM, with R^2(D̃_ GITT) =0.800.

Since the true diffusivity is known for the synthetic data, we can also define a measure of accuracy for D(c) itself:


    R^2_D(D) = 1 - ⟨(D(c) -D_ true)^2 ⟩_c/⟨(D_ true(c) - ⟨ D_ true(c)⟩_c)^2 ⟩_c
     ,


where ⟨ f(c)⟩_c denotes the mean of a function f(c) over the range of c. [It is important to note that this measure loses significance if D_ true is flat, i.e. does not exhibit variations. Nevertheless, for the cases examined in this paper, the diffusivities typically have strong dependence on concentration and R^2_D proves to be a reliable metric for assessing the differences of the inferred diffusivities from the true one.]

The metric R^2_D quantifies the proximity of the estimated diffusivity to its true value.
For the ISDM estimate, we calculate


    R^2_D(D̂) = 0.991 ,


while using the diffusivity obtained by applying the traditional GITT to the synthetic data, we get R^2(D̃_ GITT) = 0.882.

This demonstrates the self-consistency of our approach. Since the synthetic data was generated from a single particle spherical diffusion model, it is no surprise that the spherical-diffusion-based ISDM performs better than GITT, which is based on semi-infinite slab approximation.




 §.§ Validation with real data


After confirming that the ISDM approach can reproduce the true D(c) from synthetic data, we move on to analysing the real data from Ref.<cit.>. Specifically, we consider the GIT charging cycle from which we inferred D_ GITT(c) in the previous sub-section using the conventional GITT. We note that the U_ eq(c) function used in this section is also inferred from the GIT data.

After applying the partitioning and optimisation scheme outlined in section <ref>, we obtain the ISDM result for diffusivity. In Fig.<ref>, we show the inferred D̂(c) for 50
partitions.




The coefficient of determination, as defined in Eq.(<ref>), for the ISDM result is computed as


    R^2(D̂) = 0.501 ,


that is, the predictive model accounts for half of the deviations in the data from the null model. On the other hand, using D_ GITT, i.e. the diffusivity obtained by applying the Sand equation, this quantity is slightly lower, with R^2(D_ GITT) = 0.470.


Interestingly, ISDM estimates a diffusivity D̂(c) that is generally lower than the conventional GITT result D_ GITT. The mean relative deviation is
48%.
We show the plot for the percentage relative deviations for all available c values in Figure.<ref>



We see that the typical relative deviation is of order 𝒪(10)%. The largest contribution arises from the marginal values of concentration. Constraining the comparison to c∈ [0.3, 0.9], the mean relative percentage deviation drops down to

34%.
In Appendix <ref>, we showed how the GITT result suffers from errors due to approximating the spherical particle with a semi-infinite slab, which would render the ISDM result more accurate.
On the other hand, the GITT is insensitive to constant contributions to the output voltage during the pulse regime, such as internal resistance effects.



§ INFERENCE WITH CONSTANT CURRENT PROFILE


So far, our application of ISDM has been limited to GIT data, serving mainly for comparison with the traditional GITT-based methods.
On the other hand, the true strength of ISDM lies in its versatility to infer diffusivity across the range of concentration spanned by any current profile.

In this section, we consider the constant current measurements of Chen et al <cit.> and use ISDM to estimate D(c). The data contains a 10-hour measurement where the electrode is charged with a constant C/10 current.

A vital component of our inference methodology is the equilibrium potential U_ eq(c), which, until now, has been derived from GIT data. For ISDM to be used as a swift alternative to the GITT, it becomes essential to determine U_ eq(c) without resorting to the extensive, weeks-long measurements traditionally required. In this section, we perform ISDM using two different equilibrium potentials: U_ eq^ (GITT)(c) obtained from GIT data and U_ eq^ (pOCV)(c) obtained from C/20 charge/discharge data. The latter data set <cit.>
was acquired from a half cell cathode and consists of alternating charge and discharge cycles, repeated 10 times.
We estimate U_ eq^ (pOCV)(c) as outlined in Appendix <ref>, averaging over the 10 charge/discharge cycles.

For partitioning, we distribute the useful data evenly into 50
partitions. We define a useful data where the change in output voltage between two consecutive points exceeds 2× 10^-4 V.

The first stage of optimisation is performed as prescribed in Sec.<ref>, allowing us to obtain constant diffusion values θ̂_α for each partition, corresponding to the average concentrations ĉ_α. For the second stage, the quasi-global loss definition is no longer useful, since the constant current introduces nonnegligible gradients across the particles and the parameters θ̂_α are no longer localised to their corresponding partitions. We therefore use a traditional gradient descent algorithm where we use the global loss function (<ref>) and numerically estimate the gradients with respect to every parameter θ̂_α, then update all parameters simultaneously towards the minimum of the global loss until convergence.

We show the results in Figure.<ref> for inference from C/10 charging
data, compared with the corresponding GITT result, using both choices of the equilibrium potential.



The coefficients of determination for the two case are calculated as:


    R^2(D̂)    = 0.863 ,        (using  U_ eq^ (GITT)) ,
    
    R^2(D̂)    = 0.860 ,        (using  U_ eq^ (pOCV)) .


where for the computation of model predictions, we use the corresponding U_ eq(c) that was used when inferring each D̂(c). For comparison, the traditional GITT result gives  R^2(D_GITT) = 0.0017 (using U_ eq^ (GITT)).

The poor performance of D_ GITT in explaining the constant charge data underlines the inconsistency between the predictive model (single particle spherical diffusion) and the inference model (single particle slab diffusion), which becomes particularly apparent when the system is subjected to a more generalised forcing beyond the periodic pulse-relaxation signal.



§ ACCURACY COMPARISON: ISDM VS GITT


In Sections <ref> and <ref>, we have demonstrated that ISDM outperforms the traditional GITT in explaining the data. However, these assessments are primarily consistency tests that hinge on the validity of the single-particle approximation. Moreover, we have neglected any contribution from e.g. overpotentials, electrolyte. Given these constraints, the superior data fitting of ISDM compared to the GITT is expected, as ISDM is specifically designed to minimise deviations from observed data it is inferred from. Nonetheless, neither analysis provides insight into which of the inferred diffusivities is more correct.
Even in Section <ref>, where the true diffusivity was known, the synthetic data originated from a single-particle model, limiting the analysis therein to only a self-consistency test.


To fairly compare the results of ISDM and the GITT, we require an unbiased data with a known diffusivity. To this end, we used Dandeliion <cit.> for DFN model simulations of a LG M50 NMC811 cathode half-cell with various current profiles, using the parameters from Ref.<cit.>.







In these runs, we use a realistic diffusion function formulated as the exponential of a high-order polynomial, whose full form can be retrieved from the links provided in footnotes <ref>–<ref>.
























The equilibrium potential was fixed to be the fitted function presented in Ref.<cit.>.

We first generate a dataset for a cathode half-cell charging with a GIT forcing-relaxation current [ The Dandeliion run for GIT pulse-relaxation current can be accessed at:

https://simulation.dandeliion.com/legacy/simulation/?id=40751b48-307b-46fb-bcb9-5b3ab79b3b5c] (covering 9.5 days of data). We then apply the traditional GITT with Sand equation to determine diffusivity D_ GITT(c).

We then move on to applying ISDM to the simulated half-cell, which consists of two stages. In the first stage we apply charge then discharge the half-cell with a constant C/20 current [ The Dandeliion run for charge and discharge with constant C/20 current can be accessed at:

https://simulation.dandeliion.com/legacy/simulation/?id=d91db0e5-e0bb-4391-936c-50c69f267c24] (covering 40 hours of data), which allows us to determine U_ eq^ (pOCV) as described in Appendix <ref>.  In the second stage, we generate charging data with C/10 constant current [ The Dandeliion run for charging with constant C/10 current can be accessed at:

https://simulation.dandeliion.com/legacy/simulation/?id=3865ded1-60e1-439b-82b3-65a0ec8c5b0c] (covering 10 hours) and apply ISDM (with 𝒩=50 partitions) to infer diffusivity D̂(c), using the equilibrium potential U_ eq^ (pOCV).

The comparison of the inferred diffusivities along with the true one is shown in Fig.<ref>.



Coefficient of determination for diffusivity R^2_D, defined in Eq.(<ref>), which quantifies how close the inferred functions are to the true diffusivity is calculated for each case as:


    R^2_D(D̂) = 0.883 ,    
    R^2_D(D_ GITT) = 0.777 .


ISDM not only achieves a more accurate approximation of the true diffusivity compared to GITT, but it also requires only 50 hours of data collection, in contrast to the 9.5 days necessary for GITT. This represents a 78% reduction in data collection time.



§ CONCLUSIONS



We have introduced a novel methodology for inferring the concentration-dependent solid-state diffusivity 
using experimentally measured voltage data. We refer to our method as “inference from spherical diffusion model” (ISDM). The ISDM approach offers advantages over the GITT, namely: (i) it does not depend upon the questionable semi-infinite slab assumption inherent in the Sand equation, (ii) it does not require a carefully controlled pulsed current input and lengthy relaxation intervals. The proposed approach is highly versatile, and can infer diffusivity from any current excitation which causes the concentration to pass through the desired range of values, and which is consistent with the assumptions of the SPM. This allows significant savings in experimental time.

We first showed that ISDM is capable of recovering the correct diffusivity from synthetic data (generated by solving the model equations forward with a known diffusivity). This demonstrates that ISDM is consistent, in contrast to the GITT. Second, we applied the approach to real galvanostatic intermittent titration (GIT) data from an LG M50 NMC811 cathode <cit.>. The diffusivity inferred by our method is in good agreement with the results of the GITT. Third, we applied ISDM to C/10 galvanostatic charging data, which is much faster and easier to harvest than GIT data. Once again, we were able to obtain results broadly consistent with the diffusivity inferred from the previous tests. Finally, we directly compared the performance of ISDM and the GITT at inferring diffusivity. We used the DFN model to generate synthetic GIT charge, C/20 charge-discharge and C/10 charge data using a realistic pre-specified diffusivity. Diffusivities were then predicted from this data using both the GITT and ISDM, with the latter method producing substantially more accurate results. 



The most obvious application for our method is in inferring the diffusivity of an electrode at the beginning of life in order to make predictions of the device behaviour in the future. However, owing to the speed and agility of the technique it may prove useful in other ways. It is well-known that degradation, for example in the form of microscale intra-particle cracks, causes the effective diffusivity of the insertion material to reduce with cycle number. Since the ISDM is agnostic to data-type, it could be used in the field to update the model parameters from incoming data as the devices ages, maintaining the predictive capability of the model. These parameter changes might also be used as an advanced indicator of approaching device failure or nonlinear ageing. Finally, we emphasise that although we have chosen to focus on solid-state diffusivity there are no conceptual hurdles to applying the technique elsewhere. In Li-ion batteries it could be used to infer electrolytic conductivity, transference number, or activation energies for the Arrhenius temperature dependence of the diffusivity. It may also prove useful in a variety of other energy capture and storage devices.




We are grateful to Edmund Dickinson for their invaluable comments and suggestions on an earlier version of the draft. JF was supported by the Faraday Institution Multi-Scale Modelling (MSM) project Grant number EP/S003053/1.






§ VALIDITY OF SLAB APPROXIMATION ON A SPHERE



In this appendix, we discuss the validity of the slab approximation implicit in Sand equation to describe a spherical diffusion process. Similar analyses was presented in Ref.<cit.> and in the Supplementary Information of Ref.<cit.>. In particular, we quantify the validity of D(c) obtained using the Sand equation on an actual galvanostatic intermittent titration (GIT) experiment.

We consider the spherical diffusion equation (<ref>) for constant diffusivity D, forced with a constant surface flux j_0, and starting with an initial homegenous concentration c_0. We rescale concentration c, distance r and time t with initial concentration c_0, particle radius R and diffusion timescale t_d, respectively. That is,


    c̃(t,r) = c(t,r)/c_0 ,     x = r/R ,    τ = t/t_d = D t/R^2 .


With these redefinitions, Eq.(<ref>) becomes


    ∂c̃(τ,x)/∂τ = ∂^2c̃(τ, x)/∂ x^2 +2/x ∂c̃(τ, x)/∂ x ,


subject to:


    c̃(0,x)=1 ,     .∂c̃(τ,x)/∂ x|_x=0 = 0  ,    .∂c̃(τ,x)/∂ x|_x=1 = -δ ,


where we defined


    δ = j_0 R/c_0D .


The solution to Eq.(<ref>) is known analytically <cit.>,


    c̃(τ,x) = 1 - δ [3 τ +  5 x^2 - 3/10 - 2/x ∑_n=1^∞sin(α_n x)/α_n^2 sin(α_n)  e^-α_n^2τ] ,


where α_n is the n–th positive root of the transcendental equation α = tanα.

On the surface (x=1), this solution reduces to:


    c̃(τ,1) = 1 - δ [3 τ +  1/5 -2 ∑_n=1^∞ e^-α_n^2τ/α_n^2] .



At τ=0, the sum evaluates to ∑_n=1^∞α_n^-2 = 1/10 <cit.>, which verifies the initial condition c̃(0,1)=1.

One way to determine the behaviour of the finite sum is to approximate it into an integral. We observing that for large n, the recursion approximates to α_n+1≈α_n +π. Extending this relation, we have


    α_n ≈α_1 +(n-1)π ,


which is accurate at 1% level for low n, but becomes more precise for larger values. This allows the approximation


    S = ∑_n=1^∞ e^-α_n^2τ/α_n^2≈∑_n=1^∞ s(n)  ,


where we defined


    s(n) =  e^-[α_1+(n-1)π]^2τ/[α_1+(n-1)π]^2 .


Next, we approximate the infinite sum (<ref>) with the Euler-Maclaurin formula:


    S=∑_n=1^∞ s(n) ≈∫_1^∞ s(x) dx + s(∞)+s(1)/2 + ∑_k=1^∞B_2k/(2k)! (s^(2k-1)(∞) - s^(2k-1)(1)) ,


where B_2k are the even Bernoulli numbers. Since we are only interested in small τ approximation, we only keep terms up to linear order in τ. The integral is straightforward to evaluate and we find:


    ∫_1^∞ s(x)dx =  e^-α_1^2τ/α_1π- √(τ/π)  erfc(α_1 √(τ)) = 1/α_1π - √(τ/π) +𝒪(τ)


As for the correction terms, s(n) and its odd derivatives evaluated at infinity all vanish. The only remaining contributions come from the function and its derivatives at n=1.

The infinite sum of odd derivatives at n=1 are problematic. Around k=12th term, the sum starts to grow rapidly. On the other hand, these terms only give contributions at order τ^0, τ^2 and higher. Since we are only interested in the behaviour of the τ^1/2 term, we will assume that these terms can be resummed. Considering the finite limit for τ→0, we conjecture that


    S ≈1/10-√(τ/π) +𝒪(τ) .


That is, for small τ, the solution to the spherical diffusion equation is  approximated as


    c̃(τ,1)-1/δ = -2 √(τ/π) +𝒪(τ) .


Written in dimensionful quantities, the surface concentration is:


    c_ surf = c_0 - 2 j_0 √(t/D π) + 𝒪(D t/R^2) .


Thus, for t ≪ R^2/D, the spherical diffusion equation has the same surface solution as the diffusion equation for a semi-infinite slab <cit.> which forms the basis of the Sand equation. Although Eq.(<ref>) suggests that the GITT would be applicable for short pulses, we need to determine how small D t/R^2 should be to determine the diffusivity reliably using the GITT.

Moreover, at relatively large τ, the infinite sum decays exponentially, and we are left with


    c̃(τ,1)-1/δ = -1/5-3 τ ,


or


    c_ surf = c_0 - j_0R/D(1/5+3 t D/R^2) ,
         (t∼ R^2/D) .


This expression is particularly interesting since the surface concentration depends linearly on t and the slope is independent of D. Therefore if one employs the GITT in this regime, the inferred diffusion constants will be arbitrary. Again, we need to use numerical methods to determine precisely the time τ = t D/R^2 where this behaviour becomes relevant.

In Figure <ref>, we show a comparison of the full solution and demonstrate the validity of (<ref>) and (<ref>) in the corresponding regimes.




To quantify the accuracy of the approximations, we calculate the relative error by


    Error[%] = 100×c̃-c̃_ approx/c̃-1 ,


where c̃ is the exact solution while c̃_ approx is the approximation. The time dependence of the relative errors are presented in Figure <ref>.



We see that the square-root solution (<ref>) is valid within 5% for τ<0.0032, 7.5% for τ <0.0073 and 10% for τ <0.0132.

On the other hand, we see the relative error in the late solution (<ref>) catches up with the early approximation (<ref>) at τ_ eq = 0.0402, with both errors at 17.25%. For τ >τ_ eq, the linear approximation rapidly improves and the Sand equation is no longer applicable. [A typical value for the 5% accuracy of the linear solution is τ > 1.27 <cit.>. This actually corresponds to the accuracy associated with approximating (c̃-1)/δ≈ -3 τ. Instead, including the intercept (c̃-1)/δ≈ -1/5 -3 τ, the linear behaviour becomes valid at a much earlier time, reaching 5% accuracy at τ>0.0783.]

In Figure <ref> we demonstrate the effect of approximating the sphere as a semi-infinite slab using the GIT data for cathode delithiation from Ref.<cit.>. Here τ = t D/R^2 is compute for each GITT pulse of duration t, with D inferred from the Sand equation. For this data, the value of τ is consistently above the 5% cutoff, occasionally going above the 7.5% line. The sand equation is valid within 10% accuracy in general. Even though GITT might not be accurate, it is still valid since τ stays below the critical τ_ eq line, where Sand equation breaks down for spherical diffusion.




On the other hand, given a relative error Δ c/c in concentration, Eq.(<ref>) implies that the relative error in diffusivity is


    Δ D/D = -2 Δ c/c .


Since the early solution (<ref>) consistently underestimates the drop in concentration, this leads to an overestimation of diffusivity larger by a factor of 2. For instance, a 5% error in c̃ would translate into a 10% error in D.



§ OPTIMISATION SCHEME




In this appendix, we present the pseudocode for the optimisation algorithm.






§ DETERMINATION OF EQUILIBRIUM POTENTIAL USING POCV


The efficacy of ISDM, as with any inference technique, hinges on the accurate determination of the equilibrium potential U_ eq(c). In particular, if the output voltage differs from the correct U_ eq once the system relaxed, the inference process becomes futile, yielding arbitrary outcomes.

When we are presented with a GIT data, it is straightforward to reconstruct the U_ eq(c) function from the voltages at the relaxation points. However, in order to make use of the flexibility of ISDM, it is crucial to find alternative approaches to determining U_ eq(c)  without resorting to lengthy experimental techniques.

In this appendix, we outline the determination of the U_ eq^(pOCV)(c) function used in Sections <ref> and <ref>.
In order to construct the U_ eq, we first note that


    V_ measured = U_ eq(c_ surf) ± |η| ,


where η is the reaction overpotential. The sign is positive (negative) during charging (discharging). Assuming a small enough current such that the charge is quickly dissipated across the particle, we can estimate the equilibrium potential as:


    U_ eq≈1/2 (V_ charge + V_ discharge) ,


such that the effect of the overpotential is cancelled. For the data set used in Sec.<ref>, which contains 10 charge and discharge cycles, we average over all 10 cycles, and show the estimated U_ eq^(pOCV) in Fig.<ref>, along with the U_ eq^ (GITT) obtained from a full GIT charging cycle. The mean squared error between the two functions is 7× 10^-6 within the overlapping stoichiometry range.








