Generating multi-scale Li-ion battery cathode particles with radial grain architectures using stereological generative adversarial networks

Improving Li-ion battery energy-density, power-density, and cycle-life is intricately linked to addressing the climate crisis by enabling transportation electrification and increasing storage capabilities for renewable energy sources. The tremendous progress in Li-ion battery performance, achieved over the past 30 years, has largely stemmed from fine tuning electrode microstructures and electrolyte compositions, without much change in the active material elemental composition¹. Yet, opportunities still exist for controlling cathode sub-particle microstructure to enhance performance^2,3. For example, fabricating radially oriented and elongated grains within LiNi_xMn_yCo_zO₂ (NMCxyz) cathode particles has been shown to significantly increase rate capability as compared to cathode particles that have no preferential grain orientation or internal architecture^2,4,5,6. It is also expected that by designing cathode grain microstructures, chemo-mechanically induced fracture, which perpetuates loss-of-active-material and eventually leads to battery capacity fade, can be reduced³. The present paper seeks to improve the understanding of structure-property relationships between particle-level and grain-level microstructures and observed battery performance. This is achieved by developing a tool that can generate representative microstructures from experimental data. These realistic microstructures are essential in experiment data interpretation and are required for physics-based models that seek to explore particle- and grain-level physics that ultimately determine battery performance and longevity.

There is a need in the community for a tool that can generate representative microstructures for mesoscale battery models to assist in mapping realistic microstructural information to cell-level performance. Physics-based models have been used to help accelerate Li-ion battery electrode development. For example, physics-based models can be used to optimize designs for fast-charge performance^7,8,9,10. Effective optimization using battery models requires relatively fast computation times. To achieve fast computation times, physics-based models typically abstract detailed electrode-level, particle-level, and grain microstructures to “effective” parameters that approximate these intricacies¹¹. Examples of these effective parameters include: secondary-particle diameter, specific surface area, solid-phase diffusion coefficients, and electrolyte tortuosity. While these fast electrochemical models have shown great promise, there is a disconnect between the optimized effective parameters and the underlying complex microstructures that can achieve these performance metrics.

Moreover, geometric complexities not typically included in physics-based models, such as grain orientation, grain size, and secondar particle surface area, can strongly influence cell aging^2,12. For example, secondary particle cracking has a complex relationship with capacity fade–cracks reduce solid-phase diffusion lengths and increase specific surface area¹³, but can isolate active material and introduce unprotected surfaces for side reactions^14,15,16. Additionally, grain orientation, size, and shape greatly influence the apparent solid-phase diffusion coefficient¹⁷.

To bridge the gap between optimized effective parameters, identified in physics-based models and real microstructures, researchers have developed “mesoscale” models that focus on particle- and grain-level microstructures and associated electro-chemo-mechanic physics^{15,18,19,20,21,22}. However, generating statistically representative, realistic particle and grain microstructures as inputs for these mesoscale models is challenging^23,24. Consequently, researchers often approximate geometries with idealized shape, size, and orientations^{15,19,20,21,22,25,26,27,28,29}.

In the present paper, a “multi-scale” particle model is developed to generate representative secondary NMC811 particles comprised of representative grains (primary-particles). The multi-scale approach can be subdivided into two distinct models that can generate virtual morphologies that statistically resemble secondary-particles and grain architectures. A similar two-part, multi-length-scale approach was used previously³⁰. As in this study, the outer secondary-particle shell is generated using a random field on the unit sphere³¹, and the inner grain architecture is generated using a random tessellation. Roughly speaking, common (distance-based) tessellations can be considered as parametric partitions of Euclidean space, such as the Voronoi or Laguerre diagram³². Unlike the previous study³⁰, the proposed grain architecture model is not based on Laguerre tessellations³³, but on a more general class. Furthermore, in contrast to previous approaches, a neural network is implemented to generate the parameters of the tessellations. Importantly, this allows for a (gradient descent-based) fitting procedure of the tessellation model, enabling the generation of 3D grain architectures solely by using 2D planar section information.

Several stochastic approaches for 3D (re-)constructions have been explored previously^{34,35,36,37,38}. In contrast to the present model, many of these methods are based on neural networks that rely on upsampling³⁹ or transposed convolution⁴⁰ techniques to generate discrete 3D voxel representations from 2D pixel images (i.e., slices or projections). These generated 3D voxel representations do not have any underlying constraints and thus are able to generate any morphology and spatial arrangement at the expense of interpretability and computational efficiency.

In contrast to other 2D-to-3D generation methods, the present method does not solely rely on neural networks and the generation of a 3D voxel-based representation. Instead, a neural network and a Poisson point process^41,42 to generate the parameters of a (continuous) tessellation^43,44. More precisely, the parameters of the random tessellation model, i.e., the parameters of the neural network and the point process, are fit by a generative adversarial network (GAN)-based⁴⁵ approach, allowing for the simulation of realistic 3D grain architectures. These virtual grain architectures are then sliced into several 2D images and passed to the so-called discriminator to be compared with measured grain architectures. The output of the discriminator is then used to update the parameters of the grain architecture model to generate more realistic 3D grain architectures. After training, the generated 3D structures are statistically similar to ones measured experimentally. This enables characterization of 3D features like the distribution of morphological features of grains, by using only 2D image data, bridging a major systematic barrier in characterizing particle architectures in 3D. These simulated 3D grain architectures, i.e., continuous tessellations, are especially suitable for simulation applications^{3,18,20,22,46,47}. Furthermore, the fitted parameters of the random tessellation model are interpretable, and by modifying these parameters, similar structures with novel hypothetical structural characteristics can be generated and virtually tested using numerical simulations of electrochemical performance^18,48,49 or cracking behavior³. Alongside this manuscript, a data set of simulated NMC811 particles with inner grain architecture is made publicly available.

Preprocessing of nano-CT data

We present a procedure for modeling the outer shell of NMC811 particles by a random field on the sphere. The fitting of this random field model requires preprocessing the nano-CT data first.

Initially, the 3D grayscale CT image, see Fig. 1a, undergoes segmentation and labeling. This process includes convolving the grayscale image with a discrete 3D Gaussian kernel with standard deviation σ = 1.8, followed by Otsu thresholding⁵⁰ and labeling of the resulting binary 3D image using the watershed algorithm⁵¹. Subsequently, disconnected components of the watershed-labeled image are separated, and small (volume V < 10 voxels) or non-spherical (sphericity Ψ < 0.5, see Eq. (20) for a formal definition) components are merged with their corresponding neighboring components if they exist; otherwise, they are neglected to remove artifacts. Finally, components intersecting the boundary of the sampling window are removed to minimize edge effects. This procedure yields N = 1590 labeled particles in voxel representation, as shown in Fig. 1b, c.

