Tensorial harmonic bases of arbitrary order with applications in elasticity,elastoviscoplasticity and texture-based modeling

Abstract

In continuum mechanics, one regularly encounters higher-order tensors that require tensorial bases for theoretical or numerical calculations. By using results from $\mathrm{SO}\left(3\right)$ representation theory, we present a method to derive tensorial harmonic bases which unifies existing approaches in one framework. Applying this convention to symmetric second-order tensors results in a second-order harmonic basis which, for example, simplifies the depiction and numerical handling of stiffness tensors for most common material symmetries and reduces the computational effort involved in implementations of incompressible elastoviscoplastic material laws. The same convention can be applied to texture analysis to yield tensorial texture coefficients for both polycrystals and fiber-reinforced composites. Rotations of higher-order tensors are particularly efficient in harmonic bases as both material and index symmetries can be exploited. While special attention is given here to the examples of small-strain material laws and texture analysis, the framework is entirely general and can be used to simplify calculations in other physical contexts involving tensors of arbitrary order and symmetry. A Python implementation of the harmonic basis convention defined in this work is available as a Supplementary Material.

1. Introduction

Continuum mechanics deals with tensorial quantities in all subareas, i.e., kinematics, balance equations, and constitutive theory. The tensorial quanties necessitate a higher-order tensor algebra; however, matrix–vector formulations are more accessible both numerically and didactically.

Depicting higher-order tensors as vectors is as simple as picking a set of higher-order basis tensors, but the choice is not obvious. The oldest and most common approach dates back to Voigt [1]. The Voigt notation renders symmetric second-order tensors as six-dimensional component vectors by defining six orthogonal second-order basis tensors, not all of which are normalized. Therefore, co- and contravariant tensors need to be distinguished, leading to different notations for stress and strain tensors. Although this non-normalized basis requires special care [2], it is commonly used in commercial computer codes such as Abaqus, Ansys, and LS-DYNA.

The Voigt basis can be normalized, resulting in a factor of $\sqrt{2}$ for each of the shear basis tensors. Originally proposed by Voigt [1], the normalized Voigt basis is often called Mandel basis [3]. As it is an orthonormal basis, results from elementary linear algebra can be applied to Mandel basis component vectors and matrices in a straightforward manner.

While Voigt’s original decomposition into normal and shear tensors derives naturally from the standard ${\mathit{e}}_i$ basis of the underlying vector space, the normal tensors can be further split into a spherical part and two deviator parts. The deviatoric basis is due to Kocks et al. [4], who denoted it the natural basis due to its particular usefulness in descriptions of cubic polycrystals. Kocks’s original non-orthonormalized basis was followed by orthonormalized versions by, among others, Lequeu [5] and Van Houtte [6]. A lucid description as an eigentensor basis of cubic stiffness tensors can be found in Kocks et al.’s study [7], Chapter 7.

Separately from stress and strain bases, texture description has undergone a number of evolutions. Roe [8] and Bunge [9] are generally credited with texture coefficient descriptions of the orientation distribution function based on spherical harmonics series. These scalar coefficients motivated a coordinate-independent tensorial depiction, introduced by Adams et al. [10]. In actual calculation, these tensorial descriptions still necessitate a choice of basis tensors [11]. Man and Du [12] recently proposed an orthonormal basis of harmonic tensor spaces which uses results from complexified tensor algebra.

In this work, we define a harmonic basis convention for the space of real three-dimensional vectors and related tensor spaces of various symmetries. Section 2 recapitulates relevant parts of the representation theory of $\mathrm{SO}\left(3\right)$ and defines a harmonic basis convention for arbitrary tensor spaces. Section 3 applies this convention to second-order symmetric tensors, recovering an orthonormalized Kocks-type natural basis. This basis is compared to the Mandel basis for common material symmetries in linear elasticity as well as for anisotropic elastoplasticity models. In Section 4, higher-order harmonic bases are used for texture analysis in polycrystalline and fiber-reinforced microstructures.

2. Constructing harmonic bases

2.1. The rotation group $\mathrm{SO}\left(3\right)$

We recapitulate a few concepts from the representation theory of groups with a particular focus on the rotation group $\mathrm{SO}\left(3\right)$ . Of particular interest is the action of $\mathrm{SO}\left(3\right)$ on the space of three-dimensional real vectors ${\mathbb{R}}^3$ endowed with the usual scalar product, which will be denoted as $\mathcal{R}$ to reduce index confusion. In the context of this work, an element $Q$ of $\mathrm{SO}\left(3\right)$ is a linear map on $\mathcal{R}$ such that for $\mathit{v}\in \mathcal{R}$ ,

$$Q\left(\mathit{v}\right)=\mathit{Qv}=Q_{\mathrm{ij}}{v}_j{\mathit{e}}_i.$$

(1)

In other words, the action of $Q\in \mathrm{SO}\left(3\right)$ on $\mathcal{R}$ has a specific second-order tensor representation $\mathit{Q}=Q_{\mathrm{ij}}{\mathit{e}}_i\otimes {\mathit{e}}_j$ . This representation can be used to calculate the action of $\mathrm{SO}\left(3\right)$ on an $n$ th order tensor ${\mathbb{V}}^n\in {\mathcal{R}}^{\otimes n}$ via the Rayleigh product defined as

$$\begin{array}{c}\multicolumn{1}{c}{Q\left({\mathbb{V}}^n\right)=\mathit{Q}★{\mathbb{V}}^n}\\ \multicolumn{1}{c}{=Q_{\mathrm{ai}}{Q}_{\mathrm{bj}}\dots Q_{\mathrm{ck}}{V}_{\mathrm{ij}\dots k}{\mathit{e}}_a\otimes {\mathit{e}}_b\otimes \dots \otimes {\mathit{e}}_c.}\end{array}$$

(2)

After algebraic manipulation, we find the $2n$ -order tensor representation of the action of $Q$ on the higher-order tensor space,

$$\begin{array}{c}\multicolumn{1}{c}{Q\left({\mathbb{V}}^n\right)={\mathit{Q}}^{★n}\left[{\mathbb{V}}^n\right]}\\ \multicolumn{1}{c}{=(Q_{\mathrm{ai}}\dots Q_{\mathrm{bj}}{\mathit{e}}_a\otimes \dots {\mathit{e}}_b\otimes {\mathit{e}}_i\otimes \dots {\mathit{e}}_k)\left[{\mathbb{V}}^n\right],}\end{array}$$

(3)

where we call ${\mathit{Q}}^{★n}$ the $n$ th Rayleigh power of $\mathit{Q}$ . The properties of the Rayleigh power were previously discussed in a continuum mechanics context by Zheng and Spencer [13], who denote it as a Kronecker exponent ${〈\mathit{Q}〉}_n$ . Being a linear map from ${\mathcal{R}}^{\otimes n}$ to ${\mathcal{R}}^{\otimes n}$ , the action of $Q$ implies a (right) eigenvalue problem

$${\mathit{Q}}^{★n}\left[{\mathbb{E}}^n\right]=\lambda {\mathbb{E}}^n.$$

(4)

The eigentensors ${\mathbb{E}}^n$ of any particular ${\mathit{Q}}^{★n}$ depend on that rotation’s axis. However, they can be grouped into certain sets of eigentensors which span subspaces of ${\mathcal{R}}^{\otimes n}$ which are independent of the axis chosen, forming a complete decomposition of the tensor space. Formally, we call ${\mathcal{W}}^n\subset {\mathcal{R}}^{\otimes n}$ an invariant subspace of ${\mathcal{R}}^{\otimes n}$ if

$${\mathbb{W}}^n\in {\mathcal{W}}^n\Rightarrow {\mathit{Q}}^{★n}\left[{\mathbb{W}}^n\right]\in {\mathcal{W}}^n,\forall \mathit{Q}\in \mathrm{SO}\left(3\right).$$

(5)

If an invariant subspace cannot be decomposed further into smaller invariant subspaces, we call it a harmonic subspace. Correspondingly, the decomposition into harmonic subspaces is known in continuum mechanics as the harmonic decomposition of ${\mathcal{R}}^{\otimes n}$ [14]. Harmonic subspaces of ${\mathcal{R}}^{\otimes n}$ are at most of order $2n+1$ . In fact, there is exactly one harmonic subspace of dimension $2n+1$ , sometimes called irreducible subspace in continuum mechanics. As it is traceless and symmetric, we will call it the deviatoric subspace ${\mathcal{D}}^n$ . Other harmonic subspaces have dimension $2k+1$ , $k<n$ and are each isomorphic to the $k$ th-order deviatoric subspace ${\mathcal{D}}^k$ . We call the map from the $k$ th order deviatoric subspace to an $n$ th-order harmonic subspace the inclusion ${\mathbb{P}}^{\mathrm{nk}}$ .

Using these concepts, we say that a basis of an $n$ th-order tensor space is harmonic if its basis tensors are elements of harmonic subspaces. In the following, we will calculate orthonormal harmonic bases by calculating bases for deviatoric subspaces and finding their inclusions into higher-order spaces.

2.2. Harmonic bases for deviatoric tensor spaces

By the eigenvalue problem in equation (4), a basis of a deviatoric subspace can be found using deviatoric eigentensors of an arbitrarily chosen rotation $Q$ . Since the eigentensor decomposition of $Q$ depends on the axis of rotation, the harmonic basis convention depends on a choice of axis. We choose the Z axis ${\mathit{e}}_3$ to ensure compatibility to a number of other common conventions in continuum mechanics, such as the ZXZ and ZYZ Euler angle conventions. This specific choice does not affect the generality of the approach.

While it is possible to calculate these eigentensors directly, for high-order tensors, it is cumbersome. Instead, one may use the fact that $\mathrm{SO}\left(3\right)$ is a Lie group with a corresponding Lie algebra $\mathfrak{s}\mathfrak{o}\left(3\right)$ . Eigentensors of $Q$ are eigentensors of $\mathfrak{q}\in \mathfrak{s}\mathfrak{o}\left(3\right)$ as eigentensors are preserved by the exponential map

$$Q=\exp \left(\mathfrak{q}\right).$$

(6)

Just as the action of $\mathrm{SO}\left(3\right)$ on higher-order tensor spaces is defined by the Rayleigh product ★, the action of its Lie algebra $\mathfrak{s}\mathfrak{o}\left(3\right)$ is described by the infinitesimal Rayleigh product

$$\begin{array}{c}\mathfrak{q}\left({\mathbb{V}}^n\right)=\mathit{q}⊠{\mathbb{V}}^n={\mathfrak{q}}_{\mathrm{ai}}{V}_{\mathrm{ij}\dots k}{\mathit{e}}_a\otimes {\mathit{e}}_j\dots {\mathit{e}}_k\hfill \\ +{\mathfrak{q}}_{\mathrm{bj}}{V}_{\mathrm{ij}\dots k}{\mathit{e}}_i\otimes {\mathit{e}}_b\dots {\mathit{e}}_k\hfill \\ +\dots \hfill \\ +{\mathfrak{q}}_{\mathrm{ck}}{V}_{\mathrm{ij}\dots k}{\mathit{e}}_i\otimes {\mathit{e}}_j\dots {\mathit{e}}_c.\hfill \end{array}$$

(7)

Analogously to the Rayleigh power of equation (3), the infinitesimal Rayleigh power is defined as

