Artificial generation of representative single Li-ion electrode particle architectures from microscopy data

Accurately capturing the architecture of single lithium-ion electrode particles is necessary for understanding their performance limitations and degradation mechanisms through multi-physics modeling. Information is drawn from multimodal microscopy techniques to artificially generate LiNi0.5Mn0.3Co0.2O2 particles with full sub-particle grain detail. Statistical representations of particle architectures are derived from X-ray nano-computed tomography data supporting an ‘outer shell’ model, and sub-particle grain representations are derived from focused-ion beam electron backscatter diffraction data supporting a ‘grain’ model. A random field model used to characterize and generate the outer shells, and a random tessellation model used to characterize and generate grain architectures, are combined to form a multi-scale model for the generation of virtual electrode particles with full-grain detail. This work demonstrates the possibility of generating representative single electrode particle architectures for modeling and characterization that can guide synthesis approaches of particle architectures with enhanced performance.


Introduction
To enable widespread electrification of passenger cars, lithium(Li)-ion batteries must achieve high-energy densities, high-rates, and long cycle-lives. 1 The design of electrode architectures is known to significantly influence the rate performance and energy density of electrodes 2-4 but it is not well understood how the architecture of single particles influences their performance. For example, the most widely used positive electrode in Li-ion batteries consists of LiNi 1−x−y Mn y Co x O 2 (NMC) particles of various stoichiometries. The majority of NMC electrodes consist of polycrystalline particles where each grain within the particle facilitates transport of Li along 2D planes. When each grain is oriented randomly relative to its neighbor, one can envisage a sub-particle tortuosity where Li travels along the 2D plane of one grain and changes direction as it enters and travels along the 2D plane of the next grain, and so on as the particle continues to lithiate or delithiate. 5 Anisotropic expansion or contraction of the layered crystal structure of NMC occurs when the grains lithiate or delithiate. This expansion and contraction can lead to subparticle mechanical strains, 6-8 particle cracking, 9,10 and accelerated capacity fade. 7,8 Furthermore, the morphology and orientation of grains has been shown to greatly affect the rate capability of the electrode, where particles with grains oriented with exposed edge facets that transport Li radially inwards, display superior rate performance and longer life. [11][12][13] This knowledge presents an opportunity for sub-particle architectures to be designed in such a way to improve rate performance and reduce degradation within electrodes, where the efficacy of various particle architectures could † OF and LP contributed equally. * Corresponding authors Email addresses: orkun.furat@uni-ulm.de (Orkun Furat), donal.finegan@nrel.gov (Donal P. Finegan) be evaluated through 3D multi-physics models. However, to construct an accurate model of a particle, it is necessary to retrieve or generate realistic particle architectures from experimental data. X-ray nano-and micro-computed tomography (CT) has been demonstrated to be a powerful tool for non-destructively imaging and thereafter quantifying the morphological properties of many particles simultaneously 2, 9, 14-20 but falls short of capturing the detailed morphology of individual grains within particles and is incapable of determining the orientation of the grains.
Recently, focused ion beam electron backscatter diffraction (FIB-EBSD) has emerged as an effective means to quantify the morphology and orientation of sub-particle grains in 3D, 5,21 filling the capability gap of quantifying sub-particle grain detail. Merging the respective strengths of multiple techniques to compile datasets has been shown to enhance characterization capabilities, as demonstrated by merging X-ray nano-CT with FIB-SEM data, 22 or X-ray nano-CT with ptychographic CT 19 to reconstruct images with particle detail (from X-ray nano-CT) and detail on the conductive carbon and binder (from FIB-SEM or ptychography).
Here, we leverage data from both X-ray nano-CT and FIB-EBSD to generate artificial but representative single particle architectures completed with grain morphological details by parametric stochastic modeling. Using stochastic geometry models 23 for describing the 3D nano-and microstructure of various kinds of materials has several advantages. From a model-which was calibrated to image data-it is possible to generate virtual materials, so-called digital twins, the 3D morphology of which is statistically similar to that of real materials observed in tomographic image data. 24 Moreover, systematic variation of model parameters allows for the generation of a wide range of different virtual materials with predefined morphological properties. Then, the simulated 3D nano-and microstructures can be used as geometry input for (spatially resolved) numerical simulations, i.e., for virtual materials testing to study the influence of a material's geometry on its physical properties. 25,26 Furthermore, using parameter regression it is possible to predict model parameters and thus material properties based on parameters of manufacturing processes. 27 3D FIB-EBSD data depicting the grain architecture of an NMC particle  Figure 1: Scheme for stochastic multi-scale modeling of NMC particles. Using FIB-EBSD data depicting the grain architecture of an NMC particle (top, left), a 3D grain architecture model is developed which generates digital twins of the segmented FIB-EBSD data (top, middle) . On the other hand, using nano CT data depicting a large number of NMC particles (bottom, left), a 3D model for the outer shell of NMC particles is developed which generates random particle shells (bottom, middle). By combining both models, a stochastic multi-scale model is obtained for both the outer shell and the grain architecture of NMC particles (right).
In the present paper, we combine two kinds of stochastic 3D models which represent NMC particles on different length scales to obtain a multi-scale model for both the outer shell and grain architecture of NMC particles. This workflow is visualized in Figure 1. For modeling the outer shell of NMC particles we consider a random field on the unit sphere, 28 which is fitted to particles observed in CT data. 27,29,30 The fitted model can then generate virtual, but realistic outer particle shells. To model the grain architecture depicted in the FIB-EBSD data, we consider a random tessellation driven by a random marked point process. 23 Therefore, we first determine a parametric representation of the grains in the FIB-EBSD data using so-called Laguerre tessellations which are a generalization of Voronoi tessellations. 31 Then, a so-called Matérn cluster process 23 with conditionally independent marks is fitted to the parametric Laguerre representation of FIB-EBSD data, resulting in a stochastic 3D model for the grain architecture of NMC particles. Finally, in a last step both models are combined to obtain a multi-scale model from which we can generate virtual NMC particles together with their inner grain architecture. Thus, this work demonstrates how realistic particle architectures can be generated from experimental image data and sets the basis for future modeling efforts of single-particle geometries with tailored particle and grain characteristics.

Results
In this section we describe and quantitatively validate the stochastic 3D model and its individual components used in the present paper for modeling geometrical features of NMC particles at different length scales. Therefore, we first describe the imaging procedures and the image processing steps which were deployed for capturing a statistically relevant number of both outer particle shells (using nano CT) and sub-particle grains (using FIB-EBSD). Then, the stochastic models are described from which random outer particle shells and random grain architectures can be generated, respectively. Both models are combined into a multi-scale model for the outer shell and grain architecture of NMC particles. Finally, the stochastic geometry models described in this section are then calibrated to image data and validated.