These voxel-based particles are assumed to be star-shaped such that their so-called star point coincides with their centers of masses. Thereby, a star-shaped domain is one where every line segment from the star point to a boundary point lies entirely inside the domain, thus this assumption is much more general than the assumption of ellipsoidal/spherical or even convex particles. A star-shaped particle with star point (cin {{mathbb{R}}}^{3}) can be represented by a radius function (P:{S}^{2}to {{mathbb{R}}}_{+}=[0,infty )) on the unit sphere ({S}^{2}={xin {{mathbb{R}}}^{3}:| x| =1}) given by

where (xi (x)in {{mathbb{Z}}}^{3}) denotes the grid point that is closest to some point (xin {{mathbb{R}}}^{3}). Note that the inner region of a star-shaped representation P(u) of a particle is, by definition, always a subset of its segmented representation. On average, the star-shaped representation accurately captures 96% of the volume of segmented particles. For a visual validation of the star-shape assumption, see Fig. 1d.

To stochastically model the star-shaped particles, the particles radius is modeled as a function by means of random fields on S². Realizations of these models can be considered to be radius functions that describe the outer shell of particles. Therefore, we use so-called (real-valued) spherical harmonic functions ({Y}_{ell m}:{S}^{2}to {mathbb{R}}) for (ell in {{mathbb{N}}}_{0}={0,1,ldots }) and (min {mathbb{Z}}={ldots ,-1,0,1,ldots }) with ∣m∣ ≤ ℓ, which are given by

where ({P}_{ell }^{m}:[0,1]to [0,1]) are the Legendre polynomials and θ(u) ∈ [0, π), ϕ(u) ∈ [0, 2π) are the polar and the azimuthal angles of u ∈ S², respectively. These functions form an orthonormal basis of the Hilbert space L₂(S²) of square-integrable functions on the unit sphere S²⁵². Thus, loosely speaking, the surface of each star-shaped particle can be expressed as a linear combination of spherical harmonic functions. As the order of these spherical harmonic functions increases (i.e., the value of ℓ), they contribute to capturing more detailed and rougher surface features. Since the considered CT data has a finite resolution, it is not meaningful to consider significantly high orders. Therefore, in the following only spherical harmonic basis functions with ℓ ≤ L = 9 are considered. This results in a spherical-harmonics representation (widetilde{P}:{S}^{2}to {{mathbb{R}}}_{+}) of the radius function P, given by

where a_ℓm are the (L + 1)² = 100 coefficients of this basis representation. The computation of these coefficients for a given star-shaped particle’s radius function P, and therefore its approximation in terms of the spherical harmonics basis up to order L, is done by least squares regression. More precisely, the least squares fit between (widetilde{P}(Ru)) and P(u) is computed regarding 642 equidistant evaluation points u ∈ S² on the unit sphere, resulting in 642 equations with the 100 variables a_ℓm. Here, (Rin {{mathbb{R}}}^{3times 3}) denotes a rotation matrix, uniformly chosen at random, which helps to mitigate any anisotropy in the data arising from the voxel-based representation of the CT data or the segmentation procedure. For a visual inspection of the quality of this fit, see Fig. 1e.

Morphology of 2D grain architectures

As mentioned in “SEM and EBSD imaging”, the text data produced by the EBSD measurement was first converted into image data of labeled grains, as shown in Fig. 2. A few pixel values are missing due to overlapping Kikuchi patterns in EBSD imaging. However, these few missing pixels do not significantly affect the grain shape or orientation. These 2D image data resulting from EBSD measurement are from now on referred to as 2D EBSD (planar section) data. The grain architectures observable in these data show a non-uniform morphological orientation, i.e., the grains morphologies exhibit a preferred orientation towards the center of the particle planar section, which is identified as the origin (o=(0,0)in {{mathbb{R}}}^{2}) of the coordinate system. Note that morphological orientation described here is distinct from crystallographic orientation. The proposed model should be able to reflect this behavior, thus, the orientation of individual 2D grains has to be captured quantitatively. For this purpose, the principal components ({v}_{1},{v}_{2}in {{mathbb{R}}}^{2}) of the pixel positions associated with a grain, along with their corresponding eigenvalues e₁≥ e₂ > 0 are utilized. The first principal component v₁ is a direction that maximizes the (empirical) variance of the pixel positions⁵³. In other words, it describes the primary orientation of the grain. The orientation (alpha in [0,frac{pi }{2}]) of a grain is thus given by the smallest angle between the grain’s center c (in the planar section) and v₁, i.e.,

see Fig. 2. For grains with α ≈ 0 the principal component v₁ points towards the origin, whereas an orientation α = π/2 corresponds to principal components that are orthogonal to the direction pointing to the origin. Figure 2 shows the grain orientation distribution per 2D EBSD image. Furthermore, exemplary planar sections are presented in Fig. 2. Notably, for some planar sections, the corresponding EBSD data (outlined with gray boxes in Fig. 2) exhibits no preferred grain orientation. It is conjectured that these images correspond to planar sections located farther away from the particle center as compared to those highlighted with a black box, i.e., the assumption that the particle center is located at the origin of the planar section observed in the EBSD image is not valid. However, the stereological model fitting approach for the inner grain architecture model, introduced in “Stereological tessellation model calibration”, relies on image data of planar sections passing through the particle center, since in such images, grains in each radial distance from the particle center can be observed, a property not present in planar sections at other locations. Consequently, EBSD planar sections represented by gray curves in Fig. 2 are excluded from the subsequent model fitting procedure.

In this section, the stochastic models for both the outer shell and the inner grain architecture, along with their corresponding fitting procedures, are presented. The choice of a two-scale modeling approach is driven by the presumed minimal correlation between the outer shell of NMC811 particles and their inner grain architecture. Initially, the outer shell structure is modeled on a coarser length scale (by means of nano-CT data), before proceeding to the modeling of the inner 3D grain architecture, which is observable on a finer length scale (by means of 2D EBSD data).

Stochastic outer shell model