$$\begin{array}{c}\multicolumn{1}{c}{\left({\mathit{q}}_i^{⊠n}\right)={\mathfrak{q}}_{\mathrm{ai}}\dots {\delta }_{\mathrm{bj}}{\mathit{e}}_a\otimes \dots \otimes {\mathit{e}}_b\otimes {\mathit{e}}_i\otimes \dots \otimes {\mathit{e}}_j}\\ \multicolumn{1}{c}{+\dots }\\ \multicolumn{1}{c}{+{\delta }_{\mathrm{ai}}\dots {\mathfrak{q}}_{\mathrm{bj}}{\mathit{e}}_a\otimes \dots \otimes {\mathit{e}}_b\otimes {\mathit{e}}_i\otimes \dots \otimes {\mathit{e}}_j.}\end{array}$$

(8)

The following eigentensor calculation scheme is based on the angular momentum operator approach known from quantum mechanics [15]. Using the Levi-Civita or permutation tensor $ϵ$ and the imaginary number $\mathfrak{i}$ , we define angular momentum operators around the first-order axes ${\mathit{e}}_i$ as

$${\mathit{J}}_i=\mathfrak{i}ϵ\left[{\mathit{e}}_i\right].$$

(9)

These momentum operators form a basis of the tensor representation of the action of $\mathfrak{s}\mathfrak{o}\left(3\right)$ on $\mathcal{R}$ . The scalar factor $\mathfrak{i}$ is an arbitrarily chosen convention in quantum mechanics; we use it for our results to be compatible to pre-existing results and implementations.

As ${\mathit{J}}_3$ is a second-order tensor, calculating its eigenvectors to use as a basis for $\mathcal{R}$ is straightforward. However, to define bases for $n$ th-order tensor spaces, eigentensors of $2n$ th-order momentum operators ${\mathit{J}}_i^{⊠n}$ are required. Based on the eigentensors of ${\mathit{J}}_i$ , the $n$ th-order eigentensor associated with the eigenvalue $-n\mathfrak{i}$ is given by

$${\mathbb{E}}^n(-n\mathfrak{i})={\left(\frac{1}{\sqrt{2}}({\mathit{e}}_1+\mathfrak{i}{\mathit{e}}_2)\right)}^{\otimes n}.$$

(10)

To calculate other eigentensors based on ${\mathbb{E}}^n(-n\mathfrak{i})$ , we use the Z-axis raising operator given by Edmonds [15] among others as

$${\mathit{J}}_{+}^{⊠n}={\mathit{J}}_1^{⊠n}+{\mathit{J}}_2^n.$$

(11)

This operator maps a given eigentensor ${\mathbb{E}}^n\left(w\right)$ of ${\mathit{J}}_3^{⊠n}$ with eigenvalue $w$ to that with $w+\mathfrak{i}$ , but introduces a scalar factor. To calculate all $2n+1$ eigentensors, we recursively apply this operator and renormalize, resulting in

$${\mathbb{E}}^n\left((k+1)\mathfrak{i}\right)=\frac{1}{\sqrt{(n-k)(n+k+1)}}{\mathit{J}}_{+}^n\left[{\mathbb{E}}^n\left(k\mathfrak{i}\right)\right].$$

(12)

In quantum mechanics, this complexified orthonormal eigentensor basis $\{{\mathbb{E}}^n(-k\mathfrak{i}),\dots {\mathbb{E}}^n\left(k\mathfrak{i}\right)\}$ is used directly. This approach is also used by Man and Du [12] in a texture description context. We note that, respective to a complexified basis, even fully real tensors may have complex components. As classical continuum mechanics deals with fully real tensors, it is simpler to use a real basis instead. We define the nth-order real deviatoric basis ${\mathbb{D}}_i^n$ as

$${\mathbb{D}}_{2(n-k)+1}^n=\frac{1}{\sqrt{2}}({\mathbb{E}}^n\left(k\mathfrak{i}\right)+{(-1)}^k{\mathbb{E}}^n(-k\mathfrak{i})),k\in [1,n]$$

(13)

$${\mathbb{D}}_{2(n-k)+2}^n=\frac{-\mathfrak{i}}{\sqrt{2}}({\mathbb{E}}^n\left(k\mathfrak{i}\right)+{(-1)}^{k+1}{\mathbb{E}}^n(-k\mathfrak{i})),k\in [1,n]$$

(14)

$${\mathbb{D}}_{2n+1}^n=\mathbb{E}n\left(0\right).$$

(15)

These tensors are real]-valued because the eigentensors form complex conjugate pairs

$${\mathbb{E}}^n(-k\mathfrak{i})={\overline{\mathbb{E}}}^n(-k\mathfrak{i}),k\in [1,n].$$

(16)

Furthermore, they are orthonormal because the ${\mathbb{E}}^n\left(k\mathfrak{i}\right)$ are orthonormal. Finally, they form a basis because of orthonormality and the fact that there are $2n+1$ tensors ${\mathbb{D}}_i^n$ for any given $n$ .

As the above definition is linear, there exists a transformation matrix ${\underset{¯}{\underset{¯}{T}}}^n$ that maps ${\mathbb{E}}_j^n$ to ${\mathbb{D}}_i^n$ via

$${\mathbb{D}}_i^n=T_{\mathrm{ij}}^n{\mathbb{E}}^n(j\mathfrak{i}-n\mathfrak{i}).$$

(17)

With this particular convention for ${\mathbb{D}}_i^n$ ,

$${\mathbb{D}}_i^1={\mathit{e}}_i.$$

(18)

Because these deviatoric bases are fully real and orthonormalized, it is not necessary to distinguish between primal and dual bases or co- and contravariant indices.

2.3. Harmonic bases for arbitrary tensor spaces

To compose arbitrary tensor spaces from deviatoric subspaces, we use the harmonic decomposition. Unlike the usual depictions of that decomposition in continuum mechanics, e.g., Boehler et al. [16] or Forte and Vianello [14], we construct the decomposition of higher-order tensor spaces directly from that of lower-order ones. For brevity of notation, we introduce various operations on tensor spaces which can be defined by applying tensorial operations to their respective basis tensors but are nonetheless basis-independent. The dyadic product of tensor spaces contains all linear combinations of pairwise dyadic products of the respective basis tensors, reading

$${\mathcal{V}}^n\otimes {\mathcal{V}}^m=\{K_{\mathrm{ij}}{\mathbb{B}}_i^{\mathcal{V}n}\otimes {\mathbb{B}}_j^{\mathcal{V}m}:K_{\mathrm{ij}}\in \mathbb{R}\}.$$

(19)

It can be shown that the space being spanned is independent of the precise choice of bases ${\mathbb{B}}^{\mathcal{V}n}$ and ${\mathbb{B}}^{\mathcal{V}m}$ .

An analogous definition of products between tensors and tensor spaces leads to element-wise application of the product, e.g.

$$\mathbb{W}\left[\mathcal{V}\right]=\{\mathbb{W}\left[\mathbb{V}\right]:\mathbb{V}\in \mathcal{V}\}.$$

(20)

When writing the sum of two tensor spaces $\mathcal{U}$ and $\mathcal{W}$ , we presume that they form an orthogonal decomposition of a tensor space $\mathcal{V}$ such that

$$\mathcal{U}+\mathcal{W}=\{\mathbb{U}+\mathbb{W}:\mathbb{U}\in \mathcal{U}\mathrm{and}\mathbb{W}\in \mathcal{W}\}=\mathcal{V}.$$

(21)

As the fundamental units of the harmonic decomposition are the deviatoric tensor spaces, we first calculate how the dyadic product of two deviatoric tensor spaces decomposes harmonically. Fundamentally, the space decomposes harmonically into

$${\mathcal{D}}^m\otimes {\mathcal{D}}^n=\sum _{l=|m-n|}^{m+n}\mathcal{H}(m,n,l)$$

(22)

with one harmonic subspace $\mathcal{H}(m,n,l)$ per deviatoric order $l$ . For example,

$$\mathcal{H}(m,n,m+n)={\mathcal{D}}^{m+n},$$

(23)

while $\mathcal{H}(m,n,0)$ , which exists only for $m=n$ , is isomorphic to the scalars and therefore an isotropic subspace. Note that because each $\mathcal{H}(m,n,l)$ is of different dimension, this decomposition is unique. As the harmonic subspaces $\mathcal{H}(m,n,l)$ are orthogonal to each other and their dimensions add up to $\mathrm{mn}$ , the decomposition is complete.

To calculate $\mathcal{H}$ , consider the maps portrayed in Figure 1. First, there exists a unique projection from ${\mathcal{D}}^m\otimes {\mathcal{D}}^n$ to $\mathcal{H}(m,n,l)$ which can be given as a $2m+2n$ -order tensor of rank $2l+1$ that has projector properties. Its pseudoinverse is the inclusion form $\mathcal{H}(m,n,l)$ into ${\mathcal{D}}^m\otimes {\mathcal{D}}^n$ . Second, the isomorphism between $\mathcal{H}(m,n,l)$ and ${\mathcal{D}}^l$ is again a tensor, this one of order $m+n+l$ and rank $2l+1$ . This tensor, being an isomorphism, is bijective when considered as a map from $\mathcal{H}(m,n,l)$ . Being more interested in its inverse, we call this isomorphism ${\left({\mathbb{C}}^{\mathrm{mnl}}\right)}^{-1}$ . By inserting

$$\mathcal{H}(m,n,l)={\mathbb{C}}^{\mathrm{mnl}}\left[{\mathcal{D}}^l\right]$$

(24)

into equation (22), we retrieve

$${\mathcal{D}}^m\otimes {\mathcal{D}}^n=\sum _{l=|m-n|}^{m+n}{\mathbb{C}}^{\mathrm{mnl}}\left[{\mathcal{D}}^l\right],$$

(25)

which is an explicit harmonic decomposition of the dyadic product of two deviatoric tensor spaces. The superscripts $m$ , $n$ , and $l$ which indicate tensor order are not indices and the summation convention does not apply to them.

Figure 1.

Maps between harmonic and deviator subspaces.

To calculate the tensor ${\mathbb{C}}^{\mathrm{mnl}}$ , we borrow established calculation schemes from quantum mechanics. There, the deviatoric eigentensors $\mathbb{E}\left(k\mathfrak{i}\right)$ are discussed in the context of angular momentum as Eigenstates of the angular momentum operator $\mathit{J}$ . Correspondingly, ${\mathbb{C}}^{\mathrm{mnl}}$ is discussed in the context of coupled angular momenta and is known by its tensor components

$${\left(C^{\mathrm{mnl}}\right)}_{o\prime p\prime }^{q\prime }={\mathbb{C}}^{\mathrm{mnl}}·({\mathit{e}}^m(q\prime \mathfrak{i})\otimes {\mathit{e}}^n(o\prime \mathfrak{i})\otimes {\overline{\mathit{e}}}^l(p\prime \mathfrak{i})),$$

(26)

which are called Clebsch–Gordan coefficients, see e.g. Griffiths [17]. Here, the indices $q\prime$ , $o\prime$ , and $p\prime$ start at $-m$ , $-n$ , and $-l$ respectively, as is common in quantum mechanics, and the complex conjugate $\overline{e}$ is a consequence of the complexified tensor algebra involved. Edmonds [15] gives explicit recursive calculation schemes for Clebsch–Gordan coefficients but uses the Bra–Ket notation, which compares to the one above via

$${\left(C^{\mathrm{mnl}}\right)}_{o\prime p\prime }^{q\prime }=〈\mathrm{mno}\prime p\prime |\mathrm{mnlq}\prime 〉.$$

(27)

The Bra–Ket notation emphasizes the Clebsch–Gordan coefficients being linear maps from component vectors of dimension $2l+1$ to component matrices of size $2(m+1)\times 2(n+1)$ , while the notation used above emphasizes their tensorial nature. In either case, the Clebsch–Gordan coefficients, being tensor components, admit basis transformations, in particular the transformation from equation (17), which yields

