# Inferring the probability distribution over strain tensors in polycrystals from diffraction based measurements

Axel Henningsson, Adrian G. Wills, Stephen A. Hall, Johannes Hendriks, Jonathan P. Wright, Thomas B. Schön, Henning F. Poulsen

► To cite this version:

Axel Henningsson, Adrian G. Wills, Stephen A. Hall, Johannes Hendriks, Jonathan P. Wright, et al. Inferring the probability distribution over strain tensors in polycrystals from diffraction based measurements. Computer Methods in Applied Mechanics and Engineering, 2023, 417, pp.116417. 10.1016/j.cma.2023.116417. hal-05080031

HAL Id: hal-05080031

https://hal.science/hal-05080031v1

Submitted on 22 May 2025

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 Attribution 4.0 International License
# Inferring the probability distribution over strain tensors in polycrystals from diffraction based measurements☆

Axel Henningsson a,*, Adrian G. Wills b, Stephen A. Hall a, Johannes Hendriks c,d, Jonathan P. Wright c,d, Thomas B. Schön e, Henning F. Poulsen f

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

Check for
updates

## GRAPHICAL ABSTRACT

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

## ARTICLE INFO

**Keywords:**

Diffraction

Estimation

Strain tensor

X-rays

## ABSTRACT

Polycrystals illuminated by high-energy X-rays or neutrons produce diffraction patterns in which the measured diffraction peaks encode the individual single crystal strain states. While state of the art X-ray and neutron diffraction approaches can be used to routinely recover per grain mean strain tensors, less work has been produced on the recovery of higher order statistics of the strain

<footnote>a Lund University, Division of Solid Mechanics, Ole Römersväg 1, Lund, 221 00, Sweden</footnote>

<footnote>b School of Engineering, University of Newcastle, Newcastle, Australia</footnote>

<footnote>c Australian National University, Canberra, Australia</footnote>

<footnote>d ESRF-The European Synchrotron, Grenoble, France</footnote>

<footnote>e Department of Information Technology, Uppsala University, Uppsala, Sweden</footnote>

<footnote>f Department of Physics, Technical University of Denmark, Kongens Lyngby, Denmark</footnote>

<footnote>☆ This work was funded by Vetenskapsrådet - Röntgen Ångström Cluster, project no. 2017-06719.</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: axel.hennningsson@solid.lth.se (A. Henningsson).</footnote>
Polycrystals

Probability distributions

distributions across the individual grains. In the setting of small deformations, we consider the problem of estimating the crystal elastic strain tensor probability distribution from diffraction data. For the special case of multivariate Gaussian strain tensor probability distributions, we show that while the mean of the distribution is well defined from measurements, the covariance of strain has a null-space. We show that there exist exactly 6 orthogonal perturbations to this covariance matrix under which the measured strain signal is invariant. In particular, we provide analytical parametrisations of these perturbations together with the set of possible maximum-likelihood estimates for a multivariate Gaussian fit to data. The parametric description of the null-space provides insights into the strain PDF modes that cannot be accurately estimated from the diffraction data. Understanding these modes prevents erroneous conclusions from being drawn based on the data. Beyond Gaussian strain tensor probability densities, we derive an iterative radial basis regression scheme in which the strain tensor probability density is estimated by a sparse finite basis expansion. This is made possible by showing that the operator mapping the strain tensor probability density onto the measured histograms of directional strain is linear, without approximation. The utility of the proposed algorithm is demonstrated by numerical simulations in the setting of single crystal monochromatic X-ray scattering. The proposed regression methods were found to robustly reject outliers and accurately predict the strain tensor probability distributions in the presence of Gaussian measurement noise.

# 1. Introduction

Neutrons and high energy X-rays are key tools for non-destructive testing in the materials research community (e.g., [1,2]). Today, large scale synchrotron, reactor and spallation facilities offer a wide range of material characterisation possibilities where in-situ 3D volume measurements in the bulk of dense materials can be achieved, e.g., [3,4]. Coupled with modern statistical regression tools and high performance computing, a plethora of physical quantities, embedded into first principle mathematical material models, can be characterised, e.g., [5,6]. Such estimations are important for successful calibration and validation of advanced numerical materials models. Importantly, the non-destructive nature of neutrons and X-rays enables measurements on the same test specimens to be repeated, decreasing uncertainty in the estimation and unlocking the possibility for monitoring of changes during in-situ testing, e.g., loading [7] or heating [8].

Three Dimensional X-ray Diffraction Microscopy (3DXRD) [9] and the High Energy Diffraction Microscopy (HEDM) [10] are becoming established as workhorses for the diffraction microscopy imaging community. Using hard X-rays (~10–100 keV), these microscopy techniques enable per-grain characterisation of polycrystalline aggregates featuring thousands of grains [11–14]. With a penetration depth in the range of mm for dense materials, 3DXRD/HEDM has developed into a popular tool for non destructive testing of metals, alloys, rocks and granular media. The technique relies on the measurement of 2D mono-channel digital images of spotty diffraction patterns to estimate per-grain position, volume, average strain [15], average stress [16,17] and orientation [18].

Beyond the average grain state, accessible by standard 3DXRD and in pursuit of intragranular variation, Scanning Three Dimensional X-ray Diffraction Microscopy (s3DXRD) [19] and point focused High Energy Diffraction Microscopy (pf-HEDM) [20] have developed under the 3DXRD and HEDM umbrellas. The key difference in s3DXRD is the focusing of the X-ray beam such that only sub-volumes of the crystals are illuminated at any given time. By carefully scanning the sample across the beam, bundles of 2D diffraction pattern images are acquired. From these images, tomography-like reconstruction methods can be used to estimate the grain shapes as well as the intragranular strain tensor fields [21–24]. In addition to achieving sub-grain resolution, narrowing the beam helps alleviate the issue of diffraction spot overlap due to abundant scattering [25]. Unfortunately, despite recent progress in beam brilliance and detector technology, the scanning type methods are time consuming, with state of the art experiments often requiring tens of hours for a full 3D scan of a sample, e.g., [26]. On the other hand, classical (full-field) 3DXRD/HEDM microscopy, using a wide beam, offers a 10–100 times speedup.

For polycrystal samples with a moderate number of grains (100–1000), with minimal diffraction spot overlap, efficient algorithms for mapping the individual crystal orientations exist [27]. The main drawback of full-field 3DXRD/HEDM, in relation to its scanning mode, is the loss of information on the spatial origin of diffraction within the grains. This makes reconstruction of orientation and strain maps over the grain volume very challenging, especially for highly deformed samples [28]. However, valuable intragranular information can still be attained in terms of probability density distributions over the underlying spatial fields. For example, in the work of Poulsen et al. [29] and Hansen et al. [30] the probability density of orientation states within the individual crystals where shown to be recoverable. In similar spirit [31], inspired by Behnken [32], showed how the mean strain tensor state for each intragranular orientation state can be efficiently baked into the inverse problem, enriching the capabilities of the 3DXRD/HEDM microscope. The latter of these methods recovers a lattice strain density function (LSDF), sometimes also refereed to as a strain orientation density function (SODF). As explained by Bernier and Miller [33], the variation in the micromechanical state is not fully captured by the SODF. To illustrate we may consider a near-perfect crystal, such that the orientation distribution function approaches a Dirac delta, the resulting SODF will then only provide information on the mean strain tensor in the grain. None the less, the diffracting grains can possess spatially varying elastic strain fields, even if the grain orientation is close to constant. More generally, considering an orientation spread within the grain, each individual orientation state in the grain can exist on a spatially extended sub-volume of the grain over which the strain tensor is unlikely to be constant.
Just as the grain orientation variation unlocks a new type of material parameter estimation [34], higher order statistics on the distribution of the elastic strain tensors can be used to calibrate micro-mechanical material models. Motivated by this, we decouple the problem of orientation reconstruction from that of strain and develop a regression framework capable of estimating the strain tensor probability distribution (strain PDF) inside the individual crystals within a polycrystalline sample. The aim is to recover the 6-dimensional strain PDF that describes the scalar likelihood of encountering any one particular strain tensor at any randomly (uniformly) selected point within the given grain. We emphasise that our method recovers the strain PDF in a single crystal diffraction setting.

Our method is primarily aimed at upgrading the full-field mode of the 3DXRD/HEDM microscope, where peak strain broadening is a known phenomena (c.f [35]). However, our results can be cast abstractly as a regression from sets of histograms of directional strains attached to an unknown spatial volume. Other neutron and X-ray microscopy methods where such data can reliably be made available will, therefore, be applicable for a similar upgrade. For instance [36] recently showed that the second moment of such directional strain histograms can be reliably recovered from neutron Bragg edge transmission measurements on polycrystals.

The structure of the paper is as follows. In Section 2, we describe how the considered diffraction measurements relate to the strain state of the crystals and introduce the concept of histograms of directional strain. In Section 3 we define our target reconstruction quantity – the strain PDF – and present the inverse problem we are solving. In Section 4, we discuss some of the benefits and limitations of recovering the strain PDF. In Section 5, we derive the mapping between measured and reconstructed quantities. Using the established measurement model we describe our solution to the inversion problem in Section 6. In Section 7, we discuss the observability of the inverse problem for the special case of a multivariate Gaussian strain PDF and provide a closed form solution for its estimation. To verify our mathematical results we devote Section 8 to a numerical simulation study of a spherical iron (Fe) grain with a cubic crystal structure in mechanical self-equilibrium. We summarise our findings in Section 9, and give an outlook for future research in Section 10. Finally, we conclude our paper in Section 11.

## 2. Crystal strains in diffraction measurements

The conditions for diffraction to occur from a crystal can be described by the Laue equations. Under elastic scattering conditions, letting the incident wavevector, $k$, be aligned with the incident photon propagation direction and denoting the scattered wavevector by, $k'$, the Laue equations state that