After the preprocessing of nano-CT data explained in “Preprocessing of nano-CT data”, the next step towards stochastic outer shell modeling is to describe the distribution of the coefficients a_ℓm appearing in Eq. (3) by appropriately chosen random variables Z_ℓm. Due to the assumed isotropy of the underlying image data, the random variables Z_ℓm and ({Z}_{ell m{prime} }) are supposed to be identically distributed for any (m,m^{prime}). For ((ell ,m)ne (ell ^{prime} ,m^{prime} )) with (ell ,ell ^{prime} ne , 0), the random variables Z_ℓm and ({Z}_{ell ^{prime} m^{prime} }) can be assumed to be independent, see⁵⁴. Furthermore, to justify the assumption of independence of Z₀₀ and Z_ℓm for ℓ, m ≠ 0, the empirical distance correlation coefficient (EDCC)⁵⁵ of the observations of a₀₀ and a_ℓm is investigated. Note that the distance correlation coefficient of two random variables is zero if and only if the random variables are independent. Thus, the small values of the EDCC, shown in Fig. 3, justify the assumption of independence. Note that the coefficient a₀₀ predominantly determines the size of a particle, while the coefficients a_ℓm for large values of ℓ primarily capture finer surface features. Thus, the observed increase in EDCC values with increasing values of ℓ corresponds to the trend that, due to the finite resolution of image data, larger particles exhibit progressively finer surface features. Therefore, the model assumption of independent spherical harmonics coefficients might no longer hold when modeling very fine surface features (ℓ ≫ 9), which are not considered in this work.

The coefficients a₀₀ are described by a log-normally distributed random variable Z₀₀ with parameters μ = 0.4, σ = − 29.83, whose probability density ({f}_{0}:{{mathbb{R}}}_{+}to {{mathbb{R}}}_{+}) is given by

On the other hand, for each ℓ > 0, m ∈ { − ℓ, …, ℓ}, the coefficients a_ℓm are described by a univariate scaled Student’s t-distribution, i.e., for some ν, τ > 0, the probability density ({f}_{ell }:{{mathbb{R}}}_{+}to {{mathbb{R}}}_{+}) of the random variable τZ_ℓm is given by

where Γ and (ln) denote the gamma function and the natural logarithm, respectively. The fitted values of the parameters ν and τ are given in Table 1, see also Fig. 4 for a visual validation of the fitted densities ({f}_{ell }:{{mathbb{R}}}_{+}to {{mathbb{R}}}_{+}) for ℓ = 0, 1, …, 9. Both, the selection of the parametric distribution families, from the set of normal, log-normal and Student’s t distributions, and their best fitting parameters are determined using maximum likelihood estimation. This was performed based on the (deterministic) fitted spherical harmonics coefficients a_ℓm that are obtained when fitting the representation given in Eq. (3) to segmented particles observed in data.

The resulting outer shell model {X(u), u ∈ S²} is given by

where Y_ℓm denotes the spherical harmonic function introduced in Eq. (2). Realizations of this model are shown in Fig. 1f. Note that the model has been calibrated on a voxel grid, see Eq. (1). Thus, the generated outer shells are scaled by the voxel side length of 128 nm afterwards.

Stochastic grain architecture model

This section introduces the model used to describe the inner grain architecture of NMC811 particles. The model is based on a so-called generalized balanced power diagram (GBPD)⁵⁶, which is a generalization of the well-known Voronoi³² and Laguerre³³ tessellations. GBPDs are a powerful tool to describe Euclidean space filling structures, such as grains, in a low parametric way⁴³. This low-parametric representation allows for data compression, fast computation, and a simpler and more interpretable grain architecture modeling as compared to voxel-based approaches.

In general, for any integer (din {mathbb{N}}={1,2,ldots }), a GBPD is a decomposition of the d-dimensional Euclidean space ({{mathbb{R}}}^{d}) into non-overlapping subsets (or grains) ({G}_{1},…,{G}_{nu }subset {{mathbb{R}}}^{d}), where for each i ∈ {1, …, ν} it holds that

given (nu in {mathbb{N}}) distinct seed points ({s}_{1},…,{s}_{nu }subset {{mathbb{R}}}^{d}) and markers (M₁, r₁), . . . , (M_ν, r_ν). Here, ({M}_{i}in {{mathbb{R}}}^{dtimes d}) is a positive definite matrix, ({r}_{i}in {mathbb{R}}) is an offset parameter and (| x-{s}_{i}{| }_{{M}_{i}}^{2}={(x-{s}_{i})}^{top }{M}_{i}(x-{s}_{i})) is the metric induced by M_i. Later on, for d ∈ {2, 3}, each set G_i will represent a single grain and consequently will be referred to as such.

The tessellations used within the proposed grain architecture model are further restricted to GBPDs whose matrices M_i can be written as ({M}_{i}={Q}_{i}{D}_{i}{Q}_{i}^{intercal}). Here, Q_i is a basis change matrix to a given basis whose first basis vector is s_i/∣s_i∣ and D_i is a diagonal matrix. This restriction reduces the marker dimensionality from 6 + 1 to 3 + 1 and allows deriving a non-stationary but radial symmetric model, the realizations of which resemble the grain architectures observed in the image data considered in “Morphology of 2D grain architectures”. To derive such a radial symmetric model, the seed points s₁, . . . , s_ν as well as the markers (D₁, r₁), . . . , (D_ν, r_ν) of the tessellations are modeled stochastically. Thereby, the modeling of the seed points is done by a Poisson point process^41,42, i.e., a random point pattern, whereas the modeling of the markers is accomplished by a fully connected neural network as explained in “Stereological tessellation model calibration” below. Note that the neural network only determines markers (D₁, r₁), . . . , (D_ν, r_ν) of seed points s₁, …, s_ν, meaning that it defines the local shape and size of grains but not their spatial distribution or clustering behavior. Thus, this combined tessellation generation approach is not able to exactly reproduce regions of the data it was trained on—instead it is able to reproduce the distribution of geometrical statistics observed in the training data. Additionally, this neural-network-based approach allows for a stereological model fitting procedure. This is crucial since only 2D EBSD data is available to investigate the grain architecture of NMC811 particles.

The following section describes this stereological procedure in more detail.

Stereological tessellation model calibration

The neural network used for marker generation is calibrated by means of a GAN framework⁴⁵. Note that GANs represent a game-theory inspired zero-sum game framework, where two key neural networks comprise the GAN framework: The generator G and the discriminator D. These neural networks compete with each other to achieve a so-called Nash equilibrium⁵⁷, i.e., the generator continually refines its ability to produce data resembling measured samples, while the discriminator iteratively enhances its capability to distinguish between measured and simulated data. Typically, this is realized by an image generating neural network (G:{{mathcal{N}}}to {{mathcal{I}}}), i.e., a neural network that maps some noise space ({{mathcal{N}}}) to some set of images ({{mathcal{I}}}), and a discriminator (D:{{mathcal{I}}}to {mathbb{R}}) that is trained to distinguish between measured and simulated images, see refs. ^34,45 for details. In the present paper, the discriminator tries to assign measured data a label of 1 and simulated ones a label of 0. Using the least squares loss⁵⁸, this results in the following min-max problem:

