Learning heterogeneous reaction kinetics from X-ray videos pixel by pixel – Nature

Experimental methods

Micron-sized LiFePO₄ platelets were synthesized using a solvothermal method^1,26. LiFePO₄ particles were mixed with sucrose at a ratio of 5:1 and heated to 600 °C for 5 h in a tube furnace under flowing Ar to yield the carbon-coated LiFePO₄. The STXM experiment was conducted at beamlines 11.0.2.2 and 5.3.2.1 of the Advanced Light Source (ALS). LFP particles were dispersed on the working electrode of the liquid STXM imaging platform. X-ray absorption at Fe(II) and Fe(III) L₃ edge was measured in situ and used to obtain the depth-averaged Li fraction map for each particle.

The experimental uncertainty of Li fraction c is estimated by inferring c based on the X-ray absorption at several energy levels with bootstrapping, which returns a standard deviation of c. The average of the standard deviation over all pixels and particles used in the uncertainty analysis is 0.072. A previous estimate of the standard deviation of Li fraction based on reference images is 0.06 (ref. ¹). Therefore, we estimate that the standard deviation of the error of c is σ_ϵ ≈ 0.07.

More details on the synthesis, structure and electrochemical characterization, STXM and AEM imaging, data processing and analysis of experimental uncertainty can be found in Supplementary Information sections 1.1, 1.2 and 6.2.

Chemomechanics

The presence of elastic energy in the free-energy functional arises from the lattice deformation occurring as Li intercalates into the LFP crystal. The elastic constants C_ijkl can be found in Supplementary Information section 2.2. The constitutive law for the chemomechanical coupling relates the misfit strain tensor (or the deformation of the crystal lattice when it is stress free) to the local lithium concentration ε⁰(c). The elastic stress is \({\sigma }_{ij}={C}_{ijkl}({\varepsilon }_{kl}-{\varepsilon }_{kl}^{0}(c))\), in which ε_kl is the total strain.

The misfit strain of LFP is anisotropic—insertion of Li expands the lattice in the a direction and contracts in the c direction. To minimize the elastic energy, the interface between Li-rich and Li-poor phases at equilibrium tends to be aligned in the [101] direction^20,25 (Supplementary Information section 2.2). Therefore, it is crucial to determine the principal directions (a and c) for each imaged particle. By comparing images of relaxed particles and the simulation result based on the same particle geometry, we find that the major axis aligns well with the c axis, and the largest deviation between the two among all examples is found to be 10° (Supplementary Information section 3). Hence we align the c-axis direction to the major axis in our simulations.

Parameterization

For the PDE-constrained optimization, we need to parameterize the unknown homogeneous free energy g_h(c) and exchange current j₀(c). Because the LFP solubility δ is smaller than the experimental error σ_ϵ = 0.07, δ and the unknown constitutive laws cannot be determined from the images. Instead, we use previous knowledge in the literature to inform the binodal composition and the single-phase region (c < δ and c > 1 − δ), that is, we represent μ_h(c) = dg_h/dc as the sum of the ideal mixing entropy \({\rm{ln}}\frac{c}{1-c}\) and Legendre polynomials while subject to the constraint of the experimentally measured values of the binodal compositions (c = δ and 1 − δ)^31,32. lnj₀(c) is parameterized by Legendre polynomials to ensure that j₀(c) is positive. See Supplementary Information sections 2.3 and 6.8 for details of the parameterizations and prior of μ_h(c) and j₀(c).

Because of the presence of both intraparticle and interparticle heterogeneity, we assume that the spatial heterogeneity \(\psi ({\bf{x}})\equiv {\rm{l}}{\rm{n}}k({\bf{x}})\) follows a prior distribution of a Gaussian random field with spatial covariance C(x₁ – x₂) with an offset ψ₀ that follows a normal distribution independently, where x = (x,y) is the spatial coordinate. We set the correlation length l of the Gaussian random field to be the width of a pixel, 50 nm, because spatial variation at the subpixel scale cannot be identified. With this prior distribution, the overall variance of the heterogeneity in space is \({\mathbb{E}}\,[\int {(\psi ({\bf{x}})-{\mathbb{E}}[\bar{\psi }])}^{2}{\rm{d}}{\bf{x}}]\,/\,\int {\rm{d}}{\bf{x}}={\sigma }_{\psi }^{2}+{\sigma }_{{\psi }_{0}}^{2}\), in which we define the mean \(\bar{\psi }=\int \psi ({\bf{x}}){\rm{d}}{\bf{x}}/\,\int {\rm{d}}{\bf{x}}\) and the interparticle variation of the mean is \({\mathbb{V}}{\rm{ar}}\left[\bar{\psi }\right]\approx {\sigma }_{{\psi }_{0}}^{2}\) when l is much smaller than the size of the particle.

Using Karhunen–Loève expansion, ψ(x) can be parameterized as

in which the basis functions ϕ_i(x) are eigenfunctions of the covariance function and λ_i is the corresponding eigenvalue,