$$c_{\mathrm{ijk}}^{\mathrm{mnl}}={\mathfrak{i}}^{m+n+l}{T}_{\mathrm{io}}^m{T}_{\mathrm{jp}}^n\overline{T_{\mathrm{kq}}^l}{\left(C^{\mathrm{mnl}}\right)}_{o\prime p\prime }^{q\prime }.$$

(28)

The indices without primes are defined as $o=o\prime +m+1$ such that $o\in [1,2m+1]$ , and respectively for $q$ and $n$ and $p$ and $l$ . The components $c_{\mathrm{ijk}}^{\mathrm{mnl}}$ are given respective to real-valued deviatoric bases, such that the Clebsch–Gordan tensor can be written

$${\mathbb{C}}^{\mathrm{mnl}}=c_{\mathrm{ijk}}^{\mathrm{mnl}}{\mathbb{D}}_i^m\otimes {\mathbb{D}}_j^n\otimes {\mathbb{D}}_k^l.$$

(29)

Because the coefficients and basis tensors in this expression are fully real, ${\mathbb{C}}^{\mathrm{mnl}}$ is a fully real tensor. Further properties of $\mathbb{C}$ are elaborated on in Appendix 1.

Using the Clebsch–Gordan tensors, the harmonic basis for ${\mathcal{D}}^m\otimes {\mathcal{D}}^n$ has orthonormal basis tensors

$${\mathbb{H}}_{i(l,j)}^{m+n}={\mathbb{C}}^{\mathrm{mnl}}[{\mathbb{D}}_j^l],l\in [|m-n|,m+n],j\in [0,l],$$

(30)

with an appropriate enumeration convention $i(l,j)$ . By convention, we sort the harmonic subspaces starting from the least-dimensional, which leads to

$$i(l,j)=j+\sum _{k=|m-n|}^{l}2(k+1).$$

(31)

In the general case, we are interested in the harmonic decomposition of an arbitrary tensor space which is constructed as the dyadic product of two other tensor spaces ${\mathcal{V}}^m$ and ${\mathcal{V}}^n$ . In the interest of building an iterative procedure, we assume that ${\mathcal{V}}^m$ is already decomposed into harmonic subspaces of dimension $2l+1$ , of which there are $a(m,l)$ . Instead of giving a formula for $a(m,l)$ , we treat it as the number of inclusions ${\mathbb{P}}_i^{\mathrm{mp}}$ which map the deviatoric space ${\mathcal{D}}^p$ to the various $2l+1$ -dimensional harmonic subspaces and give a formula for those inclusions. Then, ${\mathcal{V}}^m$ is given as

$${\mathcal{V}}^m=\sum _{p=0}^m\sum _{i=1}^{a(m,p)}{\mathbb{P}}_i^{\mathrm{mp}}\left[{\mathcal{D}}^p\right],$$

(32)

and ${\mathcal{V}}^n$ is given as

$${\mathcal{V}}^n=\sum _{q=0}^n\sum _{j=1}^{a(n,q)}{\mathbb{P}}_j^{\mathrm{nq}}\left[{\mathcal{D}}^q\right].$$

(33)

Brevity here causes an imprecision of notation. Even if $m=n$ , ${\mathcal{V}}^m$ is not necessarily equal ${\mathcal{V}}^n$ , for example, if different index symmetries apply to these spaces. Therefore, the inclusions ${\mathbb{P}}_i^{\mathrm{mp}}$ depend on the space considered. The same is true for the number of inclusions $a(m,l)$ , which generally increases with $m$ but may be as small as $0$ depending on the index symmetries which apply to ${\mathcal{V}}^m$ , in which case, there are no summands in the respective sum. These imprecisions of notation notwithstanding, the calculation shown here is general.

The space ${\mathcal{V}}^m\otimes {\mathcal{V}}^n$ is given as

$${\mathcal{V}}^m\otimes {\mathcal{V}}^n=\sum _{p=0}^m\sum _{q=0}^n\sum _{i=1}^{a(m,p)}\sum _{j=1}^{a(n,q)}{\mathbb{P}}_i^{\mathrm{mp}}\left[{\mathcal{D}}^p\right]\otimes {\mathbb{P}}_j^{\mathrm{nq}}\left[{\mathcal{D}}^q\right].$$

(34)

We define a $p,q$ -order tensorial Kronecker product ${\times }^{\mathrm{pq}}$ such that

$$\left({\mathbb{P}}^{\mathrm{mp}}{\times }^{\mathrm{pq}}{\mathbb{P}}^{\mathrm{nq}}\right)[{\mathcal{D}}^p\otimes {\mathcal{D}}^q]={\mathbb{P}}^{\mathrm{mp}}\left[{\mathcal{D}}^p\right]\otimes {\mathbb{P}}^{\mathrm{nq}}\left[{\mathcal{D}}^q\right].$$

(35)

Using the Clebsch–Gordan tensors, we find

$${\mathbb{P}}_i^{mp}\left[{\mathcal{D}}^p\right]\otimes {\mathbb{P}}_j^{nq}\left[{\mathcal{D}}^q\right]={\displaystyle \sum }_{l=\left|p-q\right|}^{p+q}\left({\mathbb{P}}_i^{mp}\stackrel{pq}{\times }{\mathbb{P}}_j^{nq}\right)\left[{\mathbb{C}}^{pql}\left[{\mathcal{D}}^l\right]\right].$$

(36)

The sums can be rearranged with the deviatoric order $l$ and a new index variable $k$ such that a tensor product of arbitrary tensor spaces ${\mathcal{V}}^m\otimes {\mathcal{V}}^n$ can be harmonically decomposed into

$${\mathcal{V}}^m\otimes {\mathcal{V}}^n=\sum _{l=0}^{m+n}\sum _{k=0}^{a(m+n,l)}{\mathbb{P}}_k^{(m+n)l}\left[{\mathcal{D}}^l\right].$$

(37)

With the $n$ -fold tensor contraction $̣n$ , the new inclusions can be written as

$${\mathbb{P}}_k^{(m+n)l}={\left({\mathbb{P}}_i^{\mathrm{mp}}{\times }^{\mathrm{pq}}{\mathbb{P}}_j^{\mathrm{nq}}\right)}^{p\underset{.}{+}q}{\mathbb{C}}^{\mathrm{pql}}.$$

(38)

Here $k$ enumerates all possible combinations of $p$ and $q$ for which ${\mathbb{C}}^{\mathrm{pql}}$ exists, i.e. with $l\in [|p-q|,p+q]$ . Within each of these $p,q$ combinations, all combinations of $i,j$ need to be enumerated, of which there are $a(m,p)a(n,p)$ . We choose an enumeration convention for $k$ by sorting in the order $p$ , $q$ , $i$ , $j$ , each ascending. In all concrete examples considered in this paper, the number of harmonic subspaces is limited due to various index symmetries, and the enumeration of the various projectors proves quite simple.

Relative to a harmonic basis, an $n$ th-order tensor ${\mathbb{A}}^n$ is depicted as

$${\mathbb{A}}^n=\sum _{k=0}^n\sum _{i=0}^{a(n,k)}\sum _{j=1}^{2k+1}{A}_{h(k,i,j)}{\mathbb{P}}_i^{\mathrm{nk}}\left[{\mathbb{D}}_j^k\right]$$

(39)

with the harmonic basis components $A_{h(k,i,j)}$ and basis tensors ${\mathbb{P}}_i^{\mathrm{nk}}\left[{\mathbb{D}}_j^k\right]$ . In the form

$${\mathbb{A}}^n=\sum _{k=0}^n\sum _{i=0}^{a(n,k)}{\mathbb{P}}_i^{\mathrm{nk}}\left[\sum _{j=1}^{2k+1}{A}_{h(k,i,j)}{\mathbb{D}}_j^k\right],$$

(40)

it is more clear that the harmonic components $A_{h(k,i,j)}$ are organized into groups of $2k+1$ components with the same $k$ and $i$ which belong to the same harmonic subspace. Therefore, the components themselves represent a harmonic decomposition of ${\mathbb{A}}^n$ .

The construction of inclusions from Clebsch–Gordan tensors defined above allows the calculation of a harmonic decomposition from the ground up by repeatedly applying dyadic products. By using the ${\mathbb{D}}^p$ basis for the deviatoric spaces ${\mathcal{D}}^p$ , we can thus construct harmonic bases, albeit with some effort. Thankfully, in continuum mechanics, a few specific harmonic bases are most interesting, and it is usually not necessary to construct others. With that in mind, the following sections are concerned with specific bases and their practical applications.

3. Second-order harmonic bases

3.1. Definition

To calculate the harmonic basis equivalent of the Voigt basis, we harmonically decompose the space of second-order tensors ${\mathcal{R}}^{\otimes 2}$ . First, we note that both the zeroth-order tensor space $\mathbb{R}$ and the first-order tensor space $\mathcal{R}$ are fully deviatoric, i.e.,

$$\mathcal{R}={\mathcal{D}}^1.$$

(41)

The harmonic basis of ${\mathcal{R}}^{\otimes 2}$ fundamentally results from the definition of invariant subspaces (equation (5)). Working with this definition, the isotropic (or spherical) subspace

$$V^{°}=\frac{\mathrm{sp}\left(\mathit{V}\right)}{3}\mathit{I}$$

(42)

is harmonic. The traceless part of ${\mathcal{R}}^{\otimes 2}$ remains. Since rotations only map symmetric tensors to symmetric tensors, the remaining traceless tensors can be further split into the skew-symmetric and symmetric deviatoric parts

$$\mathrm{skw}\left(\mathit{V}\right)=\frac{1}{2}(\mathit{V}-{\mathit{V}}^{\top }),$$

(43)

$$\mathit{V}\prime =\mathit{V}-{\mathit{V}}^{°}-\mathrm{skw}\left(\mathit{V}\right).$$

(44)

We thus arrive at the harmonic decomposition of ${\mathcal{R}}^{\otimes 2}$ , which is exactly the decomposition into spherical, symmetric-deviatoric, and skew-symmetric parts. The same result can also be reached by direct calculation. We define ${\mathcal{R}}^{\otimes 2}$ via the product

$${\mathcal{R}}^{\otimes 2}=\mathcal{R}\otimes \mathcal{R}=D^1\otimes D^1$$

(45)

and use equation (25) to calculate the harmonic decomposition in the form

$${\mathcal{D}}^{\otimes 2}=c^{110}\left[D^0\right]+{\mathbb{C}}^{111}\left[D^1\right]+{\mathbb{C}}^{112}\left[D^2\right].$$

(46)

The relevant Clebsch–Gordan tensors can be calculated via equation (29) as

$${\mathit{c}}^{110}=\frac{1}{\sqrt{3}}\mathit{I},$$

(47)

$${\mathbb{C}}^{111}=\frac{1}{\sqrt{2}}ϵ,$$

(48)

$${\mathbb{C}}^{112}=\mathbb{P}\prime .$$

(49)

Using equation (39), a second-order tensor can therefore be depicted as

$$\mathit{A}=A_{h(0,1,1)}{\mathit{P}}_1^{20}{D}^0+A_{h(1,1,i)}{\mathbb{P}}_1^{21}\left[{\mathit{D}}_i^1\right]+A_{h(2,1,j)}{\mathbb{P}}_1^{22}\left[{\mathit{D}}_j^2\right].$$

(50)