where (Xin {{mathcal{I}}}) is a random vector following the distribution of measured data, and (Zin {{mathcal{N}}}) is some random vector which serves as input for the generator network G (where Z is sometimes chosen as white noise⁵⁸). The terms ({max }_{G}) and ({min }_{D}) in Eq. (9) refer to the maximum and minimum taken over all possible parameter values of the neural networks G and D, respectively, given that their architectures are fixed. However, in the GAN approach proposed in the present paper, the output of the generator, i.e., a 3D tessellation-based grain architecture, can not directly be compared to measured data, i.e., pixel-based 2D image data, and, therefore, the loss function given by Eq. (9) has to be adapted³⁴. Thus, we introduce two (random) functions, T_2D and T_3D, which compute (random) discretized 2D cutouts of 2D images drawn from X, and 2D cutouts of simulated 3D tessellations drawn from G(Z), respectively. This leads to the optimization problem

In the following, the precise definitions of G, T_2D, T_3D, and D will be explained in more detail. To ensure an efficient GAN training, we assume the superimposed mapping D(T_3D(G)) to be differentiable almost everywhere with respect to the parameters of G. Since G and D are defined by neural networks and thus almost everywhere differentiable, only T_3D(G) has to be analyzed in more detail. To achieve differentiability of T_3D(G), the tessellation representation given in Eq. (8) is not sufficient. Thus, a more general differentiable softmax representation is introduced instead.

For a given sample (s₁, n₁, γ₁), …, (s_ν, n_ν, γ_ν) ∈ W = [−128, 128]³ × [0, 1]² × [0, π] of seed points s₁, …, s_ν and noise (n₁, γ₁), …, (n_ν, γ_ν), the generator G computes markers (M₁, r₁), . . . , (M_ν, r_ν) and a differentiable softmax representation ({{mathcal{T}}}:{{mathbb{R}}}^{3}to {[0,1]}^{nu }) of a tessellation, see refs. ^43,59. This means that the values of ({{mathcal{T}}}) are given by

The vector-valued mapping ({{mathcal{T}}}=({{{mathcal{T}}}}_{1}ldots ,{{{mathcal{T}}}}_{nu })) directly implies a tessellation G₁, …, G_ν, as defined in Eq. (8), by considering the index i of the component ({{{mathcal{T}}}}_{i}(x)) with maximal value, i.e., ({G}_{i}={xin {{mathbb{R}}}^{3}:{{{mathcal{T}}}}_{i}(x)={max }_{j}{{{mathcal{T}}}}_{j}(x)}). Furthermore, note that the differentiable softmax representation ({{mathcal{T}}}) given in Eq. (11) for 3D tessellations can be analogously defined for planar 2D tessellations.

The generation of the sequence (s₁, n₁, γ₁), …, (s_ν, n_ν, γ_ν), i.e., the seed points s₁, …, s_ν and corresponding noise (n₁, γ₁), …, (n_ν, γ_ν) necessary for random marker creation, is done by drawing samples from a homogeneous Poisson point process with some intensity μ > 0, being restricted to the set W = [−128, 128]³ × [0, 1]² × [0, π]^41,42. However, the corresponding markers (M₁, r₁), . . . , (M_ν, r_ν) are generated by a fully connected (FC) neural network, which is denoted as G_net. Since the grain shape and the size is only influenced by nearby grains, we assume that the marker of a seed point only depends on its neighboring seed points. More precisely, for the generation of the marker (M_i, r_i) only information is considered derived from the κ = 12 nearest seed points ({s}_{{i}_{1}},ldots ,{s}_{{i}_{kappa }},{i}_{j}in {1,ldots ,nu }setminus {i}) with respect to the Euclidean distance. Furthermore, to ensure a statistically radial symmetric output of the generator network G_net, its input is modified to get rid of anisotropic information, i.e., the direction of s_i. This is done by rotating the coordinate system in such a way that s_i is always aligned in the direction of the unit vector e₁ = (1, 0, 0). That is, for a given seed point s_i we determine a rotation matrix Q_i given by

where R_s,γ denotes the rotation matrix that describes a rotation around (sin {{mathbb{R}}}^{3}) by an angle of γ ∈ [0, π], ({widetilde{s}}_{i}={s}_{i}/| {s}_{i}|) is the normalized version of s_i, and γ_i is some angle determined by the Poisson point process. Finally, the vectors ({v}_{{i}_{1}},ldots ,{v}_{{i}_{nu }}), which are used as input of G_net for marker generation, are given by rotating the displacements ({s}_{{i}_{j}}-{s}_{i}) by Q_i, i.e., ({v}_{{i}_{j}}={Q}_{i}({s}_{{i}_{j}}-{s}_{i})) for j = 1, …, κ. Note that Q_i is a basis change as described in “Stochastic grain architecture model”.

To further allow the generator network G_net to find a non-deterministic grain architecture for given seed points, and therefore being a suitable stochastic model, seed point specific noise ({n}_{{i}_{j}}in {[0,1]}^{2}) is added. Thus, the input of G_net is given by (({v}_{{i}_{1}},{n}_{{i}_{1}},ldots {v}_{{i}_{kappa }},{n}_{{i}_{kappa }})). Figure 5a shows the fully connected architecture of G_net. It uses ReLu activation functions (f(t)=max {0,t},,tin {mathbb{R}}) in the hidden layers for neuron activations and utilizes scaled sigmoid functions (phi {(t)}_{a,b}=a+(b-a)/(1+{e}^{-t}),,tin {mathbb{R}}) as activation functions in the output layers to ensure that the generated markers are in a reasonable range. This limits the elongation of the resulting grains, which makes them less likely to be disconnected (something that can occur for GBPDs), this helps to prevent unintended artifacts, and improves the generator-discriminator training. In particular, the values of (a,bin {mathbb{R}}) of the channels within the output layer are chosen such that the value of r_i and the diagonal entries of D_i are restricted to [1, 10] and [0, 50], respectively. A more detailed pseudocode representation of this procedure is shown in Algorithm 1.

Algorithm 1

Tessellation Generation

1: procedure (G , ( S=(({s}_{1},{n}_{1},{gamma }_{1}),ldots ,({s}_{nu },{n}_{nu },{gamma }_{nu })) subset {{mathbb{R}}}^{3}times {[0,1]}^{2} times [0,pi ] ))

2: for iin 1, …, ∣S∣ do

3:  ({Q}_{i}leftarrow {R}_{{e}_{1},{gamma }_{i}}cdot {R}_{{tilde{s}}_{i}+{e}_{1},pi })             ⊳ See Eq. (12)

