# Stochastic inference of effective material properties from crystallographic textures via Gaussian process regression

Bingqian Li, Piotr Breitkopf, Ludovic Cauvin

► To cite this version:

Bingqian Li, Piotr Breitkopf, Ludovic Cauvin. Stochastic inference of effective material properties from crystallographic textures via Gaussian process regression. Computer Methods in Applied Mechanics and Engineering, 2026, 453, 10.1016/j.cma.2026.118825. hal-05521213

HAL Id: hal-05521213

https://utc.hal.science/hal-05521213v1

Submitted on 20 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Stochastic inference of effective material properties from crystallographic textures via Gaussian process regression

Bingqian Li ID, Piotr Breitkopf ID, Ludovic Cauvin ID*

Université de Technologie de Compiègne, Alliance Sorbonne Université, Roberval, Centre de recherche Royallieu, CS 60319, Compiègne, France

## ARTICLE INFO

**Keywords:**

Stochastic homogenization
Crystalline texture
Gaussian process regression
Uncertainty propagation
EBSD

## ABSTRACT

This work investigates the impact of measurement uncertainty in polycrystalline texture data on homogenized material properties through full-field homogenization. Crystallographic texture, represented by Pole Density Functions, is modeled through Gaussian Process (GP) regression, in which measurement error is explicitly introduced via the noise term of the GP. Two uncertainty propagation approaches are examined. In the first, deterministic descriptors derived from the GP posterior distribution, namely the posterior mean and the bounds of the 95% confidence interval, are used as texture inputs for homogenization. In the second, Monte Carlo Pole Density Function instances are sampled directly from the GP posterior distribution and propagated through homogenization. Numerical experiments on a two-dimensional thermal conduction problem use a synthetic texture. The effective properties and their variability obtained from both approaches are compared with deterministic predictions based on the error-free reference texture. The results show that sampling-based homogenization provides a more robust estimation of effective physical properties and their associated uncertainty.

## 1. Introduction

Microscopic imaging techniques now provide increasingly detailed microstructural information, among which polycrystalline texture is a key descriptor governing anisotropic behavior. The objective of this work is to showcase a general framework for the characterization of the stochastic distribution of effective material properties corresponding to Electron Backscatter Diffraction (EBSD)-derived textures [1,2] by explicitly modeling measurement errors [3,4] in the homogenization.

Full crystalline orientations measured by EBSD, expressed as triplets of Euler angles, describe the density of orientations in Euler space via Orientation Distribution Function (ODF). In homogenization, discretization strategies employing ODF have been proposed, such as using a finite set of orientations [5] or discretizing Rodrigues space [6]. Specifically, for a given slip system, the orientation of slip planes can be represented by their normal vectors, referred to as poles, and visualized as points on the unit sphere, defining the Pole Density Function. Alternatively, the poles can be characterized experimentally using X-ray or neutron diffraction [7].

Texture-dependent mechanical responses such as Young's modulus [8] and plasticity [9,10], as well as physical phenomena such as thermal conductivity [11] and spin-orbit torque [12,13] are widely documented.

Homogenization mean-field models rely on statistical averaging and do not resolve the explicit microstructural geometry [14,15], whereas full-field approaches account for grain-scale heterogeneity by explicitly constructing a Representative Volume Element (RVE) [16] sufficiently large [17] at the expense of computational cost. Kanit et al. [18] proposed a statistical criterion for the RVE size,

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: bingqian.li@utc.fr (B. Li), piotr.breitkopf@utc.fr (P. Breitkopf), ludovic.cauvin@utc.fr (L. Cauvin).</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. From texture to Pole figure for a given slip system.</figcaption></figure>

based on the convergence of apparent properties with respect to size, the number of realizations, and the desired precision. Under such conditions, effective properties can be estimated by averaging the response of several stochastic RVEs of smaller size [19]. Combined with digital material reconstruction techniques for RVE generation [20,21], homogenization strategies are widely employed to predict effective properties from microstructural data.

The Pole Density Function, defined as a spherical function on the unit sphere, is commonly represented using spherical harmonic expansions [22,23]. When Pole Density Functions corresponding to multiple slip systems are available, they may be inverted to re-construct the ODF [24]. Since the ODF is defined in Euler-angle space, it is typically expressed using generalized spherical harmonics. The quantification of measurement uncertainty has been investigated in [25,26], where a maximal coefficient of variation of approximately 15% in pole density measurements was reported. However, sphere harmonics *per se* do not provide a consistent framework for error propagation.

Gaussian Processes (GP) provide a flexible stochastic framework for modeling noisy observations through the specification of a prior mean function and a covariance kernel [27]. They have been widely employed to construct surrogate models of experimental data [28], as well as to establish direct statistical links between microstructural descriptors and macroscopic responses, such as plastic deformation behavior [29] or effective mechanical properties [30]. These surrogate-based approaches typically rely on dimensionality-reduction techniques to represent the microstructure, sometimes also the response in a reduced space [31]. While effective in lowering computational costs, such strategies may suffer from limited interpretability and reduced generalizability outside the training domain.

In contrast to GP surrogate models, based on the customized Gaussian Process Regression (GPR) framework proposed in [32], present approach employs GPR solely as a stochastic representation of crystallographic texture, with measurement uncertainty explicitly introduced at the level of Pole Density Functions by a noise parameter. The resulting GP posterior distribution defines a stochastic description of admissible Pole Density Functions and enables the sampling of random instances [33], which can be propagated through physics-based homogenization or integrated into other stochastic learning frameworks [34].

In the present study, this GPR-based description of texture is coupled with a stochastic homogenization strategy to quantify the impact of measurement uncertainty on effective material properties. Pole Density Functions derived from the GP posterior, namely the posterior mean, the upper and lower bounds of the 95% confidence interval, and randomly sampled instances, are used to assign crystalline orientations within RVEs. Homogenization performed on these RVEs, combined with Monte Carlo simulations over the sampled Pole Density Functions, yields a stochastic distribution of effective properties induced by texture measurement uncertainty. All results are evaluated relative to a reference homogenization performed directly on the synthetic texture data.

The paper is organized as follows. Section 2 describes the methodology and formulates the three homogenization strategies based on the GPR representation of the Pole Density Function as observed, GPR mean and bounds and Monte Carlo sampling. Section 3 applies the finite element homogenization of thermal conductivities to a synthetic texture with the three approaches. Numerical results are then presented and discussed, focusing on the resulting distribution of effective conductivities under varying measurement error levels. Section 4 concludes the study.

## 2. Methodology

The goal of present work is accounting for texture information in inferring polycrystalline materials' properties. Texture information may be acquired by Electron Backscatter Diffraction (EBSD), based on Scanning Electron Microscopy (SEM), measuring three-dimensional orientations of crystals by analyzing the sample surface.

In Fig. 1a, the cells represent grains, and the hexagonal prisms indicate crystal lattices individually oriented according to Euler angles {$\phi_1, \Phi, \phi_2$} [7].

For a given grain orientation and slip system h, the slip plane normal vector is projected as a single point $y = \{az, el\}$ onto the unit sphere surface $S^2$ (Fig. 1b),

$$y \in S^2, \text{ where } S^2 = \{x \in \mathbb{R}^3 \mid \|x\| = 1\}.$$

By projecting the slip plane orientations of all the grains in the sample onto $S^2$, a discrete set of points is obtained, each associated with a measured density of grains sharing the same orientation and visualized as the discrete *Observed Pole Figure (PF)* (Fig. 1c). To facilitate interpolation, continuous spherical function $P_h(y)$ can be approximated from the Observed PF, referred to as a *Pole Density*
<figure><img src="image_4.png" /><figcaption>Fig. 2. FE homogenization based on an individual PF.</figcaption></figure>

Function satisfying

$$ \frac{dV_y}{V} = \frac{1}{4\pi} P_h(y) dy, \qquad (1) $$

where $dV_y$ denotes the volume of grains oriented towards y, and V is the total volume of the measured region. As a direct consequence, $P_h(y)$ is non-negative, and its integral over the unit sphere is equal to $4\pi$.

Homogenization techniques constitute a well-established framework for evaluating the macroscopic response of polycrystalline materials. The central concept is the Representative Volume Element (RVE) Ω reproducing the relevant microstructural features of the material, such as the texture. Once the size of the RVE is properly defined, the material and the RVE are expected to exhibit equivalent overall responses and properties at the macroscopic scale. The macroscopic response $f^{hom}$ can then be computed analytically or numerically using the Finite Element Method (FEM) as the volume average of the local response field $f(x)$, with $x \in \Omega$:

$$ f^{hom} = \langle f(x) \rangle, \quad \langle \cdot \rangle = \frac{1}{|\Omega|} \int_{\Omega} \langle \cdot \rangle d\Omega, \qquad (2) $$

where $|\Omega|$ denotes the volume of the RVE. In finite element homogenization, the computational complexity of the polycrystalline RVE primarily depends on the number of grains as well as the level of mesh refinement of each grain.

## 2.1. RVE built from pole density function