$$C^T(k' - k) = 2\pi \begin{bmatrix} h \\ k \\ l \end{bmatrix}, \quad k^T k = k'^T k' = \frac{4\pi^2}{\lambda^2}, \qquad (1)$$

where $h, k$ and $l$ are Miller indices and the columns of the matrix $C$ are the crystal unit cell vectors in direct space. Eq. (1) represents a geometrical relation between the X-ray wavelength, $\lambda$, and the lattice plane separation, $d_{hkl}$, resulting in alignment of the outgoing scattered wave crests. Reducing the Laue equations (1) to 2D results in the well known Bragg's law,

$$2d_{hkl} \sin(\theta) = m\lambda, \quad m = 1, 2, 3, \dots \qquad (2)$$

where $\theta$ is the Bragg angle formed between $k$ and the diffracting crystal planes. Considering a small scalar strain, $s$, acting on the lattice planes we find from Bragg's law (2) that

$$2(1 + s)d_{hkl} \sin(\theta) = m\lambda \quad \rightarrow \quad \sin(\theta) = \frac{m\lambda}{2(1 + s)d_{hkl}}, \qquad (3)$$

i.e, the scattering angle, $\theta$, is sensitive to directional strain. Given some reference angle of scattering, $\theta_0$, corresponding to a strain free crystal, we find that the directional strain, $s$, existing along the diffraction vector, $G = k' - k$, can be expressed as

$$s = \frac{d_{hkl} - d_{hkl}^{(0)}}{d_{hkl}^{(0)}} = \frac{\sin(\theta_0) - \sin(\theta)}{\sin(\theta)}, \qquad (4)$$

where $d_{hkl}$ and $d_{hkl}^{(0)}$ are the strained and strain free lattice plane separations respectively. The geometrical context of the above discussed quantities is presented in Fig. 1.

For a more in depth explanation of the fundamental concepts of X-ray scattering see, for example, Als-Nielsen and McMorrow [37].

### 2.1. The 3DXRD/HEDM framework

The experimental 3DXRD/HEDM techniques exploit the Laue equations (1) to find the crystal lattices of the individual grains in a diffracting polycrystalline aggregate. The typical experimental setup is depicted in Fig. 2 where the sample, in this case, rotates about a single (vertical) axis throughout the measurement sequence. When the angle between the lattice planes and $k$ fulfil Bragg's law (2) diffraction spots are measured with the 2D detector placed behind the sample. The Laue equations (1) constrain $k'$ to fall on cones of increasingly larger opening angles. As a result, the diffraction pattern seen in Fig. 2 is characterised by non-zero intensity falling on rings centred around the detector and beam intersection point.

Commonly in the 3DXRD/HEDM type of diffraction applications, the unit cell matrix, $C$, is given in a multiplicative form as

$$C = U C_c, \qquad (5)$$
<figure><img src="image_4.png" /><figcaption>Fig. 1. Conceptual drawing of scattering from a set of crystal lattice planes. The difference between incident, $k$, and scattered, $k'$, wavevector forms the diffraction vector $G$, which is orthogonal to the scattering planes. A perturbation to the lattice plane distance $d_{hkl}$ will be observed as a change in Bragg angle $\theta$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. A polycrystalline aggregate is rotated about a single axis, $\hat{z}_s = \hat{z}_l$, while fully illuminated by a beam of high energy X-rays with wavevector along $\hat{x}$. The recorded digital images of spotty diffraction data (black dots) can be used to infer the grain positions and orientations. The Cartesian sample coordinates follow the turntable while the Cartesian laboratory frame is fixed with respect to sample motion. A diffraction image is recorded at each turntable rotation $(\omega_0, \omega_1, \omega_2...)$.</figcaption></figure>

where the crystal rotation matrix, $U$, acts to rotate the unit cell vector, now contained in crystal coordinates in $C_c$, into the sample coordinate system. This multiplicative split is motivated by the fact that, while the matrix $C_c$ is fairly constant between crystals sharing a phase in the sample, the orientation of the crystals, $U$, can vary drastically between grains. The primary task in 3DXRD/HEDM is to find a set of crystal orientations, $U$, such that the measured data can be explained. This procedure is denoted as grain mapping or grain indexing.

Introducing $B = 2\pi C_c^{-T}$ and $G_{hkl} = [h \ k \ l]^T$, we write

$$G_s = U B G_{hkl}, \qquad (6)$$

where subscript $s$ has been used to denote that $G_s$ is given in sample coordinates. As the sample is rotated by an angle $\omega$ around $\hat{z}_s$ the corresponding rotation matrix, $\Omega$, provides the map between laboratory and sample coordinates. The central equation for 3DXRD/HEDM is, therefore,

$$G_l = \Omega U B G_{hkl}, \qquad (7)$$

where $G_l$ can be computed from observations of $k'_l$. Thus, given a set of unit cell parameters (defining $B$), a selected lattice plane (defining $G_{hkl}$), a crystal orientation ($U$) and a sample rotation ($\Omega$), Eq. (7), can be used to predict $G_l = k'_l - k_l$. Given experimental parameters, such as detector size and position, these equations provide a model connecting the measured diffraction patterns and the properties of the polycrystalline aggregate. Indeed, sophisticated algorithms have been developed, and refined over years, to
reconstruct the set of **U** and **B** matrices that correspond to the individual grains within the aggregate, e.g., [11,12,27]. Once estimated, any proceeding analysis can operate on a per grain basis, separating out the parts of the diffraction images that are associated to any given crystal.

## 2.2. Histograms of directional strain

The recording on the detector of diffraction peaks generated by a crystal with known orientation, **U**, are digital images representing a discrete rendering of an underlying 2D photon intensity distribution. The formation of this intensity distribution is a nontrivial process involving, among other things; the crystal shape, the crystal deformation, the beam coherence, the energy bandwidth of beam, the detector point spread function and the detector pixel size [38,39]. To recover the strain tensor probability density associated to the given crystal, we are interested in modelling the impact of the crystal strain to the peak shape. To this end we shall denote the spatial coordinate with $x \in \mathbb{R}^3$ and let the deformation gradient tensor $F \in \mathbb{R}^{3 \times 3}$ map from reference, $x_0$, to current deformed configuration as

$$ F x_0 = x. \qquad (8) $$

We denote the undeformed reference unit cell, $C_0$, and the current deformed unit cell matrix as $C_{cur}$. Let $G[\cdot]$ be the diffraction operator that, given an input unit cell matrix, returns the corresponding set of diffraction vectors, $\{G^1, G^2, \dots, G^m\}$, observed in 3DXRD/HEDM. Clearly, the measured set of diffraction vectors originates from the deformed configuration as

$$ G[C_{cur}] = \{G^1, G^2, \dots, G^m\}. \qquad (9) $$

However, there exist modes of deformation, which we, here, denote by $F_p$, to which the diffraction set is immutable, such that

$$ G[F_p C_0] = G[C_0]. \qquad (10) $$

For instance, dislocation glide that results in a constant-volume plastic deformation and leaves the crystal lattice unchanged, has no measurable effect on the diffraction patterns in 3DXRD/HEDM. This motivates the use of a multiplicative split of the deformation gradient tensor similar to that used in crystal plasticity [40–44],

$$ F = F_e F_p, \qquad (11) $$

where the elastic deformation, $F_e$, is defined, here, as any deformation to the crystal that will cause a nonzero perturbation in the measured set of diffraction vectors,

$$ G[F_e C_0] \neq G[C_0]. \qquad (12) $$

The plastic contribution, $F_p$, here maps from the reference state to an intermediate stress-free configuration while $F_e$ maps from the intermediate configuration to the current deformed configuration. See Appendix A Fig. 13 for a summary of the discussed configurations.

Building on the above arguments, we define a fictive configuration, **C**, that would result from the application of the elastic deformation to the reference unit cell,

$$ C = F_e C_0, \qquad (13) $$

with the special property,

$$ G[C] = G[C_{cur}]; \qquad (14) $$

i.e the fictive unit cell matrix, **C**, is diffraction equivalent to the true deformed unit cell matrix, $C_{cur}$. To relate the elastic, diffraction observable, deformation, $F_e$, to the crystal strain tensor field, we define the total Green–Lagrange strain in the intermediate stress free configuration [43] as

$$ e = \frac{1}{2} F_p^{-T} (F^T F - I) F_p^{-1} = \underbrace{\frac{1}{2} (F_e^T F_e - I)}_{e_e} + \underbrace{\frac{1}{2} (I - F_p^{-T} F_p^{-1})}_{e_p}, \qquad (15) $$

where the total strain, $e$, is decomposed into an elastic part, $e_e$, and a plastic part, $e_p$. Working in the limit of small strains, we are now ready to define our ultimate goal that is to reconstruct the probability density distribution of the elastic strain tensor field, $e_e$. By Taylor expansion we find

$$ e_e = e_e(x) = \begin{bmatrix} e_{11}(x) & e_{12}(x) & e_{13}(x) \\ e_{12}(x) & e_{22}(x) & e_{23}(x) \\ e_{13}(x) & e_{23}(x) & e_{33}(x) \end{bmatrix} \approx \frac{1}{2} (F_e^T + F_e) - I, \qquad (16) $$

(see e.g., [45] for an introduction to continuum-mechanics). Since all unit cell matrices are invertible, by definition, it follows from Eq. (13) that

$$ F_e = C C_0^{-1}. \qquad (17) $$
Using the Laue Equations (1), we find a reference diffraction vector, $G_0$, and a measured diffraction vector, $G$, as

$$G_0 = 2\pi C_0^{-T} G_{hkl}, \quad G = 2\pi C^{-T} G_{hkl}. \qquad (18)$$

Combining Eqs. (17) and (18) we find that

$$G = F_e^{-T} G_0. \qquad (19)$$

Multiplying Eq. (16) from left and right by $G$ and making use of Eq. (19), we find

$$G^T e_e G = \frac{1}{2} (G^T G_0 + G_0^T G) - G^T G, \qquad (20)$$

Normalising Eq. (20) by $G^T G$, we now arrive at a scalar directional strain measure, $s_{hkl} = \hat{\kappa}^T e_e \hat{\kappa}$, as

$$s_{hkl} = \hat{\kappa}^T e_e \hat{\kappa} = \frac{G^T G_0}{G^T G} - 1, \qquad (21)$$

where

$$\hat{\kappa} = \frac{G}{\sqrt{G^T G}}. \qquad (22)$$

Note that for small elastic deformations the approximation

$$\hat{\kappa} = \frac{G}{\sqrt{G^T G}} \approx \frac{G_0}{\sqrt{G_0^T G_0}}, \qquad (23)$$

is also valid. Likewise, from Eq. (21), it follows that the directional strain, $s_{hkl}$, is only weakly dependent on lattice rotations, as $G^T G = B^T U^T U B = B^T B$ and $G_0^T G = B_0^T U_0^T U B$ where $U_0^T U \approx 1$ for moderate deformations. This motivates the decoupling of strain reconstruction from that of lattice rotations in the presence of small deformations. While it may be possible to extend the derivations presented above, and in the following, to a large deformation setting we limit the scope of this paper to that of small elastic strains.

Considering that each of the detector pixels composing a single diffraction peak is associated to a distinct diffraction vector, $G$, we see from Eq. (21) that also $s_{hkl}$ will vary over the diffraction peak. As the pixel intensity values of the diffraction peak are proportional to the diffracting sub volumes of the grain we have access to a histogram, $h(s_{hkl})$, over directional strains in the crystal. To form this histogram each normalised pixel intensity weight is to be added to the corresponding bin-count for the associated directional strain value, $s_{hkl}$. Note that this procedure reduces the information contained in the diffraction peak to a one dimensional function in $s_{hkl}$. Importantly, these histograms of directional strain are directly linked to the underlying spatial strain tensor field as

$$h(s_{hkl}) = \frac{1}{V_{\Omega}} \int_{s_{hkl}-w/2}^{s_{hkl}+w/2} \int_{\mathbb{R}^3} I(\hat{\kappa}^T e_e(\mathbf{x})\hat{\kappa} - s')dx ds', \quad V_{\Omega} = \int_{\Omega} dx, \qquad (24)$$

where $I$ is an indicator function supported on the crystal grain domain $\Omega$, $V_{\Omega}$ is the volume of the grain and $w$ is the histogram bin width. To simplify our following analysis we introduce a flattened vector notation similar to Henningsson and Hendriks [24]. Utilising the symmetry of the strain tensor we write

$$E = \begin{bmatrix} e_{11} \\ e_{22} \\ e_{33} \\ e_{12} \\ e_{13} \\ e_{23} \end{bmatrix}, \quad \bar{\kappa} = \begin{bmatrix} \kappa_x^2 \\ \kappa_y^2 \\ \kappa_z^2 \\ 2\kappa_x \kappa_y \\ 2\kappa_x \kappa_z \\ 2\kappa_y \kappa_z \end{bmatrix}, \qquad (25)$$

and find that the inner product in Eqs. (21) and (24) become

$$s = \hat{\kappa}^T e_e(\mathbf{x}) \hat{\kappa} = \bar{\kappa}^T E. \qquad (26)$$

With the introduction of the histograms of directional strains we are now ready to formally define our sought quantity; the strain PDF. Once this additional definition is made we can state our full inverse problem formulation.

### 3. Problem formulation

We seek to estimate the probability density function over strain tensors (strain PDF), $p_E(\epsilon)$, that describes the scalar probability of encountering any one particular strain tensor, $\epsilon$, at a randomly (uniformly) selected point, $\mathbf{x}$, in the grain. Formally, we define the strain PDF, $p_E: \mathbb{R}^6 \to \mathbb{R}^1$, associated to a spatial strain tensor field, $E(\mathbf{x})$, as

$$p_E(\epsilon) = \frac{1}{V_{\Omega}} \int_{\mathbb{R}^3} I(E(\mathbf{x}) - \epsilon)d\mathbf{x}, \quad V_{\Omega} = \int_{\Omega} d\mathbf{x}, \qquad (27)$$
<figure><img src="image_6.png" /><figcaption>Fig. 3. The mapping of a spatial field, E, supported on the unit disc, centred at x = [0.5,0.5], results in a new spatial field, E'(x) = E(R[x]), such that E(x) ≠ E'(x). Nevertheless, the two corresponding strain PDFs remain identical, p_E(ε) = p_E'(ε), as long as the map, R[x], is area preserving.</figcaption></figure>

and state the inverse problem: Given a set of histograms, h₁(s_hkl), h₂(s_hkl), ..., hₘ(s_hkl), all of which are associated with a single grain, we seek an estimation of the strain PDF, p_E(ε), of the given grain, such that the measured histograms can be plausibly explained.

We emphasise that each histogram, hᵢ(s_hkl), corresponds to a single diffraction peak and measures directional strain in a known distinct direction, $\bar{\kappa}_i$, in the indexed crystal. The task is therefore to find a parametric expression for p_E(ε) such that the measured histograms, h₁(s_hkl), h₂(s_hkl), ..., hₘ(s_hkl), are well modelled. To achieve this we will need to establish the transformation from strain PDF to directional histograms of strains. Before pursuing this task, however, we shall discuss some general properties and limitations of the strain PDF. Especially, we emphasise how the strain PDF provides richer information about the grain deformation state compared to the grain average strain tensor.

## 4. Properties of the strain probability density function

We note that while the strain PDF exists for any one spatial strain field E(x), the mapping is not necessarily unique. As an example, let R: $R^3 \to R^3$ be an area preserving function such that det(R') = 1, ∀x, where R' is the Jacobian of R. Making use of this map we define a modified spatial strain field as E'(x) = E(R[x]) with y = R[x]. From the definition of the strain PDF it now follows that

$$ 
\begin{aligned}
p_E(\epsilon) &= \frac{1}{V_\Omega} \int_{\mathbb{R}^3} I(E(y) - \epsilon) dy = \frac{1}{V_\Omega} \int_{\mathbb{R}^3} I(E(R[x]) - \epsilon) dR[x] = \\
&= \frac{1}{V_\Omega} \int_{\mathbb{R}^3} I(E'(x) - \epsilon) \det(R') dx = \frac{1}{V_\Omega} \int_{\mathbb{R}^3} I(E'(x) - \epsilon) dx = p_{E'}(\epsilon).
\end{aligned}
\quad (28) $$

Simple examples of R on $R^3$ includes combinations of rigid body rotations and translations. In Fig. 3 we give a nonlinear example on $R^2$, mapping the unit disc centred at x = [0.5, 0.5] through the area preserving map

$$ R[x] = \left[ \frac{\sqrt{x_1}}{2x_2\sqrt{x_1}} \right], \quad x > 0. \quad (29) $$

The strain PDF provides weaker information on the deformation of the domain compared to the spatial strain tensor field itself. Nevertheless, in relation to only knowing the mean value of a spatial strain tensor field, the strain PDF is a rich quantity. We provide an illustration of the concept of a strain PDF in Fig. 4. In this example a Finite Element solver was used to solve for the equilibrium strain field given a symmetric (across y) geometry and boundary load. The body was considered linear elastic with width and height both set to 1 µm, while the out of x - y thickness was set to 100 µm such that plane strain conditions prevailed. Owing to the symmetry, a perturbation in Young's modulus, E, will correspond to a simple scaling of the strain field. As illustrated in the projections of the corresponding strain PDFs in Fig. 4 (left column), the mean of the distribution is immutable to such a scaling of the spatial strain field (the cerise dot marking the mean does not move). The topology of the strain PDF does, however, change as a result of the perturbation in material stiffness, E. This motivates the value of retaining deeper statistics on the strain PDF, beyond the commonly derived mean value. Similar examples can be constructed with respect to perturbations in boundary load, object topology, Poisson ratio, etc.

## 5. Measurement model

To solve the inverse problem defined in Section 3 we establish a measurement model that can map the strain PDF to a set of histograms of directional strain. Let us therefore consider a single histogram, hⱼ[sⱼ], and let us inquire for the bin count at the bin centre, sⱼ = sⱼ^*, such that hⱼ^* = hⱼ[sⱼ^*], where sⱼ^* is the directional strain at the bin centre, sⱼ^* = $\bar{\kappa}^T E^*$. The subdomain on strain space holding the strain tensors, E^*, that fulfil this equality is then defined by

$$ \frac{1}{\|\bar{\kappa}\|_2} \bar{\kappa}^T E^* = \frac{s}{\|\bar{\kappa}\|_2}. \quad (30) $$

Eq. (30) defines a hyperplane on $R^6$ with unit normal $\hat{n} = \bar{\kappa}/\|\bar{\kappa}\|_2$ and scalar origin offset t = s/\|\bar{\kappa}\|_2. It follows that the histogram bin count, hⱼ^*, is given by the integral across all planes with normal $\hat{n}$ and origin of-sets t ∈ [sⱼ^* - w/2, sⱼ^* + w/2] where w is the
<figure><img src="image_7.png" /><figcaption>Fig. 4. Example of the topological shape of a strain PDF (left) revealed in part by its projection onto one of the 15 Cartesian planes in strain space (zoomed in around the mean strain which is marked out by a cerise dot and a dashed cross). The spatial strain field (right) is the linear elastic Finite Element solution of a symmetric plane strain boundary value problem using 4944 triangular constant strain elements. Changing the stiffness of the material (top row vs. bottom row) results in a scaled FEM solution. While the mean of the strain PDF is invariant to a change in stiffness, the full topology of the strain PDF mutates with the change in material stiffness, E (compare left top and left bottom). This illustrates that the strain PDF is a richer quantity than the mean of the spatial strain tensor field.</figcaption></figure>

histogram bin width. Thus,

$$h_j^*[s_j^*] = \int_{s_j^*-w/2}^{s_j^*+w/2} \int_{\epsilon \in \Pi(t, \hat{n}_j)} p_E(\epsilon) d\epsilon dt, \quad (31)$$

where $\Pi(t, \hat{n}_j)$ is the hyperplane with normal $\hat{n}_j$ and scalar origin offset $t = s/\|\bar{\kappa}\|_2$. This represents a linear measurement model and requires the computation of a single volume integral per measured histogram bin. Let us now define the linear operator $M_j$ which maps any function, $f(\epsilon)$, to k histogram bins modelling all histogram bin counts for the j:th histogram as

$$M_j[f(\epsilon)] = \begin{bmatrix} \int_{s_j^*-w/2}^{s_j^*+w/2} \int_{\epsilon \in \Pi(t, \hat{n}_j)} p_E(\epsilon) d\epsilon dt \\ \vdots \\ \int_{s_j^*-w/2}^{s_j^*+w/2} \int_{\epsilon \in \Pi(t, \hat{n}_j)} p_E(\epsilon) d\epsilon dt \end{bmatrix} = \begin{bmatrix} h_j[s_j^*] \\ \vdots \\ h_j[s_j^k] \end{bmatrix}, \quad (32)$$

where $s_j^k$ denotes the k:th bin centre of the j:th histogram. A closed form expression for the inverse operator, $M_j^{-1}$, is not known. However, by selecting a parametric expression for the strain PDF and a cost function, measuring the difference between model and measured histograms of directional strain, gradient based optimisation techniques can be used to recover the parameters defining the strain PDF. We therefore proceed to introduce a finite basis expansion of the strain PDF.

# 6. Reconstruction approach

## 6.1. Finite basis strain PDF expansion

Decomposing the strain PDF into a linear combination of n basis functions, $\phi_i$, as

$$p_E(\epsilon) = \sum_{i=1}^{n} c_i \phi_i(\epsilon), \quad c_i, \phi_i \in \mathbb{R}^1, \quad (33)$$
we find the maximum-likelihood estimator in the presence of Gaussian noise as

$$ \underset{c_i}{\operatorname{argmin}} \sum_{j}^{m} \|h_j - \mathcal{M}_j[p_E]\|_2^2, \quad \text{s.t.} \quad p_E = \sum_{i}^{n} c_i \varphi_i(\epsilon) > 0. \qquad (34) $$

where $c_i$ are the basis coefficients to be optimised with respect to the scalar cost function. Since $\mathcal{M}_j$ is a linear operator it suffices to compute $\mathcal{M}_j[\varphi_i(\epsilon)]$ to evaluate the cost function and its gradient with respect to $c_i$. To achieve this we shall adopt a radial basis function as

$$ \varphi_i(\epsilon) = \exp\left(-\frac{(\epsilon - p_i)^T (\epsilon - p_i)}{2\sigma^2}\right). \qquad (35) $$

The benefit of this selection is twofold. Firstly, using the error function, we will be able to derive a closed form solution to Eq. (31) as a sum of error functions. Retaining analytical expressions for the forward model is paramount as the alternative of numerical integration scales poorly with the number of dimensions (which is $\mathbb{R}^6$ for our application). Secondly, the selection in (35) simplifies the requirement $p_E > 0$ to a linear constraint in the unknowns $c_i > 0$. The optimisation problem given in (34) is therefore rendered convex.

## 6.2. Basis function integration

To execute our forward model, $\mathcal{M}_j[\varphi_i]$, we seek the solution of the integral of the form

$$ \int_{s_j^* - w/2}^{s_j^* + w/2} \int_{\epsilon \in \Pi(\hat{\mathbf{n}}_j)} \exp\left(-\frac{(\epsilon - p_i)^T (\epsilon - p_i)}{2\sigma^2}\right) d\epsilon dt. \qquad (36) $$

Analysis of Eq. (36) (see Appendix C) leads to a bounded Gaussian integral which can be solved using the error function, erf, as

$$ I_t = 4\pi^3 \sigma^6 \left| \operatorname{erf}\left(\frac{b}{\sigma\sqrt{2}}\right) - \operatorname{erf}\left(\frac{a}{\sigma\sqrt{2}}\right) \right|, \qquad (37) $$

where $a = (s_j^* - w/2)p_i^T \hat{\mathbf{n}}_j$ and $b = (s_j^* + w/2)p_i^T \hat{\mathbf{n}}_j$ respectively.

## 6.3. Maximum-likelihood estimate

In the presence of additive Gaussian noise to the measured histogram bin counts we find the maximum-likelihood estimator of $c_i$ as

$$ \underset{c_i}{\operatorname{argmin}} \quad \|h - Lc\|_2^2, \quad \text{s.t} \quad c > 0, \qquad (38) $$

where $c \in \mathbb{R}^n$ is a column vector holding the basis coefficients, $c_i$, the matrix $L \in \mathbb{R}^{M \times n}$ contains the model integral over each individual basis function and $h \in \mathbb{R}^M$ is a column vector stacking all histogram bin values for $j = 1, 2, \dots, m$ in accordance with $L$. Specifically,

$$ L = \begin{bmatrix} \mathcal{M}_1[\varphi_1] & \dots & \mathcal{M}_1[\varphi_n] \\ \mathcal{M}_2[\varphi_1] & \dots & \mathcal{M}_2[\varphi_n] \\ \vdots & \ddots & \vdots \\ \mathcal{M}_m[\varphi_1] & \dots & \mathcal{M}_m[\varphi_n] \end{bmatrix}, \quad c = \begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_n \end{bmatrix}, \quad h = \begin{bmatrix} h_1(s^1) \\ h_1(s^2) \\ \vdots \\ h_1(s^k) \\ \vdots \\ \vdots \\ h_m(s^1) \\ h_m(s^2) \\ \vdots \\ h_m(s^k) \end{bmatrix}. \qquad (39) $$

To solve for $c$ in Eq. (38) we use the active set method described by Lawson and Hanson [46] and implemented in `scipy. optimize.nnls` [47].

## 6.4. Basis placement

Selecting basis function positions, $p_i$, over a uniform rectangular grid in $\mathbb{R}^6$ and letting $n \to \infty$ while $\sigma \to 0$ there exist $c_i > 0$ such that any positive function $p_E$ can be approximated. While this is technically true, for any practical application we have to use a finite number of basis functions, $n$, each with some finite, nonzero, radius $\sigma$. One could imagine that a grid type approach to basis placement will maintain a good performance as long as $n$ is some large finite number. However, owing to the number of dimensions of the problem, this approach scales poorly. To illustrate, consider a target strain PDF which is believed to hold relevant
<figure><img src="image_8.png" /><figcaption>Fig. 5. Illustration of increasingly good basis approximations retrieved from the iteration of Algorithm 1. The current model function and the residual field are plotted at the end of each batch. The black circles correspond to basis function locations. The radius of the basis function was relaxed as $\sigma = [0.25, 0.1, 0.06, 0.03]$ with $B = 20$ and $J = 20$.</figcaption></figure>

information in a neighbourhood of $\pm 25 \times 10^{-4}$ units of strain around some known mean value. Targeting an approximate resolution of $\pm 1 \times 10^{-4}$ units of strain and producing an equidistant grid over strain space centred around the mean strain, we find that we need a total of $\sim 10^8$ basis functions. The corresponding computation of $L \in \mathbb{R}^{M \times 10^8}$ is impractical, if not altogether unfeasible. Moreover, considering the underlying physics, $p_E$ may be narrowly supported on $\mathbb{R}^6$. For instance, in the case when the underlying spatial strain tensor field is in a state of plane strain, as is the case in Fig. 4, $p_E$ will be fully contained by a 3-dimensional hyperplane on $\mathbb{R}^6$. In a more general case, the mechanical equilibrium equations constrain the support of the strain PDF. This means that the placement of basis functions on a rectangular grid runs the risk of wasting a large portion of the basis set. As an example, in the case of plane strain, only $\approx 25^3 (\approx 10^4)$ out of the $10^8$ previously mentioned basis function would be active in approximating $p_E$, effectively wasting $\approx 99.9936\%$ of the basis function set. Thus, we are presented with the challenge of distributing a sparse set of basis functions on $\mathbb{R}^6$ such that the final approximation of $p_E$ can explain the observed histograms.

### 6.4.1. Sequential sampling approach

To address the above described sampling issue, we present a class of stochastic basis placement algorithms that are heuristically expected to perform well for mapping out the support of smooth, simply connected, positive functions. The central idea evolves around producing a sparse sequence of basis sets such that the explanatory value of the basis sets with respect to measured data increases with iteration. The approach is, in part, similar to sequential Monte Carlo methods [48] as the proposal distribution, from which new basis locations are drawn, is allowed to mutate between iterations. However, to the best of the authors knowledge, the method we use has not been reported on previously. We summarise our approach in Algorithm 1. The procedure can be summarised as follows. First a single basis function is placed at a location, $\epsilon$, at which the target function $p_E$ is believed to be supported. For instance, we have found that putting the initial basis function at the mean strain of the grain performs well. Next a basis coefficient $c_1$ for our single basis function is selected such that $p_E$ becomes a probability distribution. From this probability distribution $B$ samples, $\epsilon_1, \epsilon_2, \dots, \epsilon_B$, are randomly drawn. At each sampled location, $\epsilon_i$, a basis function is now placed such that the approximation of $p_E$ is expanded to contain $B+1$ basis functions. The coefficients of the basis set is now recomputed such that the likelihood of observing the strain histogram data is maximised, i.e a least squares problem is solved. Any basis functions that appears to have a small explanatory value, i.e a small fitted coefficient, $c_i \le \xi$, are now removed from the basis expansion. The value of $\xi$ is selected arbitrarily with higher values resulting in a sparser number of basis functions and less freedom for the algorithm to match the measured data, while a smaller value of $\xi$ will result in a more dense basis set and more freedom for the algorithm to match the measured data. For instance, we found that setting $\xi = 0$, and thus keeping all basis functions that are on the support of $p_E$, performed well in our strain estimation problem. After the basis set has been pruned of any unwanted basis functions the resulting function $p_E$ is again normalised to represent a probability distribution. A new sample of $B$ locations, $\epsilon_1, \epsilon_2, \dots, \epsilon_B$, is drawn from the updated distribution and the procedure of basis coefficient fitting is repeated. This iteration is allowed to continue for a fixed number of redraw iterations $J$. For large selections of $J$ the root mean squared error between of the maximum likelihood fit is seen to plateau, i.e a large selection of $J$ will in general perform better than a small selection of $J$, with diminishing returns as $J$ grows. So far we have summarised the inner most sub-iteration of Algorithm 1 which we denote as a batch run. After this sub-iteration/batch comes to an end the standard deviation (or radius) of the radial basis functions, $\sigma$, is decreased and the entire procedure is repeated, carrying forward all basis functions accumulated in the previous sub-iteration. Note that while the radius, $\sigma$, varies globally between iteration it does not vary between the individual basis functions in the same basis expansion set. Basis locations are now again randomly drawn using the new, smaller, $\sigma$ value in the basis expansion of $p_E$ and both old and new basis functions are marked for pruning as previously. After
K sub-iterations of the algorithm, each featuring J basis location draws of size B, the algorithm terminates at the final smallest length-scale, $\sigma_K$. By decreasing $\sigma$ slowly, in many fine steps, better results were found in general, at the expense of computational complexity. For instance, in the result presented in 8 we used $(\sigma_1, \sigma_2, \sigma_3, \sigma_4, \sigma_5) = (10.0, 8.4, 7.1, 5.9, 5.0)$ and $B = J = 25$. As $\sigma$ is gradually diminished the algorithm refines the approximation of $p_E$ and is able to make use of an increasingly larger set of basis functions.

<caption>**Algorithm 1:** A Stochastic Radial Basis Placement Strategy</caption>

Result: $n, \sigma, p_i, \forall i = 1, \dots, n$,
Select batch size, $B$, and number of batch iterations, $J$;
Select a relaxation sequence $[\sigma_1, \dots, \sigma_K]$, with $\sigma_{k+1} < \sigma_k$;
Select $p_1$ on the support of $p_E$;
Set $p_E^{(1)}(\epsilon) = c_1\varphi_1(\epsilon)$;
Fit $c_1$ with respect to data (solve Equation (38));
for $\sigma_k \in [\sigma_1, \dots, \sigma_K]$ do
    Set $\sigma = \sigma_k$;
    for $j = 1, \dots, J$ do
        Draw $B$ new basis locations, $p_i$, from $p_E^{(kj)}$;
        Let $n \leftarrow n + B$;
        Extend expansion $p_E^{(kj)}(\epsilon) = \sum_{i=1}^{i=n} c_i\varphi_i(\epsilon)$;
        Fit $c_i$ with respect to data (solve Equation (38));
        Set prune count, $p = 0$;
        for $i = 1, \dots, n$ do
            if $c_i \le \xi$ then
                Mark $\varphi_i$ for pruning;
                Increment prune count $p \leftarrow p + 1$;
            end
        end
        Prune marked $\varphi_i$ from expansion $p_E^{(kj)}(\epsilon) = \sum_{i=1}^{i=n} c_i\varphi_i(\epsilon)$;
        Decrement basis set count $n \leftarrow n - p$;
    end
end
Accept $p_E = p_E^{(KJ)} = \sum_{i=1}^{i=n} c_i\varphi_i(\epsilon)$;

We illustrate the evolving basis expansion produced by our Algorithm 1 in Fig. 5. For illustrative purposes the example case takes place on $\mathbb{R}^2$ with the target data selected as a direct observation of the underlying target function (rightmost plot in Fig. 5). The residual between the current model function and the target function is seen to decrease as the algorithm produces a sequence of customised basis sets. The algorithm parameters where set to $B = 20$, $J = 20$ and $\sigma = [0.25, 0.1, 0.06, 0.03]$. This means that the algorithm can theoretically attain a maximum of $20 \times 20 \times 4 = 1600$ basis functions over the four batch runs. We note that the final number of accumulated basis functions produced by our algorithm is much sparser ($n = 30$). This means that although we have to solve the inverse problem of fitting $c_i$ 80 times throughout the run, each inversion is fast owing to the low number of basis functions. This property is especially beneficial when the algorithm used to solve for $c_i$ scales super-linearly with $n$. A concise python implementation of Algorithm 1 in the setting of Fig. 5 is openly available at https://gist.github.com/AxelHenningsson/8a9179f859751634859eb7e051ebd804.

## 7. Observability and reduced models

It is well known that the mean of the strain PDF can be recovered from 6 measured diffraction peaks [15]. This does however not indicate whether or not the strain PDF itself is well defined from measurements of directional strain. To answer this question one needs to investigate under what perturbations the strain PDF generates identical histograms of directional strains through the model, $M[p_E]$. For an arbitrary functional form of the strain PDF a complete description of this null-space is challenging. However, a feasible generalisation of a mean value analysis is to let the strain PDF be a multivariate Gaussian.

In the following we shall derive a class of perturbations to multivariate Gaussian strain PDFs that leave the measured histograms unchanged. From this analysis we will also find a closed form estimator of the strain PDF. We start our analysis by reiterating why the mean value of the strain PDF is well defined.

### 7.1. First moments

By definition, the expected value of the strain PDF is

$$ \mathbb{E}[E] = \int E p_E(\epsilon) d\epsilon, \qquad (40) $$
where $\mathbb{E}[\cdot]$ denotes the expectation operator. Considering the expected value of any one single measured directional strain histogram, $h_j(s_j)$, and using the linearity of expectation, we have

$$ \mathbb{E}[s_j] = \mathbb{E}[\bar{\kappa}_j^T \bar{E}] = \bar{\kappa}_j^T \mathbb{E}[E], \quad (41) $$

where the index $j$ denotes association to a fixed direction, $\bar{\kappa}_j$. We note that (41) holds true regardless of the functional form of $p_E(\epsilon)$. Denoting $\mu_E = \mathbb{E}[E]$ and collecting a set of $\bar{\kappa}_j$ and $\mathbb{E}[s_j]$, with $j = 1, \dots, m$, in the matrices $K$ and $\mu_s$, as

$$ K = \begin{bmatrix} \bar{\kappa}_1^T \\ \bar{\kappa}_2^T \\ \vdots \\ \bar{\kappa}_m^T \end{bmatrix}, \quad \mu_s = \begin{bmatrix} \mathbb{E}[s_1] \\ \mathbb{E}[s_2] \\ \vdots \\ \mathbb{E}[s_m] \end{bmatrix}, \quad (42) $$

we may produce a linear system of equations,

$$ K \mu_E = \mu_s. \quad (43) $$

From Eq. (43) we conclude that $\mu_E$ is well defined if and only if $K$ is full column rank. This requires the collection of a minimum of 6 linearly independent $\bar{\kappa}$ vectors. As an illustration, the choice

$$ \begin{align} \hat{\kappa}_1 &= [1 \quad 0 \quad 0], & \hat{\kappa}_2 &= [0 \quad 1 \quad 0], & \hat{\kappa}_3 &= [1 \quad 1 \quad 1] / \sqrt{3} \notag \\ \hat{\kappa}_4 &= [1 \quad 1 \quad 0] / \sqrt{2}, & \hat{\kappa}_5 &= [1 \quad 0 \quad 1] / \sqrt{2}, & \hat{\kappa}_6 &= [0 \quad 1 \quad 1] / \sqrt{2}, \tag{44} \end{align} $$

produces the fully ranked $K$ matrix as

$$ K = \frac{1}{3} \begin{bmatrix} 3 & 0 & 0 & 0 & 0 & 0 \\ 0 & 3 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 2 & 2 & 2 \\ 3/2 & 3/2 & 0 & 3 & 0 & 0 \\ 3/2 & 0 & 3/2 & 0 & 3 & 0 \\ 0 & 3/2 & 3/2 & 0 & 0 & 3 \end{bmatrix}. \quad (45) $$

We note that in this example no $\bar{\kappa}$ vector is closer than $45^\circ$ to the $\hat{z}_s = [0, 0, 1]$ axis. This explains why it is possible to reconstruct mean crystal grains strains in 3DXRD/HEDM using only a single axis of rotation. We note that in practice, to reduce experimental error, it is common to solve an overdetermined system, using more than 6 observations [15].

## 7.2. The multivariate Gaussian case

Let the strain PDF be a multivariate Gaussian as

$$ p_E(\bar{\epsilon}) = n \exp \left( -\frac{1}{2} (\bar{\epsilon} - \mu_{\bar{\epsilon}})^T \Sigma^{-1} (\bar{\epsilon} - \mu_{\bar{\epsilon}}) \right), \quad n = \frac{1}{\sqrt{(2\pi)^6 \det[\Sigma]}}, \quad \Sigma = \Sigma^T \in \mathbb{R}^{6 \times 6}, \quad \bar{\epsilon} \in \mathbb{R}^{6 \times 1}, \quad (46) $$

such that $p_E(\bar{\epsilon})$ describes the scalar probability of encountering any one particular strain tensor, $\bar{\epsilon}$. Considering $E \sim p_E(\bar{\epsilon})$ as a stochastic variable, we can then form a set of scalar stochastic variables (directional strains) as

$$ S_j = \bar{\kappa}_j^T E, \quad j = 1, 2, 3, \dots, m. \quad (47) $$

Since the transformation from $E$ to $S_j$ is linear, it follows that $S_j \sim p_{S_j}(s)$ is Gaussian with PDF

$$ p_{S_j}(s) = a \exp \left( -\frac{1}{2} (s - \bar{\kappa}_j^T \mu_{\bar{\epsilon}}) (\bar{\kappa}_j^T \Sigma \bar{\kappa}_j)^{-1} (s - \bar{\kappa}_j^T \mu_{\bar{\epsilon}}) \right), \quad a = \frac{1}{\sqrt{2\pi \det[\bar{\kappa}_j^T \Sigma \bar{\kappa}_j]}}, \quad j = 1, 2, 3, \dots, m. \quad (48) $$

Now, considering that the measurements, $h_j(s_j)$, are the histograms formed through the binning of the PDFs, $p_{S_j}(s)$, it follows that, if $p_{S_j}(s)$ is invariant under a perturbation of $p_E(\bar{\epsilon})$, so are the measurements, $h_j(s_j)$. The task of characterising the null-space is, therefore, reduced to that of finding perturbations to $p_E(\bar{\epsilon})$ such that $p_{S_j}(s)$ remains unaltered.

### 7.2.1. Measurement invariant perturbations

To permute $p_E(\bar{\epsilon})$ inside the class of multivariate Gaussians, we can change either the mean or the covariance of the distribution. Our preceding analysis has already shown that the mean of $p_E(\bar{\epsilon})$ is well defined, regardless of the functional form. Thus, we can immediately turn our analysis to perturbations in $\Sigma$. To this end, we introduce $\Sigma^* = \Sigma + N$, where $N$ is some perturbation that preserves symmetry and positive definiteness, i.e.,

$$ N = N^T, \quad z^T(\Sigma + N)z \ge 0, \quad \forall z \ne 0. \quad (49) $$

We find the modified PDF over directional strains as

$$ p_{S_j}^*(s) = a \exp \left( -\frac{1}{2} (s - \bar{\kappa}_j^T \mu_{\bar{\epsilon}}) (\bar{\kappa}_j^T (\Sigma + N) \bar{\kappa}_j)^{-1} (s - \bar{\kappa}_j^T \mu_{\bar{\epsilon}}) \right), \quad a = \frac{1}{\sqrt{2\pi \det[\bar{\kappa}_j^T (\Sigma + N) \bar{\kappa}_j]}}, \quad (50) $$
Clearly

$$p_{S_j}^*(s) = p_{S_j}(s) \quad \text{iff} \quad \bar{\kappa}_j^T (\Sigma + N) \bar{\kappa}_j = (\bar{\kappa}_j^T \Sigma \bar{\kappa}_j), \qquad (51)$$

and it follows that the null-space requires the existence of N such that

$$\bar{\kappa}_j^T N \bar{\kappa}_j = 0, \quad \forall j. \qquad (52)$$

A thorough analysis of Eq. (52) will reveal 6 possible selections of N. The somewhat lengthy derivations of these 6 covariance perturbations are presented in Appendix B and we proceed here immediately to summarise the result of these calculations.

### 7.2.2. Parametric null-space

In conclusion, if the strain PDF is a Gaussian on $\mathbb{R}^6$, with positive definite covariance $\Sigma$, the null-space of the strain PDF on this function class is described by the addition of a non-positive-semi-definite symmetric matrix to the strain PDF Gaussian covariance matrix

$$\Sigma^* = \Sigma + N(\alpha, \beta, \gamma, \xi, \eta, \rho), \quad \text{s.t} \quad z^T \Sigma z + z^T N z > 0, \quad \forall z \in \mathbb{R}^6, \quad z \neq 0, \qquad (53)$$

where the class of matrices $N = N(\alpha, \beta, \gamma, \xi, \eta, \rho)$ is parameterised as

$$N = \begin{bmatrix} 0 & 2\alpha & 2\beta & 0 & 0 & 2\gamma \\ 2\alpha & 0 & 2\eta & 0 & 2\xi & 0 \\ 2\beta & 2\eta & 0 & 2\rho & 0 & 0 \\ 0 & 0 & 2\rho & -\alpha & -\gamma & -\xi \\ 0 & 2\xi & 0 & -\gamma & -\beta & -\rho \\ 2\gamma & 0 & 0 & -\xi & -\rho & -\eta \end{bmatrix}, \quad \alpha, \beta, \gamma, \xi, \eta, \rho \in \mathbb{R}. \qquad (54)$$

Since the scalars $\alpha, \beta, \gamma, \xi, \eta, \rho$ can be chosen arbitrarily small it will, in general, be possible to construct such perturbations while maintaining positive definite $\Sigma^*$.

### 7.2.3. Null-space verification

To verify that the null-space leaves the measured histograms unaltered, in Fig. 6 we show the distribution of the directional strain in the $\bar{\kappa}^T = \sqrt{50/27}[1/2, 1/2, 1/5]$ direction for a series of perturbations to an underlying simulated strain tensor probability density.

Each histogram in Fig. 6 was constructed by drawing 20 000 000 strain tensors at random from a multivariate Gaussian, where the unperturbed covariance matrix was selected as

$$\Sigma = \begin{bmatrix} 63 & 0 & 0 & 0 & 0 & 0 \\ 0 & 34 & 0 & 0 & 0 & 0 \\ 0 & 0 & 3 & 0 & 0 & 0 \\ 0 & 0 & 0 & 49 & 0 & 0 \\ 0 & 0 & 0 & 0 & 25 & 0 \\ 0 & 0 & 0 & 0 & 0 & 25 \end{bmatrix} \times 10^{-7}. \qquad (55)$$

As expected, the two histograms in Fig. 6 that are in the null-space described in Section 7.2.2 look identical (triangles and circles). Similarly, we see in Fig. 6 that perturbations to the strain distribution that are not in the null-space leads to changes in the histogram (squares and stars).

### 7.2.4. Maximum-likelihood estimator for a Gaussian strain PDF

The relationship between the first moment of the strain PDF and the histograms of directional strain was established in Eq. (43). To fully describe the Gaussian strain PDF we are, therefore, left with the estimation of the covariance matrix $\Sigma$. Letting $\mathbb{V}[\cdot]$ be the variance operator and denoting the (scalar) variance associated to a single histogram of directional strain by $\sigma_j^2$ we find with the use of Eq. (47) that

$$\sigma_j^2 = \mathbb{V}[S_j] = \mathbb{V}[\bar{\kappa}_j^T E] = \bar{\kappa}_j^T \mathbb{V}[E] \bar{\kappa}_j = \bar{\kappa}_j^T \Sigma \bar{\kappa}_j. \qquad (56)$$

Thus, by collecting the variance of each histogram in a column vector, $\sigma_s^2 = [\sigma_1^2, \sigma_2^2, \dots]^T$, a system of equations linear in the components of $\Sigma$ can be established. To express $\Sigma$ as a function of $\sigma_s^2$ in a least squares sense we introduce in Appendix B a flattening of Eq. (56) (in analogy with $\bar{\kappa}$ and $E$ in Eq. (26)). Letting the column vector $\bar{\Sigma}^* \in \mathbb{R}^{21 \times 1}$ hold the unique components of the strain PDF covariance matrix $\Sigma^* \in \mathbb{R}^{6 \times 6}$ and letting the matrix $\mathbb{V} \in \mathbb{R}^{m \times 21}$ operate the outer product maps defined in Eq. (56) we find that

$$\mu_E = (\mathbf{K}^T \mathbf{K})^{-1} \mathbf{K}^T \mu_s, \\ \bar{\Sigma}^* = (\mathbb{V}^T \mathbb{V})^{-1} \mathbb{V}^T \sigma_s^2 + \mathbf{u}(\alpha, \beta, \gamma, \eta, \xi, \rho), \qquad (57)$$

where the null-space of the covariance is described by the vector $\mathbf{u} \in \mathbb{R}^{21 \times 1}$ which is the flattening of $N$ in Eq. (54). The rows of $\mathbb{V}$ are explicitly described in Eq. (80) in Appendix B. Without any further information we set $\mathbf{u} = \mathbf{u}(0, 0, 0, 0, 0, 0) = \mathbf{0}$ and find in
<figure><img src="image_9.png" /><figcaption>Fig. 6. Zoom in on histograms of directional strain in the $K^T = [1/2, 1/2, 1/5]/\sqrt{27/50}$ direction corresponding to four realisations of each 20 000 000 strain tensors drawn from four distinct Gaussian strain PDFs. The covariance components of the four Gaussian strain PDFs, $\Sigma_j^*$, were modified according to the legend. It is clear to see that when the modification to the covariance is in the null-space (triangles) the original (circles) histogram is recovered. Other perturbations of the covariance (squares and stars) result in observable changes in the histogram.</figcaption></figure>

conclusion that

$$ 
\begin{aligned}
\mu_E &= (K^T K)^{-1} K^T \mu_s, \\
\Sigma^* &= (V^T V)^{-1} V^T \sigma_s^2.
\end{aligned} 
\qquad (58) $$

Eq. (58) can be thought of as an approximation to the strain PDF which, in general, does not need to be a Gaussian. The main benefits of the Gaussian approximation compared to the finite basis expansion solution presented in Section 6 is its simplicity and low computational complexity.

## 8. Simulation study

To test the derived framework we reconstruct the strain PDF for a synthetically generated data set. Importantly, the simulation framework described in this section operates directly on the spatial domain ($x \in \mathbb{R}^3$) for data generation. This ensures that we are not committing the inverse crime of using one and the same mapping for both regression and data generation. Additionally, as described in Section 8.4, we introduce challenging levels of noise and outliers to test the robustness of the proposed methods.

### 8.1. Strain field generation

Since the underlying spatial strain tensor field is expected to be in equilibrium on the unknown domain $\Omega$, we generate a strain tensor field that satisfies the point-wise equilibrium equations,

$$ 
\begin{aligned}
\frac{\partial \sigma_{11}}{\partial x} + \frac{\partial \sigma_{21}}{\partial y} + \frac{\partial \sigma_{31}}{\partial z} &= 0, & \sigma &= \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \end{bmatrix}, \\
\frac{\partial \sigma_{12}}{\partial x} + \frac{\partial \sigma_{22}}{\partial y} + \frac{\partial \sigma_{32}}{\partial z} &= 0, & \bar{\sigma} &= \begin{bmatrix} \sigma_{12} \\ \sigma_{13} \\ \sigma_{23} \end{bmatrix} = DE, \\
\frac{\partial \sigma_{31}}{\partial x} + \frac{\partial \sigma_{32}}{\partial y} + \frac{\partial \sigma_{33}}{\partial z} &= 0 & &
\end{aligned} 
\qquad (59) $$

where D is a 6 × 6 isotropic material stiffness matrix,

$$ D = \frac{E}{(1 + \nu)(1 - 2\nu)} \begin{bmatrix} 1 - \nu & \nu & \nu & 0 & 0 & 0 \\ \nu & 1 - \nu & \nu & 0 & 0 & 0 \\ \nu & \nu & 1 - \nu & 0 & 0 & 0 \\ 0 & 0 & 0 & (1 - 2\nu) & 0 & 0 \\ 0 & 0 & 0 & 0 & (1 - 2\nu) & 0 \\ 0 & 0 & 0 & 0 & 0 & (1 - 2\nu) \end{bmatrix}. \qquad (60) $$

To generate a random strain field satisfying (60), we constructed a zero mean Gaussian process, similar to Hendriks et al. [23], using a squared exponential kernel that we mapped through the partial differential operator that encodes the Maxwell equilibrium
<figure><img src="image_10.png" /><figcaption>Fig. 7. Volume rendering of spatial strain tensor field in units of [10⁻⁴] sampled over a spherical domain using an equilibrium encoded Gaussian process. The strain tensor field is used to construct synthetic histograms of directional strains considering a cubic crystal structure (space group Im-3m).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. A total of $j = 1, ..., 134$ directions, $\hat{k}_j$, generated from a cubic iron (Fe) crystal structure are displayed by projecting the unit ball onto the $k_x - k_y$-plane (sample coordinate system). The sampling is seen to be both non-uniform as well as incomplete as a result of the crystal structure and experimental limitations.</figcaption></figure>

stress solutions [49]. We used iron (Fe) like materials parameters, $E = 200$ GPa, $\nu = 0.3$, and sampled a strain field at random on a spherical domain with radius 5 µm. The length-scale of the kernel was set to the radius of the sphere and the amplitude of the kernel was set to 60. An interpolation of the resulting strain tensor field can be visualised in Fig. 7. For strain tensor sampling ~50 000 spatial coordinates, $\mathbf{x}$, where selected over the spherical domain on an equidistant grid.

## 8.2. Crystal generation

With samples from the spatial strain tensor field defined, we are left with the generation of $\hat{k}$ before histograms of directional strains can be formed. To mimic a 3DXRD/HEDM experiment we generate $\hat{k}$ based on a crystal lattice structure. Adopting the cubic symmetry of iron (Fe) with unit cell parameters,

$$a = 2.8665 \AA, \quad b = 2.8665 \AA, \quad c = 2.8665 \AA, \quad \alpha = 90^{\circ}, \quad \beta = 90^{\circ}, \quad \gamma = 90^{\circ}, \qquad (61)$$

and a wavelength of 0.19 Å, we consider all reflections that fall on a 2048 × 2048 flat area detector with a 50 µm pixel size. The resulting diffraction events where confined to fall in the range $\theta = [5^{\circ}, 15^{\circ}]$. Without loss of generality, we set the crystal orientation as $U = I$ and solve the Laue equations (1) for $G$. The corresponding $\hat{k}$ directions where computed by normalising the $G$ vectors, which, by definition, are orthogonal to the diffracting lattice planes. To simulate the limitations imposed by the finite energy of the X-rays and the single axis of rotation we remove all $\hat{k}$ that fall within 10° of the rotation ±2-axis. This procedure generated a total $j = 1, ..., 132$ directions, $\hat{k}_j$. The spread of the generated directions on the unit ball can be viewed in Fig. 8.

## 8.3. Histogram formation

To form the final histograms of directional strains we must select a histogram bin size. The practical achievable resolution in terms of bin size depends on the instrument, material, geometry and detector parameters. In this work we select to use 6 bins spaced equidistantly throughout histogram formation. We select the range of the bins to cover fully the directional strains generated by multiplying through the sample strain tensors with the generated $\hat{k}$ directions. Additionally, for visual purposes, we pad the histogram with an additional bin on each side of this range to illustrate that no strain is falling outside our interval. Note that the proposed regression method can equally be deployed in the presence of fewer or more histogram bins. Naturally, the number of accessible bins will regulate what resolution can be expected of the recovered strain PDF.

Given the spatial strain tensor sample, $E(\mathbf{x}_1), E(\mathbf{x}_2), ..., E(\mathbf{x}_m)$, we generate samples of directional strains as

$$s_j = (s_1, s_2, ..., s_m) = (\bar{k}_j^T E(\mathbf{x}_1), \bar{k}_j^T E(\mathbf{x}_2), ..., \bar{k}_j^T E(\mathbf{x}_m)), \quad \forall j. \qquad (62)$$
<figure><img src="image_12.png" /><figcaption>Fig. 9. Synthetic histograms used as input to the derived regression framework. Each column in the image is a histogram of directional strain in the direction of a fixed $\hat{k}_j$. The pixel colour values represent the (noisy) probability that the directional strain will fall into the respective bin.</figcaption></figure>

For each j the corresponding sample of $s_j$ was used to form a histogram. The histograms were normalised to represent probability densities.

## 8.4. Noise model

Before deploying our regression framework, we introduce noise on the measured histogram bin values, $h_{ij}$, where $h_{ij}$ denotes the i:th bin in the j:th histogram. We add independent, identically distributed, Gaussian noise to all bins as

$$h_{ij} \leftarrow h_{ij} + e_{ij}, \quad e_{ij} \sim \mathcal{N}(0, \sigma_{ij}^{\text{noise}}), \quad \forall i, j, \qquad (63)$$

where

$$\sigma_{ij}^{\text{noise}} = 10^{-4} + h_{ij} / \text{SNR}_1. \qquad (64)$$

The signal to noise ratio was set to $SNR_1 = 50$ representing a few percent of noise per bin with the constant $10^{-4}$ ensuring that noise is added to bins featuring zero counts. We simulate outliers in the data by adding an additional larger Gaussian noise to a randomly (uniformly) selected subset as

$$h_{ij} \leftarrow h_{ij} + e_{ij}^{\text{outlier}}, \quad e_{ij}^{\text{outlier}} \sim \mathcal{N}(0, \sigma_{ij}^{\text{outlier}}), \forall i, j \in \mathcal{O}, \qquad (65)$$

where

$$\sigma_{ij}^{\text{outlier}} = 10\sigma_{ij}^{\text{noise}} \qquad (66)$$

and ~20% of the data histogram bins were selected for outlier corruption. The final step of our noise model is a re-normalisation, ensuring that the histogram bin data are positive. We subtract the minimal bin-count,

$$h_{ij} \leftarrow h_{ij} - \min_{i}(h_{ij}), \quad \forall j, \qquad (67)$$

and normalise,

$$h_{ij} \leftarrow h_{ij} \left( \sum_{i} w h_{ij} \right)^{-1} \quad \forall j. \qquad (68)$$

Collecting all generated histograms, one for each direction $\hat{k}_j$, the total available data can be visualised as an image. In Fig. 9 the noisy histogram data used for regression is displayed. Each pixel corresponds to a single histogram bin and each column to a single histogram. The colourmap reveals the bin-count as the probability of encountering a particular directional strain in the crystal.

## 8.5. Results

We deployed the generic radial basis regression framework derived in Section 6 to the synthetic histogram data described in Section 8. Additionally, using the closed form solution given in Eq. (58), we provide estimates based on a sparsely parameterised multivariate Gaussian strain PDF model, as described in Section 7.2.4.

### 8.5.1. Data prediction quality

To assess the quality of reconstruction we present a comparison between predicted and measured histograms of directional strains in Fig. 10. A good regression model should reject noise while at the same time predict directional strain histogram data in agreement with the measurements. The root mean squared error between predicted and measured data was defined as

$$\text{rmse} = \sqrt{\frac{\sum_i \sum_j (h_{ij}^{\text{predicted}} - h_{ij}^{\text{true}})^2}{8m}}, \qquad (69)$$

where $8m$ is the total number of measured bin values across all $m = 132$ histograms.
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Predicted histograms using a sparsely parameterised multivariate Gaussian model (right) and a generic radial basis regression model (left) compared to the ground truth noisy histograms (central column). The respective residual fields are shown in the outermost columns. Histogram bins that were randomly selected to be corrupted into outliers are marked in the residual plots with a black border.</figcaption></figure>

### 8.5.2. Convergence

The parameters for Algorithm 1 were selected as $B = J = 25$ and $(\sigma_1, \sigma_2, \sigma_3, \sigma_4, \sigma_5) = (10.0, 8.4, 7.1, 5.9, 5.0)$. The convergence of the algorithm is shown in Fig. 11.

### 8.5.3. Strain PDF reconstruction quality

To further assess our two models, in Fig. 12 we present a comparison between the predicted strain PDFs and the ground truth strain PDF (originating from the strain field in Fig. 7). Since the strain PDF is supported on 6-dimensions we present 2-dimensional projections onto the 15 unique Cartesian planes that can be formed by combining components of the strain tensor. As the directional strain histogram data resolution is limited, we must expect arbitrary oscillations in the reconstructed strain PDF to occur at length-scales smaller than the directional strain histogram bin width. Any meaningful comparison should, therefore, be made on a grid with similar resolution. To this end we have binned the projected strain PDFs into histograms in Fig. 12.

## 9. Discussion

Comparing the multivariate Gaussian prediction fit in Fig. 10 with the generic radial basis expansion model we observe a reduction in the residual error using the generic model. This is expected considering that the underlying strain PDF visualised in Fig. 12 is non-Gaussian. In fact, already from visual inspection of the histograms of directional strains, we may expect that the strain PDF is non-Gaussian. Nevertheless, despite the strain PDF being non-Gaussian, we observe from Fig. 12 that the Gaussian model captures multiple features of the strain PDF. Further inspection of Fig. 12 confirms that the full radial basis model offers additional improvements in terms of predicting the strain PDF shape. Considering the noise introduced into the data, as described in Section 8.4, these results imply that both regression models are robust in the presence of multiple outliers, limited and non-uniform $\hat{\kappa}$ sampling and elevated levels of Gaussian measurement noise.

In terms of compute time, using a non-optimised python implementation, the estimation for the sparsely parameterised multivariate Gaussian model was in the range of ~20 ms, while the stochastic radial basis placement algorithm finished in ~6 s. These results were achieved with a Dell XPS 15 9560 laptop using a single Intel(R) Core(TM) i7-7700HQ CPU @ 2.80 GHz. Low
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Convergence of Algorithm 1 applied to the data generated through the procedure described in Section 8. Each batch represents a fixed radial basis radius, $\sigma$. The number of basis functions used in the expansion of the strain PDF (right axis) is seen to increase with iteration. Likewise, the root mean squared error (rmse) is seen to decrease (left axis) with iteration. The minimal obtainable rmse is bounded by the standard deviation of the noise model described in Section 8.4.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Binned projections of the strain PDF onto the 15 possible unique Cartesian planes on $\mathbb{R}^6$ (the $E_2-E_1$-plane, the $E_3-E_1$-plane, ... etc.). The ground truth strain PDF (central column) originated from the spatial strain tensor field in Fig. 7. The reconstructed strain PDFs, using both the sparsely parameterised multivariate Gaussian model (right) and the generic radial basis regression method (left) are to be compared to the ground truth strain PDF (central column). The respective residual fields are shown in the outermost columns.</figcaption></figure>
compute times are important features of the proposed algorithms, effectively enabling the analysis of polycrystals with large number of grains (~1000) without the need for high performance computing resources. On the other hand, if an optimised multi-threaded implementation is pursued, these benchmarks are reasonable to permit online reconstruction at the same pace as data are being collected.

The final number of free parameters used by the radial basis regression is observed to be 94 × 6 = 564 from Fig. 11. We note that the total number of nonzero data bins in the measured histograms are 6 × 132 = 792, which is greater than the 564 parameters, avoiding over parametrisation of the model. In contrast to the 564 parameters used by the radial basis model, the sparsely parameterised multivariate Gaussian fit uses a mere total of 6 + 21 = 27 model parameters (estimating the mean and symmetric covariance respectively). This could motivate the use of the Gaussian model in situations when the available data are scarce.

## 10. Outlook

The convergence properties of the proposed Algorithm 1 are an ongoing research topic and future research should focus on parameter selection for B, J, σ and ξ. We believe the algorithm outlined in 1 has potential applications for a wider class of estimation problems, going beyond the estimation of strain PDFs, where a smooth, positive, simply connected scalar function needs to be approximated using a sparse basis. One such possible application is (X-ray or neutron) tomographic reconstruction, where a positive scalar attenuation function is to be estimated from noisy lower dimensional projections.

In terms of designing pilot experiments for the verification of the proposed regression methods, we suggest to perform time consuming s3DXRD in conjunction with fast full-field 3DXRD. Using the s3DXRD data the spatial strain tensor field can be estimated. By transforming the reconstructed spatial strain tensor field into a strain PDF and using the methods presented in this paper to recover the strain PDF from the full-field 3DXRD data independently, experimental validation could be achieved.

The research presented in this paper shows that approximation of the strain PDF, beyond a single mean value, is possible even with noisy distributions of directional strains. As such, we have provided a catalyst for pursuing algorithms that can optimally transform raw diffraction images into strain histograms. Since the number of achievable bins in such histograms is strongly related to the resolution of the detector, we anticipate the histogram resolution will increase with time, following the general trend in 3DXRD/HEDM type experiments. Alternatively, keeping the pixel density unchanged, the detector may be moved further away from the sample, leading to diffraction spot magnification. This is similar to the high angular resolution 3DXRD method, developed by Jakobsen et al. [50] (see also [51,52]).

The framework developed within this paper is not strictly limited to 3DXRD/HEDM applications. Other diffraction based techniques capable of producing histograms of directional strains, such that at least 2 bins can be resolved, will benefit from our results. As such we foresee a broader application of our findings in the context of non-destructive diffraction based materials research using both neutrons and X-rays.

## 11. Conclusions

In the context of high energy X-ray diffraction measurement from polycrystals, we have discussed the estimation of the strain tensor probability density function from noisy histograms of directional strains. Working in the limit of small strains, we derive a generic, iterative framework for strain PDF regression using a radial basis expansion as well as a sparsely parameterised multivariate Gaussian fit to which we provide a closed form solution. For the latter Gaussian class of strain tensor probability density functions, we give an exact description to the null-space of the inverse problem. The parametric description of the null-space provides insights into the strain PDF modes that cannot be accurately estimated from the diffraction data. Understanding these modes prevents erroneous conclusions from being drawn based on the data. The radial basis expansion method was shown to provide more accurate strain PDF estimations, at the expense of being mathematically involved, compared to the multivariate Gaussian fit which, on the other hand, represents a conceptually simple first order approximation to the strain PDF. We demonstrate our findings by numerical simulation considering a cubic iron (Fe) crystal in self-equilibrium. The simulation results imply that the proposed method performs well in the presence of Gaussian noise, outliers and non-uniform strain sampling. Overall, our findings represent an upgrade to the full-field 3DXRD/HEDM microscope methodology with the potential of broader application to other diffraction based methods.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

Data will be made available on request.

## Acknowledgements

This work was funded by Vetenskapsrådet - Röntgen Ångström Cluster, project no. 2017-06719. The authors are grateful for the grants provided by The Fund of the Walter Gyllenberg Foundation, Sweden (no. 43096) as well as the Travel Grants for Young Researchers (no. 43000) granted by the Royal Physiographic Society of Lund, Sweden. These grants made necessary international collaboration between the authors possible.
<figure><img src="image_16.png" /><figcaption>Fig. 13. Connection between single crystal deformation configurations associated with a multiplicative decomposition of the deformation gradient tensor, $F = F_e F_p$. The physical process that yields diffraction from a crystal configuration is denoted by the operator $\mathcal{G}$. The fact that the diffraction data in 3DXRD/HEDM is independent of the plastic deformation, $F_p$, is illustrated to the right and left.</figcaption></figure>

## Appendix A. Deformation configurations

In Fig. 13 we illustrate that the deformed unit cell matrix is reached from the undeformed state by applying the deformation gradient tensor as

$$C_{cur} = F C_0 = F_e F_p C_0. \qquad (70)$$

The set of diffraction vectors measured from the current deformed state is denoted $\mathcal{G}[C_{cur}] = \{G^1, G^2, \dots, G^m\}$ and could equally have originated from a crystal with unit cell matrix

$$C = F_e C_0, \qquad (71)$$

as indicated to the right of Fig. 13. Similarly the unit cell matrices $C_0$ and $C_p = F_p C_0$ are diffraction equivalent, which is illustrated to the left in Fig. 13. We note that while we call $F_e$ an elastic deformation and $F_p$ a plastic deformation the technical definition is here that $F_e$ is any deformation such that

$$\mathcal{G}[F_e C_0] \neq \mathcal{G}[C_0], \qquad (72)$$

while $F_p$ is any deformation such that

$$\mathcal{G}[F_p C_0] = \mathcal{G}[C_0]. \qquad (73)$$

This means that, in our formulation, $F_e$ can embody deformations that persist in the crystal when the boundary tractions are revoked.

## Appendix B. Covariance perturbations

We seek solutions to Eq. (52) which we reiterate here as

$$\bar{\kappa}_j^T N \bar{\kappa}_j = 0, \quad \forall j. \qquad (74)$$

We divide our analysis into two cases; one where $N$ is symmetric-positive-definite and one where $N$ is symmetric-non-positive-definite.

### B.1. Symmetric positive-definite covariance perturbations

Let us start by analysing the case when $z^T N z \ge 0$, which can be parameterised by the matrix outer product $N = AA^T$, $A \in \mathbb{R}^{6 \times k}$. Adopting our previous matrix format we have

$$K A A^T K = 0, \qquad (75)$$

and it follows from the quadratic format that we must require

$$K A = 0, \qquad (76)$$

i.e., we must require $K$ to have a null-space. This is in contradiction with our previous results in Section 7.1 were we showed that $K$ can be fully ranked in general. Thus, no positive semi-definite perturbations to the covariance of the strain PDF exists in general, such that the measurements remain invariant.
It follows that if we can give example of a rank 21 − 6 = 15 matrix, R , such that

$$
(\mathbf{R} + \mathbf{U}^T)\mathbf{v}_j = \mathbf{R}^T\mathbf{v}_j, \qquad (83)
$$

the column space of U depletes the sought null-space. Numerous examples of such R can be generated with computer aid. One
possible selection, with all kappa forming angles >26° to the z-axis, is, for example,

$$
\begin{align}
\hat{\mathbf{k}}_1 &= \frac{1}{\sqrt{3}} \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}, & \hat{\mathbf{k}}_2 &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix}, & \hat{\mathbf{k}}_3 &= \frac{1}{\sqrt{3}} \begin{bmatrix} 1 \\ 1 \\ -1 \end{bmatrix}, & \hat{\mathbf{k}}_4 &= \frac{1}{\sqrt{6}} \begin{bmatrix} 1 \\ 1 \\ 2 \end{bmatrix}, & \hat{\mathbf{k}}_5 &= \frac{1}{\sqrt{6}} \begin{bmatrix} 1 \\ 1 \\ -2 \end{bmatrix}, \nonumber \\
\hat{\mathbf{k}}_6 &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 0 \\ 1 \end{bmatrix}, & \hat{\mathbf{k}}_7 &= \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, & \hat{\mathbf{k}}_8 &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 0 \\ -1 \end{bmatrix}, & \hat{\mathbf{k}}_9 &= \frac{1}{\sqrt{5}} \begin{bmatrix} 1 \\ 0 \\ 2 \end{bmatrix}, & \hat{\mathbf{k}}_{10} &= \frac{1}{\sqrt{3}} \begin{bmatrix} 1 \\ -1 \\ 1 \end{bmatrix}, \tag{84} \\
\hat{\mathbf{k}}_{11} &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -1 \\ 0 \end{bmatrix}, & \hat{\mathbf{k}}_{12} &= \frac{1}{\sqrt{3}} \begin{bmatrix} 1 \\ -1 \\ -1 \end{bmatrix}, & \hat{\mathbf{k}}_{13} &= \frac{1}{\sqrt{6}} \begin{bmatrix} 1 \\ 2 \\ 1 \end{bmatrix}, & \hat{\mathbf{k}}_{14} &= \frac{1}{\sqrt{5}} \begin{bmatrix} 1 \\ 2 \\ 0 \end{bmatrix}, & \hat{\mathbf{k}}_{15} &= \frac{1}{\sqrt{6}} \begin{bmatrix} 1 \\ -2 \\ 1 \end{bmatrix}. \nonumber
\end{align}
$$

Appendix C. Basis integration

We seek to solve the integral given in Eq. (36) which we reiterate here as

$$
\int_{s_j^* - w/2}^{s_j^* + w/2} \int_{\epsilon \in \Pi(t, \hat{n}_j)} \exp\left(-\frac{(\epsilon - p_i)^T (\epsilon - p_i)}{2\sigma^2}\right) d\epsilon dt. \quad (85)
$$

We denote the inner integral I_ε and the outer integral I_l and initiate our analysis with I_ε. To solve the integral we start by parameterising all points on Π by the range of P ∈ ℝ^{6×5} such that the columns of P forms a basis for Π. For reasons that will become apparent later we select the parametrisation

$$
\epsilon = \epsilon(\mathbf{y}) = \mathbf{P}\mathbf{y} + t\hat{\mathbf{n}}_j + \mathbf{P}\mathbf{P}^T p_i \in \Pi(t, \hat{\mathbf{n}}_j), \quad \mathbf{y} \in \mathbb{R}^5. \tag{86}
$$

This simplifies our problem to an unbounded integral over ℝ⁵ as

$$
I_{\epsilon} = \int_{\mathbb{R}^5} \exp\left(-\frac{(\epsilon(\mathbf{y}) - p_i)^T (\epsilon(\mathbf{y}) - p_i)}{2\sigma^2}\right) d\mathbf{y}. \qquad (87)
$$

Without loss of generality, we now select the columns of P to be orthonormal, such that P^T P = I. Expanding the square in (87) we find that

$$
I_e = \int_{\mathbb{R}^5} \exp\left(-\frac{(y^T y - 2t p_i^T \hat{n}_j + t^2 - p_i^T P P^T p_i + p_i^T p_i)}{2\sigma^2}\right) dy, \quad (88)
$$

where it was used that P ˆn_j = 0. Introducing the scalar minimum squared distance, d^2, between Π and p_i as

$$
d^2 = (p_i - \mathbf{P}\mathbf{P}^T p_i - t\hat{n}_j)^T (p_i - \mathbf{P}\mathbf{P}^T p_i - t\hat{n}_j) = p_i^T p_i - 2t p_i^T \hat{n}_j + t^2 - p_i^T \mathbf{P}\mathbf{P}^T p_i, \quad (89)
$$

we find by insertion into (87) that

$$
I_e = \int_{\mathbb{R}^5} \exp\left(-\frac{(y^T y + d^2)}{2\sigma^2}\right) dy. \qquad (90)
$$

The benefits of the specific parametric selection of Π is now evident from the cancelled terms. Factoring out the part that is
independent of y, we have

$$
I_e = \exp\left(-\frac{d^2}{2\sigma^2}\right) \int_{\mathbb{R}^5} \exp\left(-\frac{\mathbf{y}^T \mathbf{y}}{2\sigma^2}\right) d\mathbf{y}, \qquad (91)
$$

where the rightmost integral features an non-normalised Gaussian with covariance Iσ². We solve this Gaussian integral by
normalisation as

$$
\int_{\mathbb{R}^5} \exp\left(-\frac{\mathbf{y}^T \mathbf{y}}{2\sigma^2}\right) d\mathbf{y} = \sigma^5 \sqrt{(2\pi)^5} \int_{\mathbb{R}^5} \frac{1}{\sigma^5 \sqrt{(2\pi)^5}} \exp\left(-\frac{\mathbf{y}^T \mathbf{y}}{2\sigma^2}\right) d\mathbf{y} = \sigma^5 \sqrt{(2\pi)^5}. \quad (92)
$$

Inserting (92) into (91) we find

$$
I_e = \sigma^5 \sqrt{(2\pi)^5} \exp\left(-\frac{d^2}{2\sigma^2}\right). \qquad (93)
$$

We proceed to solve for the outer integral in (36). Using the result in (93) we have

$$
I_t = \int_{s_j^* - w/2}^{s_j^* + w/2} I_e dt = \sigma^5 \sqrt{(2\pi)^5} \int_{s_j^* - w/2}^{s_j^* + w/2} \exp\left(-\frac{d^2}{2\sigma^2}\right) dt. \quad (94)
$$
Since the distance *d* is linear in *t* and the increments, *dd* = *dt*, are one to one, we now find a bounded Gaussian integral that can be solved using the error function, erf, as

$$I_t = \sigma^5 \sqrt{(2\pi)^5} \left| \int_a^b \exp\left(-\frac{d^2}{2\sigma^2}\right) dd \right| = 4\pi^3 \sigma^6 \left| \operatorname{erf}\left(\frac{b}{\sigma\sqrt{2}}\right) - \operatorname{erf}\left(\frac{a}{\sigma\sqrt{2}}\right) \right|, \quad (95)$$

where the transformed bounds, *a* and *b*, can be computed as $a = (s_j^* - w/2)p_i^T \hat{n}_j$ and $b = (s_j^* + w/2)p_i^T \hat{n}_j$, respectively. This concludes our integration analysis.

## References

[1] H. Poulsen, Multi scale hard x-ray microscopy, Curr. Opin. Solid State Mater. Sci. 24 (2) (2020) 100820.

[2] V. Santoro, K.H. Andersen, M. Bernasconi, M. Bertelsen, Y.B. Ier, D. Campi, V. Czamler, D.D.D. Julio, E. Diane, K. Dunne, P. Fierlinger, A. Gaye, G. Gorini, C. Happe, T. Kittelmann, E.B. Klinkby, Z. Kokai, R. Kolevatov, B. Lauritzen, R. Linander, J.I.M. Damian, B. Meirose, F. Mezei, D. Milstead, G. Muhrer, K. Ramic, B. Rataj, N. Rizzi, S. Samothrakitis, J.R. Selknaes, S. Silverstein, M. Strobl, M. Strothmann, A. Takibayev, R. Wagner, P. Willendrup, S.C. Yiu, L. Zanini, O. Zimmer, Development of a high intensity neutron source at the European Spallation source: The HighNESS project, 2022.

[3] A. Cereser, M. Strobl, S.A. Hall, A. Steuwer, R. Kiyanagi, A.S. Tremsin, E.B. Knudsen, T. Shinohara, P.K. Willendrup, A.B. da Silva Fanta, S. Iyengar, P.M. Larsen, T. Hanashima, T. Moyoshi, P.M. Kadletz, P. Krooß, T. Niendorf, M. Sales, W.W. Schmahl, S. Schmidt, Time-of-flight three dimensional neutron diffraction in transmission mode for mapping crystal grain structures, Sci. Rep. 7 (1) (2017) 9561.

[4] Y. Hayashi, D. Setoyama, Y. Hirose, T. Yoshida, H. Kimura, Intragranular three-dimensional stress tensor fields in plastically deformed polycrystals, Science 366 (6472) (2019) 1492-1496.

[5] E. Demir, J.-S. Park, M. Miller, P. Dawson, A computational framework for evaluating residual stress distributions from diffraction-based lattice strain data, Comput. Methods Appl. Mech. Engrg. 265 (2013) 120-135.

[6] R. Quey, L. Renversade, Optimal polyhedral description of 3D polycrystals: Method and application to statistical and synchrotron X-ray diffraction data, Comput. Methods Appl. Mech. Engrg. 330 (2018) 308-333.

[7] M.M. Thakur, N.A. Henningsson, J. Engqvist, P.-O. Autran, J.P. Wright, R.C. Hurley, On mesoscale modeling of concrete: Role of heterogeneities on local stresses, strains, and representative volume element, Cem. Concr. Res. 163 (2023) 107031.

[8] J.M. Dake, J. Oddershede, H.O. Sørensen, T. Werz, J.C. Shatto, K. Uesugi, S. Schmidt, C.E. Krill, Direct observation of grain rotations during coarsening of a semisolid Al-Cu alloy, Proc. Natl. Acad. Sci. 113 (41) (2016) E5998-E6006.

[9] H. Poulsen, 3DXRD - A New Probe for Materials Science (Ph.D. thesis), Risø National Laboratory, 2004.

[10] J.V. Bernier, R.M. Suter, A.D. Rollett, J.D. Almer, High-energy X-Ray diffraction microscopy in materials science, Annu. Rev. Mater. Res. 50 (1) (2020) 395-436.

[11] E.M. Lauridsen, S. Schmidt, R.M. Suter, H.F. Poulsen, Tracking: a method for structural characterization of grains in powders or polycrystals, J. Appl. Crystallogr. 34 (6) (2001) 744-750.

[12] H. Sørensen, S. Schmidt, J. Wright, G. Vaughan, S. Techert, E. Garman, J. Oddershede, J. Davaasambuu, K. Paithankar, C. Gundlach, H. Poulsen, Multigrain crystallography, Z. Kristallogr.-Cryst. Mater. 227 (1) (2012) 63-78.

[13] H. Sharma, R.M. Huizenga, S.E. Offerman, A fast methodology to determine the characteristics of thousands of grains using three-dimensional X-ray diffraction. I. Overlapping diffraction peaks and parameters of the experimental setup, J. Appl. Crystallogr. 45 (4) (2012) 693-704.

[14] J.M. Wozniak, H. Sharma, T.G. Armstrong, M. Wilde, J.D. Almer, I. Foster, Big data staging with MPI-IO for interactive X-ray science, in: 2014 IEEE/ACM International Symposium on Big Data Computing, 2014, pp. 26-34.

[15] L. Margulies, T. Lorentzen, H. Poulsen, T. Leffers, Strain tensor development in a single grain in the bulk of a polycrystal under loading, Acta Mater. 50 (7) (2002) 1771-1779.

[16] H.F. Poulsen, S.F. Nielsen, E.M. Lauridsen, S. Schmidt, R.M. Suter, U. Lienert, L. Margulies, T. Lorentzen, D. Juul Jensen, Three-dimensional maps of grain boundaries and the stress state of individual grains in polycrystals and powders, J. Appl. Crystallogr. 34 (6) (2001) 751-756.

[17] J. Oddershede, S. Schmidt, H.F. Poulsen, H.O. Sørensen, J. Wright, W. Reimers, Determining grain resolved stresses in polycrystalline materials using three-dimensional X-ray diffraction, J. Appl. Crystallogr. 43 (3) (2010) 539-549.

[18] J.V. Bernier, N.R. Barton, U. Lienert, M.P. Miller, Far-field high-energy diffraction microscopy: a tool for intergranular orientation and strain analysis, J. Strain Anal. Eng. Des. 46 (7) (2011) 527-547.

[19] Y. Hayashi, D. Setoyama, Y. Seno, Scanning three-dimensional X-Ray diffraction microscopy with a high-energy microbeam at spring-8, Mater. Sci. Forum 905 (2017) 157-164.

[20] W. Li, H. Sharma, K. Peter, S. Ravi, H. Sehitoglu, A. Bucsek, Resolving intragranular stress fields in plastically deformed titanium using point-focused high-energy diffraction microscopy, J. Mater. Res. 38 (1) (2023) 165-178.

[21] J. Hendriks, C. Wensrich, A. Wills, V. Luzin, A. Gregg, Robust inference of two-dimensional strain fields from diffraction-based measurements, Nucl. Instrum. Methods Phys. Res. B 444 (2019) 80-90.

[22] N.A. Henningsson, S.A. Hall, J.P. Wright, J. Hektor, Reconstructing intragranular strain fields in polycrystalline materials from scanning 3Drd data, J. Appl. Crystallogr. 53 (2) (2020) 314-325.

[23] J.N. Hendriks, C.M. Wensrich, A. Wills, A Bayesian approach to triaxial strain tomography from high-energy X-ray diffraction, Strain 56 (3) (2020) e12341, e12341 10.1111/str.12341.

[24] A. Henningsson, J. Hendriks, Intragranular strain estimation in far-field scanning X-ray diffraction using a Gaussian process, J. Appl. Crystallogr. 54 (4) (2021) 1057-1070.

[25] Y. Hayashi, Y. Hirose, Y. Seno, Polycrystal orientation mapping using scanning three-dimensional X-ray diffraction microscopy, J. Appl. Crystallogr. 48 (4) (2015) 1094-1101.

[26] J. Hektor, S.A. Hall, N.A. Henningsson, J. Engqvist, M. Ristinmaa, F. Lenrick, J.P. Wright, Scanning 3Drd measurement of grain growth, stress, and formation of Cu6Sn5 around a tin whisker during heat treatment, Materials 12 (3) (2019).

[27] S. Schmidt, GrainSpotter: a fast and robust polycrystalline indexing algorithm, J. Appl. Crystallogr. 47 (1) (2014) 276-284.

[28] S. West, S. Schmidt, H. Sørensen, G. Winther, H. Poulsen, L. Margulies, C. Gundlach, D. Juul Jensen, Direct non-destructive observation of bulk nucleation in 30% deformed aluminum, Scr. Mater. 61 (9) (2009) 875-878.

[29] H.F. Poulsen, U. Lienert, W. Pantleon, Characterisation of orientation distributions of individual grains within deformed metals, Mater. Sci. Technol. 21 (12) (2005) 1397-1400.

[30] P.C. Hansen, H.O. Sørensen, Z. Sükösd, H.F. Poulsen, Reconstruction of single-grain orientation distribution functions for crystalline materials, SIAM J. Imaging Sci. 2 (2) (2009) 593-613.

[31] N.R. Barton, J.V. Bernier, A method for intragranular orientation and lattice strain distribution determination, J. Appl. Crystallogr. 45 (6) (2012) 1145-1155.

[32] H. Behnken, Strain-function method for the direct evaluation of intergranular strains and stresses, Phys. Status Solidi (A) 177 (2) (2000) 401-418.
[33] J.V. Bernier, M.P. Miller, A direct method for the determination of the mean orientation-dependent elastic strains and stresses in polycrystalline materials from strain pole figures, J. Appl. Crystallogr. 39 (3) (2006) 358-368.

[34] D. Boyce, P. Shade, W. Musinski, M. Obstalecki, D. Pagan, J. Bernier, T. Turner, Estimation of anisotropic elastic moduli from high energy X-Ray data and finite element simulations, Materialia 12 (2020) 100795.

[35] W. Pantleon, H. Poulsen, J. Almer, U. Lienert, In situ X-ray peak shape analysis of embedded individual grains during plastic deformation of metals, Mater. Sci. Eng. A 387-389 (2004) 339-342, 13th International Conference on the Strength of Materials.

[36] K. Fogarty, E. Ametova, G. Burca, A.M. Korsunsky, S. Schmidt, P.J. Withers, W.R.B. Lionheart, Recovering the second moment of the strain distribution from neutron bragg edge data, Appl. Phys. Lett. 120 (16) (2022) 164102.

[37] J. Als-Nielsen, D. McMorrow, Elements of Modern X-ray Physics, John Wiley & Sons, Ltd, 2011.

[38] S.L. Wong, J.-S. Park, M.P. Miller, P.R. Dawson, A framework for generating synthetic diffraction images from deforming polycrystals using crystal-based finite element formulations, Comput. Mater. Sci. 77 (2013) 456-466.

[39] P.R. Dawson, M.P. Miller, A virtual diffractometer for creating synthetic HEDM images of tessellated and meshed finite element polycrystals, 2023.

[40] E. Kröner, Allgemeine kontinuumstheorie der versetzungen und eigenspannungen, Arch. Ration. Mech. Anal. 4 (1959) 273-334.

[41] E.H. Lee, D.T. Liu, Finite-strain elastic—Plastic theory with application to plane-wave analysis, J. Appl. Phys. 38 (1) (1967) 19-27.

[42] R.J. Clifton, On the equivalence of FeFp and FpFe, J. Appl. Mech. 39 (1) (1972) 287-289.

[43] J. Clayton, D. McDowell, A multiscale multiplicative decomposition for elastoplasticity of polycrystals, Int. J. Plast. 19 (9) (2003) 1401-1444.

[44] Y. Jiao, J. Fish, On the equivalence between the multiplicative hyper-elasto-plasticity and the additive hypo-elasto-plasticity based on the modified kinetic logarithmic stress rate, Comput. Methods Appl. Mech. Engrg. 340 (2018) 824-863.

[45] N.S. Ottosen, M. Ristinmaa, The Mechanics of Constitutive Modeling, Elsevier, 2005.

[46] C.L. Lawson, R.J. Hanson, Solving Least Squares Problems, SIAM, 1995.

[47] P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S.J. van der Walt, M. Brett, J. Wilson, K.J. Millman, N. Mayorov, A.R.J. Nelson, E. Jones, R. Kern, E. Larson, C.J. Carey, İ. Polat, Y. Feng, E.W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E.A. Quintero, C.R. Harris, A.M. Archibald, A.H. Ribeiro, F. Pedregosa, P. van Mulbregt, SciPy 1.0 Contributors, Scipy 1.0: Fundamental algorithms for scientific computing in python, Nature Methods 17 (2020) 261-272.

[48] A.G. Wills, T.B. Schön, Sequential Monte Carlo: A unified review, Annu. Rev. Control Robot. Auton. Syst. 6 (1) (2023) null.

[49] R. Rostamian, The completeness of Maxwell's stress function representation, J. Elasticity 9 (4) (1979) 349-356.

[50] B. Jakobsen, H.F. Poulsen, U. Lienert, J. Almer, S.D. Shastri, H.O. Sørensen, C. Gundlach, W. Pantleon, Formation and subdivision of deformation structures during plastic deformation, Science 312 (5775) (2006) 889-892.

[51] T. Ungár, G. Ribárik, L. Balogh, A.A. Salem, S.L. Semiatin, G.B. Vaughan, Burgers vector population, dislocation types and dislocation densities in single grains extracted from a polycrystalline commercial-purity Ti specimen by X-ray line-profile analysis, Scr. Mater. 63 (1) (2010) 69-72.

[52] C. Nisr, G. Ribárik, T. Ungár, G.B.M. Vaughan, P. Cordier, S. Merkel, High resolution three-dimensional X-ray diffraction study of dislocations in grains of MgGeO3 post-perovskite at 90 GPa, J. Geophys. Res.: Solid Earth 117 (B3) (2012).