Structure tensors

Structure tensors ${\bf a}^{(k)}$ are defined as the mean of the $k$-th repeated outer product of a crystallographic axis with itself, say $\bf n$:

\[ {\bf a}^{(k)} := \langle {\bf n}^k \rangle. \]

In the case of a discrete ensemble of $N$ grains with slip plane normals ${\bf n}_i$ for $i=1,\cdots,N$, the structure tensors are

\[ {\bf a}^{(k)} = \dfrac{\sum_i w_i {\bf n}_i^k}{\sum_i w_i}, \]

where $w_i$ is the weight attributed to the $i$-th grain, typically its size.

Alternatively, if the distribution function $n(\theta,\phi)$ of ${\bf n}$-axes is known, the structure tensors are

\[ {\bf a}^{(k)} = \dfrac{\int_{S^2} n \hat{{\bf r}}^k \, \mathrm{d}\Omega}{\int_{S^2} n \, \mathrm{d}\Omega} , \]

where $\mathrm{d}\Omega = \sin(\theta) \mathrm{d}\theta \mathrm{d}\phi$ is the infinitesimal solid angle and $\hat{{\bf r}}(\theta,\phi)$ is the radial unit vector.

Principal frame

Since $n(\theta,\phi)$ and $b(\theta,\phi)$ are antipodally symmetric, odd moments (odd $k$) vanish identically. Hence, ${\bf a}^{(2)}(n)$ and ${\bf a}^{(2)}(b)$ measure the variance of $\bf n$- and $\bf b$-axes, respectively, around the three coordinate axes.

Posing ${\bf a}^{(2)}$ in its principal frame

\[ {\bf a}^{(2)} = \left[\begin{matrix} \lambda_1 & 0 & 0\\ 0 & \lambda_2 & 0\\ 0 & 0 & \lambda_3\\ \end{matrix}\right] \]

therefore has a similar interpretation as in PCA: the first principal component (eigenvector ${\bf m}_1$) is the direction that maximizes the variance (eigenvalue $\lambda_1$) of the projected data (red curve), the second component is the direction orthogonal to the first component that maximizes the variance of the projected data, and so on with the third component.

Convert to spectral

Converting between spectral and tensorial representations is a linear problem in the sense that

\[ {\bf a}^{(k)} = {\bf f}(\hat{n}_2^{m}, \hat{n}_4^{m}, \cdots, \hat{n}_k^{m}) , \qquad\text{(for all $m$)} \]

where ${\bf f}$ is linear in its arguments and $\hat{n}_l^m = n_l^m/n_0^0$.

In the case of ${\bf a}^{(2)}$, the relation is simple:

\[ {\bf a}^{(2)} = \frac{{\bf I}}{3} + \sqrt{\frac{2}{15}} \left[\begin{matrix} \operatorname{Re}[\hat{n}_2^2] - \dfrac{1}{2}\sqrt{\dfrac{2}{3}} \hat{n}_2^0 & -\operatorname{Im}[\hat{n}_2^2] & -\operatorname{Re}[\hat{n}_2^1] \\ & -\operatorname{Re}[\hat{n}_2^2] - \dfrac{1}{2}\sqrt{\dfrac{2}{3}} \hat{n}_2^0 & \operatorname{Im}[\hat{n}_2^1] \\ \mathrm{sym.} & & \sqrt{\dfrac{2}{3}} \hat{n}_2^0 \end{matrix}\right] , \]

but for higher-order structure tensors the expressions are long (not shown). The above applies to $b(\theta,\phi)$ as well.

The following code example shows how to convert between the two CPO representations:

import numpy as np
from specfabpy import specfab as sf
lm, nlm_len = sf.init(8)
nlm = np.zeros((nlm_len), dtype=np.complex64) # array of expansion coefficients

### a2 to nlm
a2 = np.diag([0.0,0.25,0.75]) # arbitrary second-order structure tensor
nlm[:sf.L2len] = sf.a2_to_nlm(a2) # l<=2 expansion coefficients
a2 = sf.a2(nlm) # nlm back to a2
print('a2 is: ', a2)

### a4 to nlm
p = np.array([0,0,1]) # unidirectional CPO
a4 = np.einsum('i,j,k,l', p,p,p,p) # a4 for ODF = deltafunc(r-p) 
nlm[:sf.L4len] = sf.a4_to_nlm(a4) # l<=4 expansion coefficients
a4 = sf.a4(nlm) # nlm back to a4 
print('a4 is: ', a4)

### a6 to nlm
p = np.array([0,0,1]) # unidirectional CPO
a6 = np.einsum('i,j,k,l,m,n', p,p,p,p,p,p) # a6 for ODF = deltafunc(r-p) 
nlm[:sf.L6len] = sf.a6_to_nlm(a6) # l<=6 expansion coefficients
a6 = sf.a6(nlm) # nlm back to a6
print('a6 is: ', a6)

Construct from measurements

The harmonic expansion coefficients of any CPO can be determined from discrete measurements of its crystallographic axes. This requires constructing the corresponding structure tensors (of each crystallographic axis), from which the expansion coefficients may be derived.

In the case of glacier ice where ${\bf n} = {\bf c}$, this can be done as follows:

import numpy as np
from specfabpy import specfab as sf
lm, nlm_len = sf.init(8) 

### Replace with your own array/list of measured c-axes
# caxes = [[c1x,c1y,c1z], [c2x,c2y,c2z], ...] 

### Determine sixth-order structure tensor, a6
a6 = np.zeros((3,3,3,3,3,3))
for c in caxes
    a6 += np.einsum('i,j,k,l,m,n', c,c,c,c,c,c) # sixth outer product of c-axis with itself
a6 /= len(caxes) # normalize by number of c-axes (grains) if grain sizes are not known

### Determine spectral expansion coefficients
nlm = np.zeros((nlm_len), dtype=np.complex64) # array of expansion coefficients
nlm[:sf.L6len] = sf.a6_to_nlm(a6) # l<=6 expansion coefficients

Note that constructing ${\bf a}^{(6)}$ is to be preferred over ${\bf a}^{(4)}$ and ${\bf a}^{(2)}$, since it contains more information on the fine-scale structure of the distribution: ${\bf a}^{(6)}$ encodes information about the expansion coefficients $l\leq 6$, ${\bf a}^{(4)}$ about $l\leq 4$, and ${\bf a}^{(2)}$ about $l\leq 2$.