Because ${\mathcal{R}}^2$ is the product of two deviatoric spaces, the inclusions are Clebsch–Gordan tensors. With $D^0=1$ and ${\mathit{D}}_i^1={\mathit{e}}_i$ ,

$$\mathit{A}=A_{h(0,1,1)}{\mathit{c}}^{110}+A_{h(1,1,i)}{\mathit{c}}^{111}\left[{\mathit{e}}_i\right]+A_{h(2,1,j)}{\mathit{c}}^{112}\left[{\mathit{D}}_j^2\right].$$

(51)

To find the basis of ${\mathcal{D}}^2$ , we apply the eigentensor calculation scheme of subsection 2.2, finding as eigentensors

$${\mathit{E}}^2\left(0\right)=\frac{1}{\sqrt{6}}(3{\mathit{e}}_3\otimes {\mathit{e}}_3-\mathit{I}),$$

(52)

$${\mathit{E}}^2\left(\mathfrak{i}\right)=-\mathrm{sym}({\mathit{e}}_1\otimes {\mathit{e}}_3)+\mathfrak{i}\mathrm{sym}({\mathit{e}}_2\otimes {\mathit{e}}_3),$$

(53)

$${\mathit{E}}^2(-\mathfrak{i})=\mathrm{sym}({\mathit{e}}_1\otimes {\mathit{e}}_3)+\mathfrak{i}\mathrm{sym}({\mathit{e}}_2\otimes {\mathit{e}}_3),$$

(54)

$${\mathit{E}}^2\left(2\mathfrak{i}\right)=\frac{1}{2}({\mathit{e}}_1\otimes {\mathit{e}}_1-{\mathit{e}}_2\otimes {\mathit{e}}_2)-\mathfrak{i}\mathrm{sym}({\mathit{E}}_2\otimes {\mathit{e}}_1),$$

(55)

$${\mathit{E}}^2(-2\mathfrak{i})=\frac{1}{2}({\mathit{e}}_1\otimes {\mathit{e}}_1-{\mathit{e}}_2\otimes {\mathit{e}}_2)+\mathfrak{i}\mathrm{sym}({\mathit{E}}_2\otimes {\mathit{e}}_1),$$

(56)

from which ${\mathit{D}}_i^2$ can be calculated via equation (17). After combining the deviatoric bases and the inclusion tensors, the second-order harmonic basis tensors in the Z-axis-convention read

$${\mathit{H}}_1^2=\frac{1}{\sqrt{3}}\mathit{I},$$

(57)

$${\mathit{H}}_2^2=\frac{1}{\sqrt{2}}\mathrm{skw}({\mathit{e}}_2\otimes {\mathit{e}}_3),$$

(58)

$${\mathit{H}}_3^2=\frac{1}{\sqrt{2}}\mathrm{skw}({\mathit{e}}_3\otimes {\mathit{e}}_1),$$

(59)

$${\mathit{H}}_4^2=\frac{1}{\sqrt{2}}\mathrm{skw}({\mathit{e}}_1\otimes {\mathit{e}}_2),$$

(60)

$${\mathit{H}}_5^2={\mathit{D}}_2^2=\frac{1}{\sqrt{2}}({\mathit{e}}_1\otimes {\mathit{e}}_1-{\mathit{e}}_2\otimes {\mathit{e}}_2),$$

(61)

$${\mathit{H}}_6^2={\mathit{D}}_1^2=\sqrt{2}\mathrm{sym}({\mathit{e}}_1\otimes {\mathit{e}}_2),$$

(62)

$${\mathit{H}}_7^2={\mathit{D}}_3^2=\sqrt{2}\mathrm{sym}({\mathit{e}}_1\otimes {\mathit{e}}_3),$$

(63)

$${\mathit{H}}_8^2={\mathit{D}}_4^2=\sqrt{2}\mathrm{sym}({\mathit{e}}_2\otimes {\mathit{e}}_3),$$

(64)

$${\mathit{H}}_9^2={\mathit{D}}_5^2=\frac{1}{\sqrt{6}}\mathrm{sym}(3{\mathit{e}}_3\otimes {\mathit{e}}_3-\mathit{I}).$$

(65)

3.2. Application 1: small-strain elasticity

Small-strain elasticity is described by Hooke’s Law,

$$\mathit{\sigma }=\mathbb{C}\left[\mathit{\epsilon }\right],$$

(66)

a linear constitutive material law which links the second-order Cauchy stress $\mathit{\sigma }$ and the second-order strain tensor $\mathit{\epsilon }$ via the fourth-order stiffness tensor $\mathbb{C}$ . In the full second-order harmonic basis of equations (57)–(65), Hooke’s Law is a nine-dimensional matrix-vector equation. Since both stresses and strains are symmetric, the skew-symmetric parts of the second-order harmonic basis can be neglected, and stresses and strains may be written as vector components in a six-dimensional space,

$${\underset{¯}{\epsilon }}_{\mathrm{H}}=\left(\begin{array}{c}\frac{1}{\sqrt{3}}({\epsilon }_{11}+{\epsilon }_{22}+{\epsilon }_{33})\hfill \\ \multicolumn{1}{c}{\frac{1}{\sqrt{2}}({\epsilon }_{11}-{\epsilon }_{22})}\\ \multicolumn{1}{c}{\sqrt{2}{\epsilon }_{12}}\\ \multicolumn{1}{c}{\sqrt{2}{\epsilon }_{13}}\\ \multicolumn{1}{c}{\sqrt{2}{\epsilon }_{23}}\\ \multicolumn{1}{c}{\frac{1}{\sqrt{6}}(2{\epsilon }_{33}-{\epsilon }_{11}-{\epsilon }_{22})}\end{array}\right),$$

(67)

where we denote this symmetrized harmonic basis $\left({\mathit{H}}^2\right)$ as $\mathrm{H}$ for brevity. We note here that Zheng and Spencer [13] show how the full harmonic basis including skew-symmetric tensors is useful in micropolar elasticity. Their canonical basis is similar in origin to our harmonic basis, but differs in convention. As a result, our symmetric second-order harmonic basis tensors are associated with elementary strain modes, as shown in Figure 2.

Figure 2.

Harmonic second-order basis tensors visualized as elementary strain modes. (a) ${\mathit{H}}_1^2$ ; (b) ${\mathit{H}}_2^2$ ; (c) ${\mathit{H}}_6^2$ ; (d) ${\mathit{H}}_3^2$ ; (e) ${\mathit{H}}_4^2$ ; and (f) ${\mathit{H}}_5^2$ .

When purely deviatoric tensors are considered, the resulting five-dimensional ${\mathit{D}}^2$ -basis is equivalent to the various Kocks-type bases, for example those of [5]

$${\underset{¯}{\epsilon }}_{\mathrm{Lequeu}}=\left(\begin{array}{c}\frac{1}{\sqrt{2}}({\epsilon }_{22}-{\epsilon }_{11})\hfill \\ \multicolumn{1}{c}{\sqrt{\frac{3}{2}}({\epsilon }_{11}+{\epsilon }_{22})}\\ \multicolumn{1}{c}{\sqrt{2}{\epsilon }_{23}}\\ \multicolumn{1}{c}{\sqrt{2}{\epsilon }_{13}}\\ \multicolumn{1}{c}{\sqrt{2}{\epsilon }_{12}}\end{array}\right)$$

(68)

and [6]

$${\underset{¯}{\epsilon }}_{\mathrm{vanHoutte}}=\left(\begin{array}{c}\frac{(\sqrt{3}+1){\epsilon }_{22}+(\sqrt{3}-1){\epsilon }_{33}}{2}\hfill \\ \multicolumn{1}{c}{\frac{(\sqrt{3}-1){\epsilon }_{22}+(\sqrt{3}+1){\epsilon }_{33}}{2}}\\ \multicolumn{1}{c}{\sqrt{2}{\epsilon }_{23}}\\ \multicolumn{1}{c}{\sqrt{2}{\epsilon }_{13}}\\ \multicolumn{1}{c}{\sqrt{2}{\epsilon }_{12}}\end{array}\right).$$

(69)

A full discussion of those bases can be found in Manik’s study [18].

The basis most commonly used in commercial software is the non-normalized basis of Voigt [1]. Also common is its normalized equivalent, the Mandel basis described by e.g. Mandel [3]. Conversions between the harmonic basis and the Mandel basis can be written as

$${\underset{¯}{\epsilon }}_{\mathrm{H}}=\underset{¯}{\underset{¯}{M}}{\underset{¯}{\epsilon }}_{\mathrm{M}}$$

(70)

with

$$\underset{¯}{\underset{¯}{M}}=\left(\begin{array}{cccccc}\frac{1}{\sqrt{3}}\hfill & \frac{1}{\sqrt{3}}\hfill & \frac{1}{\sqrt{3}}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{1}{\sqrt{2}}}& \frac{-1}{\sqrt{2}}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{-1}{\sqrt{6}}}& \frac{-1}{\sqrt{6}}\hfill & \frac{2}{\sqrt{6}}\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right).$$

(71)

For conversions from the non-normalized Voigt basis, rows 3–5 acquire a factor of $\sqrt{2}$ or its inverse depending on whether the strain or stress Voigt basis is used. As some commercial software uses different index conventions, the shear components in rows 3–5 may need to be reordered.

We now compare the harmonic and the Mandel basis for depicting the stiffness $\mathbb{C}$ given common material symmetries. First, for the isotropic case $\mathbb{C}={\mathbb{C}}^{\mathrm{I}}$ ,

$${\underset{¯}{\underset{¯}{C}}}_{\mathrm{H}}^{\mathrm{I}}=\left(\begin{array}{cccccc}3K\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 2G\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 2G\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 2G\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 2G\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 2G\hfill \end{array}\right)$$

(72)

$$=\mathrm{diag}(3K,2G,2G,2G,2G,2G),$$

(73)

with the compression modulus $K$ and the shear modulus $G$ . In the Mandel notation,

$${\underset{¯}{\underset{¯}{C}}}_{\mathrm{M}}^{\mathrm{I}}=\left(\begin{array}{cccccc}\lambda +2G\hfill & \lambda \hfill & \lambda \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\lambda }& \lambda +2G\hfill & \lambda \hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{\lambda }& \lambda \hfill & \lambda +2G\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 2G\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 2G\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 2G\hfill \end{array}\right),$$

(74)

with $\lambda =K-2G/3$ . The diagonal structure of ${\underset{¯}{\underset{¯}{C}}}_{\mathrm{H}}^{\mathrm{I}}$ arises naturally from the fact that the symmetric isotropic projectors ${\mathbb{P}}^{°}$ and $\mathbb{P}\prime$ , which are the projectors on spherical and deviatoric second-order tensors respectively, are partial identity matrices in a basis which derives from that decomposition, i.e.

$${\underset{¯}{\underset{¯}{P}}}_{\mathrm{H}}^{°}=\mathrm{diag}(1,0,0,0,0,0),$$

(75)

$${{\underset{¯}{\underset{¯}{P}}}^{\prime }}_{\mathrm{H}}=\mathrm{diag}(0,1,1,1,1,1).$$

(76)

Diagonal matrices have many numerical advantages to non-diagonal matrices, in particular when inverting

$$\left({\underset{¯}{\underset{¯}{C}}}_{\mathrm{H}}^{\mathrm{I}}\right)=\mathrm{diag}{\left(\frac{1}{C_i}\right)}^{-1},$$

(77)

or applying non-linear functions $\mathbb{F}$ with a scalar equivalent $f$ ,

$$\underset{¯}{\underset{¯}{F}}{\left(\mathbb{C}\right)}_{\mathrm{H}}=\mathrm{diag}\left(f\left(C_i\right)\right),$$

