Polyadic random fields

The paper considers mean-square continuous, wide-sense homogeneous, and isotropic random fields taking values in a linear space of polyadics. We find a set of such fields whose values are symmetric and positive-definite dyadics, and outline a strategy for their simulation.


Introduction
Stochastic continuum physics (mechanics, conductivity, electromagnetism, and coupled field) problems hinge on tensor-valued random fields (TRFs). Two types of TRFs are of particular interest: (i) fields of dependent quantities (displacement, velocity, deformation, rotation, stress,. . . ) and (ii) fields of constitutive responses (conductivity, stiffness, permeability. . . ); see also [16]. In the first case, one can work out restrictions on scalar correlation functions which enter the full representations of the tensor-valued correlation structure.
In the second case, there is a need to construct TRFs having a positive-definite property everywhere. This need is especially pronounced in two areas: stochastic partial differential equations and stochastic finite elements. One wants to be able to simulate TRFs as input into stochastic boundary value problems. While a theory of second-order homogeneous and isotropic TRFs has recently been formulated in [15,17], the present study provides a theoretical framework of polyads to meet this goal.
The paper starts by introducing basic concepts of polyads and polyadics and then discusses the symmetry classes for a range of different problems. Armed with this formulation, we develop secondorder TRFs with homogeneous and isotropic correlation polyadics. The paper culminates in an outline of simulation strategy for rank 2 TRFs which are positive-definite. See also [28].
x i y i , (1) and the norm of a vector x ∈ R d is The symbol e i , 1 ≤ i ≤ d, denotes the element of R d whose jth component is equal to the Kronecker delta, δ ij . We have where the first equality means that the set { e i : 1 ≤ i ≤ d } is a basis of the real linear space R d , and the second equality means that the above basis is orthonormal. The symbol r denotes a nonnegative integer. The symbol (R d ) r denotes the linear space R 1 if r = 0 and the Cartesian product of r copies of the set R d otherwise.

Definition 1.
A function r T : (R d ) r → R is called the polyadic (monadic if r = 1, dyadic if r = 2, triadic if r = 3, and so on), if either r = 0 and 0 T is linear or r ≥ 1 and r T is r-linear, that is, for all integers i with 1 ≤ i ≤ d, for all x 1 , x 2 , . . . , x i−1 , x, y, x i+1 , . . . , x r in R d , and for all α, β ∈ R we have r T(x 1 , x 2 , . . . , x i−1 , αx + βy, x i+1 , . . . , x r ) = α r T(x 1 , x 2 , . . . , x i−1 , x, x i+1 , . . . , x r ) The number r is called the rank of a polyadic. Following [3], we denote a polyadic with a capital sans serif letter with a superscript preceding it to indicate its rank, for example, 2 T is a dyadic.
It is easy to see that the set of all polyadics is a real linear space of dimension d r . Indeed, if r = 0, then (R d ) r = R 1 , d r = 1, and a linear functional f : R 1 → R has the form f (x 1 ) = a 1 x 1 with a 1 ∈ R 1 . Otherwise, if r ≥ 1, we define a special polyadic as follows.
It is easy to see that the d r polyads e i1 e i2 · · · e ir with 1 ≤ i j ≤ d for all j with 1 ≤ j ≤ r form a basis in the linear space of all polyadics.
We equip this space by such a scalar product that the above basis becomes orthonormal, and call it the r-dot product. In particular, if r = 1, then the 1-dot product of two monads a and b coincides with the scalar product (1). The founder of polyadic calculus Josiah Willard Gibbs denotes the 2-dot product of two dyads by and the 3-dot product of two triads by He introduced these notions in the paper Elements of Vector Analysis privately printed in New Haven in two volumes in 1881 and 1884 and reprinted in [5, pp. 17-90].