3D Imaging of NMC particles and their grain architecture
The electrode particles analyzed in this work consisted of TODA LiNi 0.5 Mn 0.2 Co 0.2 O 2 (NMC532) particles in a 70 µm thick electrode coating of 90 wt% NMC532, 5 wt% C45 Timcal conductive carbon, and 5 wt% PVdF binder, on a 20 µm thick aluminum current collector. The electrode sheet was manufactured at Argonne National Laboratory's Cell Analysis, Modeling, and Prototyping (CAMP) facility. The TODA NMC532 particles had a size distribution of D10, D50, and D90 or 7.1, 9.3, and 12.1 µm, respectively.
For imaging the outer shell of such NMC particles, a sample was laser-milled into a pillar of approximately 80 µm diameter using a technique described previously. 32 The pillar of electrode material was then imaged using X-ray nano-computed tomography (CT) with a quasi-monochromatic beam of 5.4 keV in absorption contrast model (Zeiss Ultra 810, Zeiss Xradia XRM, Pleasanton, CA, USA) with a voxel size of 128 nm. The image was reconstructed using the Zeiss XRM Reconstructor software package that is based on filtered back projection. The reconstructed greyscale data was then used for segmenting the NMC particles and characterizing their morphology. A resource for micro-and nano-CT data of battery electrode architectures is the NREL Battery Microstructure Library. 33 For the forthcoming analyses, we consider the raw reconstructed CT data to be a map I CT : W CT → R, where W CT = {1, . . . , 499}×{1, . . . , 508}×{1, . . . , 1616}, see Figure 2a. Then, I CT (x) corresponds to the grayscale value of I CT of the voxel x ∈ W CT . For simplicity's sake, we extend the domain of I CT to Z 3 , i.e., we consider I CT : In order to quantitatively analyze the image data and stochastically model the outer shape of NMC particles, we have to compute a particle-wise segmentation map from I CT . Therefore, we first crop the image data I CT , followed by denoising using a Gaussian filter with standard deviation 4, see Figure 2b. 34 Then, the image is binarized using Otsu's method, 35 i.e., we obtain an image B CT : 1, if x corresponds to an NMC particle, 0, if x corresponds to the background.
Then, the binary image (cf. Figure 2c) is segmented into individual particles using the procedure described in Spettl et al. (2015). 36 The resulting segmentation of the CT data is visualized in Figure 2d. Since some particles are cut off by the cuboid sampling window, we do not consider them in our subsequent analysis-resulting in n = 239 a b c d e f g fully imaged particles as a basis for modeling the outer shape. From here on, we denote the segmented particles by P CT 1 , . . . , P CT n : Z 3 → {0, 1} which are given by for each i = 1, . . . , n. Note that this representation of particles allows us to easily check, whether some point x on the grid Z 3 belongs to a particle P CT i . However, we also have to make such inquiries for points in-between gridpoints, i.e., for any x ∈ R 3 . To do so, we extend the maps P CT i to R 3 using nearest neighbor interpolation, i.e., putting where x i denotes the closest integer to x i with x i = x i − 0.5 if 2x i is an odd integer. The resulting particles P CT 1 , . . . , P CT n will be used for modeling the outer shape of NMC particles. In addition to the CT data, further measurements were performed for mapping the inner grain architecture of NMC particles. For that purpose, NMC particles were removed from the electrode sheet by partially dissolving the electrode binder with N-Methyl-2-pyrrolidone (NMP) solvent, wiping the partially dissolved region with a rubber spatula, and smearing the electrode from the spatula onto a copper sheet in order to disperse individual electrode particles for analysis. Single particles with good electrical connection were identified on the copper sheet using Scanning Electron Microscopy (SEM) and were targeted for FIB-EBSD, which was conducted using an FEI Helios NanoLab 600i equipped with an EDAX-EBSD detector. A 30 kV 2.5 nA Ga beam was used for FIB-slicing of ca. 50 nm width and the EBSD measurements were taken at 50 nm steps in the x-and y-directions. Diffraction patterns were fit to a trigonal crystal system (space group R-3m) with a = b = 2.87 ¦ A, and c = 14.26 ¦ A to obtain orientation data from the measurements. The software produced text files containing spatially resolved confidence index, image quality (IQ), and Bunge-Euler angle data.
For modeling the grain architecture we utilize the segmentation of the FIB-EBSD data obtained in Furat et al. (2021). 21 Therefore, we give a short overview of the segmentation procedure which was used to obtain a grain-wise segmentation of the partially imaged NMC particle. 21 First, the stack of IQ images was aligned to obtain a 3D image, see Figures 2e and f. Then, the grain boundaries were manually labeled in four different slices of the 3D image. To enhance the grain boundaries in the entire 3D image, a convolutional neural network, namely a 3D U-Net 37, 38 with a modified loss function, was trained to predict grain boundaries from the IQ image. Finally, a grain wise segmentation of the 3D image was achieved by applying a marker-based watershed algorithm 36,39 to the neural network's grain boundary predictions, see Figure 2g.