4:  Compute i_j, j = 1, …, κ       ⊳ Nearest neighbor indices

5:  ({v}_{{i}_{j}}leftarrow {Q}_{i}({s}_{{i}_{j}}-{s}_{i}),,j=1,ldots ,kappa )       ⊳ Relative directions

6:  ({D}_{i},{r}_{i}leftarrow {G}_{{{rm{net}}}}({v}_{{i}_{1}},{n}_{{i}_{1}},ldots ,{v}_{{i}_{kappa }},{n}_{{i}_{kappa }}))

7:  ({M}_{i}leftarrow {Q}_{i}^{top }{D}_{i}{Q}_{i})                 ⊳ Markers

8: end for

9: return ({{mathcal{T}}})                   ⊳ See Eq. (11)

10: end procedure

Recall that the sequence (s₁, n₁, γ₁), …, (s_ν, n_ν, γ_ν) ∈ W = [−128, 128]³ × [0, 1]² × [0, π], consisting of seed points s_i, noise n_i, and uniformly sampled angles γ_i, is generated by a homogeneous Poisson point process with some intensity μ > 0. This point process corresponds to the noise generation random variable Z considered in Eq. (10). A Poisson-distributed random variable is used to determine the number ν of sampling points. Then, the locations of the respective points are independently and uniformly sampled within the window W = [−128, 128]³ × [0, 1]² × [0, π]. Note that this kind of point process has only one parameter, the intensity μ, i.e., the expected number of points per unit volume. Since each of these points correspond to a grain, the parameter μ has to be fitted such that the average numbers of visible grains of generated and measured data in planar 2D cutouts match. However, this cannot be accomplished a priori, since the number of visible grains in a 2D planar section of a 3D tessellation is influenced by the number of seed points as well as their markers, as can be seen in Fig. 6. To overcome this problem, the intensity μ is adopted iteratively while training, i.e., the intensity μ_e+1 of the Poisson point process at training epoch (e+1,,ein {mathbb{N}}) is defined as follows:

where (N{G}^{{{rm{gt}}}},N{G}_{e}^{{{rm{gen}}}} , > ,0) correspond to the average number of observable grains per unit area in planar 2D cutouts of ground truth data and generated data at training epoch e, respectively. To avoid large jumps of μ_e, which would hinder the learning process of both generator and discriminator, the intensity μ_e+1 is chosen as a convex combination of the previous intensity μ_e and the desired intensity ({mu }_{e},(N{G}^{{{rm{gt}}}}/N{G}_{e}^{{{rm{gen}}}})). The initial value of this procedure is set heuristically to μ₁ = 0.0018.

As mentioned above, the output of the generator G, i.e., ({{mathcal{T}}}) as defined in Eq. (11), can not be directly compared to measured EBSD data. One reason for this is the circumstance that there are several pixel positions in the EBSD data to which no grain is assigned at all. Moreover, measured data is only present as planar 2D sections, therefore planar sections have to be extracted from generated 3D grain architectures to ensure comparability.

To address the issue of unassigned pixels, a differentiable representation ({{mathcal{T}}}) (see Eq. (11)) of the restricted GBPD with d = 2 is fitted to the EBSD data using the method described in⁴³. This fitted 2D tessellation is referred to as “ground truth data”. Figure 7 shows an exemplary fit. The approach of considering a 2D tessellation fit to the experimental EBSD data as ground truth data has at least three advantages. First, it provides a meaningful prescription for missing measurement information (inner white pixels in Fig. 7, left), while maintaining the overall shape of grains. Second, it ensures that the planar sections with which the generator attempts to imitate can indeed be generated, i.e., the ground truth data can be represented by GBPDs. Without assurance that the ground truth data can be represented by GBPDs, it is possible that the discriminator differentiates between measured and simulated data solely based on a difference of the representation, which would offer no insightful feedback to guide the generator’s improvement. Third, it facilitates effective data augmentation during training, i.e., by modifying the seed points and markers of the fitted 2D tessellation slightly, similar but new 2D grain architectures can be generated as additional training data. Note that the distribution of the random variable X considered in Eq. (10) describes the distribution of 2D training data, i.e., the data that arise from the ground truth data by this kind of augmentation.

Despite the fact that now both ground truth and simulated data are represented as differentiable tessellations, their dimensionalities are not yet consistent (being equal to 2 and 3, respectively). To overcome this problem and to generate input that is feasible for the discriminator D, i.e., a matrix of fixed size, randomly located (2D) pixel-based cutouts of these tessellations are computed.

For this purpose, a 32 × 32 grid of equidistant points is chosen, say ({1,ldots ,32}times {1,ldots ,32}subset {{mathbb{R}}}^{2}). The tessellation ({{mathcal{T}}}) of the restricted GBPD with d = 2 is then evaluated at these grid points, resulting in a matrix I with 32 × 32 × ν entries. Such a matrix is often referred to as a multichannel 2D image, with dimensions representing the x-axis, y-axis, and channels, respectively. For each (x, y, i) ∈ {1, …, 32}² × {1, …, ν}, the entry I_x,y,i of I corresponds to the probability that the grid point (x, y) belongs to grain i. To reduce the number of channels of I to a (sufficiently small) fixed value, the set of channels is truncated to those channels which belong to grains being most present among the considered grid points. In the following, for the case of a 2D tessellation ({{mathcal{T}}}), this procedure is described in more detail. First, the originally chosen 32 × 32 grid {1, …, 32}² of equidistant points in ({{mathbb{R}}}^{2}) is shifted and rotated at random, i.e., a (random) grid of equidistant points is determined by drawing a radial distance p and a polar angle β uniformly from the intervals (0, 80] and [0, 2π], respectively. Then, for each (x, y) ∈ {1, …, 32}², the corresponding grid point (({x}^{{prime} },{y}^{{prime} })) of the transformed grid is given by

Evaluating the tessellation ({{mathcal{T}}}) at these points gives the (non-truncated) image I = (I_x,y,i), where the entry I_x,y,i of the matrix I at (x, y, i) ∈ {1, …, 32}² × {1, …, ν} is given by