Remark 1.
A reader familiar with tensor calculus may identify a polyad a 1 a 2 · · · a r with the tensor product a 1 ⊗ a 2 ⊗ · · · ⊗ a r and the linear space of all polyadics with the tensor product (R d ) ⊗r of r copies of the space R d . In what follows, we do use the notation (R d ) ⊗r but keep the Gibbs' notation a 1 a 2 · · · a r for polyads and extend it to include the case of polyadics: r T 1 Theorem 1. For r ≥ 1, any r-adics can be represented as a sum of not more than d r−1 r-ads.
Proof. We use mathematical induction. Induction base. Any monadics a is a monad. The statement of the theorem takes the form of a trivial identity a = a. Induction hypothesis. Assume that r ≥ 2 and any (r − 1)-adics has the form and the r-adics r T 1 , . . . , r T d by The sum of all constructed r-adics acts by By the induction hypothesis, we have The terms of this sum are r-ads.

Remark 2.
The representation of Theorem 1 is not unique. Later we will use other presentations.

Symmetry classes
Let I be the d × d identity matrix. The corresponding dyadic is consists of all d × d matrices g with real-valued entries such that gg = I. This group acts on (R d ) ⊗r by g · α = α for r = 0 and g · α i1···ir e i1 · · · e ir = α i1···ir (ge i1 ) · · · (ge ir ). (2) In the above definition of a symmetry class, the group O(d) can be replaced with an arbitrary compact Lie group G and the space U by a finite-dimensional real linear space where an orthogonal representation of G acts. It is well known that the number of symmetry classes is finite, see [2,18,19]. For a particular case of a representation of O(3) in a linear space consisting of tensors, this result is known to physicists as the Hermann theorem [8]. An algorithm for an effective enumeration of symmetry classes for a given finite-dimensional real orthogonal representation of the group O(3) was proposed in [20]. In what follows, we put d = 3.
We introduce a partial order relation on the set of the strata: Σ G2 Σ G1 if and only if G 1 is conjugate to a subgroup in G 2 . It turns out that the partially ordered set of strata has the greatest element, call it the generic stratum, and denote by Σ G0 . It is an open and dense subset of U . The least element of the set of strata is called the minimum stratum. We denote it by Σ GN−1 . The linear space U is partitioned into strata: This partition is called the isotropic stratification of U . For proofs, see [2].
The action (2) on R 3 ⊗ R 3 becomes g · 2 T = g 2 Tg −1 . The subspace U = S 2 (R 3 ) of symmetric dyadics is invariant under that action. The classification of symmetry classes easily follows from elementary linear algebra. If all three eigenvalues of a dyadic 2 T are different, then its symmetry class is [D 2 × Z c 2 ], where D 2 is the dihedral group generated by the two 3 × 3 diagonal matrices with elements −1, −1, 1, and 1, −1, −1, and Z c 2 is the group of order 2 generated by the matrix −I. The corresponding orbit is a prism manifold, see [10,25,26] and modern description in [12,22,27] [20].

Definition 7.
The fixed point set of a group G is the set of all tensors r T ∈ U which are fixed by G: In contrast with a stratum, which is a manifold, a fixed point set is a linear space. Moreover, U G = { r 0} if and only if [G] is a symmetry class. In that case, U G is the biggest linear subspace of U , where G acts trivially: g · r T = r T for all r T ∈ U G and for all g ∈ G.
Consider a group H, which satisfies the following conditions: it is a subgroup of O(3), G is a subgroup of H, and U G is invariant under the action of H. It turns out that these conditions are equivalent to the following: G is a subgroup of H, and H is a subgroup of the normalizer of G: In what follows, we fix the following data: r, U , [G], U G , and H.
A possible physical interpretation of the above data may be as follows, see [1,29]. A continuous medium occupies a compact subset D of the linear space R 3 . A microstructure is attached to any material point x ∈ D. The group G is the group of symmetries of the above microstructure; for simplicity we assume that the medium is homogeneous and this group is the same for all material points. At the macroscopic scale, all the details of the microstructure are lost, all what remains is the material symmetry group G.
On the other hand, the physical properties of the media are encoded by the physical symmetry group, the group H. The Curie principle states that G ⊆ H.
The properties of a medium are encoded in a function r T : D → U G . At the macroscopic scale, this function is deterministic and is usually a solution of a boundary value problem for a system of partial differential equations. At the microscopic scale, this function fails to be deterministic any more and becomes stochastic. It turns out that a relevant description for such a function is a random field.