Stochastic 3D model for the outer particle shell
In this section, we describe the stochastic model which is used to generate random outer particle shells being statistically similar to those of the segmented particles P CT i . Therefore, NMC particles are assumed to be pore-free with the constraint that each line segment between any point within the particle and another fixed point (e.g. particle centroid) is contained in the particle. This assumption limits the modeling to so-called star-shaped particles which reduces the complexity for modeling the 3D outer particle shell while still providing a reasonably flexible family of outer shell shapes and sizes. 27, 29 Thus, for each particle P CT i observed in the CT data, we assume that there is some point s i ∈ R 3 with P CT i (s i ) = 1 such that for any other point s ∈ P CT i the segment from s i to s is fully contained in the particle, i.e., the equation holds for all t ∈ [0, 1]. Furthermore, we assume that for each particle P CT i the so-called star point s i can be chosen as the particle's centroid x i which is given by where #{x ∈ Z 3 : P CT i (x) = 1} denotes the number of voxels associated with P CT i . Star-shaped particles P CT i have the advantage that they can be represented by a radius function r CT i : S 2 → [0, ∞) which is easier to model stochastically, where S 2 = {x ∈ R 3 : x = 1} denotes the unit sphere and · is the Euclidean norm. More precisely, the radius function r CT i assigns to any orientation vector u ∈ S 2 the distance between the star point s i and the particle's boundary into direction u. More precisely, the radius function r CT i is given by Thus, we replace the volumetric (discretized) representation P CT i by the surface representation r CT i . In practice, we will consider a discretized surface representation of r CT i , i.e., we will consider evaluation vectors (r CT i (u 1 ), . . . , r CT i (u m )) of the radius function for some gridpoints u 1 , . . . , u m ∈ S 2 . In other words, we transfer the data from a grid in space to a 2D grid-resulting in a complexity reduction which is advantageous for modeling purposes.
The goal of the stochastic model described below is to define a random function r : S 2 → [0, ∞) on the unit sphere S 2 , which statistically behaves similarly as the radius functions r CT 1 , . . . , r CT n computed from tomographic image data. Note that, from a representation as a radius function r we can determine the corresponding particle P : R 3 → {0, 1} uniquely up to the position of its star point by where o = (0, 0, 0) ∈ R 3 . The stochastic model which we use for describing the outer particle shell is a random field {X(u) : u ∈ S 2 } on the unit sphere S 2 . 28 Note that X(u) is a real-valued random variable for each orientation vector u ∈ S 2 . Since we are solely interested in modeling the shape of the particle shell and not in its location nor its orientation in space, we assume our random field to be isotropic, i.e., for any rotation matrix R ∈ SO(3) the random fields {X(u) : u ∈ S 2 } and {X(Ru) : u ∈ S 2 } have the same distribution. Specifically, we will fit a mixture of Gaussian random fields to reassemble the radius functions r CT 1 , . . . , r CT n computed from tomographic image data. Therefore, to make our paper more self-contained, we first explain the notion of an isotropic Gaussian random field on the unit sphere S 2 .
A random field {Z(u) : u ∈ S 2 } is called Gaussian if for any m > 0 and pairwise distinct orientation vectors u 1 , . . . , u m ∈ S 2 the random vector (Z(u 1 ), . . . , Z(u m )) follows a multivariate normal distribution. 40 Then, the random field {Z(u) : u ∈ S 2 } is uniquely defined by its mean function µ : S 2 → R and its covariance function σ : S 2 × S 2 → R which are given by and where E[·] denotes the expected value. For isotropic Gaussian random fields the mean function µ is constant-thus, it is described by a single numerical parameter which we denote by µ ∈ R. Moreover, in the isotropic case the covariance function of a Gaussian random field has the representation where P 0 , P 1 , . . . denote the Legendre polynomials and u 1 · u 2 is the inner product of u 1 and u 2 . 40,41 The coefficients a 0 , a 1 , . . . ≥ 0 are called the angular power spectrum and they, together with the mean value µ, uniquely determine the isotropic Gaussian random field {Z(u) : u ∈ S 2 }. In other words, the model parameters of an isotropic Gaussian random field on the unit sphere are the mean value µ and the angular power spectrum (a ) ∞ =0 . For computational reasons, we only consider Gaussian random fields such that a = 0 for all > L, where L > 0 is some fixed integer. Then, the angular power spectrum of these random fields is characterized by a finite-dimensional vector a = (a 0 , . . . , a L ). Consequently, the covariance function σ of {Z(u) : u ∈ S 2 } is given by Note that for the fitting procedure described below it is necessary to consider the finite-dimensional distributions of an isotropic Gaussian random field {Z(u) : u ∈ S 2 } with mean value µ ∈ R and angular power spectrum a = (a 0 , a 1 , . . . , a L ) ∈ R L+1 . More precisely, for any fixed integer m > 0, let u 1 , . . . , u m be pairwise distinct points on S 2 . Then, the random vector (Z(u 1 ), . . . , Z(u m )) has a multivariate normal distribution with mean vector µ = (µ, . . . , µ) ∈ R m and covariance matrix Specifically, assuming that Σ is not singular, i.e. det(Σ) = 0, the probability density function f µ,a : R m → [0, ∞) of (Z(u 1 ), . . . , Z(u m )) is given by Thus, for simulating such a Gaussian random field at some (arbitrary, but fixed) evaluation points u 1 , . . . , u m ∈ S 2 we have to draw samples from a multivariate normal distribution. Note that in particular, according to Eq. (13), the random variable Z(u) is normal distributed with mean value µ and variance σ(u, u) for each u ∈ S 2 . However, the histograms of the values r CT 1 (u), . . . , r CT n (u) of the radius functions extracted from the segmented image data for any fixed direction u ∈ S 2 can not be fitted by the density function of a normal distribution, see Figure 3. Specifically, the histogram shown in Figure 3 indicates a multimodal distribution which motivates using mixtures of Gaussian random fields to model the outer shell of particles.
More precisely, we randomly select among K > 0 Gaussian random fields and generate an outer shell by drawing a sample from the selected model. Therefore, let {Z (1) (u) : u ∈ S 2 }, . . . , {Z (K) (u) : u ∈ S 2 } be independent, isotropic Gaussian random fields with mean values µ (1) , . . . , µ (K) and angular power spectra a (1) = (a (1) ) L =0 , . . . , a (K) = (a (K) ) L =0 , respectively. We call the random fields {Z (1) (u) : u ∈ S 2 }, . . . , {Z (K) (u) : u ∈ S 2 } the components of the mixture. In order to model the random selection of a component, we introduce further parameters p 1 , . . . , p K ≥ 0 with p 1 + · · · + p K = 1 which correspond to the probability of selecting a component. Then, the mixture {X(u) : u ∈ S 2 } of K Gaussian random fields is given by with probability p 1 , . . . with the vector θ = (L, K, p 1 , . . . , p K , µ (1) , . . . , µ (K) , a (1) , . . . , a (K) ) of model parameters. Furthermore, the finitedimensional probability density function f θ : R m → [0, ∞) of the random field {X(u) : u ∈ S 2 } for some evaluation points u 1 , . . . , u m ∈ S 2 is given by To generate a realization of {X(u) : u ∈ S 2 } at the evaluation points u 1 , . . . , u m we just have to randomly pick a component according to the probabilities p 1 , . . . , p K and then draw a sample from the corresponding multivariate normal distribution. Now that we have defined the mixture {X(u) : u ∈ S 2 } of Gaussian random fields, we describe how its parameter vector θ = (L, K, p 1 , . . . , p K , µ (1) , . . . , µ (K) , a (1) , . . . , a (K) ), for any fixed pair L, K > 0, can be fitted to data. Recall that previously in this section we described the procedure for representing particles extracted from tomographic image data as radius functions r CT 1 , . . . , r CT n : S 2 → R. Now, we describe how to fit the parameter vector θ such that the realizations of the model {X(u) : u ∈ S 2 } statistically reassemble the radius functions r CT 1 , . . . , r CT n : S 2 → R. Therefore, let u 1 , . . . , u m ∈ S 2 be some pairwise distinct evaluation points. Then, the random field {X(u) : u ∈ S 2 } is fitted by determining an optimal parameter vector θ for which the evaluations r CT 1 = r CT 1 (u 1 ), . . . , r CT 1 (u m ) , . . . , r CT n = r CT n (u 1 ), . . . , r CT n (u m ) are most likely with respect to the distribution of the random vector (X(u 1 ), . . . , X(u m )). More precisely, the optimal parameter θ is obtained by maximizing the (log-)likelihood function, i.e., by where f θ denotes the probability density of (X(u 1 ), . . . , X(u m )) given in Eq. (15) and Θ L,K is the set of admissible model parameters for fixed L, K > 0. Since the probability density f θ is differentiable with respect to the remaining model parameters in θ, the optimization problem given in Eq. (16) can be solved using gradient descent algorithms. 42 However, since the number of model parameters can be-depending on the number K of components and the number L + 1 of angular power spectrum coefficients per component-relatively large, we require a good initial parameter constellation for applying gradient descent algorithms. Therefore, we present a method for determining a good initial parameter vector θ init for performing optimization. Instead of directly fitting the mixture {X(u) : u ∈ S 2 } of Gaussian random fields to the evaluations r CT 1 , . . . , r CT n by solving the optimization problem given in (16), we first fit a mixture of finite-dimensional Gaussian random vectors to r CT 1 , . . . , r CT n . Then, from this fit we derive an initial parameter constellation θ init ∈ Θ L.K to solve Eq. (16). Therefore, we shortly describe the construction of mixtures of finite-dimensional Gaussian random vectors which is similar to the construction of mixtures of Gaussian random fields described in Eq. (14). Let Z 1 , . . . , Z K be pairwise independent, normally distributed random vectors taking values in R m , with mean vectors µ (1) , . . . , µ (K) ∈ R m and non-singular covariance matrices Σ (1) , . . . , Σ (K) ∈ R m×m . Thus, for k = 1, . . . , K, their probability densities f k are given by for any x ∈ R m . Furthermore, let p 1 , . . . , p K ≥ 0 be some probabilities with p 1 + · · · + p K = 1. Then, the mixture X of the finite-dimensional Gaussian random vectors Z 1 , . . . , Z K is given by The parameters of the distribution of X-namely the mixing probabilities p 1 , . . . , p K , the mean vectors µ (1) , . . . , µ (K) and the covariance matrices Σ (1) , . . . , Σ (K) -can be fitted to the evaluation vectors r CT . . , n, using the so-called expectation-maximization algorithm. 43 Note that the distribution of the random vector X, fitted in this manner, does not yet describe the outer particle shells extracted from CT data sufficiently well. First of all, there can be some k ∈ {1, . . . , K} such that the components of the mean vector µ (k) are not identical, which is in contradiction with the assumption of anisotropy. Furthermore, the random vector X does not describe the behaviour of the outer shell at evaluation points u ∈ S 2 such that u = u i for all i = 1, . . . , m. Therefore, in the next step, we describe how to determine a parameter constellation θ init ∈ Θ of a mixture {X(u) : u ∈ S 2 } of isotropic Gaussian random fields which is in some sense similar to the random vector X at the evaluation points u 1 , . . . , u m , i.e., for which it holds that (X(u 1 ), . . . , X(u m )) ≈ X.
First, for each k = 1, . . . , K we determine the parameters of an isotropic Gaussian random field The angular power spectrum a (k) = (a (12), in some sense coincides with the covariance matrix Σ (k) of the corresponding multivariate normal distribution, i.e., we have to solve a system of linear equations with unknowns a where the equations are given by for each i, j = 1, . . . , m. Furthermore, we have the constraints a L ≥ 0. This system of linear equations does not necessarily have an exact solution. However, using an iterative algorithm 44 we can determine a good solution a (k) in a least-squares sense. Then, the initial value θ init for solving the optimization problem in Eq. (16) is given by Stochastic 3D model for the grain architecture using Laguerre tessellations In this section we explain our stochastic modeling approach, which is used to simulate the 3D grain architecture inside NMC particles, and its calibration to the FIB-EBSD data described above. For this, we consider a certain sampling window W ⊂ R 3 , being a cuboid W = [a 1 , for some a 1 , a 2 , b 1 , b 2 , c 1 , c 2 ∈ R with a 1 < a 2 , b 1 < b 2 , c 1 < c 2 , and employ (random) tessellations to randomly partition the sampling window W into non-overlapping cells, each of which corresponds to one of the grains in the FIB-EBSD data. In combination with the stochastic model of the outer shells described in the previous section, we are then able to present a holistic stochastic model for complete NMC particles, along with their inner 3D grain architecture.
To parametrically subdivide the sampling window W ⊂ R 3 into cells, Laguerre tessellations-generalizations of the well-known Voronoi tessellation 31, 45 -are used. For this purpose, for any s, x ∈ R 3 and w ∈ R, the power distance pow(x, (s, w)) (also called Laguerre distance function) is considered, which is defined by A Laguerre tessellation T of W with n T cells is fully characterized by a sequence of weighted generating points with seed points s i ∈ R 3 and weights w i ∈ R for i = 1, . . . , n T . The i-th cell C L i of the Laguerre tessellation T is then defined as For simplicity, we use the notation to indicate whether a point x ∈ W belongs to the i-th Laguerre cell. It can be shown that Laguerre cells are convex polyhedra. Note that for modeling non-convex grains other tessellation models, e.g. Johnson-Mehl tessellations 46,47 can be deployed analogously. Furthermore, note that the seed points roughly determine where the corresponding cell manifests (though depending on the weights, the seed point might actually not belong to the cell, or the cell might even be empty), and that high additive weights compared to those of neighboring generating points lead to larger cells. A big advantage of Laguerre tessellations is that they provide much more flexibility compared to simpler Voronoi tessellations through their additional weights. For more details on this topic, see e.g. Lautensack and Zuyev (2008). 48 In summary, Laguerre tessellations allow us to parametrically describe a system of convex grains. For stochastic modeling it is thus sufficient to come up with a random sequence of seed points and weights. With these, the spatial extends of the simulated grains can then be easily computed. However, before we can discuss the stochastic model for the 3D grain architecture inside NMC particles and how to calibrate it to experimental image data, it is necessary to find a representation of the grains in the segmented FIB-EBSD data (Figure 4a) in terms of a Laguerre tessellation, see Figure 4b. We will then use the extracted seed points and weights as basis for the modeling.
Analogous to the representation of particles {P CT i } n i=1 extracted from segmented CT data which has been introduced above, we consider the i-th grain of the experimental FIB-EBSD data as a map C E i : Z 3 → {0, 1} given by with i = 1, . . . , n E and n E the number of grains in the sampling window W ⊂ R 3 . Again, nearest neighbor interpolation can be used to extend the domain of C E i to the continuous Euclidean space. Furthermore, let X F = {x F j } nF j=1 ∈ Z 3×nF be the coordinates of voxels that belong to the foreground of the image data, i.e., for each At its heart, fitting a Laguerre tessellation to image data is an optimization problem of finding the best generators such that the discrepancy of the tessellation and the image data is minimized. The ground truth for that is given by the grain labels We employ a volume-based discrepancy measure which counts all voxels where the experimental FIB-EBSD image data and the discretized Laguerre tessellation have the same label. To be more precise, the value of the objective function E : where the cells C L 1 , . . . , C L nE of the Laguerre tessellation T depend on the choice of the generators in G such that n T = n E . The corresponding fitting problem is thus to determine an optimal sequence of generators G opt defined as Eq. (27) is probably the most natural choice to define the fitting problem, but in practice it is computationally expensive. In the following, we address this issue by slightly altering the original fitting problem.
It is easy to see that C L j (x F ) with x F ∈ X F can be reformulated as where argmin * j is the j-th component of the n T -dimensional vector-valued argmin function, i.e., In cases where the minimum is not unique, i.e., there are indices j 1 , j 2 ∈ {1, . . . , n T } with j 1 = j 2 and z j1 = z j2 , only the component with the smaller index is set equal to 1. The function argmax * is defined analogously. Even though the Laguerre distance function pow is differentiable (with respect to the generators), the fact that argmin * in Eq. (28) is not makes the objective function E defined in Eq. (26) non-differentiable. This leaves us only with derivative-free optimization algorithms to solve Eq. (27), which in most cases converge slower than gradient-descent methods. 49 In order to increase efficiency, we slightly deviate from the original Laguerre tessellation formulation by replacing the argmin * function in Eq. (28) with a 'softmin * ' function, i.e., a softmax * function with a negative argument Here, the n T -dimensional softmax * function is a smooth version of the argmax * function. So instead of returning a vector whose components are either 0 or 1, the softmax * function is a vector-valued map where each component is a continuous function with values between 0 and 1. In fact its output vector softmax * z for some argument z ∈ R n T defines a discrete probability measure (i.e., the values of all components are between 0 and 1 and their sum is equal to 1), which assigns the highest probability to the index j if z j ≥ z i for all i ∈ {1, . . . , n T }. The softmax * function is often used for multi-class classification problems in machine learning, 50 which show many similarities with our tessellation fitting. With this, the value of a component of the vector ( C L 1 (x F ), . . . , C L n T (x F )) is thus the largest one if the corresponding generator has the shortest power distance to the given evaluation point x F . Additionally, since argmax * is a composition of differentiable functions and therefore differentiable, the function C L j is differentiable, too, for each j ∈ {1, . . . , n T }. Because it is now no longer possible with C L j to assign each voxel a unique label, it is necessary to also adapt the objective function. Thus, instead of E defined in Eq. (26), we consider the function E : with the (binary) cross-entropy loss function : [0, 1] 2 → [0, ∞) given by (ỹ, y) = −y logỹ − (1 − y) log(1 −ỹ). Note that the cross-entropy loss is often used in machine learning 50 to compare the output of a classifier with ground truth data-basically the same purpose it serves here: If an evaluation point x F belongs to the i-th grain (i.e., ) also needs to be close to 1 in order to maximize and vice versa. Note that the modified objective function E is differentiable with respect to the Laguerre generators. The fitted sequence of generators G opt can then be obtained by solving In summary, we modified the original fitting problem, formulated in Eq. (27), to get a differentiable version of it in Eq. (33). This allows us to employ fast, gradient-based optimization algorithms, such as in our case, Adam. 51 The initial guess for the generators is computed using the procedure described in Spettl et al. (2016), 52 which is based on a global optimization algorithm minimizing a fast, but less accurate objective function. Note that for the discretization of a fitted Laguerre tessellation, we compute the (unique) cell labels using the classical definition in Eq. (28).
Once a sequence of generators To establish a suitable parametric model for the sequence of random seed points , we employ so-called Matérn cluster processes 23 to describe the spatial clustering of seed points extracted from experimental FIB-EBSD data. Matérn cluster processes are comprised of two modeling components: First, points are drawn from a certain parent point process, which play the role of cluster centers. Then, an individual cluster of descendant points is created by a different (independent) point process at each of the previously generated locations of parent points. The parent points themselves are disregarded for the final point pattern. More specifically, a homogeneous Poisson point process with intensity λ M > 0 is used to draw the cluster centers. For each cluster, the number of descendant points is drawn from a Poisson distribution with rate parameter µ M > 0, and the descendant points are uniformly distributed in a sphere with radius r M > 0 surrounding the cluster center. Thus, the sequence of random seed points is characterized by three model parameters: λ M , µ M and r M ; The intensity of the parent point process λ M and the cluster size r M are estimated using the so-called minimum contrast method with respect to the pair correlation function, 23 a functional characteristic which statistically describes the distances between pairs of seed points, where a peak at a given distance indicates that inter-point distances of that length occur frequently. In this context, the minimum contrast method minimizes the discrepancy (with respect to the L 2 -norm) between the theoretical pair correlation function of the Matérn cluster process and the one estimated from the sequence of seed points {s i } n T i=1 . Finally, the third model parameter µ M can be estimated in the following way. Since a closed formula for the overall intensity of Matérn cluster processes is known, 23 it can be set equal to n T /|W|, i.e., the number of seed points of the ground truth divided by the volume of the sampling window |W|. The resulting equation can then be solved for the mean number of points in each cluster µ M . The next step is to model the additive weights {w i } n T i=1 of the Laguerre tessellation fitted in Eq. (33). The easiest way to do this would be to estimate the probability distribution of additive weights in the fitted Laguerre tessellation and to independently sample from that distribution when the stochastic grain architecture model is simulated. However, it turns out that this approach is too simple and does not capture the interdependencies of the fitted generators well enough. We therefore choose a different approach: For each random seed point S M i , we draw the corresponding additive weight W M i conditionally on the random distance D M i to its nearest neighboring seed point. By doing so, it is possible to give seed points which are grouped together in a spatial cluster different weights than those points that are isolated. This leaves us with the question of how to model the joint probability distribution of the additive weights and the nearest neighbor distances, from which the conditional probability distribution can be computed directly. One possible answer can be given by means of copulas. 53 A two-dimensional copula C : [0, 1] 2 → [0, 1] is a cumulative distribution function with [0, 1]-uniformly distributed marginal distributions. Sklar's theorem 53 states that for any two (real-valued) random variables X 1 and X 2 with the cumulative distribution functions F X1 and F X2 defined by F X1 (x 1 ) = P(X 1 ≤ x 1 ) and F X2 (x 2 ) = P(X 2 ≤ x 2 ) for all x 1 , x 2 ∈ R, respectively, there is a copula C such that their joint cumulative distribution function F (X1,X2) , defined by F (X1,X2) (x 1 , x 2 ) = P(X 1 ≤ x 1 , X 2 ≤ x 2 ) for all x 1 , x 2 ∈ R, can be written as for all x 1 , x 2 ∈ R.
In practice it is often beneficial to consider probability densities instead of cumulative distribution functions. So assuming that C, F X1 , and F X2 are differentiable, by applying the chain rule to Eq. (34) the joint probability density f (X1,X2) of (X 1 , X 2 ) can be obtained as with f X1 and f X2 the probability densities of X 1 and X 2 , respectively, and c the probability density of C, i.e., The conditional probability density f (X1|X2) (·|x 2 ) of X 1 provided that X 2 = x 2 for some x 2 ∈ R is then given by see, e.g., Aas et al. (2009). 54 In our case, we consider the nearest neighbor distance D 0 of the typical seed point-a characteristic which describes the Euclidean distance from a seed point, selected at random, to its nearest neighbor among the remaining seed points. 23 Furthermore, we consider the additive weight W 0 of the typical seed point. Thus, X 1 corresponds to the weight W 0 of the typical seed point, and X 2 to its nearest neighbor distance D 0 . This allows us to model the univariate marginal distributions of the random vector (W 0 , D 0 ) with parametric distributions first, and then the dependence structure of (W 0 , D 0 ) using (parametric) copulas. This drastically facilitates the modeling of higherdimensional random vectors compared to deriving their joint distributions directly.
Motivated by the statistical analysis performed below regarding model validation, for the sequence of generators only via the random distance D M i of S M i to its nearest neighboring seed point, (iii) the probability density of the two-dimensional random vector (W 0 , D 0 ) of additive weight and nearest neighbor distance of the typical seed point is parametrically modeled as follows, using the representation formula given in Eq. (35). In particular, we assume that the first component W 0 can be modeled by a (shifted) inverse-gamma distribution. That is, its density f W0 is given by with shape parameter α M > 0, scale parameter β M > 0 and location parameter θ M ∈ R. Here, Γ denotes the gamma function. The name of this distribution originates from the fact that the reciprocal of a gamma-distributed random variable is inverse-gamma-distributed (and vice versa). Furthermore, we use the fact that for the cumulative distribution function F D0 of the nearest neighbor distance D 0 of the typical seed point from a Matérn cluster process the following formula holds: 55 where V (d, r M , x) is the volume of the intersection of two balls with radii d and r M such that the distance between their centers is equal to x, and F : R → [0, 1] is the cumulative distribution function of the so-called contact distance, which is given by For modeling the joint distribution of the random vector (W 0 , D 0 ) we employ a special case of the Tawn coupla, 56 whose asymmetric behavior matches that observed in the data. The general definition of this copula is given by where u 1 , u 2 ∈ (0, 1) and (41) with the parameters ψ M 1 , ψ M 2 ∈ [0, 1] and κ M > 1. In our case, we set ψ M 1 = 1. To estimate the parameters α M , β M , θ M , ψ M 2 , and κ M of the joint distribution of (W 0 , D 0 ), we use the representation formula of the conditional probability density f (W0|D0) given in (37), i.e., maximizing the likelihood function where {d i } n T i=1 is the sequence of nearest neighbor distances corresponding to the sequence of seed points {s i } n T i=1 . In summary, the simulation of the 3D grain architecture is comprised of the following steps: First, the Matérn cluster process is simulated in the sampling window W ⊂ Stochastic multi-scale model for the outer particle shell and grain architecture In this section we describe our approach to combine the outer shell model {X(u) : u ∈ S 2 } and the grain architecture , introduced in the previous sections, resulting in a multi-scale model for the outer shell and grain architecture of NMC particles, see Figure 5. Note that the outer shell model {X(u) : u ∈ S 2 } is fitted to reassemble the shape and size of NMC particles observed in nano CT data, while the parameters of the random Laguerre are attuned to statistically represent the grain architecture of a single particle which was partially imaged with FIB-EBSD. Since the grain architecture of NMC particles might depend on their size, in the following we will consider a conditional version of the fitted outer shell model {X(u) : u ∈ S 2 } to only generate outer shells with a similar size as the NMC particle imaged using FIB-EBSD. Since the latter was only partially imaged its true size is unknown-therefore, in a first step, we describe an approach to estimate the particle size. We start by computing a binary image B : Z 3 → {0, 1} depicting the (partially visible) outer shell of the NMC particle imaged with FIB-EBSD by x is neighbor of a grain but belongs to no grain, 0, else.
Then we fit the parameters of a sphere, i.e., the center x c ∈ R 3 and the diameter d > 0, to the outer shell image B by minimizing the discrepancy where p xc,d (x) denotes the point on a sphere with center x c and diameter d which is closest to x. Then, the fitted sphere's diameter d opt = 7.44 µm is an estimate for the size of the particle which was partially imaged with FIB-EBSD. Using this size estimate, we generate virtual multi-scale particle morphologies according to the following scheme: (i) For some m > 0 and pairwise distinct orientation vectors u 1 , . . . , u m ∈ S 2 , generate realizations r sim of the random vector (X(u 1 ), . . . , X(u m )) until the volume-equivalent diameter d(r sim ), given by belongs to the interval [0.98 d opt , 1.02 d opt ], see Figure 5a. Details on the computation of the volume V (r sim ) are given below. Figure 5b.
for which the centroids, see Eq. (5), are located inside the outer shell defined by r sim . We denote these cells by c I 1 , . . . , c I nI . (iv) The discretized multi-scale particle morphology P M : Z 3 → R is then given by Figure 5c visualizes a realization of the multi-scale particle model. Note that due to this construction scheme the outer shell of a realization P M differs from the outer shell which was generated in step (i). More precisely, this approach introduces surface roughness which similarly can be observed in FIB-EBSD data. 21

Model calibration
In this section we describe the procedure for calibrating the model parameters of the stochastic geometry models for the outer shell and the grain architecture. The stochastic outer shell model was fitted using the Adam stochastic gradient descent method. 51 First, we manually tuned the truncation parameter L in the series expansion given in Eq. (11) and the number K of Gaussian components, putting L = 15 and K = 7. For fitting the remaining model parameters, we chose m = 258 orientation vectors u 1 , . . . , u m ∈ S 2 which provide a triangular mesh on the unit sphere S 2 . 57 The fitted model parameters p 1 , . . . , p K , µ (1) , . . . , µ (K) , a (1) , . . . , a (K) are given in Table 1. Note that even though the model parameters of {X(u) : u ∈ S 2 } were fitted using m = 258 evaluation points, the realizations of the model can be evaluated at arbitrarily many sampling points.  The parameters of the grain architecture model were fitted using the procedure described above. With this, we get the following parameter values for the Matérn cluster process of seed points: λ M = 0.00108, µ M = 0.418, and r M = 2.58. Using the aforementioned maximum likelihood approach results in parameters for the distribution of the additive weights α M = 16.8, β M = 1247, θ M = −76.4, and for the copula κ M = 1.25, ψ M 2 = 0.372. Note that the model calibration was performed based on the FIB-EBSD image data and uses its voxel length scale.
In order to evaluate the goodness-of-fit of the grain architecture model with respect to the Laguerre tessellation fitted to FIB-EBSD data, we first consider various characteristics of the Laguerre generators which have been used for model fitting. The goodness-of-fit with respect to further characteristics, not used for model fitting, will be analyzed in the next section.
In Figure 6a the pair correlation functions of the seed points of the fitted Laguerre tessellation and those of a realization of the stochastic grain architecture model are shown. Additionally, the theoretical pair correlation function of the fitted Matérn cluster process is depicted in Figure 6a. It becomes clear that the pattern of seed points of the fitted Laguerre tessellation is far from the complete spatial randomness assumption (indicated as dashed line in Figure 6a at the y-axis value 1). In fact the high values of the pair correlation function for small distances prove that the seed points are spatially clustered. The fitted Matérn cluster process captures this behavior quite well.
When considering the probability densities of the nearest neighbor distances in Figure 6b, it is important to keep in mind that the theoretical one is fully specified by the Matérn cluster process and its parameters. While the theoretical density follows the general shape of the density, which has been estimated for the nearest neighbor distances of the seed points fitted to FIB-EBSD data (i.e., both densities are concentrated in roughly the same interval and the modes are at similar positions), the finer details are not fully captured. On the other hand, the (small) differences between the theoretical density and the density estimated for the simulated nearest neighbor simulated data d Figure 6: Calibration of the grain architecture model. Comparison of characteristics of the Laguerre tessellation fitted to FIB-EBSD data and the stochastic grain architecture model. Pair correlation function of seed points (a), probability densities of nearest neighbor distances (b) and additive weights (c), and two-dimensional (joint) density of additive weights and nearest neighbor distances of seed points (d). All distances are given in terms of the voxel length scale of FIB-EBSD data.
distances are due to the finite sample size and the smoothing effect of the kernel density estimator. The densities estimated for the additive weights of the Laguerre generators fitted to FIB-EBSD data and the simulated weights shown in Figure 6c indicate a good fit. The small discrepancies can be explained by the fact that the latter are not drawn directly from the marginal distribution, but conditionally on the nearest neighbor distances of the corresponding seed points.
Finally, by considering a bivariate kernel density estimator for the joint density of nearest neighbor distances and additive weights, see Figure 6d, a generally good correspondence can be observed. Most of the slight mismatches can be explained by the differences mentioned above between the theoretical and estimated densities of the nearest neighbor distances, cf.

Model validation
Having fitted the models for the outer shell and the grain architecture to the experimental data, model validation is performed by comparing probability distributions of size and shape characteristics of particles and grains extracted from tomographic image data to those of virtual particles and grains generated by the stochastic models. Note that these characteristics have not been used for parameter fitting.
For a quantitative validation of the fitted outer shell model {X(u) : u ∈ S 2 } we first generate n discretized realizations, i.e., we independently draw n = 239 realizations r sim 1 , . . . , r sim n of the random vector (X(u 1 ), . . . , X(u m )). Then, we compare size and shape characteristics of the corresponding particles to those extracted from CT data. More precisely, for each outer particle shell, represented by an evaluation vector r P = (r P 1 , . . . , r P m ) ∈ R m , we compute the volume-equivalent diameter d(r P ) given by where V (r P ) denotes the volume of a particle corresponding to the evaluation vector r P . Since the entry r P i of the latter represent the radius of the corresponding particle into the direction of the orientation vector u i , it is impossible to determine the true volume of the particle based on a finite number of directional radii r P 1 , . . . , r P m . Therefore, we now specify the procedure for computing an estimate V (r P ) for the particle's true volume.
Recall that the orientation vectors u 1 , . . . , u m chosen in the previous section form a triangular mesh of the unit sphere. More precisely, the surface of the unit sphere is approximated by triangles, the vertices of which are given by certain triplets (u i1 , u i2 , u i3 ). We denote this set of triangle vertices by V ⊂ {u 1 , . . . , u m } 3 . Then, by scaling the vertices in V with the corresponding radii we obtain a triangular mesh V P of the outer particle shell, i.e., From this representation of a particle we can easily compute an estimate for its volume by where V tetra (v 1 , v 2 , v 3 , o) denotes the volume of a tetrahedron with vertices v 1 , v 2 , v 3 and the origin o = (0, 0, 0) ∈ R 3 . Note that the representation of a particle's surface as a mesh V P allows for the estimation of its surface area. More precisely, by summing up the areas of the triangles represented by the triplets of vertices in V P we obtain an estimate A(r P ) for the surface area of a particle represented by the evaluation vector r P . Inserting the estimate V (r P ) of the particle volume given in Eq. (51) into Eq. (49) we computed the volumeequivalent diameters d(r CT 1 ), . . . , d(r CT n ) and d(r sim 1 ), . . . , d(r sim n )) of particles observed in CT and of model realizations, respectively. Figure 7a indicates a good fit between kernel density estimates of both samples.
To validate the shape of model realizations we investigate the sphericity s(r P ) of the outer particle shell represented by the evaluation vectors r P , which is given by The corresponding kernel density estimates of the sphericity values s(r CT 1 ), . . . , s(r CT n ) and s(r sim 1 ), . . . , s(r sim n ) are visualized in Figure 7b which again indicates a good fit.
Moreover, we computed bivariate kernel density estimates for pairs of volume-equivalent diameter and sphericity of outer particle shells. More precisely, the bivariate probability density corresponding to pairs of characteristics (d(r CT 1 ), s(r CT 1 )) . . . , (d(r CT n ), s(r CT n )) computed from CT data is visualized in Figure 7c, whereas the bivariate density computed from the characteristics (d(r sim 1 ), s(r sim 1 )) . . . , (d(r sim n ), s(r sim n )) of realizations of the outer shell model is depicted in Figure 7d. A visual comparison indicates a good match between both bivariate probability densities. The bivariate probability densities indicate that small particles are preferentially more spherical though with a higher variance in sphericity values, while larger particles have a narrower sphericity distribution. Such distributions are material-specific with different trends for other electrode materials as shown in a previous work. 58 Recall that for the multi-scale model we generated solely outer shells from the model {X(u) : u ∈ S 2 } for which the corresponding particles have a volume-equivalent diameter close to d opt = 7.44 µm. To validate the sphericity distribution of particles whose outer shells were generated with this specific volume-equivalent diameter we computed the conditional probability densities of the sphericity from the bivariate probability densities depicted in Thus, the outer shell ∂P M does not necessarily coincide with the realization of the outer shell model {X(u) : u ∈ S 2 } which was used for generating P M . Therefore, we perform further validation-by comparing the shape of particles extracted from CT data to the shape of outer shells which were computed from realizations of the multi-scale model.
To do so, we generated (discretized) realizations P M 1 , . . . , P M n of the multi-scale model. For better comparability to the outer shells extracted from CT image data, we first apply the same degree of smoothing to the realizations P M 1 , . . . , P M n as to the CT data during image processing described above. A smoothed realization of P M is visualized in Figure 8. Then, we computed sphericity values of the smoothed versions of P M 1 , . . . , P M n and their corresponding kernel density estimate, see Figure 7e, which indicates a slight underestimation of sphericity values for the multi-scale model.
After the validation for the outer shell model, we now consider cell characteristics computed from realizations of the stochastic grain architecture model, which have not been used for model fitting, and compare them to those of the grains in the experimental FIB-EBSD data. These comparisons serve on the one hand as (further) measures for the goodness-of-fit of the fitted Laguerre tessellation with respect to the (segmented) FIB-EBSD image data, but also as validation for the stochastic model since grain-based characteristics have not been used for model fitting. A more comprehensive analysis of the FIB-EBSD image data was provided by Furat et al. (2021). 21 The probability density of volume-equivalent diameters of the fitted Laguerre cells, visualized in Figure 9a, closely matches that computed from the experimental FIB-EBSD data, whereas the stochastic grain architecture model tends to produce more medium sized grains and fewer very large grains. A similar behavior can be asserted for the densities of grain/cell surface areas, see Figure 9b. Regarding the sphericity of cells, the densities computed from a realization of the stochastic model and from its ground truth, the fitted Laguerre tessellation, indicate a good match, where, however, some slight mismatches to the density computed from the FIB-EBSD image data are visible, see Figure 9c. Note that for all three characteristics, the corresponding mean values are all in close proximity, and most discrepancies come from a slightly different skewness and tail-behavior of the probability densities. The histogram of the numbers of cell neighbors for the simulated data, visualized in Figure 9d, is shifted a little bit to the right, compared to the similarly distributed numbers of cell/grain neighbors in the fitted Laguerre and segmented FIB-EBSD data. In summary, together with a visual inspection in Figure 4, a very good match of the fitted Laguerre tessellation and the (segmented) experimental FIB-EBSD data can be observed. Furthermore, the stochastic grain architecture model captures all considered structural characteristics reasonably well.

Discussion
As shown in the previous section the proposed multi-scale model for the outer shell and grain architecture of NMC particles represents the particle(s) mapped with CT and FIB-EBSD sufficiently well. However, in principle, modifications and further extensions of the outer shell model and the grain architecture model are possible, although their fitting and validation with respect to tomographic image data might then be computationally more expensive.
Regarding the outer shell model the following modifications are possible: Recall that according to Eqs. (13) and (15)  distributions increases with the number of considered orientation vectors, i.e., when the realization is to be highly resolved, we have to draw from a high-dimensional normal distribution which can be inefficient. For such inquiries it can be useful to utilize an alternative representation of Gaussian random fields {Z(u) : u ∈ S 2 } considering a series expansion with respect to so-called spherical harmonics functions. More precisely, an isotropic Gaussian random field {Z(u) : u ∈ S 2 } on the unit sphere with mean value µ and angular power spectrum a = (a 0 , a 1 , . . . , a L ) adheres to the representation Re(a ,k )Re(Y ,k (u)) + Im(a ,k )Im(Y ,k (u)) for each u ∈ S 2 , where the a ,k are independent random variables with a ,0 ∼ N (µ, a ) for ≥ 0 and Re(a ,k ), Im(a ,k ) ∼ N (0, a /2) for , k > 0, see. 29,40 Note that the spherical harmonics functions Y ,k : S 2 → C are given by where i denotes the imaginary unit and P ,k are the associated Legendre functions. 59 The representation of {Z(u) : u ∈ S 2 } given in Eq. (54) has the advantage that we can first generate realizations of the random coefficients a ,k , the number of which does no depend on the number of evaluation points u 1 , . . . , u m . Then, we can compute the values of the corresponding realization of the random field {Z(u) : u ∈ S 2 } at arbitrarily many orientation vectors, i.e., we can obtain highly resolved realizations of {Z(u) : u ∈ S 2 } without having to simulate high-dimensional normally distributed random vectors. The spherical harmonics representation given in Eq. (54) can also be utilized for simulating the fitted Gaussian mixture model {X(u) : u ∈ S 2 }. More precisely, we have to randomly pick a component according to the probabilities p 1 , . . . , p K introduced in Eq. (14) and then generate realizations of the spherical harmonics coefficients of the corresponding Gaussian model. The validation of the outer shell model {X(u) : u ∈ S 2 } indicates a good match between distributions of size/shape characteristics of model realizations and particles extracted from CT data, see Figure 7. However, Figure 7e indicates that the sphericity distribution of the outer shell of realizations P M of the multi-scale model is slightly underestimated. Therefore, further modifications to the outer shell model could be taken into account which can ensure an even better match between the sphericity distributions of particles extracted from CT data and model realizations. More precisely, rejection sampling 60 could be applied: Let f data denote the (size-conditioned) probability density of the sphericity values of particles extracted from CT data. Furthermore, let f model be the corresponding probability density of some outer shell model {X(u) : u ∈ S 2 }, being, e.g., a Gaussian random field or a mixture of Gaussian random fields. Note that f model can be determined by extensive simulation, e.g., by generating a large number of outer shells from which the sphericity distribution can be estimated. For rejection sampling to be applicable the inequality f data (x) ≤ M f model (x) has to hold for each x ∈ [0, 1], where M > 1 is some constant. Then, the rejection sampling is performed as follows: (i) Generate a realization r sim of (X(u 1 ), . . . , X(u m )) and a random number u ∈ [0, 1] drawn from the uniform distribution on the unit interval.
(ii) If u < f data (s(r sim )) M f model (s(r sim )) , the vector r sim is accepted as a realization of the rejection sampling scheme, otherwise go back to step (i).
Note that the sphericity distribution of outer shells generated by rejection sampling matches f data perfectly. However, since this modeling approach solely considers the sphericity value it does not necessarily lead to a good match with distributions of other shape characteristics. On the other hand, the model {X(u) : u ∈ S 2 } for the outer shell described in the present paper was not fitted directly to size and shape characteristics of particles extracted from CT data, but to the representation of outer shells as radius functions which were evaluated on a mesh of the unit sphere's surface. In other words, the model {X(u) : u ∈ S 2 } was fitted directly to a (approximative) representation of the outer shell which is more exhaustive than the aggregated sphericity value. This issue of the rejection sampling scheme described above could be remedied by considering multivariate probability densities which incorporate further particle shape characteristics. In this manner more informative vectors of such characteristics can be used to represent outer shells. However, this can lead to an increase in the value M which makes the simulation procedure computationally unfeasible. This is due to the fact, that, on average, step (i) in the rejection sampling scheme has to be performed M times until a single realization is accepted.
Furthermore, note that due to the choice of the parametric covariance functions (see Eq. (20)) for the individual components of the mixture {X(u) : u ∈ S 2 } of Gaussian random fields, the realizations of the outer shell model exhibit smooth surfaces-this can also be seen by the representation of the corresponding isotropic Gaussian random fields given in Eq. (54). However, in the context of the present paper, the smoothness of realizations of the outer shell model is no issue, since the realizations reassemble the particles depicted in CT data reasonably well. For Gaussian random fields on the sphere, the realizations of which exhibit a rough surface, we refer the reader to, 30 where various parametric covariance functions are considered to generate fractal surfaces. Alternatively, surface roughness can be introduced by additional subsequent modeling steps, e.g., the outer shell of realizations of the presented multi-scale model exhibit faceted surfaces, see Figure 8c.
As described at the end of the Results section, the general fit of the grain architecture model is reasonably good. Nevertheless, also with respect to this component of our multi-scale model, modifications and extensions are possible: Note that most of the aforementioned discrepancies can be explained by the differences regarding the seed point pattern of the fitted Laguerre generators and the ones from the simulated data, and the thus resulting differences regarding the nearest neighbor distributions. An obvious solution is to fit the parameters of the Matérn cluster process by the minimum contrast method with respect to the nearest neighbor distance distribution instead of the pair correlation function. While this improves the quality-of-fit of the nearest neighbor distance distribution, it degrades the one for the pair correlation function and pronounces the trend of the simulated cells being too large. So the general fit of the simulated grain architecture would be worse compared to the results obtained by the approach proposed in the present paper.
Another way to further improve the fit would be to employ a more sophisticated (seed) point process model. For example, instead of the Matérn cluster process, another type of random point processes could be considered, such as a suitably chosen Gibbs point process, 23 that could capture both the pair correlation and the nearest neighbor distances of the fitted Laguerre generators better. However, fitting the parameters of Gibbs point processes is computationally much more expensive than fitting the parameters of a Matérn cluster process. Additionally, formulas for the nearest neighbor distance distribution for more sophisticated types of point processes (including Gibbs processes) are hard to derive theoretically, and, since they are necessary to model the conditional probability density of additive weights, must therefore be approximated, e.g., through Monte-Carlo techniques. This would introduce a further source of error and complexity.
In a forthcoming research, it will be investigated whether modeling both, the outer shells and the grain architectures, conditionally on each other (e.g. by having smaller grains near the particle surfaces) improves the quality of the fit. For this, however, further FIB-EBSD measurements are necessary. The presented technique overcomes systematic limitations of standalone techniques, such as the inability of X-ray nano-CT to capture grain boundaries, and the challenge of imaging large volumes and multiple particles with FIB-EBSD, by leveraging the respective strengths of both techniques to generate equivalent high-resolution and numerous virtual particle architectures. This represents a significant step forward to acquiring 3D single particle architectures with effective resolution in the range of 10's of nanometers. With this technique, electrode manufactures can build representative 3D images of the particles that they synthesize using two accessible lab-based techniques-X-ray nano-CT and FIB-EBSD. Input to the particle generator can be tuned to produce particles of different grain or particle sizes or shapes. This capability of generating representative particles facilitates exploration of different particle architectures through multi-physics modeling with respect to power density and the propensity of the particle to crack due to intra-grain mechanical strain during operation, thus potentially elucidating favorable particle architectures for specific operating conditions, e.g., architectures with minimal mechanical strain and maximum rate of lithium transport. Therefore, modeling various particle architectures will be the focus of future work in order to clarify the influence of grain properties on cell performance and also to analyze the influence of interfacial transport resistances and grain orientations on the performance and degradation of particles.

Computational details
A random field model used to characterize and generate the outer shells, and a random tessellation model used to characterize and generate grain architectures, are combined to form a multi-scale model for the generation of virtual electrode particles with full grain detail. This methodological approach demonstrates the possibility of generating representative single electrode particle architectures for modeling and characterization that can guide synthesis approaches of particle architectures with enhanced performance. In particular, this was achieved by drawing statistical representations of particle surfaces (outer shells) from X-ray nano-CT data, and sub-particle grain detail from FIB-EBSD data which has higher resolution and much greater sensitivity to grain boundaries. The random field model for generating particle surfaces, or outer shells, was calibrated using the gradient descent implementation of the Python package TensorFlow 61 and showed a good fit to the experimental X-ray nano-CT data. Additionally, TensorFlow allowed us to fit Laguerre tessellations to experimental FIB-EBSD data in a highly parallel or even GPU-accelerated manner. The resulting fits showed a good agreement and were then used to calibrate a random tessellation model for the grain architecture. The two models were combined to generate full NMC particles with statistically representative shapes and sizes of particle boundaries and grain architectures.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding authors on reasonable request.