In Fig. 7c the three steps of this cutout computation procedure are illustrated, where two different coordinate systems are used, one for the original pixel positions (x, y) ∈ {1, …, 32}² and one for the grid points ((x^{prime} ,y^{prime} )) introduced in Eq. (14). The procedure starts with a cutout, whose center lies on the x-axis, i.e., p = β = 0. Grains located in this cutout that are strongly orientated towards the center of the planar section (i.e., (alpha < < frac{pi }{2}), see Eq. (4)) are elongated in the direction of the x-axis. In a second step, the cutout is shifted in the direction of x-axis, according to the selected value of p (where p = 70 in Fig. 7). As before, grains within the cutout which have an orientation (alpha < < frac{pi }{2}) show elongations in the direction of the x-axis of the cutout. Finally, in a third step, the cutout is rotated by an angle of β ∈ [0, 2π]. This does not only displace the cutout but also changes its orientation. Consequently, grains that are located in this cutout having an orientation of (alpha < < frac{pi }{2}), are orientated in the direction of the (new) ({x}^{{prime} })-axis of the cutout. This procedure is performed to generate consistent cutouts which are feasible for the discriminator D and, simultaneously, to be able to use a wide range of cutout positions within the ground truth data for training.

To achieve a fixed (small) number of channels necessary for the discriminator input, the ν channels of I are sorted by the sum of their entries in descending order. Then, all but the first 12 channels are omitted, i.e, only the channels corresponding to the 12 most prominent grains regarding the grid points are retained, in order to retain the channels that contain the most information of grain morphologies. This (random) procedure will be denoted by T_2D in the following.

A similar procedure, analogous to that described above for T_2D, is used to determine 32 × 32 × 12 matrices which represent (random) 2D cutouts of the simulated 3D tessellations for a given number ν of seed points. More precisely, for a uniformly drawn radial distance p ∈ (0, 80], a matrix I = (I_x,y,i) with 32 × 32 × ν entries is computed, where

Due to the rotational invariance of both the Poisson point process of seed points and the marker generation, and the ability to create infinity many simulated grain architectures, no rotation matrix is needed in Eq. (16) for the generation of grain locations, unlike the procedure for T_2D described above. The resulting matrix I (a multichannel 2D image) is then channelwise sorted by the sum of their entries and cropped to maintain a constant number of channels, using a procedure similar to the one already described. This (random) procedure of computing 32 × 32 × 12 multichannel 2D images from artificially generated 3D tessellations is denoted by T_3D. Recall that both, T_2D and T_3D pay attention to differentiability and to preserve a consistent orientation of the cutout with respect to the origin of the tessellation. This is crucial to enable the discriminator to capture preferred grain orientations with respect to the particle center.

The discriminator D, which is trained to decide whether a 32 × 32 × 12 cutout image originates from ground truth or artificially generated data, is structured as a convolutional neural network (CNN). Note that CNNs are a neural networks that can handle spatial dependencies in their input and are therefore especially useful in computer vision tasks such as image classification⁶⁰. To avoid overfitting issues, the discriminator D considered in the present paper has a rather simple architecture. Namely, the discriminator consists of stacked convolutional layers, batch normalization layers, and ReLU activation functions. These components are responsible for feature extraction, learning acceleration, as well as neuron activation and deactivation. Additionally, max-pooling layers are used for dimension reduction, and dense layers are used to process the extracted features to a single output. Figure 5b shows a detailed representation of the network, and the number of features per layer.

The GAN-based training process involves iteratively training the generator G and discriminator D through gradient-descent and a non-gradient-descent based optimization of the point process intensity μ, as described in Eq. (13). The gradient-descent based training uses an Adam optimizer with a learning rate of 10⁻⁴ and gradient normalization⁶¹. The whole training procedure is done over 200 epochs with 100 steps per epoch and a batch size of 64 and 128 for generator and discriminator training, respectively. While training the discriminator the pixel values of its input images (see Eqs. (15) and (16)) are rounded to the closest value in {0, 1} to suppress features that arise from the dimensions of the underlying artificially generated and ground truth tessellations, respectively. To address the issue of overfitting the discriminator, the update of the discriminator’s weights is skipped during training, provided that its current mean squared error given by Eq. (10) is below 0.3². Furthermore, training data augmentation is achieved through marker modification, i.e., the diagonal entries of the matrices D_i introduced in “Stochastic grain architecture model”, and the additive marks r_i of the fitted 2D tessellation are modified through uniform augmentations up to 20%, i.e., (frac{| x-x{prime} | }{| x| } < 0.2), where x is equal to r_i or a diagonal entry of D_i, and (x{prime}) is its augmented version. This augmentation of training data results in more diverse grain architectures that still adhere to the constraints of the restricted GBPD representation, see Fig. 7. However, this is contrary to augmentations achievable through conventional techniques of image data processing⁶². Figure 8 and Algorithm 2 provide an overview of the training procedure described above, where the pseudocode is presented in such a way as to improve readability and should not be considered computationally efficient. For a more efficient implementation, the seed point generation and, thus, the 3D grain architecture generation, as well as the 2D grain architecture augmentation can be restricted to a small sampling window that contains the observed planar section.

Algorithm 2

Training Procedure

1: procedure Train

2:   lr ← 0.0001                         ⊳ Learning rate

3:   bs ← 128                           ⊳ Batch size

4:   μ₁ ← 0.0018                     ⊳ Point process intensity

5:   for e in 1, …, 200 do                       ⊳ Epochs

6:   for step in 1, …, 100 do

7:    Draw S₁, …, S_bs/2 from Z with intensity μ_e      ⊳ Seeds, noises and angles

8:    ({{rm{loss}}}leftarrow frac{2}{bs}mathop{sum }_{i}^{bs/2}{(D({T}_{{{rm{3D}}}}(G({S}_{i}))))}^{2})

9:    Maximize loss with respect to G_net and lr  ⊳ One step of Adam gradient descent

10:    Draw S₁, …, S_bs from Z with intensity μ_e       ⊳ Seeds, noises and angles

11:    Draw x₁, …, x_bs from X        ⊳ Augmented measured planar sections

12:    ({{rm{loss}}}leftarrow frac{1}{bs}mathop{sum }_{i}^{bs}{(D(xi ({T}_{{{rm{2D}}}}({x}_{i})))-1)}^{2} +D{(xi ({T}_{{{rm{3D}}}}(G({S}_{i}))))}^{2})          ⊳ (xi (a) = mathop{{arg},{min}}_{bin {{mathbb{N}}}^{32 times 32 times 12}}| a-b|) (rounding)

13:     if loss > 0.09 then             ⊳ Avoid discriminator overfitting

14:     Minimize loss with respect to D and lr    ⊳ One step of Adam gradient descent

15:     end if

16:     ({mu }_{e + 1}leftarrow {mu }_{e}(0.5 + 0.5frac{N{G}^{{{rm{gt}}}}}{N{G}_{e}^{{{rm{gen}}}}}),)      ⊳ Iterative intensity adaption, see Eq. (13)