For a given texture, the RVE should be adjusted to produce a unique macroscopic response independently of random orientation mapping. The overall process is schematically illustrated in Fig. 2, where the Observed PF is used as an example to construct the RVE and perform homogenization, leading to a single material response. The material response shown on the right can be any type of linear/non-linear response related to the crystallographic texture: strain-stress, thermal conductivity, etc., obtained through the chosen homogenization framework.

The creation of the RVE involves two major steps: the generation of the microstructural geometry and the assignment of material properties. In the case of the polycrystalline material studied in the present work, the microstructure is represented by a collection of grains. Since the material properties at the grain level are assumed orientation-dependent, the second step reduces to assigning orientations to the grains according to the given texture, described by the Pole Density Function.

In the continuous form, for a given $P_h(y)$ and assuming that the n individual grains of the RVE have uniform size, the assignment of orientations consists of selecting a set $y^* = \{y_i^* | i = 1, ..., n\}$ such that

$$ \frac{1}{n} \sum_{i=1}^{n} \delta(y_i^*, y) \approx \frac{1}{4\pi} P_h(y), \quad y \in S^2, \qquad (3) $$

where $\delta$ denotes the Kronecker delta. According to Eq. (1), the normalization

$$ \int_{S^2} \frac{1}{4\pi} P_h(y) dy = 1 $$

allows interpreting $\frac{1}{4\pi} P_h(y)$ as a probability density function.

In the discrete form, as in the case of the Observed PF, consider a set of discrete slip plane orientations $y = \{y_1, y_2, ..., y_N\}$ with a corresponding measured Pole Density Vector $d = \langle d_1, d_2, ..., d_N \rangle^T$, satisfying $\sum_{i=1}^{N} d_i = 1$. Then, Eq. (3) becomes

$$ \frac{1}{n} \sum_{i=1}^{n} \delta(y_i^*, y_j) \approx d_j. $$

In practice, an acceptance-rejection algorithm is employed for the orientation sampling. The algorithm for the discrete case is detailed in Algorithm 1.

## 2.2. RVE built from GP mean and bounds

Approximation of Pole Density Functions typically relies on spherical harmonic expansions [22], neglecting measurement uncertainty, and potentially yielding non-physical negative density values.
<caption>**Algorithm 1: Random sampling of crystalline orientations following a Pole Density Vector.**</caption>

Input: **y**, **d**, *n* (number of grains in RVE)
Output: **y**\* = {y\*<sub>1</sub>, y\*<sub>2</sub> ..., y\*<sub>n</sub>}
1 *i* = 1;
2 **while** *i* < *n* + 1 **do**
3     *j* ← Random integer with 1 ≤ *j* ≤ *N*;
4     *p* ← Random scalar with 0 < *p* < 1;
5     **if** *p* < *d*<sub>*j*</sub> **then**
6         *y*\*<sub>*i*</sub> = *y*<sub>*j*</sub>;
7         *i* = *i* + 1;
8 **return** **y**\* = {y\*<sub>1</sub>, y\*<sub>2</sub> ..., y\*<sub>n</sub>}

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** GPR-based $P_h(y)$ estimate, posterior mean ($M(P_h)$), lower and upper bounds of 95% CI (LB/U $B_{0.95}(P_h)$).</figcaption></figure>

To overcome these limitations, a customized Gaussian Process (GP) framework was introduced in previous work [32].