(78)

or in matrix-vector multiplications, which are of complexity order $O\left(n\right)$ for diagonal matrices, instead of $O\left(n^2\right)$ for non-diagonal matrices. Therefore, isotropic linear elastic material laws, when evaluated in the harmonic basis, are less expensive than their Mandel counterparts.

For crystals aligned with the axes ${\mathit{e}}_1$ , ${\mathit{e}}_2$ , and ${\mathit{e}}_3$ , cubic stiffness tensors take the form

$${\underset{¯}{\underset{¯}{C}}}_{\mathrm{H}}^{\mathrm{C}}=\mathrm{diag}({\lambda }_1^{\mathrm{C}},{\lambda }_2^{\mathrm{C}},{\lambda }_3^{\mathrm{C}},{\lambda }_3^{\mathrm{C}},{\lambda }_3^{\mathrm{C}},{\lambda }_2^{\mathrm{C}})$$

(79)

with

$${\lambda }_1^{\mathrm{C}}=C_{1111}+2C_{1122},$$

(80)

$${\lambda }_2^{\mathrm{C}}=C_{1111}-C_{1122},$$

(81)

$${\lambda }_3^{\mathrm{C}}=2C_{2323}.$$

(82)

Similarly to the isotropic case, the diagonality of $\mathbb{C}$ can be explained using the symmetric cubic projectors ${\mathbb{P}}_1^{\mathrm{C}}$ , ${\mathbb{P}}_2^{\mathrm{C}}$ , and ${\mathbb{P}}_3^{\mathrm{C}}$ , which are partial identity matrices in the crystal-aligned harmonic basis, reading

$${\underset{¯}{\underset{¯}{P}}}_{1\mathrm{H}}^{\mathrm{C}}=\mathrm{diag}(1,0,0,0,0,0),$$

(83)

$${\underset{¯}{\underset{¯}{P}}}_{1\mathrm{H}}^{\mathrm{C}}=\mathrm{diag}(0,1,0,0,0,1),$$

(84)

$${\underset{¯}{\underset{¯}{P}}}_{1\mathrm{H}}^{\mathrm{C}}=\mathrm{diag}(0,0,1,1,1,0).$$

(85)

For transversal-isotropic or hexagonal stiffness tensors whose symmetry axis is aligned with ${\mathit{e}}_3$ ,

$${\underset{¯}{\underset{¯}{C}}}_{\mathrm{H}}^{\mathrm{H}}=\left(\begin{array}{cccccc}{C}_{\mathrm{H}11}^{\mathrm{H}}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & C_{\mathrm{H}16}^{\mathrm{H}}\hfill \\ \multicolumn{1}{c}{0}& C_{\mathrm{H}22}^{\mathrm{H}}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & C_{\mathrm{H}22}^{\mathrm{H}}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 2C_{2323}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 2C_{2323}\hfill & 0\hfill \\ \multicolumn{1}{c}{C_{\mathrm{H}16}^{\mathrm{H}}}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & C_{\mathrm{H}66}^{\mathrm{H}}\hfill \end{array}\right).$$

(86)

with

$$C_{\mathrm{H}11}^{\mathrm{H}}=\frac{1}{3}(2C_{1111}+2C_{1122}+4C_{1133}+C_{3333}),$$

(87)

$$C_{\mathrm{H}22}^{\mathrm{H}}=C_{1111}-C_{1122}$$

(88)

$$C_{\mathrm{H}66}^{\mathrm{H}}=\frac{\sqrt{2}}{3}(-C_{1111}-C_{1122}+C_{1133}+C_{3333}),$$

(89)

$$C_{\mathrm{H}16}^{\mathrm{H}}=\frac{1}{3}(C_{1111}+C_{1122}-4C_{1133}+2C_{3333}).$$

(90)

Note that ${\underset{¯}{\underset{¯}{C}}}_{\mathrm{H}}^{\mathrm{H}}$ is not diagonal as the eigentensors of a hexagonal stiffness tensor depend on its exact values; the matrix exhibits spherical-deviatoric coupling. If necessary, a full diagonalization can be achieved by diagonalizing the remaining 2x2 block matrix. By comparison, the hexagonal Mandel stiffness matrix reads

$${\underset{¯}{\underset{¯}{C}}}_{\mathrm{M}}^{\mathrm{H}}=\left(\begin{array}{cccccc}{C}_{1111}\hfill & C_{1122}\hfill & C_{1133}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{C_{1122}}& C_{1111}\hfill & C_{1133}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{C_{1133}}& C_{1133}\hfill & C_{3333}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 2C_{2323}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 2C_{2323}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & C_{\mathrm{H}22}^{\mathrm{H}}\hfill \end{array}\right).$$

(91)

In the harmonic basis, trigonal, tetragonal, and monoclinic stiffnesses are also simplified, as shown in Appendix 2. Orthotropic and fully triclinic stiffnesses are of a similar degree of complexity as in the Mandel notation. We conclude that choosing the harmonic basis instead of the Mandel basis, which is purely a matter of convention, simplifies calculations and representations for a large class of material symmetries.

Anisotropic quantities are not always aligned to the coordinate axes as in the cases above. A non-aligned stiffness tensor $\mathbb{C}$ can be transformed from an aligned stiffness tensor ${\mathbb{C}}^A$ via an active rotation

$$\mathbb{C}=\mathit{Q}★{\mathbb{C}}^A.$$

(92)

Respective to the harmonic basis, ${\mathbb{C}}^A$ is depicted as a component matrix ${\underset{¯}{\underset{¯}{C}}}_{\mathrm{H}}^A$ , such that

$${\underset{¯}{\underset{¯}{C}}}_{\mathrm{H}}={\underset{¯}{\underset{¯}{R}}}_{\mathrm{H}}{\underset{¯}{\underset{¯}{C}}}_{\mathrm{H}}^A{{\underset{¯}{\underset{¯}{R}}}_{\mathrm{H}}}^{\top }$$

(93)

with $6\times 6$ rotation matrices given by

$$R_{\mathrm{H}\mathrm{ij}}={\mathit{Q}}^{★2}·({\mathit{H}}_i^2\otimes {\mathit{H}}_j^2).$$

(94)

To calculate these harmonic rotation matrices, we begin with deviatoric rotation matrices. The nth order eigentensor basis ${\mathbb{E}}^n$ from equation (12) is defined via eigentensor decomposition of Z rotations. Therefore, Z rotations act on ${\mathcal{D}}^n$ as

$${\mathit{Q}}_{\mathrm{Z}}^{★n}\left(\alpha \right)\left[{\mathcal{D}}^n\right]=\sum _{k=-n}^n{e}^{\mathrm{ik}\alpha }({\mathbb{E}}^n\left(k\mathfrak{i}\right)\otimes {\mathbb{E}}^n\left(k\mathfrak{i}\right))\left[{\mathcal{D}}^n\right].$$

(95)

After transforming into the deviatoric basis via equation (17),

$$\sum _{k=-n}^n{e}^{-\mathrm{ik}\alpha }{\mathbb{E}}_k^n\otimes {\mathbb{E}}_k^n=d_{\mathrm{Z}\mathrm{ij}}^n{\mathbb{D}}_i^n\otimes {\mathbb{D}}_j^n$$

(96)

with