17:   end for

18:   end for

19: end procedure

Multi-scale model

Combining the outer shell model and the grain architecture model follows the approach presented in ref. ³⁰. More specifically, to get a simulated NMC811 particle, an outer shell and a grain architecture are independently drawn from the models introduced in “Stochastic outer shell model” and “Stochastic grain architecture model”, respectively, and overlaid on a 3D domain. Next, grains whose centers of mass are not located inside the outer shell are removed. Figure 9 illustrates this procedure and shows 3D renderings of samples drawn from the outer shell model and the grain architecture model, as well as their combination to a virtual NMC811 particle.

Validation of the modeling approaches

In this section, both models are validated using interpretable descriptors that were not part of the training process. Note that due to the limited number of particles measured by EBSD, the training data will be used to compute these interpretable descriptors that will serve as ground truth during validation. Note that due to the choice of the models these are not able to reproduce training data and thus a validation data split is not necessary, see “Stochastic grain architecture model”.

Validation of the outer-shell model

To evaluate the outer-shell model introduced in “Stochastic outer shell model”, four different descriptors are chosen to characterize the 3D morphology of ground truth and simulated microstructures, respectively. These descriptors include the diameter ({d}_{max }), the compactness γ, the aspect ratio A_R, and the sphericity Ψ of particles⁶³. More precisely, for a particle represented by its radius function (P:{S}^{2}to {{mathbb{R}}}_{+}), the values of these descriptors are given by

where A( ⋅ ) is the surface area, V( ⋅ ) the enclosed volume, and ∂Conv( ⋅ ) the boundary of the convex hull. Note that the compactness γ measures the degree of convexity of a particle, whereas the aspect ratio A_R describes the particle’s elongation, and the sphericity Ψ measures the roundness of a particle, i.e., its similarity to a sphere. Figure 10 (a) shows histograms of these descriptors for both ground truth and simulated particles.

In the outer-shell analysis, simulated particles having a volume of less than 10 voxels, i.e., less than 10 ⋅ (0.128 μm)³, are neglected. These volumes are omitted for further analysis, since the described segmentation procedure generates only particles with larger sizes. The histograms of the descriptors ({d}_{max }), γ, A_R, and Ψ of the resulting particles show nice agreement with corresponding histograms computed from nano-CT data, see Fig. 10. It is noteworthy that the model does not produce particles with extremely low sphericities, as observed in the segmentation process (see the values for Ψ ≤ 0.6 in Fig. 10). Further investigations are required to check whether these particles can be attributed to an imperfect segmentation, and consequently, whether their absence in model realizations is important.

Recall that the relatively small EDCC values displayed in Fig. 3 provided the basis for assuming that the model’s random spherical harmonics coefficients are independent of each other. However, an increasing trend of the EDCC along with the order of the basis functions can be observed. It is assumed that this trend can partially be attributed to missing fine structures in particles with a small diameter due to the discrete voxel-based resolution. Additionally, the small expected values of ∣Z_ℓm∣ for ℓ > 4 suggest that these coefficients can be neglected, without decreasing the model quality significantly. A systematic investigation of the influence of the maximum spherical harmonics order L on the model quality will be the subject of future work.

Validation of the grain architecture model

The model introduced in “Stochastic grain architecture model” for the inner 3D grain architecture of NMC811 particles is based solely on 2D EBSD data. Since only 2D planar section data of ground truth grain architectures is available, the validation of simulated grain architectures is done by a statistical comparison of pixel-based planar sections. For this comparison, the distributions of four microstructure descriptors are evaluated: the distribution of the size A_cross of planar grain sections, the distributions of their elongation η and their orientation α, and the chord-length distribution⁶⁴ of the ensemble of planar grain sections. Note that the size A_cross of a planar grain section is given by the number of pixels occupied by the given grain section. The orientation α is defined in Eq. (4), and the elongation η is given by

where e₁≥ e₂ > 0 are the eigenvalues of the principal components of the planar section of the given grain, see “Morphology of 2D grain architectures”. The chord-length distribution is given by the distribution of the lengths c of all chords of the planar section, where a chord of a pixel-based planar section of grains is a sequence of consecutive aligned pixels that belong to the same grain and that can not be extended further without containing pixels belonging to other grains. Figure 10b shows histograms of these descriptors for both ground truth and simulated data, which are quite similar to each other. However, the elongation η shows slightly larger values for simulated data, which means that the grains observed in the planar sections of the stereologically fitted model are slightly more elongated than those in the ground truth data. A similar behavior can be observed for the grain orientation α, where simulated grain architectures show a slightly more preferred orientation as compared to the ground truth. Nevertheless, the histograms computed from ground truth and simulated data, respectively, show nice agreement, which indicates that planar 2D sections of grains observed in EBSD data are well captured by the stochastic grain architecture model introduced in “Stochastic grain architecture model”. This gives strong evidence that effective properties such as ion transport or cracking behavior are also well captured. However, the investigation of these based on numerical simulation is beyond the scope of this paper and will be subject of future work.

Computational speed

An important factor in practical applications is computational speed, particularly for generating a large number of model realizations to compute effective properties for homogenization. To accurately assess this, it is necessary to distinguish between two key components: the computation time required for model calibration and the computation time for generating virtual 3D particles with the calibrated model. Since the former step only has to be performed once per data set, the latter one is of particular interest. All computation times were measured on a personal computer equipped with an AMD Ryzen 9 5900X CPU and an NVIDIA RTX A4000 GPU. For data already available in a format suitable for training, i.e., spherical harmonics coefficients determined for segmented outer shells or parameters of a deterministic 2D tessellation fitted to 2D image data “Preprocessing of nano-CT data” and “Stereological tessellation model calibration” model calibration takes less than one minute for the outer shell model and approximately half a day for the grain architecture model. In terms of simulation speed, generating a realization of the outer shell model takes 2.37 × 10⁻³ s per particle shell, and for the grain architecture model, it takes 7.06 × 10⁻⁶ s per grain. These times are averaged over 10⁶ generated outer shells or grains, respectively. It is important to note that these times correspond to generating analytical representations of the morphologies. In practical applications, one has to account for additional computation time, necessary for transforming these analytical representations of morphologies into discretized representations, e.g., discretized 3D image or surface meshes, which are suitable for further computations, such as numerical simulations. The time required for this transformation depends on the desired resolution and can be significant.

Discussion of model assumptions