Random fields
The stochastic function r T : D → U G describes properties of homogeneous random medium. It is convenient to consider this function as the restriction to D of another stochastic function, which maps R 3 to U G . We denote it by the same symbol: r T : Note that the r-dot product induces the norm and the distance on U G in a standard way. Let B(U G ) be the σ-field of all Borel subsets of U G , that is, the minimal σ-field that contains all open subsets of U G . Let (Ω, F, P) be a probability space.
is an event, an element of F.
The random field r T(x) is completely characterized by its finite-dimensional distributions, that is, the (U G ) n -valued random variables where n is a positive integer, and where x 1 , . . . , x n are n distinct points in R 3 . Under a shift, the map Definition 9. A random field r T(x) is called strictly homogeneous if for any positive integer n, for all distinct points x 1 , . . . , x n ∈ R 3 , and for all y ∈ R 3 , the (U G ) n -valued random variables (3) and (4) are equally distributed. What happens under the action of an element g ∈ H? A point x i , 1 ≤ i ≤ n, becomes the point gx i . The random polyadic r T(x) becomes g · r T(gx i ). The finite-dimensional distribution (3) becomes (g · r T(gx 1 ), . . . , g · r T(gx n )).
(5) Definition 10. A random field r T(x) is called strictly isotropic if for any positive integer n, for all distinct points x 1 , . . . , x n ∈ R 3 , and for all g ∈ H, the (U G ) n -valued random variables (3) and (5) are equally distributed.
In what follows, we suppose that a random field where the norm is induced by the r-fold product. Moreover, we suppose that the random field It is easy to see that under a shift, the one-point correlation polyadic and the two-point correlation polyadic of a strictly homogeneous random field do not change.
Definition 11. A second-order random field is called homogeneous if its one-and two-point correlation polyadics do not change under a shift.
The one-point correlation polyadic of the random field g · r T(gx) is To calculate the two-point correlation polyadic of the above random field, define the linear operator ρ(g): There exists a unique linear operator The two-point correlation polyadic of the random field g · r T(gx) becomes This leads to the following definition.