$$d_{\mathrm{Z}\mathrm{ij}}^n=\{\begin{array}{cc}\cos \left(k\alpha \right)\hfill & i=j=2(n-k)+1\hfill \\ \multicolumn{1}{c}{\cos \left(k\alpha \right)}& i=j=2(n-k)+2\hfill \\ \multicolumn{1}{c}{-\sin \left(k\alpha \right)}& i=2(n-k)+2,j=2(n-k)+1\hfill \\ \multicolumn{1}{c}{\sin \left(k\alpha \right)}& i=2(n-k)+1,j=2(n-k)+2\hfill \\ \multicolumn{1}{c}{1}& i=j=2n+1\hfill \\ \multicolumn{1}{c}{0}& \mathrm{else}.\hfill \end{array}$$

(97)

For $n=1$ , this is the classic rotation matrix around ${\mathit{e}}_3$ . Increasing the tensorial order from $n-1$ to $n$ adds a $2\times 2$ block analogous to an $n\alpha$ rotation matrix. To further clarify this point, we may consider the involved rotations. The rotation of ${\mathcal{D}}^1$ around the Z-axis proceeds as usual: of the three basis tensors ${\mathit{D}}_i^1={\mathit{e}}_i$ , ${\mathit{e}}_3$ does not rotate while the other two complete a full rotation after $\alpha =2\pi$ . In Figure 3, the second-order tensors are visualized as strain modes applied to a cylinder which may be imagined to rotate around ${\mathit{e}}_3$ . ${\mathit{D}}_5^2={\mathit{e}}_3\otimes {\mathit{e}}_3$ does not rotate. ${\mathit{D}}_3^2$ and ${\mathit{D}}_4^2$ complete a full rotation after $\alpha =2\pi$ , while ${\mathit{D}}_1^2$ and ${\mathit{D}}_2^2$ , due to symmetry, are indistinguishable from their starting form after $\alpha =\pi$ . In other words, ${\mathit{D}}_1^2$ and ${\mathit{D}}_2^2$ rotate with $2\alpha$ .

Figure 3.

${\mathit{D}}^2$ -transformed cylinders visualize the Z-rotation of ${\mathcal{D}}^2$ .

For rotations around other axes, the $d$ -matrices are different. Rotations around Y and X are calculated in quantum mechanics via Wigner-d-matrices. These are components respective to the ${\mathbb{E}}^n$ -basis and are given by Wigner [19] as

$$d_{\mathrm{Y}*\mathrm{ij}}^n\left(\beta \right)={((n+m\prime )!(n-m\prime )!(n+m)!(n-m)!)}^{\frac{1}{2}}\sum _{s=s_{\min }}^{s_{\max }}\frac{{(-1)}^{m\prime -m+s}{\left(\cos \frac{\beta }{2}\right)}^{2n+m-m\prime -2s}{\left(\sin \frac{\beta }{2}\right)}^{m\prime -m+2s}}{(n+m-s)!s!(m\prime -m+s)!(n-m\prime -s)!}$$

(98)

and

$$d_{\mathrm{X}*\mathrm{ij}}^n\left(\beta \right)={((n+m\prime )!(n-m\prime )!(n+m)!(n-m)!)}^{\frac{1}{2}}\sum _{s=s_{\min }}^{s_{\max }}\frac{{(-1)}^s{\mathrm{ı}}^{m-m\prime }{\left(\cos \frac{\beta }{2}\right)}^{2n+m-m\prime -2s}{\left(\sin \frac{\beta }{2}\right)}^{m\prime -m+2s}}{(n+m-s)!s!(m\prime -m+s)!(n-m\prime -s)!}$$

(99)

with

$$i=m\prime +n,j=m+n,$$

(100)

and

$$s_{\min }=\max (0,m-m\prime ),$$

(101)

$$s_{\max }=\min (n+m,n-m\prime ).$$

(102)

The Wigner- $d$ -functions can be transformed into the corresponding deviatoric basis via equation (17). For ${\mathit{D}}^2$ , we retrieve rotation matrices as shown in Appendix 3, Equations (184)–(186), from which the usual ZXZ or ZYZ Euler conventions can be assembled by multplication.

The full rotation ${\underset{¯}{\underset{¯}{R}}}_{\mathrm{H}}$ of $\mathrm{Sy}{m}^2$ is

$${\underset{¯}{\underset{¯}{R}}}_{\mathrm{H}}={\left(\begin{array}{cc}{\underset{¯}{\underset{¯}{d}}}_{\mathrm{ZXZ}{D}^0}^0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& {\underset{¯}{\underset{¯}{d}}}_{\mathrm{ZXZ}{\mathit{D}}^2}^2\hfill \end{array}\right)}_{\mathrm{H}}\mathrm{with}{d}_{\mathrm{ZXZ}}^0=1.$$

(103)

Even for the full three-angle rotation, the block matrix structure of $\underset{¯}{\underset{¯}{R}}$ remains. This sparsity simplifies calculations particularly in a theoretical context. For example, the effect of a rotation on a cubic stiffness is evident in this form: the compression modulus of a cubic single-crystal is rotation-invariant.

The above discussion is limited to elastic small-strain materials, which are a particularly simple material class. Similar issues of symmetry and anisotropy also arise in kinematically sophisticated approaches such as gradient materials or generalized continua (for general references on these topics, see Gurtin et al. [20] and Neff et al. [21], respectively). These continuum theories involve higher-order tensors which can result in lengthy expressions in the anisotropic case (see, e.g., Lazar et al. [22] and the references therein) for which the harmonic basis offers a concise representation. As the inclusion of gradient materials would go beyond the scope of the manuscript, we will not go into detail here.

3.3. Application 2: computational elastoviscoplasticity

The natural basis as defined by Kocks arose from a context of computational crystal plasticity. In this context, the natural basis is advantageous because it diagonalizes cubic stiffnesses and explicitly splits second-order tensors into spherical and deviatoric parts, as discussed by Kocks et al. [7] and, more recently, Mánik et al. [23]. When implementing crystal plasticity material laws, this spherical-deviatoric split reduces the size of equations systems that need to be solved numerically. As discussed in Section 3.2, the harmonic basis has the same properties. Mánik [18] showed that basis representations with these properties also reduce the numerical effort for calculations involving more general elastoplastic models that are not cubic, particularly if plastic strains conserve volume. This result can be extended beyond associated elastoplasticity by using the generalized standard material (GSM) framework defined by Halphen and Nguyen [24] and Germain et al. [25] in a small- deformation context. A GSM is specified by a free energy density $\mathrm{\Psi }$ , which governs, e.g., elastic effects, and a dissipation potential $\varphi$ , which describes the dissipative behaviour of the material such as viscoplastic effects. As a simple example, we assume an elastoviscoplastic material for which the elastic behavior is linear, which has no defect energies and for which the plastic strain ${\mathit{\epsilon }}_{\mathrm{p}}$ is the only internal variable. In that case, the elastic energy density is given as

$$\psi (\mathit{\epsilon },{\mathit{\epsilon }}_{\mathrm{p}})=\frac{1}{2}(\mathit{\epsilon }-{\mathit{\epsilon }}_{\mathrm{p}})·\mathbb{C}[\mathit{\epsilon }-{\mathit{\epsilon }}_{\mathrm{p}}].$$

(104)

If plastic strains conserve volume, the plastic strain is fully deviatoric,

$${\mathit{\epsilon }}_{\mathrm{p}}^{\prime }={\mathit{\epsilon }}_{\mathrm{p}},$$

(105)

and the dual variable of the plastic strain is likewise the deviatoric stress, here written with the deviatoric projector $\mathbb{P}\prime$ as

$$-\frac{\partial \psi }{\partial {\mathit{\epsilon }}_{\mathrm{p}}}=\mathbb{P}\prime \mathbb{C}[\mathit{\epsilon }-{\mathit{\epsilon }}_{\mathrm{p}}]=\mathit{\sigma }\prime .$$

(106)

The dual dissipation potential ${\phi }^{*}$ then depends only on $\mathit{\sigma }\prime$ . In this example, we will make no further assumptions regarding its form. The evolution equation is given by

$${\stackrel{·}{\mathit{\epsilon }}}_{\mathrm{p}}=\frac{\partial {\phi }^{*}(\mathit{\sigma }\prime )}{\partial \mathit{\sigma }\prime }.$$

(107)

The viscoplastic strain rate given by this evolution equation directly depends on the stress and therefore describes rate-dependent, i.e. viscoplastic behaviour. Rate-independent behaviour can be described in the GSM framework as well, but requires a more complex mathematical apparatus involving non-differential dissipation potentials as described by Halphen and quoc Son [26]. For reasons of simplicity, we continue with the rate-dependent form specified above, which we discretize via an implicit Euler approach, yielding

$$\frac{\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}}{\mathrm{\Delta }t}=\frac{\partial {\phi }^{*}(\mathbb{P}\prime \mathbb{C}[\mathit{\epsilon }-{\mathit{\epsilon }}_{\mathrm{p}}^{\mathit{n}}-\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}])}{\partial \mathit{\sigma }\prime },$$

(108)

where ${\mathit{\epsilon }}_{\mathrm{p}}^n$ is the plastic strain from the previous time step. This implicit equation for $\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}$ has second-order tensors as values and parameters. In the Mandel notation, it would appear as a six-dimensional vector equation; using the harmonic basis, it is five-dimensional. When solving via Newton’s method, the harmonic basis directly reduces the computational expense of evaluating the residual function

$$\mathit{f}\left(\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}\right)=\frac{\partial {\phi }^{*}\left(\mathbb{C}{[\mathit{\epsilon }-{\mathit{\epsilon }}_{\mathrm{p}}^{\mathit{n}}-\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}]}^{\prime }\right)}{\partial \mathit{\sigma }\prime }-\frac{\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}}{\mathrm{\Delta }t}$$

(109)

by one dimension. Furthermore, the Jacobian

$$\mathbb{J}=\frac{\partial \mathit{f}\left(\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}\right)}{\partial \mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}}=\frac{{\partial }^2{\phi }^{*}\left(\mathbb{C}{[\mathit{\epsilon }-{\mathit{\epsilon }}_p^n-\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}]}^{\prime }\right)}{\partial \mathit{\sigma }\prime \partial \mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}}-\frac{1}{\mathrm{\Delta }t}\mathbb{P}\prime$$

(110)

is reducible to a $5\times 5$ matrix in the harmonic basis. The residual update equation for the iteration step $i+1$ ,

$$\mathbb{J}\left(\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}^i\right)\left[\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}^{i+1}\right]=\mathit{f}\left(\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{p}}^{\mathit{i}}\right),$$

(111)

can be solved via a direct method with computational expense of order $O\left(n^3\right)$ . For the concrete example of a Cholesky factorization as described by e.g. [27], reducing the dimension of the problem from 6 to 5 by using the harmonic basis reduces the number of operations by 35%. More importantly, when using the Mandel notation, the Jacobian is a $6\times 6$ matrix of rank 5, i.e. singular. Because of numerical imprecision, a naive numerical inversion attempt might yet succeed, yielding erroneous results; a problem which the harmonic basis avoids entirely.

The example material considered here is both anisotropic and generally nonlinear, but simple in the sense of having a small number of internal variables. We expect similar results for materials with hardening, both isotropic and anisotropic. As evidenced by e.g. the existence of the various Kocks-type bases, the approach of modelling incompressible problems in an isochoric parameter space is not a new development. The harmonic basis is useful because transforming to an isochoric parameter space is as simple as dropping the first component. Even more generally, the harmonic basis can be of use whenever a continuum-mechanical problem can be separated into compressive and deviatoric parts.

4. Texture descriptions using the harmonic basis

4.1. Tensorial texture descriptions

The properties of locally anisotropic microstructures generally depend on local orientation, such as grain orientations in a polycrystal. One method of describing such microstructures is the orientation distribution function (ODF), a one-point probability distribution of orientation, which is defined as a probability density on the orientation space

$$f:\mathrm{SO}\left(3\right)\to {\mathbb{R}}_{\ge 0},$$

(112)

with the normalization property

$${\int }_{\mathrm{SO}\left(3\right)}f\left(Q\right)\mathrm{d}V\left(Q\right)=1.$$

(113)

The ODF is directly linked to anisotropy in various material quantities, such as effective stiffnesses (e.g. Böhlke and Bertram [28]) and plasticity (e.g. Man and Huang [29]). The calculation of effective material quantities from the ODF generally requires additional micromechanical assumptions.

We assume the ODF to be square-integrable. By choosing a basis of the space of square-integrable functions, the ODF can be described with a series of coefficients relative to this basis, cf. Bunge [9]. A coordinate-independent series expansion is possible by using tensorial texture coefficients coefficients (Adams et al. [10]). Lobos Fernández and Böhlke [11] define these tensors as general deviatoric basis tensors ${\mathit{H}}_{〈n〉}$ equivalent to our ${\mathbb{D}}^n$ , but do not choose a specific basis. Building on those results, Man and Du [12] have chosen the complexified eigentensor bases which we refer to in this work as ${\mathbb{E}}^n$ , which leads to complexified texture coefficients. We instead remain in a real tensor space by using the deviatoric bases ${\mathbb{D}}^n$ , and write for arbitrary crystal symmetry the $i$ th texture coefficient of order $n$ as

$${\mathbb{V}}_i^n={\int }_{\mathrm{SO}\left(3\right)}f\left(Q\right)Q\left({\mathbb{D}}_i^n\right)\mathrm{d}V\left(Q\right).$$

(114)

The series description of the ODF follows as

$$f\left(Q\right)=\sum _{n=0}^{\infty }\sum _{i=0}^{2n+1}(2n+1){\mathbb{V}}_i^n·Q\left({\mathbb{D}}_i^n\right).$$

(115)

In micromechanics, a common expression is the expectation value of an orientation-dependent tensor $Q\left({\mathbb{A}}^n\right)$ , which may be computed as an $\mathrm{SO}\left(3\right)$ average weighted by $f$ . Using the harmonic basis,

$${\int }_{\mathrm{SO}\left(3\right)}f\left(Q\right)Q\left({\mathbb{A}}^n\right)\mathrm{dV}\left(Q\right){\int }_{\mathrm{SO}\left(3\right)}f\left(Q\right)\mathit{Q}★\sum _{k=0}^n\sum _i\sum _{j=1}^{2k+1}{A}_{h(k,i,j)}{\mathbb{P}}_i^{\mathrm{nk}}\left[{\mathbb{D}}_j^k\right]\mathrm{d}V\left(Q\right).$$

(116)

After rearranging the integral and using the isotropy of the inclusion tensors, equation (114) can be used to arrive at

$${\int }_{\mathrm{SO}\left(3\right)}f\left(Q\right)Q\left({\mathbb{A}}^n\right)\mathrm{dV}\left(Q\right)=\sum _{k=0}^n\sum _i^{a(n,k)}\sum _{j=1}^{2k+1}{A}_{h(k,i,j)}{\mathbb{P}}_i^{\mathrm{nk}}\left[{\mathbb{V}}_j^k\right].$$

(117)

We note that the harmonic basis components $A_{h(k,i,j)}$ are directly reflected in the calculation of orientation averages.

4.2. Application 3: orientation averages of elasticity tensors

As an example, we calculate the orientation average of a major and minor symmetric fourth-order tensor $\mathbb{A}$ . This operation is useful in calculating effective stiffnesses of arbitrarily anisotropic textured polycrystals even beyond simple averages. For example, Böhlke and Lobos Fernández [30] have calculated Hashin–Shtrikman bounds for cubic polycrystals in this way, with an extension to bounds of polycrystals of arbitrary material symmetry in Lobos Fernández and Böhlke [11]. We call the space of major and minor symmetric fourth-order tensors $\mathcal{C}$ . To determine a full harmonic basis of $\mathcal{C}$ , we need to calculate the harmonic projectors, or equivalently, the harmonic decomposition. The harmonic decomposition of $\mathbb{A}$ is well-known in the literature as