$$P_h(y) \sim GP(m(y), \sigma_j^2 k(\text{dist}(y, y'), \ell)) \quad (4)$$

where the optimal configuration was identified as a zero prior mean $m(y) = 0$, a Matérn covariance kernel with $\nu = 5/2$,

$$k(\text{dist}(y, y'), \ell) = \left(1 + \frac{\sqrt{3} \text{dist}(y, y')}{\ell} + \frac{5 \text{dist}(y, y')^2}{3\ell^2}\right) \exp\left(-\frac{\sqrt{5} \text{dist}(y, y')}{\ell}\right). \quad (5)$$

and a spherically periodic geodesic distance $\text{dist}(y, y')$, geometrically defined as the shortest arc on the unit sphere connecting orientations $y$ and $y'$, ensuring spherical periodicity. The hyperparameters {$\sigma_f, \ell$}-representing the prior variance and characteristic length scale-are optimized by maximizing the marginal likelihood.

The non-negativity of the estimated density is ensured by the log-linear data transformation. The original density $d$ is first transformed to $d_{\text{lll}}$ (Eq. (6)),

$$d_{\text{lll}} = f_{\text{lll}}(d) = \begin{cases} \ln(d) & \text{if } d \le d_0, \\ \frac{1}{d_0} d + \ln(d_0) - 1 & \text{if } d > d_0, \end{cases} \quad (6)$$

on which the GP model is applied. The posterior values are then mapped back to the original space (Eq. (7)).

$$d = f_{\text{lll}}^{-1}(d_{\text{lll}}) = \begin{cases} \exp(d_{\text{lll}}) & \text{if } d_{\text{lll}} \le \ln(d_0), \\ d_0 \cdot (d_{\text{lll}} - \ln(d_0) + 1) & \text{if } d_{\text{lll}} > \ln(d_0). \end{cases} \quad (7)$$

When applied to a vector, $f_{\text{lll}}$ and $f_{\text{lll}}^{-1}$ operate elementwise, i.e., the transformation is performed independently on each component.

Finally, the GP-based function reconstruction is scaled such that its integral over the unit sphere equals $4\pi$. When the function is discretized to a density vector **d**, the normalization condition is simplified to $\sum_i d_i = 1$.

Gaussian Process Regression (GPR) takes the posterior mean function as the regression estimate. This deterministic reconstruction of the Pole Density Function is referred to as the *Mean Pole Density Function*, denoted by $M(P_h)$ in this paper. The posterior mean density vector $\tilde{\mathbf{d}}$ evaluated at the measurement points $y = \{y_1, y_2, \dots, y_N\}$ is

$$\begin{aligned} \tilde{\mathbf{d}} &= f_{\text{lll}}^{-1}(\tilde{\mathbf{d}}_{\text{lll}}) \\ &= f_{\text{lll}}^{-1}(\mathbf{C}(\mathbf{C} + \sigma_n^2 \mathbf{I}))^{-1} f_{\text{lll}}(\mathbf{d}), \end{aligned} \quad (8)$$

where $\mathbf{C}$ is the $N \times N$ covariance matrix,

$$C_{ij} = \sigma_j^2 k(\text{dist}(y_i, y_j), \ell), \quad i, j = 1, \dots, N, \quad (9)$$

and $\sigma_n$ is the noise level modeling measurement uncertainty, expressed as a fraction $\alpha$ of the maximum observed density:

$$\sigma_n = \alpha \cdot d_{\text{max}}, \quad \text{with} \quad d_{\text{max}} = \max_{i=1, \dots, N} f_{\text{lll}}(d_i). \quad (10)$$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Material response estimation based on deterministic GP-based Pole Density Functions’ mean and CI.</figcaption></figure>

Eq. (8) shows that the error is assumed to act on the measured density values, not on the orientations. Note that for $\alpha \neq 0$, $\tilde{d} \neq d$, meaning that the posterior mean is not interpolant-it does not reproduce exactly the rough measured densities.

Confidence Intervals (CI) of the GP model can also be derived. For a 95% CI, the lower and upper bounds, denoted by $LB/UB_{0.95}(P_h)$, correspond approximately to $\pm 1.96\sigma$, with $\sigma$ the posterior standard deviation. $LB/UB_{0.95}(P_h)$ can thus be written as

$$CI_{0.95} = [LB_{0.95}(P_h), UB_{0.95}(P_h)] \\ = [f_{\text{II}}^{-1}(\tilde{d}_{\text{II}} - 1.96 \tilde{\sigma}_{\text{II}}), f_{\text{II}}^{-1}(\tilde{d}_{\text{II}} + 1.96 \tilde{\sigma}_{\text{II}})], \quad (11)$$

where $\tilde{\sigma}_{\text{II}}$ is the vector of posterior standard deviations at evaluation points. For a given $y \in S^2$,

$$\tilde{\sigma}_{\text{II}}(y) = \sqrt{\sigma_j^2 k(y, y) - c^\top (C + \sigma_n^2 I)^{-1} c}, \quad (12)$$

with $c$ the $N \times 1$ vector of covariance between $y$ and the measurement points:

$$c_i = \sigma_j^2 k(y, y_i), \quad i = 1, \dots, N. \quad (13)$$

Because of the log-linear data transformation (Eqs. (6)-(7)), the posterior CI is also strictly non-negative:

$$CI \subset [0, +\infty). \quad (14)$$

Building on the expressions above, the complete procedure for obtaining the deterministic reconstructions of $M(P_h)$ and $LB/UB_{0.95}(P_h)$ illustrated in Fig. 3 is summarized in Algorithm 2.

<table><caption>Algorithm 2: GP model reconstruction for $M(P_h)$ and CI bounds.</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>Input:</td><td>y, d, α</td><td></td></tr><tr><td>Output:</td><td>tilde d (M(P<sub>h</sub>)), LB/UB<sub>0.95</sub>(P<sub>h</sub>)</td><td></td></tr><tr><td>1</td><td>D<sub>N×N</sub> ← D<sub>ij</sub> = dist(y<sub>i</sub>, y<sub>j</sub>);</td><td>/* Distance matrix */</td></tr><tr><td>2</td><td>{σ<sub>f</sub><sup>*</sup>, ε<sup>*</sup>} ← arg max<sub>{σ<sub>f</sub><sup>*</sup>, ε<sup>*</sup>}</sub> L(<b>D</b>, d, σ<sub>f</sub><sup>2</sup> k(dist(y, y'), ε'));</td><td>/* Hyperparameters optimization by maximizing marginal likelihood */</td></tr><tr><td>3</td><td>C ← C<sub>ij</sub> = σ<sub>f</sub><sup>*</sup> k(D<sub>ij</sub>, ε<sup>*</sup>);</td><td>/* Covariance matrix */</td></tr><tr><td>4</td><td>d<sub>II</sub> ← d;</td><td>/* log-linear data transformation Eq. (6) */</td></tr><tr><td>5</td><td>σ<sub>n</sub> ← α, d<sub>II</sub>;</td><td>/* Noise level hyperparameter Eq. (10) */</td></tr><tr><td>6</td><td>d̃<sub>II</sub> ← C, σ<sub>n</sub>, d<sub>II</sub>;</td><td>/* Posterior mean Eq. (8) */</td></tr><tr><td>7</td><td>d̃ ← d̃<sub>II</sub>;</td><td>/* M(P<sub>h</sub>), Eq. (7) */</td></tr><tr><td>8</td><td>σ̃<sub>II</sub> ← σ<sub>f</sub><sup>*</sup>, σ<sub>n</sub>, k(dist(y, y'), ε), y;</td><td>/* Eq. (12) */</td></tr><tr><td>9</td><td>LB/UB<sub>0.95</sub>(P<sub>h</sub>) ← d̃<sub>II</sub>, σ̃<sub>II</sub>;</td><td>/* Eq. (11) */</td></tr><tr><td>10</td><td>return d̃ (M(P<sub>h</sub>)), LB/UB<sub>0.95</sub>(P<sub>h</sub>)</td><td></td></tr></tbody></table>

Fig. 4 illustrates the homogenization procedure based on the GP mean and CI, which requires three simulations, one for each approximated $P_h(y)$. Although the three responses depicted in this schematic may appear to form a distribution of the material response, this interpretation is misleading. The Gaussian nature of the Pole Density Functions in the transformed space is not necessarily preserved after inversion and homogenization, and therefore the resulting material responses do not, in general, constitute a statistically meaningful distribution. Besides, although $M(P_h)$ and $LB/UB(P_h)$ provide a guideline of the GP posterior distribution, neither the mean nor the bounds necessarily correspond to physically realizable Pole Density Functions. From a stochastic perspective, the probability of obtaining a deterministic estimate of $P_h(y)$ is negligible. Consequently, while the corresponding homogenized responses may still offer illustrative insights, the inherent stochasticity is not represented, and limited quantitative assessment of confidence can be made for the resulting estimates. This limitation motivates the investigation of stochastic sampling from the posterior distribution of $P_h(y)$, which is addressed in the next section.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Random sampling of admissible posterior GP instances.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Monte Carlo simulations based on $S(P_h)$ revealing the variability of the material response.</figcaption></figure>

### 2.3. RVE sampling from posterior GP distribution

Monte Carlo samples generated from the posterior GP distribution respect the covariance structure defined by the optimized covariance function, and correspond to admissible Pole Density Functions capturing the uncertainty of the texture. They are referred to as *Sampled Pole Density Functions*, denoted by $S(P_h)$ in this work (Fig. 5).

GP instances are sampled discretely from a multivariate normal distribution defined by the posterior mean $\tilde{d}_{ll}$ and the posterior covariance matrix $\Sigma$ at the evaluation points,

$$f_{ll}(S(P_h)) \sim \mathcal{N}(\tilde{d}_{ll}, \Sigma), \qquad (15)$$

where

$$\Sigma = \frac{1}{\sigma_f^2} \text{diag}(\tilde{\sigma}_{ll}) C \text{diag}(\tilde{\sigma}_{ll}). \qquad (16)$$

To generate one $S^{(i)}(P_h)$,

$$S^{(i)}(P_h) = f_{ll}^{-1}(\tilde{d}_{ll} + Lx^{(i)}), \qquad (17)$$

where $L$ is the lower-triangular matrix from the Cholesky decomposition

$$\Sigma = LL^T,$$

and $x^{(i)}$ is a sample from the standard multivariate normal distribution of matching dimension,

$$x \sim \mathcal{N}(0, I).$$

Overall sampling process is given in Algorithm 3.

<table><caption>Algorithm 3: $S(P_h)$ sampling from the GP model.</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>Input:</td><td>y</td><td>d</td><td>α</td></tr><tr><td>Output:</td><td>$S^{(i)}(P_h)$</td><td></td><td></td></tr><tr><td>1</td><td>$\tilde{d}_{ll}$, $\tilde{\sigma}_{ll}$, $\sigma_f^*$, $C$ ← y, d, α;</td><td></td><td>/* Algorithm 2 */</td></tr><tr><td>2</td><td>$\Sigma \leftarrow \sigma_f^*$, $\tilde{\sigma}_{ll}$, C;</td><td></td><td>/* Cov. matrix Eq. (16) */</td></tr><tr><td>3</td><td>$L \leftarrow \Sigma = LL^T$;</td><td></td><td>/* Cholesky decomposition */</td></tr><tr><td>4</td><td>$x^{(i)} \leftarrow x \sim \mathcal{N}(0_{N \times 1}, I_{N \times N});$</td><td></td><td>/* Random sampling */</td></tr><tr><td>5</td><td>$S^{(i)}(P_h) \leftarrow x^{(i)}, L, \tilde{d}_{ll}$;</td><td></td><td>/* Eq. (17) */</td></tr><tr><td>6</td><td>$\text{return } S^{(i)}(P_h)$</td><td></td><td></td></tr></tbody></table>
It is important to note that, in the present work, since the analysis is restricted to the measurement points and the hyperparameters of the GP model are kept fixed at their optimized values, the posterior variability of the Sampled Pole Density Functions $S(P_h)$ is primarily governed by the prescribed noise level $\sigma_n$ (Eq. (12)).

With a potentially infinite number of admissible $S(P_h)$ available, a Monte Carlo simulation (MCS) can be employed to quantify the variability arising from texture uncertainty and measurement error. Each $S(P_h)$ is used to assign crystalline orientations within the RVE, whose homogenization yields a corresponding material response (Fig. 6). The ensemble of responses is then analyzed statistically.

Determining the number $m$ of samples is essential before proceeding. One might expect that, since Eq. (15) defines the $S(P_h)$ as following a multivariate normal distribution with a known posterior mean, $m$ could be determined analytically. However, this is valid only in the transformed log-linear space. The inverse transformation $f_{\text{III}}^{-1}(x)$ is a non-linear and convex function in its logarithmic region-particularly for densities that were initially close to zero. According to Jensen's inequality [35],

$$ \forall y, \quad f_{\text{III}}^{-1}(\tilde{d}_{\text{III}}(y)) = M(P_h, y) \leq \mathbb{E}[S^{(i)}(P_h, y)], \quad i = 1, \dots, \infty. \quad (18) $$

This means that $M(P_h)$ no longer coincides with the mean of the $S(P_h)$ in the original space. Consequently, the multivariate normal assumption becomes invalid, and $m$ cannot be estimated analytically. Therefore, it is proposed to determine $m$ numerically based on the convergence of the resulting effective responses (Algorithm 4).

<caption>**Algorithm 4: Monte Carlo on $S(P_h)$.**</caption>

Input: $y, d, \alpha, \xi$ (tolerance)
Output: $m$ (number of $S(P_h)$ required for convergence)
1 $S^{(1)}(P_h) \leftarrow y, d, \alpha$; /* A first $S(P_h)$, Algorithm (3) */
2 $\langle f(x) \rangle^{(1)} \leftarrow S^{(1)}(P_h)$; /* Homogenized response */
3 $\tilde{f} = \langle f(x) \rangle^{(1)}$; /* Mean response initialization */
4 $a = 1$; /* Initial relative difference (100%) */
5 $i = 2$;
6 while $a > \xi$ do
7 $S^{(i)}(P_h) \leftarrow y, d, \alpha$; /* Sample a new random $S(P_h)$ */
8 $\langle f(x) \rangle^{(i)} \leftarrow S^{(i)}(P_h)$;
9 $\tilde{f}' = ((i-1)\tilde{f} + \langle f(x) \rangle^{(i)}) / i$; /* New mean response */
10 $a = |\tilde{f}' - \tilde{f}| / |\tilde{f}|$; /* Difference between the two means */
11 $\tilde{f} = \tilde{f}'$; /* Update running mean */
12 $i = i + 1$;
13 $m = i - 1$;
14 return $m$

## 3. Numerical experiments

The proposed methodology is now applied to a synthetic example for quantitative analysis. In this section, developed methodology is illustrated on a two-dimensional thermal conduction problem.

### 3.1. Synthetic data

The in silico texture measurement introduced in [32] is used in the present study. A single measurement dataset is considered, consisting of 1500 distinct orientations randomly distributed over the surface of the unit sphere. These orientations are interpreted as an Observed PF (Fig. 7), with associated density values assumed to be free from error.

In real experiments, measurement uncertainty can be quantified through repeated tests followed by uncertainty analyses. The error related to instrument precision are typically specified by the equipment manufacturer. In the present in silico example, however, the error level, denoted by $\alpha$ in Eq. (10), is arbitrarily prescribed with I.I.D. Gaussian noise $\sigma_n$ for this single dataset, aiming to investigate its influence on the homogenized material response.

Applying GPR to the Observed PF, the hyperparameters $\{\sigma_f, \ell\}$ in Eq. (4) are optimized by maximizing the marginal likelihood and reported in Table 1 (see [32] for details). These hyperparameters are reevaluated following the log-linear transformation and kept fixed throughout all numerical experiments, with only the error level $\alpha$ being adjusted in the GP model. Raissi and Karniadakis [36] suggest that the noise can also be considered as a hyperparameter to be optimized. However, this is of particular interest if the experimental data is error-integrated. In the present case, since the error level is artificially prescribed, it is treated as a known input to the model.

### 3.2. Homogenization of thermal conductivity

Present framework is applied on a 2D thermal conduction problem aiming to evaluate the effective conductivity, in which case, the following assumptions are adopted throughout the in silico experiments:
<figure><img src="image_9.png" /><figcaption>Fig. 7. In silico observed PF of 1500 measured orientations with density (upper view of the unit sphere).</figcaption></figure>

<caption>Table 1</caption>

Optimized hyperparameters of the GP model for the synthetic data [32].

<table><thead><tr><th colspan="2">Optimized hyperparameters for Matérn (ν = 5/2)</th></tr></thead><tbody><tr><td>σ<sub>f</sub> = 1.258</td><td>ε = 0.1467</td></tr></tbody></table>

* the material is of single-phase, hexagonal close-packed (HCP) lattice structure;

* the texture is characterized by the Pole Density Function of the basal slip plane family {0001};

* at the grain scale, thermal conductivity is three-dimensional transversely isotropic (polar anisotropic), with the axis of transverse isotropy aligned with the [0001] crystallographic direction;

* grains are of uniform size;

* the measurement error affects only the measured densities but not the orientation.

In the general case, a full ODF is required to characterize the texture of polycrystalline materials. However, for mechanical behavior, only specific crystallographic slip systems effectively contribute to the macroscopic response. In such cases, the Pole Density Functions associated with these slip systems contain the necessary and sufficient texture information to determine the homogenized response.

Although the numerical experiments in this work are conducted on a 2D problem, the local thermal conductivity, assumed to be transversely isotropic, is inherently three-dimensional. Its projection onto the homogenization plane is fully determined by an orientation $y \in S^2$. As a result, a single Pole Density Function is sufficient to compute the homogenized thermal conductivity while preserving the crystalline symmetry of the material. The methodology proposed here readily extends to more general situations by considering Pole Density Functions associated with multiple slip systems, which can be inverted to reconstruct the corresponding ODF [24] and subsequently used for orientation assignment in the RVE.

The local 3D conductivity tensor of each orientation-assigned grain is projected onto the X-Y homogenization plane (Fig. 1), and is assigned to the corresponding finite elements in the RVE.

Under the transverse isotropy assumption, the local thermal conductivity tensor in the crystallographic coordinate system is chosen as

$$ \mathbf{k}_0 = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 100 \end{bmatrix} W/(m \cdot K). \qquad (19) $$

Let $y \in S^2$ be a unit vector representing the orientation of the basal plane normal for a given grain. The thermal conductivity tensor in the global coordinate system is then

$$ \mathbf{k}(y) = \mathbf{R} \mathbf{k}_0 \mathbf{R}^{\mathrm{T}} = \begin{bmatrix} k_{XX} & k_{XY} & k_{XZ} \\ k_{YY} & k_{YY} & k_{YZ} \\ \mathrm{sym} & k_{ZZ} & k_{ZZ} \end{bmatrix}. \qquad (20) $$

R is a rotation matrix satisfying $\mathbf{Rn} = y$, with $\mathbf{n} = (0, 0, 1)^{\mathrm{T}}$ the local normal to the basal plane.

Since the focus is on 2D thermal conduction in the X-Y plane, only the in-plane components of $\mathbf{k}(y)$, denoted by $\mathbf{k}_{XY}(y)$ are retained and assigned to each grain.

Once the orientations assigned, the projected properties of each grain are determined following Eq. (20). Then, the effective properties of the RVE are computed using the Uniform Temperature Gradient (UTG) boundary condition [16]

$$ T(\mathbf{x}) = \tilde{\mathbf{G}} \cdot \mathbf{x}, \quad \forall \mathbf{x} \in \partial\Omega, \qquad (21) $$
where $\tilde{G}$ denotes the prescribed macroscopic temperature gradient applied on the boundary $\partial\Omega$ of the RVE. Two simulations are carried out with

$$ \tilde{G}_1 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \quad \tilde{G}_2 = \begin{bmatrix} 0 \\ 1 \end{bmatrix}, \qquad (22) $$

resulting in temperature fields $T^{(1)}(x)$ and $T^{(2)}(x)$, $x \in \Omega$. The slip plane's orientation of a point $x$ is given by $y(x)$. The localization tensor of the temperature gradients,

$$ A(x) = \begin{bmatrix} \frac{\partial T^{(1)}}{\partial X}(x) & \frac{\partial T^{(2)}}{\partial X}(x) \\ \frac{\partial T^{(1)}}{\partial Y}(x) & \frac{\partial T^{(2)}}{\partial Y}(x) \end{bmatrix}, \qquad (23) $$

yields the effective conductivity tensor

$$ \langle \mathbf{k} \rangle = \langle \mathbf{A}^{\mathrm{T}}(\mathbf{x}) \mathbf{k}_{XY}(y(\mathbf{x})) \mathbf{A}(\mathbf{x}) \rangle \qquad (24) $$

$$ = \begin{bmatrix} \langle k_{XX} \rangle & \langle k_{XY} \rangle \\ \langle k_{XY} \rangle & \langle k_{YY} \rangle \end{bmatrix}, \qquad (25) $$

for an arbitrary Pole Density Function $P_h(y)$.

The whole process is summarized in Algorithm 5.

<caption>Algorithm 5: Homogenization of in-plane effective conductivity using uniform temperature gradient (UTG) load cases, T3 element.</caption>

Input: $k_0$, Node coordinates, Element connectivity, $y^*$
Output: $\langle \mathbf{k} \rangle$
/* Pre-processing: material assignment */
1 foreach grain g do
2 $\mathbf{k}_{XY}(y_g^*) \leftarrow k_0, y_g^*$;
3 $\forall$ element e in g do
4 $\left[ \mathbf{k}_e \leftarrow k_{XY}(y_e^*) \right]$;
/* Element in-plane conductivity */
/* Node temperature calculated by EF under two UTG B. C. */
5 $U_1, U_2 \leftarrow$ EF simulation with $\tilde{G}_1$ and $\tilde{G}_2$;
/* Post-processing: computation of effective conductivity */
6 $\mathbf{k} \leftarrow 0_{2 \times 2}$;
/* Initialization */
7 foreach element e do
8 $\mathbf{n}_e \leftarrow$ vector $(3 \times 1)$ of node indices of element e;
9 $U_e \leftarrow [U_1(\mathbf{n}_e), U_2(\mathbf{n}_e)]_{3 \times 2}$;
10 $X_e \leftarrow$ coordinates of the 3 nodes $\mathbf{n}_e$;
11 $N_1(\mathbf{x}), N_2(\mathbf{x}), N_3(\mathbf{x}) \leftarrow X_e$;
12 $\mathbf{D}_x = \begin{bmatrix} \frac{\partial N_1(\mathbf{x})}{\partial x} & \frac{\partial N_2(\mathbf{x})}{\partial x} & \frac{\partial N_3(\mathbf{x})}{\partial x} \\ \frac{\partial N_1(\mathbf{x})}{\partial y} & \frac{\partial N_2(\mathbf{x})}{\partial y} & \frac{\partial N_3(\mathbf{x})}{\partial y} \end{bmatrix}$;
13 $\mathbf{k} \leftarrow \mathbf{k} + U_e^T \mathbf{D}_x^T \mathbf{k}_e \mathbf{D}_x U_e A_e$;
/* A_e, area of element */
14 $A(\mathrm{RVE}) \leftarrow$ area of the RVE;
15 $\langle \mathbf{k} \rangle \leftarrow \mathbf{k}/A(\mathrm{RVE})$;
16 return $\langle \mathbf{k} \rangle$

### 3.3. Deterministic homogenization with synthetic data

The orientation sampling is performed using the acceptance-rejection algorithm described in Algorithm 1.

A convergence study carried out to determine an appropriate RVE size (c.f. A.2) retained an RVE containing 2400 grains, meshed with T3 finite elements, with approximately 270 elements per grain, resulting in a total of 667,635 elements. Under this configuration, the numerical error is estimated to be within 1–2%. The microstructure was generated using the open-source software Neper [20], and the mesh was created with the integrated Gmsh module [37]. The resulting RVE geometry and mesh is shown in Fig. 8.

This RVE configuration is employed for subsequent homogenization tests, with different Pole Density Function reconstructions affecting solely the orientation of the grains.

Following Kanit et al. [18], the effective conductivities computed from the Observed PF are reported in Table 2, and denoted by $\langle \mathbf{k} \rangle_{\mathrm{obs}}$ in the following sections.
<figure><img src="image_10.png" /><figcaption>Fig. 8. Meshed 2400 grains 2D RVE with 667,635 T3 elements.</figcaption></figure>

<caption>Table 2</caption>

Effective conductivity of deterministic homogenization with the Observed PF.

<table><thead><tr><th>(k)obs (W/(m·K))</th><th>⟨k<sub>XX</sub>⟩</th><th>⟨k<sub>YY</sub>⟩</th><th>⟨k<sub>XY</sub>⟩</th></tr></thead><tbody><tr><td>Observed PF</td><td>2.658</td><td>15.766</td><td>0.039</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Fig. 9. GP-based reconstructions of $M(P_h)$ and $LB/UB_{0.95}(P_h)$ for $\alpha = 0.05$, upper view of the sphere.</figcaption></figure>

The results exhibit a pronounced anisotropy, as evidenced by the significant difference between $\langle k_{YY} \rangle$ and $\langle k_{XX} \rangle$. Under present assumptions, this outcome is consistent with the combination of a highly anisotropic local conductivity (Eq. (19)) and a texture strongly localized along the Y and Z axes.

The coupling term $\langle k_{XY} \rangle$ remains close to zero, indicating minimal correlation between the X and Y directions. This outcome is consistent with the near-symmetry of the PF with respect to the X-Z and Y-Z planes.

These values serve as the reference for subsequent comparisons discretized on the same set of 1500 measured orientations.

### 3.4. Mean pole density function and CI bounds

To integrate measurement error within the GP framework, a noise hyperparameter $\sigma_n$ characterized by a noise level coefficient $\alpha$ (Eq. (10)) is introduced.

When $\alpha = 0$, the GP model is interpolant. Consequently, reconstructions discussed in the following sections, whenever $\alpha = 0$, yield the effective conductivity of Table 2. However, when $\alpha \neq 0$, the non-interplant model provides a tolerance interval, or equivalently a CI, whose width is related to $\sigma_n$. In this way, the measurement error is explicitly integrated into the model. In this section, $M(P_h)$ and $LB/UB_{0.95}(P_h)$ are considered. Fig. 9 illustrates these three cases for a noise level of $\alpha = 0.05$.

The lower bound enhances high-density regions, producing a more localized texture, whereas the upper bound is flatter and closer to a uniform distribution. This difference stems from the CI construction (Eq. (11)): the posterior deviation is subtracted or added to the posterior mean, and the result is renormalized to maintain the sum of 1. Subtracting the deviation drives low-density regions closer to zero, amplifying the relative intensity of high-density zones after normalization. Conversely, adding the deviation elevates all densities, including zero or near-zero ones, reduces overall contrast and promoting a more uniform texture after the normalization.
<table><caption>**Table 3**
Effective $\langle k_{XX} \rangle$ and $\langle k_{YY} \rangle$ for $M(P_h)$ and $LB/UB_{0.95}(P_h)$ under varying measurement error levels. Represented by the relative difference to $\langle k \rangle_{\text{obs}}$.</caption><thead><tr><th rowspan="2">$\langle k_{XX} \rangle$<br>($W/(m \cdot K)$)</th><th colspan="6">$\alpha$</th></tr><tr><th>0% ($\langle k \rangle_{\text{obs}}$)</th><th>1%</th><th>2%</th><th>3%</th><th>4%</th><th>5%</th></tr></thead><tbody><tr><td>$M(P_h)$</td><td>2.658</td><td>+0.6%</td><td>+0.3%</td><td>+0.8%</td><td>+0.5%</td><td>+1.1%</td></tr><tr><td>$LB_{0.95}(P_h)$</td><td>/</td><td>-0.7%</td><td>-2.4%</td><td>-1.7%</td><td>-2.1%</td><td>-1.8%</td></tr><tr><td>$UB_{0.95}(P_h)$</td><td>/</td><td>+13.0%</td><td>+21.5%</td><td>+28.0%</td><td>+35.0%</td><td>+38.5%</td></tr><tr><th rowspan="2">$\langle k_{YY} \rangle$<br>($W/(m \cdot K)$)</th><th colspan="6">$\alpha$</th></tr><tr><th>0% ($\langle k \rangle_{\text{obs}}$)</th><th>1%</th><th>2%</th><th>3%</th><th>4%</th><th>5%</th></tr><tr><td>$M(P_h)$</td><td>15.766</td><td>-1.2%</td><td>+0.7%</td><td>-0.8%</td><td>-1.4%</td><td>-0.5%</td></tr><tr><td>$LB_{0.95}(P_h)$</td><td>/</td><td>+0.5%</td><td>+1.7%</td><td>+1.3%</td><td>+1.7%</td><td>+3.1%</td></tr><tr><td>$UB_{0.95}(P_h)$</td><td>/</td><td>-7.4%</td><td>-10.1%</td><td>-13.7%</td><td>-15.4%</td><td>-16.4%</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Evolution of the relative difference for $\langle k_{XX} \rangle$ and $\langle k_{YY} \rangle$ obtained from $M(P_h)$ and from $LB/UB_{0.95}(P_h)$, compared with $\langle k \rangle_{\text{obs}}$. Noise levels vary from $\alpha = 0$ to 0.05.</figcaption></figure>

<table><caption>**Table 4**
Effective $\langle k_{YY} \rangle$ in exact value for $M(P_h)$ and $LB/UB_{0.95}(P_h)$ under varying measurement error levels.</caption><thead><tr><th rowspan="2">$\langle k_{YY} \rangle$<br>($W/(m \cdot K)$)</th><th colspan="6">$\alpha$</th></tr><tr><th>0% ($\langle k \rangle_{\text{obs}}$)</th><th>1%</th><th>2%</th><th>3%</th><th>4%</th><th>5%</th></tr></thead><tbody><tr><td>$M(P_h)$</td><td>0.039</td><td>0.056</td><td>0.023</td><td>-0.035</td><td>0.003</td><td>0.002</td></tr><tr><td>$LB_{0.95}(P_h)$</td><td>/</td><td>0.064</td><td>-0.065</td><td>0.029</td><td>-0.022</td><td>0.039</td></tr><tr><td>$UB_{0.95}(P_h)$</td><td>/</td><td>0.064</td><td>0.069</td><td>0.015</td><td>0.023</td><td>-0.013</td></tr></tbody></table>

The effective values of $\langle k_{XX} \rangle$ and $\langle k_{YY} \rangle$ obtained from $M(P_h)$ and $LB/UB_{0.95}(P_h)$ are listed in Table 3 and visualized in Fig. 10 as relative difference with respect to $\langle k \rangle_{\text{obs}}$ components

$$ \epsilon = \frac{\langle k \rangle_{\text{test}} - \langle k \rangle_{\text{obs}}}{\langle k \rangle_{\text{obs}}} $$

When $\alpha = 0$, $M(P_h)$ and both $LB/UB_{0.95}(P_h)$ coincide with the Observed PF, and thus the effective conductivities reported in Table 3 are identical to those presented in Table 2. As the error level increases, $M(P_h)$-based results show negligible variation ($\le 1.4\%$), demonstrating strong robustness of the reconstruction. For $LB_{0.95}(P_h)$, $\langle k_{XX} \rangle$ decreases slightly and stabilizes beyond $\alpha = 0.02$, while $\langle k_{YY} \rangle$ increases monotonically with $\alpha$, indicating enhanced anisotropy. This trend is consistent with the corresponding PF at $\alpha = 0.05$ (Fig. 9a), where the X-direction exhibits mostly zero-density regions that cannot decrease further, leading to the saturation of $\langle k_{XX} \rangle$. In contrast, the Y-direction becomes increasingly concentrated, causing $\langle k_{YY} \rangle$ to rise. However, the variation remains moderate, with the maximum difference of $3.1\% < 5\%$ obtained for $\langle k_{YY} \rangle$ at $\alpha = 0.05$.

For $UB_{0.95}(P_h)$, the opposite trend is observed: $\langle k_{XX} \rangle$ increases while $\langle k_{YY} \rangle$ decreases with $\alpha$, indicating a progressive tendency toward isotropy. Meanwhile, the magnitude of variation is substantial: at $\alpha = 0.05$, deviations from the $\langle k \rangle_{\text{obs}}$ reach 37% for $\langle k_{XX} \rangle$ and 16% for $\langle k_{YY} \rangle$. This behavior is consistent with the uniformization behavior observed in the corresponding PF (Fig. 9c).

This can be interpreted in two ways. First, the underlying texture and the strong intrinsic anisotropy of the chosen local conductivity amplify the contrast between directions, allowing for large variations in effective properties towards the isotropy. In contrast, the $LB_{0.95}(P_h)$ is constrained by the non-negativity imposed through the log-linear transformation.

Second, as discussed in Section 2.2, although $LB/UB(P_h)$ satisfy the mathematical conditions of a Pole Density Function, they are not statistically admissible instances of the GP. They serve only as illustrative extremes of the possible variability rather than physically realizable textures. This is confirmed by the observation of the evolution.
<figure><img src="image_13.png" /><figcaption>Fig. 11. $S(P_h)$ samples for $\alpha = 0.05$, upper view of the sphere.</figcaption></figure>

<caption>**Table 5**</caption>

Effective values and corresponding 95% CI bounds $LB/UB_{0.95}(\langle k \rangle)$ of $\langle k_{XX} \rangle$ and $\langle k_{YY} \rangle$ based on log-normal distribution models fitted to the results of 600 $S(P_h)$-based homogenizations at each measurement error level $\alpha$. Exact values are reported for the modes, together with their relative differences (%) with respect to $\langle k \rangle_{\text{obs}}$. Relative differences (%) with respect to the mode are provided for $LB/UB_{0.95}(\langle k \rangle)$.

<table><thead><tr><th>⟨k<sub>XX</sub>⟩</th><th colspan="5">α</th></tr><tr><th>(W/(m·K))</th><th>0% (⟨k⟩<sub>obs</sub>)</th><th>1%</th><th>2%</th><th>3%</th><th>4%</th><th>5%</th></tr></thead><tbody><tr><td>Mode</td><td>2.658 (+0%)</td><td>2.729 (+2.7%)</td><td>2.782 (+4.7%)</td><td>2.839 (+6.8%)</td><td>2.887 (+8.6%)</td><td>2.916 (+9.7%)</td></tr><tr><td>UB<sub>0.95</sub>(⟨k⟩)</td><td>/</td><td>+1.45%</td><td>+1.86%</td><td>+2.86%</td><td>+2.86%</td><td>+2.00%</td></tr><tr><td>LB<sub>0.95</sub>(⟨k⟩)</td><td>/</td><td>-1.42%</td><td>-1.81%</td><td>-2.74%</td><td>-2.74%</td><td>-1.94%</td></tr><tr><th>⟨k<sub>YY</sub>⟩</th><th colspan="5">α</th></tr><tr><th>(W/(m·K))</th><th>0% (⟨k⟩<sub>obs</sub>)</th><th>1%</th><th>2%</th><th>3%</th><th>4%</th><th>5%</th></tr><tr><td>Mode</td><td>15.766 (+0%)</td><td>15.464 (-1.9%)</td><td>15.340 (-2.7%)</td><td>15.151 (-3.9%)</td><td>15.001 (-4.9%)</td><td>14.854 (-5.8%)</td></tr><tr><td>UB<sub>0.95</sub>(⟨k⟩)</td><td>/</td><td>+1.75%</td><td>+1.66%</td><td>+3.07%</td><td>+3.08%</td><td>+2.62%</td></tr><tr><td>LB<sub>0.95</sub>(⟨k⟩)</td><td>/</td><td>-1.71%</td><td>-1.62%</td><td>-2.93%</td><td>-2.94%</td><td>-2.52%</td></tr></tbody></table>

Finally, the effective values of $\langle k_{XY} \rangle$ reported in Table 4 remain consistently close to zero, with a maximum magnitude of 0.069, indicating a negligible coupling effect. Moreover, no systematic variation is observed with increasing measurement error. This behavior suggests that the symmetry characteristics of the texture are not affected by the introduced noise, and minor fluctuations can be attributed primarily to numerical error.

### 3.5. Sampled pole density functions

As discussed in Section 2.2, $M(P_h)$ and $LB/UB(P_h)$ are both deterministic reconstructions and are not consistent necessarily with the covariance structure defined by the covariance function. In contrast, the $S(P_h)$ random instances drawn from the GP follow the covariance matrix determined by the optimized covariance function (Eq. (15)). They thus represent statistically and physically admissible PFs incorporating measurement uncertainty. Examples of $S(P_h)$ for $\alpha = 0.05$ are presented in Fig. 11. Compared with $LB/UB(P_h)$ (Fig. 9a, c), these realizations exhibit moderate yet physically meaningful deviations to $M(P_h)$ (Fig. 9b), capturing the intrinsic variability implied by the GP model.

In the analysis of $S(P_h)$, $m = 600$ instances are sampled from the GP model at each measurement-error level and used for homogenization, yielding 600 corresponding sets of effective conductivities. $m$, the minimum number of $S(P_h)$ samples required for statistical convergence of the effective properties is determined using the Monte Carlo convergence algorithm described in Algorithm 4. The resulting distributions of the homogenized properties, which reflect the variability introduced by the random GP instances, are then characterized statistically. For $\langle k_{XX} \rangle$ and $\langle k_{YY} \rangle$, which are strictly non-negative, log-normal distributions are adopted, and their modes are taken as the most statistically representative values. For $\langle k_{XY} \rangle$, a normal distribution is fitted, of which the mode coincides with the mean. Based on these fitted continuous models, the 95% CI bounds of the distributions, denoted by $LB/UB_{0.95}(\langle k \rangle)$, are subsequently determined.

The representative values and corresponding bounds are summarized in Tables 5 and 6. In Table 5, the exact values of the modes are reported, together with their relative differences with respect to $\langle k \rangle_{\text{obs}}$. The CI bounds of the fitted distributions are expressed as relative differences with respect to the mode. For $\langle k_{XY} \rangle$, all values are given in absolute form. Because $\langle k_{XX} \rangle$ and $\langle k_{YY} \rangle$ are modeled using log-normal distributions, their CI bounds are not symmetric with respect to the mode. In contrast, for $\langle k_{XY} \rangle$, which follows a normal distribution, the CI bounds are symmetric about the mode (or equivalently, the mean).

Fig. 12 illustrates the evolution of the mode-here considered as the representative effective conductivity-along with its corresponding 95% CI. The evolution follows a trend similar to that of the GP upper bound (Fig. 10), exhibiting a gradual tendency toward isotropy. This behavior arises from the log-linear data transformation, which introduces asymmetry in the GP's CI near zero density
<caption>**Table 6**</caption>

Effective values and corresponding 95% CI bounds $LB/UB_{0.95}(\langle k \rangle)$ of $\langle k_{XY} \rangle$ based on normal distribution models fitted to the results of 600 $S(P_h)$-based homogenizations at each measurement error level $\alpha$.

<table><thead><tr><th>⟨k<sub>XY</sub>⟩<br>(W/(m·K))</th><th colspan="6">α</th></tr><tr><th></th><th>0% (⟨k⟩<sub>obs</sub>)</th><th>1%</th><th>2%</th><th>3%</th><th>4%</th><th>5%</th></tr></thead><tbody><tr><th>Mode</th><td>0.039</td><td>0.023</td><td>0.010</td><td>0.041</td><td>0.042</td><td>0.015</td></tr><tr><th>LB/UB<sub>0.95</sub>(⟨k⟩)</th><td>/</td><td>±0.068</td><td>±0.128</td><td>±0.157</td><td>±0.197</td><td>±0.187</td></tr></tbody></table>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Evolution of the mode and the 95% CI of the $S(P_h)$-based effective conductivities, obtained from 600 random $S(P_h)$ at each noise level $\alpha$, varying from 0 to 0.05.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Reconstructed log-linear distribution models for $\langle k_{XX} \rangle$ and $\langle k_{YY} \rangle$ based on $S(P_h)$. Noise levels vary from $\alpha = 0$ to 0.05.</figcaption></figure>

values. For measured points where the measured density is initially zero (orientation not detected in the material), only positive densities can be sampled, leading to a strictly positive mean at these locations. Consequently, artificial density appears in zero or near-zero regions, while the overall normalization redistributes density across the sphere. This results in a more uniform texture, explaining the observed isotropic tendency of the $S(P_h)$-based effective conductivities as the measurement error increases.

For $\langle k_{XY} \rangle$, the mode remains nearly constant as $\alpha$ increases, with a total variation of only 0.032 between $\alpha = 0.02$ and $\alpha = 0.04$. However, the width of the confidence interval generally increases with the measurement error, reaching a non-negligible maximum of $\pm 0.197$. As discussed in Section 3.4, the values of $\langle k_{XY} \rangle$ close to zero primarily result from the near-symmetric character of the PF with respect to the X and Y axes. The observed increase in dispersion may be explained by locally asymmetric distributions of individual $S(P_h)$, as illustrated in Fig. 11. Nevertheless, when averaged over multiple $S(P_h)$ realizations, the statistical mean remains stable and consistent with the overall symmetry of the Pole Density Function.

Considering the evolution of the CI of the homogenization results, the fitted log-normal distributions for $\langle k_{XX} \rangle$ and $\langle k_{YY} \rangle$ are visualized in Fig. 13. It can be observed that the width of the 95% CI does not increase monotonically with the measurement error level. Although the CI is narrower for $\alpha = 0.01$ and 0.02 compared with higher error levels, this tendency does not occur at $\alpha = 0.05$ for both $\langle k_{XX} \rangle$ and $\langle k_{YY} \rangle$. Considering this observation together with the evolution of the mode values, it can be inferred that, as the
<figure><img src="image_16.png" /><figcaption>Fig. 14. Comparison of $\langle k \rangle$ obtained using $M(P_h)$, $LB/UB_{0.95}(P_h)$, and the mode $\langle k \rangle$ from multiple $S(P_h)$. Noise levels vary from $\alpha = 0$ to 0.05.</figcaption></figure>

measurement error increases, the $S(P_h)$ exhibit a consistent shift toward isotropy across all instances, leading to similar directional changes in the mode. However, the dispersion among $S(P_h)$ at a given noise level does not show a clear dependence on the global measurement error.

### 3.6. Concluding comparison of effective properties for different reconstructions

Fig. 14 compares the effective properties obtained from all homogenization strategies for measurement error levels ranging from $\alpha = 0$ to 0.05. The reference conductivity $\langle k \rangle_{\text{obs}}$ is shown in green. The deterministic results derived from $M(P_h)$ and from $LB/UB_{0.95}(P_h)$ are displayed in blue using continuous and dotted/dashed lines, respectively. The stochastic approach based on $S(P_h)$ is shown in red, where the continuous line connects the mode of the realizations at each error level and the associated 95% CI are represented as error bars.

For each conductivity component, all curves converge at $\alpha = 0$, corresponding to $\langle k \rangle_{\text{obs}}$. As the measurement error increases, the effective conductivities obtained from the posterior mean $M(P_h)$ remain stable and close to $\langle k \rangle_{\text{obs}}$, indicating a strong robustness to uncertainty. In contrast, the GP's CI bounds exhibit significant variation, particularly $UB_{0.95}(P_h)$, which tends to promote isotropy. The physically admissible realizations $S(P_h)$ follow the same global tendency as $UB_{0.95}(P_h)$, though with a more moderate amplitude.

The coupling term $\langle k_{XY} \rangle$ remains close to zero in all cases. This arises from the near-symmetry of the example PF with respect to the X and Y axes, which leads to statistical cancellation of the XY coupling effects. All reconstruction types preserve this symmetry, and the small variations observed are likely due to numerical dispersion.

This comparison highlights the contrast between the two homogenization strategies. The deterministic approach produces either stable results or extreme variations, and therefore does not provide a realistic quantification of uncertainty propagation. As a consequence, its reference value remains limited. In contrast, the stochastic approach captures a clear evolution of the effective properties with respect to the measurement error level and maintains a variability range that remains consistent with the magnitude of the error.

## 4. Conclusion

The homogenization based on GP posterior mean of Pole Density Function yields stable effective conductivities across measurement error levels, closely matching the results obtained directly from the synthetic data. In contrast, the confidence interval bounds, leading to extreme uniform or localized textures, produce increasing variations in effective properties, excessively high relatively to the measurement error.

In contrast, Monte Carlo sampling of admissible Pole Density Functions from the GP posterior distribution leads to stable variability of the homogenized conductivities unaffected by measurement uncertainty levels.

These conclusions may be explained by the non-linearity of the microstructural texture-property relationship. The full stochastic approach, based on the random sampling from the posterior distribution appears therefore as the most relevant considering measurement uncertainty propagation.

The present work, based on *in silico* experiments, needs to be further extended to accommodate *in vitro* data. Future work considers the extension of current framework to non-linear mechanical response by prescribing crystal plasticity laws on slip systems at grain level.

### CRediT authorship contribution statement

**Bingqian Li:** Writing – original draft, Methodology, Investigation, Data curation; **Piotr Breitkopf:** Writing – review & editing, Validation, Supervision, Project administration, Methodology, Investigation, Conceptualization; **Ludovic Cauvin:** Writing – review & editing, Validation, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization.
## Data availability

Data will be made available on request.

## Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Bingqian Li reports financial support was provided by Hauts-de-France Region. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgement

The authors gratefully acknowledge the financial support of the Région Hauts-de-France, France and the SAFE AI Chair of the Heudiasyc Laboratory (UMR-CNRS 7253), France, under the co-funding of the Ph.D. project.

## Appendix A. RVE construction

### A.1. Neper commands for microstructure and mesh generation

The generation of 2D microstructure for the RVE construction is proceeded by the open-source software Neper [20]. The command to generate a 2D microstructure with chosen grain number *n* is given as:

```
neper -T -n <grain number> -reg 1 -dim 2 -o <name.tess>
```

The microstructure will be recorded in a file ".tess". If one want to impose a uniform size to all the grains, the size distribution can be chosen as a Dirac distribution by adding:

```
-morpho 'diameq:dirac(1)'
```

The T3 mesh is generated using integrated Gmsh module [37] by the following command with the parameter "-rcl" controlling the mesh refinement and tested thoroughly in the following sections:

```
neper -M -dim 2 <name.tess> -rcl <value> -elttype tri
```

Table A.7 provides the correspondence between -rcl values and the average number of elements per grain ($N_e/n$).

<caption>Table A.7</caption>

Number of elements per grain w.r.t. the -rcl setting for T3 elements.

<table><thead><tr><th></th><th>-rcl</th><th>4.5</th><th>3.5</th><th>3.0</th><th>2.5</th><th>2.0</th><th>1.8</th><th>1.5</th><th>1.0</th><th>0.9</th></tr></thead><tbody><tr><th>N<sub>e</sub>/n</th><td>4.4</td><td>7.3</td><td>9.1</td><td>11.4</td><td>15.9</td><td>18.2</td><td>24.6</td><td>53.7</td><td>62.8</td></tr><tr><th>-rcl</th><td>0.8</td><td>0.7</td><td>0.6</td><td>0.5</td><td>0.4</td><td>0.3</td><td>0.2</td><td>0.15</td><td>0.12</td></tr><tr><th>N<sub>e</sub>/n</th><td>79.5</td><td>96.6</td><td>133.8</td><td>192.1</td><td>273.4</td><td>473.9</td><td>1032.0</td><td>1811.3</td><td>2808.3</td></tr><tr><th>-rcl</th><td>0.11</td><td>0.10</td><td>0.09</td><td>0.085</td><td>0.08</td><td>0.075</td><td>0.07</td><td></td><td></td></tr><tr><th>N<sub>e</sub>/n</th><td>3365.3</td><td>4068.7</td><td>4985.3</td><td>5613.1</td><td>6359.8</td><td>7211.8</td><td>8261.3</td><td></td><td></td></tr></tbody></table>

### A.2. RVE mesh refinement and grain count choice

The finite element convergence study is carried out on RVEs with $n = 10, 50, 100$, and $1600$ grains (Fig. A.15). Grain orientations are initially randomly assigned and remain unchanged during the refinement of the element size.

<figure><img src="image_17.png" /><figcaption>Fig. A.15. Voronoi tessellations with 10, 50, 100, and 1600 grains (Neper).</figcaption></figure>

The local conductivity tensor is set to as it is in Eq. (19). In the absence of an analytical solution for effective macroscopic conductivities, the relative error in a homogenized conductivity $\langle k \rangle$ is evaluated with respect to a reference value $\langle k \rangle^*$ obtained using
<figure><img src="image_18.png" /><figcaption>Fig. A.16. Convergence of the homogenized conductivity error with respect to $N_e/n$, for RVEs with 10, 50, 100, and 1600 grains.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. A.17. Evolution of $\langle k_{XX} \rangle$, $\langle k_{YY} \rangle$ and $\langle k_{XY} \rangle$ w.r.t. $n$.</figcaption></figure>

the finest mesh:

$$ \epsilon(\langle \mathbf{k} \rangle) = \frac{\|\langle \mathbf{k} \rangle - \langle \mathbf{k} \rangle^*\|_F}{\|\langle \mathbf{k} \rangle^*\|_F} \quad (A.1) $$

where $\| \cdot \|_F$ denotes the Frobenius norm.

Fig. A.16 shows the convergence curves of homogenized conductivities with respect to $N_e/n$.

The results indicate that when $N_e/n \approx 270$, corresponding to $-rcl = 0.4$, the relative error falls below 1% for all cases. This mesh density is adapted for all subsequent simulations.

Fig. A.17 shows the evolution of the homogenized conductivity components with respect to the number of grains. The dots represent the median values across all realizations for each grain count, with dotted lines connecting the medians to indicate overall trends. The error bars of the 95% CI for $\langle k_{XX} \rangle$ and $\langle k_{YY} \rangle$ are computed based on the log-normal distribution, explaining the asymmetry of the confidence intervals.

$n = 2400$ is selected as a compromise between computational efficiency and numerical precision for subsequent simulations with the relative error $\epsilon$ lies between 1–2%.

## References

[1] J. Segurado, R.A. Lebensohn, J. Llorca, Computational homogenization of polycrystals, Advances in Applied Mechanics 51, Elsevier, 2018, 1–114. https://doi.org/10.1016/bs.aams.2018.07.001

[2] J. Jin, L. Cauvin, B. Raghavan, P. Breitkopf, S. Dutta, M. Xiao, Towards a data-driven paradigm for characterizing plastic anisotropy using principal components analysis and manifold learning, Comput. Mater. Sci. 235 (2024) 112834. https://doi.org/10.1016/j.commatsci.2024.112834

[3] M. Kamaya, Assessment of local deformation using EBSD: quantification of accuracy of measurement and definition of local gradient, Ultramicroscopy 111 (8) (2011) 1189–1199. https://doi.org/10.1016/j.ultramic.2011.02.004

[4] G. Nolze, M. Jürgens, J. Olbricht, A. Winkelmann, Improving the precision of orientation measurements from technical materials via EBSD pattern matching, Acta Mater. 159 (2018) 408–415. https://doi.org/10.1016/j.actamat.2018.08.028

[5] L. Tóth, P. Houtte, Discretization techniques for orientation distribution functions, Textures Microstruct. 19 (1992). https://doi.org/10.1155/TSM.19.229
[6] N.R. Barton, D.E. Boyce, P.R. Dawson, Pole figure inversion using finite elements over rodrigues space Texture Stress Microstruct. 35(2) (2002) 372523, https://doi.org/10.1080/073033002100000182

[7] H.J. Bunge, Texture Analysis in Materials Science: Mathematical Methods, Butterworths, London Boston, english ed edition, Butterworths, London Boston, 1982.

[8] P.J. Madrid, D. Sulsky, R.A. Lebensohn, Uncertainty quantification in prediction of the in-plane young's modulus of thin films with fiber texture, J. Microelectromech. Syst. 23 (2) (2014) 380-390. https://doi.org/10.1109/JMEMS.2013.2279500

[9] L. Cauvin, B. Raghavan, J. Jin, S. Bouvier, X. Wang, F. Meraghni, J.L. Alves, Investigating the plastic anisotropy and hardening behavior of a commercial Zn-Cu-Ti alloy: experimental & modeling approach, Mech. Mater. 164 (2022) 104103. https://doi.org/10.1016/j.mechmat.2021.104103

[10] B. Panicaud, K. Saanouni, A. Baczmański, M. François, L. Cauvin, L.L. Joncour, Theoretical modelling of ductile damage in duplex stainless steels - Comparison between two micro-mechanical elasto-plastic approaches, Comput. Mater. Sci. 50 (6) (2011) 1908-1916. https://doi.org/10.1016/j.commatsci.2011.01.038

[11] G. Yuan, X. Li, Z. Dong, A. Westwood, Z. Cui, Y. Cong, H. Du, F. Kang, Graphite blocks with preferred orientation and high thermal conductivity, Carbon 50 (1) (2012) 175-182. https://doi.org/10.1016/j.carbon.2011.08.017

[12] C. Bronkhorst, S. Kalidindi, L. Anand, Polycrystalline plasticity and the evolution of crystallographic texture in FCC metals, Philos. Trans. R. Soc. London Ser. Math. Phys. Eng. Sci. 341 (1992) 443-477. https://doi.org/10.1098/rsta.1992.0111

[13] Y. Yang, J.-P. Wang, Enhancing unconventional spin-orbit torque efficiency: numerical study on the influence of crystallographic texture and polycrystalline effects on low-Symmetry materials, 2025, IEEE Trans. Magn., 1-1. https://doi.org/10.1109/TMAG.2025.3537602

[14] R.A. Lebensohn, C.N. Tomé, A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: application to zirconium alloys, Acta Metall. Mater. 41 (9) (1993) 2611-2624. https://doi.org/10.1016/0956-7151(93)90130-K

[15] L. Cauvin, B. Raghavan, S. Bouvier, X. Wang, F. Meraghni, Multi-scale investigation of highly anisotropic zinc alloys using crystal plasticity and inverse analysis, Mater. Eng. A 729 (2018) 106-118. https://doi.org/10.1016/j.msea.2018.05.038

[16] J. Yvonnet, Computational Homogenization of Heterogeneous Materials with Finite Elements, Vol. 258 of Solid Mechanics and Its Applications, Springer International Publishing, Cham, 2019. https://doi.org/10.1007/978-3-030-18383-7

[17] C. Huet, Application of variational concepts to size effects in elastic heterogeneous bodies, J. Mech. Phys. Solids 38 (6) (1990) 813-841. https://doi.org/10.1016/0022-5096(90)90041-2

[18] T. Kanit, S. Forest, I. Galliet, V. Mounoury, D. Jeulin, Determination of the size of the representative volume element for random composites: statistical and numerical approach, Int. J. Solids Struct. 40 (13) (2003) 3647-3679. https://doi.org/10.1016/S0020-7683(03)00143-4

[19] D. Jeulin, T. Kanit, S. Forest, Representative Volume Element: A Statistical Point of View, Continuum Models and Discrete Systems, Springer, 2004. Dordrecht, 21-27 https://doi.org/10.1007/978-1-4020-2316-3_5

[20] R. Quey, P. Dawson, F. Barbe, Large-scale 3D random polycrystals for the finite element method: generation, meshing and remeshing, Comput. Methods Appl. Mech. Eng. 200, (17-20) publisher: Elsevier, 2011. https://doi.org/10.1016/j.cma.2011.01.002

[21] B. Chen, D. Li, P. Davies, R. Johnston, X. Ge, C. Li, Recent Progress of Digital Reconstruction in Polycrystalline Materials, Archives of Computational Methods in Engineering, 2025. https://doi.org/10.1007/s11831-025-10245-4

[22] F. Bachmann, R. Hielscher, H. Schaeben, Texture Analysis with MTEX-Free and Open Source Software Toolbox 160 (2010). https://doi.org/10.4028/WWW.scientific.net/SSP.160.63

[23] M. Knezevic, N.W. Landry, Procedures for reducing large datasets of crystal orientations using generalized spherical harmonics, Mech. Mater. 88 (2015) 73-86. https://doi.org/10.1016/j.mechmat.2015.04.014

[24] H.J. Bunge, Mathematische Methoden der Texturanalyse, Walter de Gruyter GmbH & Co KG, 1969, google-Books-ID: pUwzEQAAQBAJ.

[25] A. Creuziger, K. Syed, T. Gnäupel-Herold, Measurement of uncertainty in orientation distribution function calculations, Scr. Mater. 72-73 (73) (2014) 55-58. https://doi.org/10.1016/j scriptamat.2013.10.017

[26] P. Acar, V. Sundararaghavan, Uncertainty quantification of microstructural properties due to variability in measured pole figures, Acta Mater. 124 (2017) 100-108. https://doi.org/10.1016/j.actamat.2016.10.070

[27] C.E. Rasmussen, C.K.I. Williams, Gaussian Processes for Machine Learning, Adaptive Computation and Machine Learning, oCLC: ocm61285753, MIT Press, Cambridge, Mass, Cambridge, Mass, 2006, oCLC: ocm61285753.

[28] A. Madra, P. Causse, F. Trochu, J. Adrien, E. Maire, P. Breitkopf, Stochastic Characterization of Textile Reinforcements in Composites Based on X-Ray Microtomographic Scans, 224, Elsevier BV, 2019. https://doi.org/10.1016/j.compstruct.2019.111031

[29] R. Saunders, C. Butler, J. Michopoulos, D. Lagoudas, A. Emwany, A. Bagchi, Mechanical behavior predictions of additively manufactured microstructures using functional Gaussian process surrogates, Comput. Mater. 7 (2021) 81. https://doi.org/10.1038/s41524-021-00548-y

[30] A. Venkatraman, D.M. D.O. Zapiain, H. Lim, S.R. Kalidindi, Texture-sensitive prediction of micro-spring performance using Gaussian process models calibrated to finite element simulations, Mater. Des. 197 (2021) 109198. https://doi.org/10.1016/j.matdes.2020.109198

[31] H. Dorward, M. Mostafavi, D.M. Knowles, M.J. Peel, Reduced-order representations of crystallographic texture for application to surrogate modelling of austenitic stainless steel, J. Mech. Phys. Solids 208 (2026) 106444. https://doi.org/10.1016/j.jmps.2025.106444

[32] B. Li, P. Breitkopf, L. Cauvin, Customized Gaussian process for representing polycrystalline texture, Comput. Methods Appl. Mech. Eng. 440 (2025) 117934. https://doi.org/10.1016/j.cma.2025.117934

[33] M. Lè, J. Unkelbach, N. Ayache, H. Delingette, Sampling Image Segmentations for Uncertainty Quantification, 34, Elsevier BV, 2016. https://doi.org/10.1016/j.media.2016.04.005

[34] P. Chen, J. Guilleminot, C. Soize, Concurrent multiscale simulations of nonlinear random materials using probabilistic learning, Comput. Methods Appl. Mech. Eng. 422 (2024) 116837. https://doi.org/10.1016/j.cma.2024.116837

[35] J.L. W.V. Jensen, Sur les fonctions convexes et les inégalités entre les valeurs moyennes, Acta Math. 30 (1) (1906) 175-193. https://doi.org/10.1007/BF02418571

[36] M. Raissi, G.E. Karniadakis, Hidden physics models: machine learning of nonlinear partial differential equations, J. Comput. Phys. 357 (2018) 125-141. https://doi.org/10.1016/j.jcp.2017.11.039

[37] C. Geuzaine, J.-F. Remacle, GMSH: a 3-D finite element mesh generator with built-in pre- and post-processing facilities, Int. J. Numer. Methods Eng. (11), 1309-1331, 2009. https://doi.org/10.1002/nme.2579