Definition 12.
A second-order random field is called isotropic if its correlation polyadics satisfy It is easy to see that if a second-order random field is strictly homogeneous and/or strictly isotropic, it is homogeneous and/or isotropic. The converse is wrong.
How to find the correlation polyadics of a homogeneous and isotropic random field? By homogeneity, the one-point correlation polyadic is a constant, say r T ∈ U G . By isotropy, the polyadic r T satisfies g · r T = r T for all g ∈ H. That is, the one-point correlation polyadic of a homogeneous and isotropic random field is an arbitrary H-invariant tensor of the space U G .
We give a plan for calculation the two-point correlation polyadic of a homogeneous and isotropic random field. For details, see [15,17] and the references given there. The addition of complex numbers and the multiplication of a complex number by a real one turns the set C of complex numbers into a real two-dimensional linear space. Denote by cU G the complex linear space generated by the dyads z r T, with scalar-vector multiplication z 1 (z 2 r T) = (z 1 z 2 ) r T. There exists a unique r-dot product on cU G with (z 1 r T 1 ) · (z 2 r T 2 ) = z 1 z 2 r T 1 · r T 2 . The map U G → cU G , r T → 1 r T, enables to consider U G as a subset of cU G . Therefore, a random field r T(x) can be considered as taking values in the complex linear space cU G .
Recall that a real structure on a complex finite-dimensional linear space V is a map j : There exists a unique real structure j on cU G with j(z 1 r T 1 + z 2 r T 2 ) = z 1 r T 1 + z 2 r T 2 . In coordinates, the real structure j maps a vector of cU G to the vector with complex-conjugate components.
A centered cU G -valued random field r T(x) is homogeneous if and only if Here,R 3 is a copy of the space domain R 3 , which is called the wavenumber domain or the Fourier space, and F is a measure on the σ-field B(R 3 ) of Borel subsets of the Fourier space taking values in the set of Hermitian nonnegative-definite linear operators on cU G . The map that maps a polyad z 1 r T 1 z 2 r T 2 to the linear operator z r T → (z r T · z 1 r T 1 )z 2 r T 2 , may be extended by linearity to an isomorphism between (cU G ) ⊗2 and the linear space L(cU G ) of all linear operators in cU G . There exists a unique real structure on the complex linear space L(cU G ), call it j ⊗ j, that maps a linear operator z r T → (z r T · z 1 r T 1 )z 2 r T 2 to the linear operator z r T → (z r T · jz 1 r T 1 )jz 2 r T 2 . In coordinates, the real structure j ⊗ j maps a matrix with complex entries to the matrix with complexconjugate entries. A Hermitian matrix is mapping to the transposed matrix.
It is easy to prove that a cU G -valued homogeneous random field r T( where dg is the unique probabilisticH-invariant measure onH known as the Haar measure, and where dν n (λ) is an arbitrary finite measure on (R 3 /H) n . Equation (6) takes the form For g ∈ H n , we have go n = o n and f (o n , λ) =ρ(g) f (o n , λ). It follows that f (o n , λ) takes values in the convex compact set C n = C ∩Ũ n , whereŨ n is the maximal linear subspace of the spaceŨ where the representationρ of the group H n acts trivially. The next step is to calculate the inner integral in the right-hand side of Eq. (7). The function g → e i(gon,λ)·(x−y) is continuous onH and constant on the left cosets ofH with respect to H n . By the Fine Structure Theorem [9], for any integer n, 0 ≤ n ≤ N − 1, there exists a subset Irr(H) n of the set Irr(H) of equivalence classes of irreducible orthogonal representations of G, and a subset A nρ of the set {1, 2, . . . , dim ρ} such that the matrix entries { ρ ij (g): ρ ∈ Irr(H) n , 1 ≤ i ≤ dim ρ, j ∈ A nρ } form an orthogonal basis in the space of complex-valued functions onH which are constants on the left cosets of H with respect to H n and square-integrable with respect to the Haar measure. The Fourier coefficients of the above continuous function have the form The Fourier series of a continuous function converges uniformly onH. We may substitute this series to Eq. (7) and interchange integration and summation. We obtain There is a basis in the subspaceŨ , in which the representationρ breaks into irreducible components. In this basis, the integral overH of any matrix entry is equal either to 0 or to ρ ij (g) 2 . Instead of going into details, we consider an example.
We determine the structure of the representationρ. First, following [7], we describe the representatives of equivalence classes of irreducible orthogonal representations of O(3). For a nonnegative integer , the representation ρ acts in the linear space H of homogeneous harmonic (with null Laplacian) polynomials p(x) of degree in three real variables by [ρ (p)](x) = p(g −1 x), while the representation ρ * acts in the same space by [ρ (p)](x) = det gp(g −1 x). Each irreducible representation is equivalent to one of the above. Note that the representations ρ with even and ρ * with odd sent −I ∈ O(3) to the identity operator, while the remaining representations send it to the minus identity. The Clebsch-Gordan rule tells that the tensor product of the representations with indices 1 and 2 is equivalent to the direct ZAMP Polyadic random fields Page 9 of 21 204 sum of representations with indices , | 1 − 2 | ≤ ≤ 1 + 2 . Later we will see how to determine which components have a star as an upper index. The representation g → g sends −I to −I and therefore is equivalent to ρ 1 . The tensor square ρ 1 ⊗ ρ 1 and all its irreducible components send −I to −I ⊗ −I = I. The Clebsch-Gordan rule gives For complex representations, the orthonormal bases in the spaces of irreducible components are calculated with the help of the so called Clebsch-Gordan coefficients which can be found in any book on quantum mechanics. For real representations, they have been calculated in [6]; we call them the Godunov-Gordienko coefficients. If the vectors e m constitute an orthonormal basis in the space of the representation ρ , then we have We use an algorithm for their calculations proposed in [24].
The first tetradics of the above basis is 4 T 1 = 2 T 1 2 T 1 that corresponds to the component ρ 0 ⊗ ρ 0 of the representationρ.
The intersection of symmetric part of the component ρ 0 ⊗ ρ 2 ⊕ ρ 2 ⊗ ρ 0 withŨ 0 is generated by the Finally, the symmetric part of the tensor product ρ 2 ⊗ ρ 2 is given by ρ 0 ⊕ ρ 2 ⊕ ρ 4 . The component ρ 0 contributes by and the component ρ 4 by In the basis (9), a generic element of the linear spaceŨ 0 takes the form of the 6 × 6 matrix f with the following nonzero entries: The matrix f must be nonnegative-definite and have unit trace. Moreover, it is block diagonal and has 3 different blocks: f 33 0 0 f 33 , and f 55 0 0 f 55 . Accordingly, the set of extreme points of the convex compact set C 0 contains 3 connected components. The first component is an ellipse
We have Let (λ,θ,φ) be the spherical coordinates inR 3 . The representation ρ 0 acts by is the real-valued spherical harmonic. The representation ρ 2 acts by where we renamed the dyadics (9) according to the scheme 2 T 2 → 2 T 0 , 2 T 3 → 2 T −1 , 2 T 4 → 2 T 1 , 2 T 5 → 2 T −2 , 2 T 6 → 2 T 2 to be consistent with notation of [6,24]. The representation ρ 4 acts by Calculations give the following nonzero elements of the 6 × 6 symmetric matrixρ(g)A 1 that are located on and over the main diagonal: Similarly, Calculations give the following nonzero elements of the 6 × 6 symmetric matrixρ(g)A 2 that are located on and over the main diagonal: Finally, Calculations give the following nonzero elements of the 6×6 symmetric matrixρ(g)A 3 (λ) that are located on and over the main diagonal: We substitute this expansion to Eq. (11) and obtain where the functions J ijk (λ, r, θ, ϕ) are symmetric in the first two indices and their nonzero values are given in Table 1. In Eq. (12), we expressed the spherical Bessel functions in term of ordinary ones and denoted by (r, θ, ϕ) (resp. (λ,θ,φ)) the spherical coordinates of the point x − y (resp. p).
To simulate a random field 2 T(x), we apply the Rayleigh expansion (12) to the plane waves e ip·x and e −ip·y separately, see details in [14,15]. In the next section, we propose an alternative method for simulation.

A strategy for simulation
We would like to develop a strategy for simulation of a homogeneous and isotropic random field 2 T(x) that takes values in the set of symmetric and positively-definite matrices. We proceed similarly to the case of rank 0. Let 0 T(x) be a centered Gaussian homogeneous and isotropic R 1 -valued random field with unit variance. The random field 0 T(x) 0 T(x) is homogeneous, isotropic, and takes positive values. Such a field is called a field of Chi-square type, because its one-dimensional distributions are Chi-square.
Let 1 T(x) be a homogeneous and isotropic R 3 -valued random field. It is necessarily centered. Its twopoint correlation function may be calculated similarly to the case of rank 2. We give the answer. An equation similar to (11) has the form 1 T(x), 1 T(y) = e ip·(x−y)ρ (g)A k dg dΦ k (λ) (13)