$$\mathbb{A}=a_1\mathbb{P}°+a_2\mathbb{P}\prime +{\mathbb{J}}_1\left[{\mathit{A}\prime }_1]+{\mathbb{J}}_2[{\mathit{A}\prime }_2\right]+\mathbb{A}\prime ,$$

(118)

with appropriate tensors ${\mathbb{J}}_1$ and ${\mathbb{J}}_2$ as discussed, e.g., by Boehler et al. [16]. Differing from most common accounts, our $\mathbb{J}$ are

$${\left({\mathbb{J}}_1\right)}_{\mathrm{ijklmn}}={\delta }_{\mathrm{ij}}P{\prime }_{\mathrm{ijmn}}+P{\prime }_{\mathrm{ijmn}}{\delta }_{\mathrm{kl}},$$

(119)

$${\left({\mathbb{J}}_2\right)}_{\mathrm{abcdef}}=P{\prime }_{\mathrm{abij}}P{\prime }_{\mathrm{cdkl}}({\delta }_{\mathrm{ik}}{P\prime }_{\mathrm{jlef}}+{\delta }_{\mathrm{jl}}{P\prime }_{\mathrm{ikef}}+{\delta }_{\mathrm{il}}{P\prime }_{\mathrm{jkef}}+{\delta }_{\mathrm{jk}}{P\prime }_{\mathrm{ilef}}),$$

(120)

such that

$$\mathbb{J}\left[{\mathit{A}}^{°}\right]=0$$

(121)

and

$${\mathbb{J}}_1\left[\mathit{A}\right]·{\mathbb{J}}_2\left[\mathit{B}\right]=0$$

(122)

for all $\mathit{A}$ and $\mathit{B}$ . In the literature, ${\mathbb{J}}_2$ is usually defined as

$${\left({\mathbb{J}}_2\left[{\mathit{A}\prime }_2\right]\right)}_{\mathrm{ijkl}}={\delta }_{\mathrm{ik}}A{\prime }_{2\mathrm{jl}}+{\delta }_{\mathrm{jl}}A{\prime }_{2\mathrm{ik}}+{\delta }_{\mathrm{il}}A{\prime }_{2\mathrm{jk}}+{\delta }_{\mathrm{jk}}A{\prime }_{2\mathrm{il}},$$

(123)

in which case ${\mathbb{J}}_1\left[\mathit{A}\right]$ and ${\mathbb{J}}_2\left[\mathit{B}\right]$ are not orthogonal if $\mathit{A}$ and $\mathit{B}$ are not deviatoric.

Forte and Vianello [14] state that fourth-order harmonic decompositions are not unique, since tensorial constants of the same order, such as ${\mathbb{J}}_1$ and ${\mathbb{J}}_2$ , may be linearly recombined. By using the Clebsch–Gordan construction of equation (25), our decomposition nonetheless arises uniquely from the second-order harmonic decomposition, meaning that it is the canonical harmonic decomposition for major and minor symmetric tensors which are maps from $\mathrm{Sy}{m}^2$ to $\mathrm{Sy}{m}^2$ , as opposed to other possible interpretations of those tensors, e.g., as maps from $\mathcal{R}$ to ${\mathcal{R}}^{\otimes 3}$ . We therefore reconstruct this harmonic decomposition by constructing $\mathcal{C}$ as ${\mathrm{sym}}^{\mathrm{H}}(\mathrm{Sy}{m}^2\otimes \mathrm{Sy}{m}^2)$ . First, we recall that

$$\mathrm{Sy}{m}^2={\mathbb{P}}^{20}\left[\mathbb{R}\right]+{\mathbb{P}}^{22}\left[{\mathcal{D}}^2\right]$$

(124)

with ${\mathbb{P}}^{20}={\mathit{c}}^{110}=\mathit{I}/\sqrt{3}$ and ${\mathbb{P}}^{22}={\mathbb{C}}^{112}=\mathbb{P}\prime$ . The inclusions from the various deviatoric subspaces to $\mathcal{C}$ follow from equation (38) as

$${\mathbb{P}}_1^{40}=\mathrm{sy}{\mathrm{m}}^{\mathrm{H}}\left(\left({\mathbb{P}}^{20}{\times }^{00}{\mathbb{P}}^{20}\right)\underset{̣}{0}{c}^{000}\right)={\mathbb{P}}^{°},$$

(125)

$${\mathbb{P}}_1^{42}=\mathrm{sy}{\mathrm{m}}^{\mathrm{H}}\left(\left({\mathbb{P}}^{20}{\times }^{02}{\mathbb{P}}^{22}\right)\underset{̣}{2}{\mathbb{C}}^{112}\right)=\frac{1}{2\sqrt{3}}{\mathbb{J}}_1,$$

(126)

$${\mathbb{P}}_2^{40}=\mathrm{sy}{\mathrm{m}}^{\mathrm{H}}\left(\left({\mathbb{P}}^{22}{\times }^{22}{\mathbb{P}}^{22}\right)\underset{̣}{4}{\mathbb{C}}^{220}\right)=\frac{1}{\sqrt{5}}\mathbb{P}\prime ,$$

(127)

$${\mathbb{P}}^{41}=\mathrm{sy}{\mathrm{m}}^{\mathrm{H}}\left(\left({\mathbb{P}}^{22}{\times }^{22}{\mathbb{P}}^{22}\right)\underset{̣}{4}{\mathbb{C}}^{221}\right)=0,$$

(128)

$${\mathbb{P}}_2^{42}=\mathrm{sy}{\mathrm{m}}^{\mathrm{H}}\left(\left({\mathbb{P}}^{22}{\times }^{22}{\mathbb{P}}^{22}\right)\underset{̣}{4}{\mathbb{C}}^{222}\right)=-\sqrt{\frac{3}{28}}{\mathbb{J}}_2,$$

(129)

$${\mathbb{P}}^{43}=\mathrm{sy}{\mathrm{m}}^{\mathrm{H}}\left(\left({\mathbb{P}}^{22}{\times }^{22}{\mathbb{P}}^{22}\right)\underset{̣}{4}{\mathbb{C}}^{223}\right)=0,$$

(130)

$${\mathbb{P}}^{44}=\mathrm{sy}{\mathrm{m}}^{\mathrm{H}}\left(\left({\mathbb{P}}^{22}{\times }^{22}{\mathbb{P}}^{22}\right)\underset{̣}{4}{\mathbb{C}}^{224}\right)={\mathbb{I}}^{\mathrm{De}{v}_4},$$

(131)

with the last term being the eighth-order identity on $\mathrm{De}{v}_4$ . We could recover the exact form of equation (118) by choosing $a_1$ , $a_2$ , ${\mathit{A}\prime }_1$ , ${\mathit{A}\prime }_2$ , and $\mathbb{A}\prime$ to offset the scalar factors. However, the scalar factors ensure that ${\mathbb{P}}^{\mathrm{nk}}\left[{\mathbb{D}}^k\right]$ are normed basis tensors. The resulting symmetric harmonic basis tensors are

$${\mathbb{H}}_1^{\mathcal{C}}={\mathbb{P}}^{°}$$

(132)

$${\mathbb{H}}_2^{\mathcal{C}}=\frac{1}{\sqrt{5}}\mathbb{P}\prime$$

(133)

$${\mathbb{H}}_{3\dots 7}^{\mathcal{C}}=\frac{1}{2\sqrt{3}}{\mathbb{J}}_1\left[{\mathit{D}}_{1\dots 5}^2\right]$$

(134)

$${\mathbb{H}}_{8\dots 12}^{\mathcal{C}}=-\sqrt{\frac{3}{28}}{\mathbb{J}}_2\left[{\mathit{D}}_{1\dots 5}^2\right]$$

(135)

$${\mathbb{H}}_{13\dots 21}^{\mathcal{C}}={\mathit{D}}_{1\dots 9}^4.$$

(136)

According to equation (117), the orientation average of $\mathbb{A}$ resolves to

$${\int }_{\mathrm{SO}\left(3\right)}f\left(Q\right)Q\left(\mathbb{A}\right)\mathrm{dV}\left(Q\right)=\sum _{k=0}^4\sum _{l=1}^{a(4,k)}\sum _{i=1}^{2k+1}{A}_{h(k,l,i)}{\mathbb{P}}_l^{\mathrm{nk}}\left[{\mathbb{V}}_i^k\right],$$

(137)

with $l\in [1,2]$ for $k=0$ and $k=2$ , $l=1$ for $k=4$ and vanishing $l$ for odd $k$ . The projectors ${\mathbb{P}}^{\mathrm{nk}}$ and the enumeration convention $h(k,l,i)$ correspond to the ${\mathbb{H}}^{\mathcal{C}}$ -basis. Because that basis is harmonic, the components $A_h$ of a rotated and averaged $\mathbb{A}$ are precisely the ${\mathbb{H}}^{\mathcal{C}}$ -components of the unrotated $\mathbb{A}$ .

In the general anisotropic case, the $A_{h(k,l,i)}$ form a 21-dimensional vector, which can be reduced depending on the material symmetry. For example, the harmonic decomposition of a cubic fourth-order tensor ${\mathbb{A}}^{\mathrm{c}}$ aligned to the ${\mathit{e}}_i$ -axes is

$${\mathbb{A}}^{\mathrm{c}}=A_1^{\mathrm{c}}{\mathbb{P}}^{°}+A_2^{\mathrm{c}}\mathbb{P}\prime +A_3^{\mathrm{c}}(-\sqrt{\frac{25}{2}}{\mathbb{D}}_1^4-\sqrt{\frac{35}{2}}{\mathbb{D}}_9^4),$$

(138)

where

$$A_1^{\mathrm{c}}={\lambda }_1^{\mathrm{c}},$$

(139)

$$A_2^{\mathrm{c}}=\frac{2}{5}{\lambda }_2^{\mathrm{c}}+\frac{3}{5}{\lambda }_3^{\mathrm{c}},$$

(140)

$$A_3^{\mathrm{c}}=\frac{1}{5}(-{\lambda }_2^{\mathrm{c}}+{\lambda }_3^{\mathrm{c}}).$$

(141)

As noted by Böhlke and Bertram [28], $A_3^{\mathrm{c}}$ is a measure of the crystal’s anisotropy, while the other two parameters represent the expectation values of the isotropic stiffness tensor for a isotropic ODF. As there are three parameters involved in the description of a cubic single crystal, three texture coefficients are sufficient to describe the ODF of a cubic polycrystal according to equation (137). Two of these are the same scalar, and redundant with the normalization condition (equation (113)). Therefore, one fourth-order deviatoric texture tensor fully describes a cubic ODF, which resolves to nine scalar components in the deviatoric basis, as discussed by Lobos Fernández and Böhlke [11]. The deviatoric basis is therefore particularly efficient in texture-based material modeling. Texture evolution models as discussed by, e.g., Böhlke et al. [31] have long used results from tensor representation theory; the harmonic basis is particularly useful because it incorporates those results into an orthonormal basis.

4.3. Application 4: fiber orientation distributions

The results found here are equivalently true for fiber orientation modeling. In the fiber orientation community, textures are modeled over the sphere $S\left(2\right)$ instead of $\mathrm{SO}\left(3\right)$ as fibers are often assumed to be geometrically and materially transversal-isotropic. We focus on two kinds of fiber orientation distribution tensors defined by Kanatani [32]. The orientation tensors of first kind are defined as