and \({Z}_{0},\ldots ,{Z}_{{N}_{{\rm{KL}}}}\) are components of the parameter Z_j used in the optimization for particle j. Supplementary Information section 8.4 discusses in detail the regularization of these parameters given σ_ψ and \({\sigma }_{{\psi }_{0}}\).

PDE-constrained optimization

The optimization is constrained by the PDE model described in the main text and Supplementary Information section 2.1. The coupled PDEs are discretized in space using finite elements (Supplementary Information section 6.3) to obtain a set of differential-algebraic equations (DAEs). The set of DAEs is then solved using numerical differentiation formula with adaptive time stepping and variable order method. See Supplementary Information sections 6.4 and 6.5 for details on the choice of time and spatial accuracy.

The initial condition for c(x) is set as the initial frame of each half-cycle. For mechanical equilibrium, we enforce a zero-traction boundary condition. Because \(\mu ={\mu }_{{\rm{h}}}(c)-\kappa {\nabla }^{2}c-{c}_{{\rm{s}}}^{-1}{\boldsymbol{\sigma }}::{{\boldsymbol{\varepsilon }}}^{0{\prime} }(c)\) (refs. ^7,20,24), a boundary condition for c is needed. The Neumann boundary condition on the normal gradient n ⋅ ∇c, which relies on knowledge of the surface energy between the solid and the electrolyte, cannot be applied in our case. Hence, to avoid determining the gradient from the image, we impose the Dirichlet boundary condition based on boundary values from the images interpolated in time^31,32 (Supplementary Information section 6.9 shows that this results in a smaller error compared with the no-wetting boundary condition, in which n ⋅ ∇c = 0).

The reaction rate R is directly determined by the local overpotential η and concentration c. In simulations, we can specify the voltage Δϕ as a control variable. However, in experimental settings, it is challenging to precisely measure the local overpotential owing to unknown electric potential losses caused by ohmic and contact resistance, as well as the overpotential loss resulting from concentration polarization. Therefore, we adopt an alternative approach by constraining the total reaction rate of each particle, \(I(t)=\int \frac{\partial c}{\partial t}{\rm{d}}V=\frac{{\rm{d}}\bar{c}(t)}{{\rm{d}}t}\), through interpolation of the time-dependent trajectory of the average concentration \(\bar{c}(t)\) observed in the experiment. Using this approach, the local potential Δϕ in the simulation becomes an algebraic variable that satisfies the constraint of the total reaction rate. See Supplementary Information section 6.7 for the analysis of error introduced by the interpolation of \(\bar{c}(t)\). We note that the dependence of the rate on the overpotential can only be tested if a precise measurement of the local potential at the particle can be obtained.

We use the trust-region method to solve the PDE-constrained optimization, perform forward sensitivity analysis to compute the gradient of the objective function and use Gauss–Newton approximation for its Hessian^31,32,52. The optimizer updates p_global and Z_i (i = 1, …, N) concurrently. During each iteration, the forward evaluation of the model and its sensitivities for all half-cycles and particles are performed independently and in parallel. By contrast, the adjoint sensitivity analysis computes the gradient of the objective function much faster than forward sensitivity analysis, however, without the approximation for the Hessian, gradient descent converges to a suboptimal solution (Supplementary Information section 6.6).

Identifiability

We examine the identifiability of the parameters in the model based on a simulated dataset that contains three particles using two approaches. The first approach is HMC sampling of the posterior distribution. At the current level of pixel error σ_ϵ = 0.07, the 99% confidence region of μ_h(c) and j₀(c) is tightly distributed around the truth, indicating that they can be correctly identified in the neighbourhood of the minimum objective function (Supplementary Fig. 31). A higher error σ_ϵ = 0.3 causes the truth to lie outside the 99% confidence region. The second approach is to run the optimizations starting at different initial guesses. Supplementary Fig. 32 shows that, when the optimizer tolerance is sufficiently low, all of the final RMSE from the optimizer can be below the model solver error. However, if the gradient penalty coefficient κ in the chemical potential is unknown and included as a parameter in the optimization, over 80% of the final RMSE is above the solver error (Supplementary Information section 7). Therefore, we use the literature value for κ = 5.02 × 10⁻¹⁰ J m⁻¹ to avoid non-identifiability. See Supplementary Information section 7 for further details.

Regularization

In our study, the number of parameters for p_global, which represents g_h(c) and j₀(c), is substantially lower than the number of pixels (Supplementary Information section 5.2). Sensitivity analysis of g_h(c) and j₀(c) reveals that the variance of the MAP estimate of g_h(c) and j₀(c) scales inversely with the number of particles, denoted as N (Supplementary Information section 8.2 and Supplementary Fig. 34). Also, many particles consist of both charge and discharge cycles, which further decreases the uncertainty of g_h(c) and j₀(c) (ref. ³²). We performed regularization on a simulated dataset of three particles with the same level of noise at σ_ϵ = 0.07, which exhibited a negligible reduction in validation error owing to regularization on g_h(c) and j₀(c) (Supplementary Information section 8.1 and Supplementary Fig. 33). We anticipate an even smaller impact on our dataset with N = 39. Therefore, we set ρ₁ = 0 to minimize the bias of p_global and let p_global be informed by the available data.