The trained grain architecture model, which was fitted and evaluated by means of planar 2D sections, can be used to generate realistic 3D grain architectures of NMC811 particles and to investigate structural properties of these 3D morphologies. For example, Fig. 11 shows how the distribution of grain orientations observed in planar sections depends on the section’s position (i.e., its cut height). It can be observed that the grains that are visible in planar sections near the origin exhibit an increasingly preferred orientation, in contrast to more distant planar sections, which show an almost uniform distribution of grain orientations. This observation coincides with the assumption stated in “Morphology of 2D grain architectures” that some planar sections observed in EBSD data may not pass through the particle center and therefore were neglected for training purposes. However, the GAN-based training of the grain architecture model could be adapted such that planar sections that are taken farther away from the particle center can also be used. Nevertheless, such modifications would necessitate providing the discriminator with additional information about the observed data. Such additional information could include the distance of the planar section to the particle center. This adjustment would enable the discriminator to evaluate the plausibility of observed grain elongations or preferred orientations based on the position of the planar section within the particle. However, this type of information, i.e., the relative location of the planar section with respect to the particle center, would have to be acquired by laboratory measurements.

The proposed multi-scale model currently accounts only for pore-free polycrystalline particles. However, the investigation of polycrystalline particles with pores, particularly those resulting from cracking during cycling (i.e., degradation), is of significant interest. The implementation of an additional scale for the proposed model⁶⁵ to stereologically simulate internal pores and cracks caused by degradation will be explored in future work.

In the proposed multi-scale model, we assumed that the inner grain architecture can be modeled independently of the outer shell. While this may be true for large particles, it will not be true for particles with very small radii, i.e., particles that consist only of hundreds of grains. To investigate this issue further and, consequently, adapt the modeling approach, EBSD data of differently sized particles has to be acquired. Importantly, these rather small particles account for a small mass percentage of an electrode and, thus, their inclusion/exclusion is not expected to significantly influence predicted electrode performance.

The multi-scale modeling approach described here enables realistic 3D particle geometries with sub-particle grain detail. The 3D sub-particle grain morphologies were generated by only using 2D EBSD data, facilitating characterization of 3D grain morphological properties from 2D experimental information. Representative particle generation capabilities are not only required for characterizing otherwise challenging 3D morphological properties but are also required for high-fidelity cathode degradation models that require significant test geometries to mitigate stochastic effects resulting from particle geometry variations. Notably, it is significantly easier and faster to simulate chemo-mechanics on hundreds of generated geometries, as compared to relying solely on a sparse set of empirical image data. A bulk of such simulated NMC811 particles with inner grain architectures, suitable for numerical simulations, is made publicly available. Previous cathode degradation modeling efforts focused mostly on NMC532 particles with randomly distributed grains, which was primarily due to the availability of full 3D imagery³. Now, with the capability of inferring 3D geometries of grain architectures with curved grain boundaries from only 2D EBSD slices, it is possible to model NMC811 degradation and consider radially oriented grain geometries, which will be the focus of future work. Additionally, the proposed method can enable quicker development of generation methods for new and emerging electrode chemistries.

In the future, we plan to expand the model not only to capture the geometry of the grains but also to describe their crystallographic orientation as well as pores within the grain architecture that arise during cycling⁶⁵. These modifications will be the basis of future numerical-simulation-based analysis of microstructure-property relationships in NMC811 particles during battery cycling, i.e., different degradation levels. Additionally, we intend to adjust the fitting for the outer shell model such that calibration can be performed from 2D image data to achieve a fully stereological multi-scale model.

Overview of experimental approach

To supply sufficient high-quality image data on particle and grain architectures, a multi-modal approach was deployed. Electron backscatter diffraction (EBSD) was used, as previously described in ref. ⁶⁶, to facilitate grain and grain boundary segmentation, as well as populate the segmented grains with local crystallographic orientation information. Note, there is a distinction between crystallographic and morphological grain orientation; crystallographic orientation refers to the orientation of the c-axis of the NMC811 crystal, whereas morphological grain orientation refers to the direction of the major axes of the best fitting elongated spheroid shapes of the grains with respect to the particle center. Conducting EBSD on a small number of particles provides a view of hundreds of grain planar sections. However, extending EBSD to 3D via focused ion beam (FIB)-EBSD and acquiring morphological data on hundreds of full particles is not yet practical experimentally. Instead, the (secondary-) particle morphological information, or outer shell, was captured in detail for many particles by X-ray nano-computed tomography (nano-CT) at the cost of a lack of inner grain architecture information. The outer shell and sub-particle grain architectures are independent and therefore must each be quantified with a sufficient number of samples to be representative. The combination of 2D EBSD and 3D nano-CT facilitated detailed information on the sub-particle grain features and full particle outer shell morphological detail, respectively, with sufficient volume to achieve representivity in both.

Materials

A pristine sample of NMC811 electrode was used for this work. The electrode consisted of 96 wt% Targray NMC811, 2 wt% Timcal C-45, and 2 wt% Solvay 5130 PVDF binder. The Al foil was 20 μm thick, and the coating thickness was 58 μm with a porosity of 33% and areal capacity of 3.07 mAh cm⁻². These cathodes have been extensively studied in fast-charge applications for electric vehicles^2,12,67,68. The preferential radial orientation of these particles is argued to significantly increase these cathodes rate performance and cycle-life resiliency¹².

X-ray nano-computed tomography (nano-CT)

For X-ray nano-CT imaging, cylindrical pillars of ca. 90 μm diameter were prepared using a micro-machining laser ablation approach described previously in ref. ⁶⁹. The pillars were then imaged in a Zeiss Xradia Ultra 810 X-ray nano-CT system in Large Field of View (LFOV) absorption mode with binning 2, giving a voxel size of 128 nm. The field of view was 64 μm  × 64 μm. 1601 images were acquired for the reconstruction of each tomogram.

SEM and EBSD imaging

The NMC811 cathode sample was argon-milled with a JEOL CP ion beam planar section polisher (JEOL, USA). This provided a wide (ca. 2 mm) smooth planar section of the NMC811 electrode. Scanning electron microscopy (SEM) and EBSD images were taken on several planar sectioned particles using an FEI Nova NanoSEM 630 equipped with an EBSD detector (EDAX, USA). EBSD data was collected using step sizes of 50 nm rastered across the surface of particle planar sections. EBSD data was processed with OIM Analysis v.8 (EDAX, USA). Diffraction patterns were fit to a trigonal crystal system (space group R-3 m) with a = b = 2.875 Å and c = 14.248 Å to obtain the orientation of the crystal at each, where a, b, c denote the side lengths of the rhombohedral unit cell associated with the crystal’s atomic lattice. The software produced text files containing a spatially resolved confidence index, image quality (IQ), and Bunge-Euler angle data. These data were then converted into images, where individual grains are labeled separately.