$${\mathbb{N}}^n={\int }_{S\left(2\right)}f\left(\mathit{n}\right){\mathit{n}}^{\otimes n}\mathrm{d}V\left(\mathit{n}\right).$$

(142)

To use our results, we transform the integration domain to $\mathrm{SO}\left(3\right)$ , leading to

$${\mathbb{N}}^n={\int }_{\mathrm{SO}\left(3\right)}f\left(Q\right)Q\left({\mathit{e}}_3^n\right)\mathrm{d}V\left(Q\right).$$

(143)

Because ${\mathit{e}}_3^{\otimes n}$ is fully index-symmetric, there is only one inclusion ${\mathbb{P}}^{\mathrm{nk}}$ for each $n$ . In addition, ${\mathit{e}}_3^{\otimes n}$ is transversally isotropic and aligned with the harmonic basis convention axis. Since only the last basis tensor ${\mathbb{D}}_{2k+1}^k$ of ch deviatoric basis is transversally isotropic, the harmonic decomposition of ${\mathit{e}}_3^{\otimes n}$ reads

$${\mathit{e}}_3^{\otimes n}=\sum _{k=0}^{n/2}{H}_k^n{\mathbb{P}}^{n2k}\left[{\mathbb{D}}_{4k+1}^{2k}\right].$$

(144)

The basis tensors are given by

$${\mathbb{H}}_k^n={\mathbb{P}}^{\mathrm{nk}}\left[{\mathbb{D}}_{2k+1}^k\right]=\frac{\mathrm{sym}\left({\mathit{I}}^{\otimes (n/2-k)}\otimes {\mathbb{D}}_{2k+1}^k\right)}{\parallel \mathrm{sym}\left({\mathit{I}}^{\otimes (n/2-k)}\otimes {\mathbb{D}}_{2k+1}^k\right)\parallel },$$

(145)

where the norm resolves to

$$\parallel \mathrm{sym}({\mathit{I}}^{\otimes (n/2-k)}\otimes {\mathbb{D}}_{2k+1}^k)\parallel =\frac{1}{\sqrt{\left(\begin{array}{c}2n-2k\hfill \\ \multicolumn{1}{c}{2k}\end{array}\right)}}.$$

(146)

The harmonic components $H_k^n$ are given by

$$H_k^n=\sqrt{\frac{(2n-4k)!}{(2n-2k)!}}\frac{n!2^k}{(n-2k)!}.$$

(147)

Applying a rotation integral weighted with the ODF $f\left(Q\right)$ over equation (144) leads to

$${\mathbb{N}}^n=\sum _{k=0}^{n/2}{H}_k^n{\mathbb{P}}^{n2k}\left[{\mathbb{V}}_{4k+1}^{2k}\right]$$

(148)

$$=\sum _{k=0}^{n/2}\frac{n!2^k\mathrm{sym}({\mathit{I}}^{\otimes (n/2-k)}\otimes {\mathbb{V}}_{2k+1}^k)}{(n-2k)!\sqrt{\left(2k\right)!}}$$

(149)

The constants $H_k^n$ and ${\mathbb{P}}^{n2k}$ provide an explicit link between texture coefficients ${\mathbb{V}}_{4k+1}^{2k}$ and the first-kind orientation tensor ${\mathbb{N}}^n$ . The texture coefficients ${\mathbb{V}}_{4k+1}^{2k}$ are normalized versions of the orientation tensors of third kind introduced by Kanatani [30], such that

$$\parallel {\mathbb{V}}_{4k+1}^{2k}\parallel \le 1.$$

(150)

Due to the normalization condition equation (113), the zeroth-order third-kind orientation tensor is always $1$ . A suitably minimal basis of the other third-kind tensors is given by the deviatoric bases ${\mathbb{D}}^k$ . The harmonic decomposition shows that a texture is equivalently characterized by either the single first-kind orientation tensor of order $n$ or the third-kind orientation tensors of order up to and including $n$ , where the latter choice requires one fewer parameter due to the normalization condition. In the harmonic basis, this relationship is made explicit.

4.4. Computational effort of tensor rotations

In practice, texture analysis involves large numbers of rotations of high-order tensors. We therefore investigate the computational efficiency of rotations. Given a generally anisotropic $n$ th-order tensor ${\mathbb{A}}^n$ , calculating $\mathit{Q}★{\mathbb{A}}^n$ via the definition of the Rayleigh product in equation (2) requires $n$ multiplications of a $3\times 3$ matrix with a $3^{\times n}$ component array. This resolves to $3^{n-1}n$ matrix-vector products with nine floating-point multiplications each, or $3^{n+1}n$ floating-point multiplications. If we instead use the $n$ th-order harmonic basis, we write ${\mathit{Q}}^{★n}$ as a $3^n\times 3^n$ matrix. Crucially, this matrix is sparse, comprising copies of $d_{\mathrm{ZXZ}}^k$ -matrices up to order $k=n$ , where each matrix has size $(2k+1)\times (2k+1)$ . We may ignore $d$ -matrices of order 0. The total number of floating-point multiplications resolves to

$$\sum _{k=1}^{n}a(n,k){(2k+1)}^2,$$

(151)

where $a(n,k)$ is the number of $2k+1$ -dimensional harmonic subspaces of order $n$ . Because the dimensions of all harmonic subspaces must add up to the dimension of the original tensor space,

$$\sum _{k=0}^{n}a(n,k)(2k+1)=3^n.$$

(152)

Therefore, an upper bound for the total number of floating-point multiplications is

$$\begin{array}{c}\multicolumn{1}{c}{\sum _{k=1}^{n}a(n,k){(2k+1)}^2\le \sum _{k=1}^{n}a(n,k)(2k+1)(2n+1)}\\ \multicolumn{1}{c}{=(2n+1)3^n.}\end{array}$$

(153)

Notably, these are fewer than $3^{n+1}n$ for $n>1$ . When considering total computational effort, floating point additions are roughly as expensive as floating point multiplications. While the number of floating point additions increases with matrix size and is therefore potentially higher when using the harmonic approach, the harmonic approach is still less expensive in total. When the tensor to be rotated has index symmetries, the computational expense of a harmonic rotation can be reduced further, as fewer harmonic subspaces are involved, whereas index symmetries are generally difficult to incorporate into the standard definition of the Rayleigh product. A visualization of exact values is given in Figure 4. While the harmonic basis rotation is generally faster than a naive Rayleigh product, it requires more space, as $n$ matrices of size $(2k+1)\times (2k+1)$ need to be stored as opposed to one $3\times 3$ matrix. Non-harmonic tensorial bases generally require fully occupied $3^n\times 3^n$ rotation matrices, greatly increasing the computational cost compared to even the Rayleigh product.

Figure 4.

Number of additions and multiplications involved in rotations of $n$ th order tensors by using the tensorial Rayleigh product or using a harmonic basis for the spaces ${\mathcal{R}}^n$ , $\mathrm{Sy}{m}^n$ and $\mathcal{C}$ .

If the rotations in question are always applied to a given reference state, as is common for orientations of single crystals in a polycrystal, it is also possible to implement material symmetries to reduce the computational effort further. As is well known in the texture analysis community, material symmetries reduce the number of tensorial texture coefficients involved in a texture description [11]. Equivalently, one might say that material symmetries reduce the $2n+1$ -dimensional deviatoric space ${\mathcal{D}}^n$ to a $k$ -dimensional symmetric subspace. This subspace is however not invariant under arbitrary rotations; for example, the one-dimensional transversal isotropic subspace may change its axis of symmetry under rotations which are not aligned to that axis. Let ${\mathbb{T}}_j^n$ with $j=[1,k]$ denote the $n$ -th order texture coefficient of a materially symmetric single crystal in a reference orientation. Using equation (103), the rotation tensor may be implemented as a $(2n+1)\times k$ matrix given by

$$Q_{\mathrm{ij}}={\mathit{Q}}^{★n}·({\mathbb{D}}_i^n\otimes {\mathbb{T}}_j^n)=d_{\mathrm{ZXZ}\mathrm{il}}^n{\mathbb{D}}_i^n·{\mathbb{T}}_{j.}^n$$

(154)

Applying $Q_{\mathrm{ij}}$ to orient a ${\mathcal{D}}^n$ tensor given in the reference orientation involves $(2n+1)k$ multiplications and $(2n+1)(k-1)$ additions.

For example, fiber orientation tensors are both fully index symmetric and transversal isotropic. In the reference orientation aligned with ${\mathit{e}}_3$ , the only transversally isotropic tensor of the $n$ th-order deviatoric subspace is given by

$${\mathbb{T}}_0^n={\mathbb{D}}_{2n+1}^n.$$

(155)

Correspondingly, rotating a transversally isotropic $n$ th-order tensor such as a fiber material property from a reference configuration by individually rotating all third-kind orientation tensors up to $n$ th order requires no additions and only

$$\sum _{l=1}^{n}2l+1=n^2+2n$$

(156)

multiplications. We summarize that rotating high-order tensors in the harmonic basis is computationally efficient compared to a naive Rayleigh product implementation in the general case, and increasingly specialized rotation implementations can be found for both index- and material symmetries, further increasing the efficiency.

5. Conclusion

This work proposes a novel convention and calculation scheme for harmonic bases which expands and unifies existing schemes in material modelling and texture quantification. The new convention

• simplifies theoretical approaches in crystallography and micromechanics by being aligned to the eigendecompositions of common material symmetries

• reduces computational effort in numerical approaches by allowing a problem-specific choice of basis which, compared to commonly used basis conventions, exploits material symmetry of tensorial quantities to reduce the dimension of equation systems and exploit sparsity of various matrices.

These improvements are accomplished by defining a basis scheme which decomposes tensor spaces into harmonic subspaces. Mathematically, the representation theory of $\mathrm{SO}\left(3\right)$ led us to the deviatoric basis consisting of deviatoric tensors, which was then extended to harmonic bases of arbitrary tensor spaces by using inclusion tensors based on Clebsch-Gordan coefficients. In using these mathematical constructs, our convention builds on group representation theory results well-known in quantum mechanics.

For small deformation material modeling, we derived a harmonic basis of symmetric second-order tensors, which is equivalent to orthonormalized Kocks-type bases. We showed that the basis tensors correspond to elementary strain modes. In linear elasticity theory, stiffness tensors of different material symmetry are efficiently represented in the harmonic basis (Application 1). In addition, we show that for a large set of phenomenological incompressible plasticity models, the time integration problem can also be efficiently formulated and solved in a harmonic basis (Application 2).

In the context of texture analysis, the harmonic basis derived here is a fully real version of the basis recently proposed by Man and Du [12], removing the need for complexified tensor algebra. The deviatoric bases parametrize tensorial texture coefficients in an efficient manner. In addition, the harmonic extension of arbitrary tensor spaces clarifies the relationship between texture coefficients and orientation averages (Application 3). For example, specifying stiffness tensors in a 21-dimensional harmonic basis considerably simplifies the computation of orientation averages if texture coefficients are known. For the description of fiber orientation distributions, e.g., in fiber-reinforced polymers, the harmonic basis tensors are closely related to established fabric or fiber tensors of third kind, and therefore help in describing the orientation distribution (Application 4). Finally, we show how the numerical effort of rotating tensors can be reduced by proper subspace-representations using the harmonic basis.

The harmonic basis formalism described in this paper can be applied to tensors of arbitrary order and symmetry. Correspondingly, harmonic bases can be used in every tensor-based theory to find tensorial bases which are tailored to specific theories and applications.

A Python implementation of the harmonic basis convention defined in this work is available at www.itm.kit.edu/harmonic-bases.