In contrast to p_global, regularization on k(x, y) is necessary, particularly for particles with nearly uniform Li concentration field (Supplementary Information section 8.3). As the rate (R/k) increases, leading to a more uniform c(x) profile^7,14,19,53, the uniformity of c(x) can provide information about the magnitude of k(x, y). However, as c(x) becomes nearly uniform, the uncertainty in k(x, y) increases. As a result, the squared error associated with these particles shows insensitivity or lacks local minima with respect to the long-wavelength modes of ψ(x, y), which are constant and smoothly varying in space (Supplementary Information section 8.3). Hence, to constrain the long-wavelength modes, we assume a common prior distribution for ψ(x, y) across all particles. Consequently, the regularization parameter ρ₂ is shared by the Z_i of all particles. Regularization on k(x, y) is also critical to prevent overfitting and leads to a marked reduction in validation error and the error of MAP compared with the truth (Supplementary Information section 8.1 and Supplementary Fig. 33).

We also test using different regularization coefficients for the constant ψ₀ and spatially varying components of ψ(x, y) and find that minimum validation error is reached around when the regularization coefficients are equal (Supplementary Information section 8.4). Because k(x, y) multiplies with j₀(c), we impose a normalization constraint on k(x, y), that is, the average of ψ(x, y) over all particles weighted by area is 0 (Supplementary Information section 4.1).

Cross-validation

We consider two approaches for the cross-validation at different values of ρ₂: (1) allowing p_global to vary as ρ₂ changes or (2) fixing p_global at its MAP values when \({\rho }_{2}={\sigma }_{{\epsilon }}^{2}/{\sigma }_{\psi }^{2}=0.01\), in which σ_ϵ ≈ 0.7 was first estimated using finite differences (see Supplementary Information section 4.2). Studies using simulated datasets show that using a small ρ₂ reduces the bias of p_global (Supplementary Fig. 35). Results on simulated datasets show that the two approaches lead to similar optimal ρ₂ and a comparable reduction in the validation error and MAP error. In the case of experimental datasets, we also perform cross-validation using these two approaches and find that, although the first approach leads to a lower training error, the second approach has a substantially smaller validation error (Supplementary Information section 8.5 and Supplementary Fig. 40). Therefore, we adopt the latter approach, which also allows us to perform k-fold cross-validation independently for each particle, as only Z_i is updated while varying ρ₂ (Supplementary Information section 8.6).

We then performed k-fold cross-validation, training the model on k − 1 half-cycles and validated on the other one. Using the one-standard-error rule⁵⁴, we determine the optimal ρ₂ = 0.88. As mentioned in the main text, the training RMSE at ρ₂ = 0.88 is 6.8%. By contrast, the training RMSE at ρ₂ = 0.01 is 6.0%, suggesting overfitting. Conversely, the training RMSE at ρ₂ → ∞ is 10.6% (Supplementary Table 3 and Supplementary Information section 9.4), that is, the kinetic prefactor is spatially uniform (k(x, y) = 1), which indicates underfitting. At larger ρ₂, the bias in g_h(c) and j₀(c) is also expected to be large because the bias increases with ρ₂ regardless of the number of particles (Supplementary Information section 8.2 and Supplementary Fig. 35). For a visual comparison of the RMSE using different parameters, refer to Supplementary Fig. 50.

Most of the particles with three half-cycles show a strong correlation among k(x, y) trained on the basis of the three training datasets, and all three validation curves have clear minima around the optimal ρ₂. The presence of a few particles that show a lack of consistency among k(x, y) in cross-validation reflects the possibility of unmodelled effects in certain half-cycles in the heterogeneous dataset (Supplementary Information section 8.6). Therefore, to avoid bias, the reported spatial heterogeneities k(x, y) in Supplementary Fig. 57 are obtained by training on all half-cycles with equal weights using the optimal regularization coefficient.

Learning from the uniformity coefficient

Using the uniformity coefficient only as the training data and hence discarding spatial information, we use both the full 2D continuum model on a square and a simplified ‘0D’ model that excludes all spatial correlation (spatial heterogeneity is included in both models but only σ_ψ is an unknown parameter, whereas k(x, y) is a fixed and randomly generated Gaussian random field). HMC sampling shows that the posterior distribution of parameters of j₀(c) and μ_h(c) approximately lies on a manifold that is defined by the contours of the terms in the normalized autocatalytic rate s/R ((lnj₀)′ and \({j}_{0}{\mu }_{{\rm{h}}}^{{\prime} }\); see Supplementary Fig. 13), resulting in coupling between the magnitude of j₀ and the slope of μ_h(c) and between parameters of j₀(c). Despite the coupling, the persistent features that arise from the HMC sampling are the non-monotonicity of μ_h(c) and the asymmetry of j₀(c).

We also include ohmic film resistance, which alters the asymptotic relationship between rate and overpotential at high overpotential, and found that it has a broad posterior distribution over many orders of magnitude and has little correlation with other parameters. Hence, ohmic film resistance is omitted in the full-image inversion.