D Marinucci

G Peccati

Representations of SO(3) and angular polyspectra

Keywords: Group Representations, Isotropy, Polyspectra, Spherical Random Fields AMS 2000 Classification. 60G10, 60G35, 20C12, 20C35

We characterize the angular polyspectra, of arbitrary order, associated with isotropic fields defined on the sphere S 2 = (x, y, z) : x 2 + y 2 + z 2 = 1 . Our techniques rely heavily on group representation theory, and specifically on the properties of Wigner matrices and Clebsch-Gordan coefficients. The findings of the present paper constitute a basis upon which one can build formal procedures for the statistical analysis and the probabilistic modelization of the Cosmic Microwave Background radiation, which is currently a crucial topic of investigation in cosmology. We also outline an application to random data compression and "simulation" of Clebsch-Gordan coefficients.

Introduction

The connection between probability theory and group representation theory has led to a long tradition of fruitful interactions. A well-known reference is provided by [START_REF] Diaconis | Group Representations in Probability and Statistics[END_REF]; see e.g. [2, Section 40-41], [START_REF] Diaconis | The Poisson-Dirichlet law is the unique invariant distribution for uniform split-merge transformations[END_REF], [START_REF] Fulman | Convergence rates of random walk on irreducible representations of finite groups[END_REF], [START_REF] Guivarc'h | Marches Aléatoires sul les Groupes de Lie[END_REF], [START_REF] Pycke | A decomposition for invariant tests of uniformity on the sphere[END_REF], [START_REF] Raimond | Flots browniens isotropes sur la sphère[END_REF], [START_REF] Yadrenko | Spectral Theory of Random Fields[END_REF], and the references therein, for other relevant contributions. In this paper we shall focus in particular on the profound connection between the probabilistic notion of isotropy, i.e. invariance in law under the action of a group, and the representation theory of the group itself. One instance of this connection is well-known, i.e. the celebrated Peter-Weyl Theorem, which allows the construction of spectral representations for isotropic random fields on homogeneous spaces of general compact groups, see [START_REF] Peccati | Decompositions of stochastic processes based on irreducible group representations[END_REF] for a general construction and [START_REF] Marinucci | High-frequency asymptotics for subordinated stationary fields on an Abelian compact group[END_REF] , [START_REF] Marinucci | Group representations and high-resolution central limit theorems for subordinated spherical random fields[END_REF] for examples related, respectively, to the torus and the sphere. Our aim here is to use these representations in order to characterize random fields by means of a higher order spectral theory; in particular, one of our main goals will be to establish the link between the so-called polyspectra (or higher order spectra) and alternative (tensor product and direct sum) representations of the underlying isotropy group. In particular, we shall provide a general expression for higher order spectra of isotropic spherical random fields in terms of convolutions of Clebsch-Gordan or Wigner coefficients. The latter where introduced in Mathematics in the XIX century for the analysis of Algebraic Invariants; they have since then played a crucial role in the development of Quantum Physics in the XX century (see for instance [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF] for a comprehensive reference); their role in Group Representation theory will be discussed below, while more details can be found for instance in [START_REF] Vilenkin | Representation of Lie Groups and Special Functions[END_REF].

Our analysis may have an intrinsic mathematical interest, but it is also strongly motivated by applications to Physics and Cosmology. Concerning the latter, the analysis of higher order spectra for isotropic spherical random fields is currently at the core of several research efforts which are related to the analysis of Cosmic Microwave Background (CMB) radiation data, see for instance [START_REF] Dodelson | Modern Cosmology[END_REF] for a general introduction and [START_REF] Hu | The Angular trispectrum of the CMB[END_REF][START_REF] Komatsu | Acoustic Signatures in the Primary Microwave Background Bispectrum[END_REF][START_REF] Marinucci | High-resolution asymptotics for the angular bispectrum of spherical random fields[END_REF][START_REF] Marinucci | A Central Limit Theorem and Higher Order Results for the Angular Bispectrum[END_REF] for some references on the bi-and trispectrum. A general characterization of the theoretical properties of higher order angular power spectra can yield several insights into the statistical analysis of the massive datasets that are or will be made available by satellite experiments such as WMAP or Planck. For instance, the current understanding of the behaviour of the bispectrum for some simple physical models has already led to many applications ( [START_REF] Cabella | The integrated bispectrum as a test of cosmic microwave background non-Gaussianity: detection power and limits on f N L with WMAP data[END_REF], [START_REF] Yadav | Detection of Primordial non-Gaussianity (f N L ) in the WMAP 3-Year Data at above 99.5% Confidence[END_REF], [START_REF] Yadav | Fast Estimator of Primordial Non-Gaussianity from Temperature and Polarization Anisotropies in the Cosmic Microwave Background II: Partial Sky Coverage and Inhomogeneous Noise[END_REF]), aiming at obtaining constraints on nonlinearity parameters of utmost physical significance; needless to say, a proper understanding of higher order spectra can lead to more efficient statistical procedures and better constraints, which may help to solve some of the important scientific issues at stake in CMB analysis (primarily a proper understanding of the Big Bang inflationary dynamics, which is tightly linked with the CMB nonlinear structure, see [START_REF] Dodelson | Modern Cosmology[END_REF], [START_REF] Babich | The Shape of non-Gaussianities[END_REF], [START_REF] Bartolo | Non-Gaussianity from Inflation: Theory and Observations[END_REF], [START_REF] Maldacena | Non-Gaussian features of primordial fluctuations in single field inflationary models[END_REF]).

The relevance of the current results need not be limited to cosmological applications. Indeed, the analysis of spherical random fields has currently led to remarkable developments in the Geophysical and Planetary Sciences, and even in Medical Imaging, see for instance ( [START_REF] Chung | Tensor-based cortical surface morphometry via weighted spherical harmonics representation[END_REF], [START_REF] Simons | Spatiospectral concentration on a sphere[END_REF], [START_REF] Wieczorek | Minimum-variance multitaper spectral estimation on the sphere[END_REF]). Moreover, we shall show below how the relationship which we establish leads very naturally to some numerical algorithms for the estimation of Clebsch-Gordan and Wigner coefficients. The latter represent probability amplitudes of quantum interactions and as such a rich literature in Mathematical Physics has been concerned with recipes for their numerical estimation: our procedure lends itself to easy implementation and can be simply extended to very general compact groups, although in this paper we focus solely on SO [START_REF] Babich | The Shape of non-Gaussianities[END_REF].

The plan of this paper is as follows: in Section 2 we introduce our general probabilistic setting and provide some preliminary notation and background material. In Section 3 we present some background material on representation theory, while in Section 4 and Section 5 we obtain our main results, including the explicit characterization of polyspectra. These results are applied in Section 6 to derive explicit expressions in some important cases (such as χ 2 random fields). Section 7 is devoted to further issues that we see as the seed for future research: they concern, in particular, the connection with the representation theory for the symmetric group and the Monte Carlo estimation of Clebsch-Gordan coefficients.

In the subsequent sections, every random element is defined on an appropriate probability space (Ω, F , P ).

General setting

In this paper, we focus on real-valued, centered, square-integrable and isotropic random fields on the sphere S 2 = (x, y, z) ∈ R 3 : x 2 + y 2 + z 2 = 1 . A centered and square integrable random field T on S 2 is just a collection of random variables of the type T = T (x) : x ∈ S 2 such that, for every x ∈ S 2 , ET (x) = 0 and ET 2 (x) < ∞. In the following, whenever we write that T is a field on S 2 , we will implicitly assume that T is real-valued, centered and square-integrable. From now on, we shall distinguish between two notions of isotropy, which we name strong isotropy and weak isotropy of order n (n ≥ 2).

Strong isotropy -The field T is said to be strongly isotropic if, for every k ∈ N, every x 1 , ..., x k ∈ S 2 and every g ∈ SO(3) (the group of rotations in R 3 ) we have Weak isotropy -The field T is said to be n-weakly isotropic (n ≥ 2) if E|T (x)| n < ∞ for every x ∈ S 2 , and if, for every x 1 , ..., x n ∈ S 2 and every g ∈ SO(3),

E [T (x 1 ) × • • • × T (x n )] = E [T (gx 1 ) × • • • × T (gx n )] .
The following statement, whose proof is elementary, indicates some relations between the two notions of isotropy described above.

Proposition 1

1. A strongly isotropic field with finite moments of some order n ≥ 2 is also n-weakly isotropic.

2. Suppose that the field T is n-weakly isotropic for every n ≥ 2 (in particular, E|T (x)| n < ∞ for every n ≥ 2 and every x ∈ S 2 ) and that, for every k ≥ 1 and every (x 1 , ..., x k ), the law of the vector {T (x 1 ) , ..., T (x k )} is determined by its moments. Then, T is also strongly isotropic.

Now suppose that T is a strongly isotropic field, and denote by dx the Lebesgue measure on S 2 . Since the variance ET (x)

2 is finite and independent of x (by isotropy), one deduces immediately that

E S 2 T (x) 2 dx < ∞,
from which one infers that the random path x → T (x) is a.s. square integrable with respect to the Lebesgue measure. Then, it is a standard result that the following spectral representation holds:

T (x) = ∞ l=0 l m=-l a lm Y lm (x) , where a lm S 2 T (x) Y lm (x)dx, (2.2)
and where the complex-valued functions {Y lm : l ≥ 0, m = -l, ..., l} are the so-called spherical harmonics, to be defined below. The spectral representation (2.2) must be understood in the L 2 (Ω × S 2 ) sense, i.e.

lim

L→∞ E T - L l=0 l m=-l a lm Y lm 2 L 2 (S 2 ) = 0,
where L 2 (S 2 ) is the complex Hilbert space of functions on S a lm Y lm (x) = 0 for all x ∈ S 2 , a.s.-P , see for instance [START_REF] Adler | Random Fields and Geometry[END_REF] or [START_REF] Yadrenko | Spectral Theory of Random Fields[END_REF]. The spherical harmonics {Y lm } m=-l,...,l are the eigenfunctions of the Laplace-Beltrami operator on the sphere, denoted by ∆ S 2 , satisfying the relation ∆ S 2 Y lm = -l(l + 1)Y lm . These functions can be represented by means of spherical coordinates x = (θ, ϕ) as follows:

Y lm (θ, ϕ) = 2l + 1 4π (l -m)! (l + m)! P lm (cos θ) exp(imϕ) , for m > 0 , Y lm (θ, ϕ) = (-1) m Y l,-m (θ, ϕ) , for m < 0 , 0 ≤ θ ≤ π, 0 ≤ ϕ < 2π ,
where P lm (cos θ) denotes the associated Legendre polynomial of degree l, m, i.e.

P lm (x) = (-1) m (1 -x 2 ) m/2 d m dx m P l (x) , P l (x) = 1 2 l l! d l dx l (x 2 -1) l , m = 0, 1, 2, ..., l , l = 0, 1, 2, 3, .... .
The random spherical harmonics coefficients {a lm } appearing in (2.2) form a triangular array of zero-mean and square-integrable random variables, which are complex-valued for m = 0 and such that Ea lm a l ′ m ′ = δ l′ l δ m′ m C l , the bar denoting complex conjugation. Here, and for the rest of the paper, the symbol δ a b is equal to one if a = b and zero otherwise. We also write C l = E |a lm | 2 , l ≥ 0, to indicate the angular power spectrum of T (we stress that the quantity C l does not depend on m -see e.g. [START_REF] Baldi | Some characterizations of the spherical harmonics coefficients for isotropic random fields[END_REF] for a proof of this fact). Observe that, by definition of the spherical harmonics, a lm = (-1) m a l-m . Note also that a convenient route to derive (2.2) is by means of an appropriate version of the stochastic Peter-Weyl theorem -see for instance [START_REF] Baldi | On the characterization of isotropic random fields on homogeneous spaces of compact groups[END_REF] or [START_REF] Peccati | Decompositions of stochastic processes based on irreducible group representations[END_REF], as well as Section 3.1 below.

Observe that the representation (2.2) still holds for fields {T (x)} that are not necessarily isotropic, but such that the random path x → T (x) is P -a.s. square integrable with respect to the Lebesgue measure dx. Indeed, if the last property holds, then one has that, P -almost surely, lim

L→∞ S 2 T (x) - L l=0 l m=-l a lm Y lm (x) 2 dx = 0. (2.3)
In this case, however, none of the previously stated properties on the array {a lm } holds in general. By an argument similar to those displayed above, a sufficient condition to have that x → T (x) is P -a.s. Lebesgue-square integrable is that sup

x∈S 2 ET (x) 2 < ∞.
The next result, that we record for future reference, is proved in [START_REF] Baldi | Some characterizations of the spherical harmonics coefficients for isotropic random fields[END_REF].

Proposition 2 Let T be a centered, square-integrable and strongly isotropic random field. Let the coefficients {a lm } be defined according to (2.2). Then, for every l, m, one has that E |a lm | 2 < ∞. Moreover, for every l ≥ 1, the coefficients {a l0 , ..., a ll } are independent if, and only if, they are Gaussian. If the vector {a l0 , ..., a ll } is Gaussian, one also has that ℜ (a lm ) and ℑ (a lm ) are independent and identically distributed for every fixed m = 1, ..., l (ℜ (z) and ℑ (z) stand, respectively, for the real and imaginary parts of z).

The following result formalizes the fact that, in general, one cannot deduce strong isotropy from weak isotropy. The proof makes use of Proposition 1.

Proposition 3 For every n ≥ 2, there exists a n-weakly isotropic field T such that T is not strongly isotropic.

Proof. Fix l ≥ 1, and consider a vector b m , m = -l, ..., l, of centered complex-valued random variables such that:

(i) b 0 is real, (ii) b -m = (-1)
m b m (m = 1, ..., l), (iii) the vector {b 0 , ..., b l } is not Gaussian and is composed of independent random variables, (iv) for every k = 1, ..., n, the (possibly mixed) moments of order k of the variables {b 0 , ..., b l } coincide with those of a vector {a 0 , ..., a l } of independent, centered and complex-valued Gaussian random variables with common variance C l and such that a 0 is real and, for every m = 1, ..., l, the real and imaginary parts of a m are independent and identically distributed (the existence of a vector such as {b 0 , ..., b l } is easily proved). Now define the two fields

T (x) = l m=-l b m Y lm (x) and T * (x) = l m=-l a m Y lm (x) .
By Proposition 2, T * is strongly isotropic, and also n-weakly isotropic by Proposition 1. By construction, one also has that T is n-weakly isotropic. However, T cannot be strongly isotropic, since this would violate Proposition 2 (indeed, if T was isotropic, one would have an example of an isotropic field whose harmonic coefficients {b 0 , ..., b l } are independent and non-Gaussian).

In what follows, we use the symbol A ⊗ B to indicate the Kronecker product between two matrices A and B. Given n ≥ 2, we denote by Π (n) the class of partitions of the set {1, ..., n}. Given an element π ∈ Π (n), we write π = {b 1 , ..., b k } to indicate that the sets b j ⊆ {1, ..., n}, j = 1, ..., k, are the blocks of π. The blocks of a partition are always listed according to the lexicographic order, that is: the block b 1 always contains 1, the block b 2 contains the least element of {1, ..., n} not contained in b 1 , and so on. Also the elements within each block b j are written in increasing order. For instance, if a partition π of {1, ..., 5} is composed of the blocks {1, 3} , {5, 4} and {2}, we will write π in the form π = {{1, 3} , {2} , {4, 5}} . Definition A. (A1) Let the field T admit the representation (2.2), and suppose that, for some n ≥ 2, one has that E |a lm | n < ∞ for every l, m. Then, T is said to have finite spectral moments of order n. (A2) Suppose that T has finite spectral moments of order n ≥ 2, and, for l ≥ 0, use the notation a l. = (a l-l , ..., a l0 , ..., a ll ) .

(2.4)

The polyspectrum of order n -1, associated with T , is given by the collection of vectors

S l1...ln = E [a l1. ⊗ a l2. ⊗ • • • ⊗ a ln. ] , (2.5)
where 0 ≤ l 1 , l 2 , ..., l n . Note that the vector S l1...ln appearing in (2.5) has dimension (2l 1 + 1) × • • • × (2l n + 1). (A3) Suppose that T has finite spectral moments of order n ≥ 2. The (mixed) cumulant polyspectrum of order n -1, associated with T , is given by the vectors

S c l1...ln = π={b1,...,b k }∈Π(n) (-1) k-1 (k -1)!E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ] , (2.6)
where 0 ≤ l 1 , l 2 , ..., l n , and, for every block b j = {i 1 , ..., i p }, we use the notation

E ⊗ i∈bj a li. = E a li 1 . ⊗ • • • ⊗ a li p .
(recall that we always list the elements of b j in such a way that i 1 ≤ • • • ≤ i p ). Plainly, the vector S c l1...ln in (2.6) has also dimension (2l

1 + 1) × • • • × (2l n + 1) .
Remark. Suppose that T has finite spectral moments of order n ≥ 2. Then, by selecting frequencies

l 1 = l 2 = • • • = l 3 = l ≥ 0, one obtains that S c l...l n times := S c l...l (n) = π={b1,...,b k }∈Π(n) (-1) k-1 (k -1)!E (a l. ) ⊗|b1| ⊗ • • • ⊗ E (a l.
)

⊗|b k | (2.7)
where |b j | stands for the size of the block b j , and we use the notation

(a l. ) ⊗|bj | = a l. ⊗ • • • ⊗ a l. |bj | times .
3 Preliminary material

Representation Theory for SO(3)

We start by reviewing briefly some background material on the special group of rotations SO(3), i.e. the space of 3 × 3 real matrices A such that A ′ A = I 3 (the three-dimensional identity matrix) and det(A) = 1. We first recall that each element g ∈ SO(3) can be parametrized by the set (ϕ, ϑ, ψ) of the so-called Euler angles (0 ≤ ϕ < 2π, 0 ≤ ϑ ≤ π, 0 ≤ ψ < 2π); indeed each rotation in R 3 can be realized sequentially as

A = A(g) = R(ψ, ϑ, ϕ) = R z (ψ)R x (ϑ)R z (ϕ) (3.8)
where R z (ϕ), R x (ϑ), R z (ψ) ∈ SO(3) can be expressed by means of the following general definitions, valid for every angle α,

R z (α) =   cos α -sin α 0 sin α cos α 0 0 0 1   , R x (α) =   1 0 0 0 cos α -sin α 0 sin α cos α   .
The representation (3.8) is unique except for ϑ = 0 or ϑ = π, in which case only the sum ϕ + ψ is determined. In words, the rotation is realized by rotating first by ϕ around the axis z, then rotating around the new x axis by ϑ, then rotating by ψ around the new z axis. It is clear that the first two rotations identify one point on the sphere, so the whole operation could be also interpreted as moving the North Pole to a new orientation in S 2 and then rotating by ψ the tangent plane at the new location.

In these coordinates, a complete set of irreducible matrix representations for SO(3) is provided by the Wigner's D matrices D l (ψ, ϑ, ϕ) = D l mn (ψ, ϑ, ϕ) m,n=-l,...,l , of dimensions (2l + 1) × (2l + 1) for l = 0, 1, 2, ...; we refer to classical textbooks, such as [START_REF] Vilenkin | Representation of Lie Groups and Special Functions[END_REF], [START_REF] Bump | Lie Groups[END_REF] or [START_REF] Diaconis | Group Representations in Probability and Statistics[END_REF], for any unexplained definition or result concerning group representation theory. An analytic expression for the elements of Wigner's D matrices is provided by D l mn (ψ, ϑ, ϕ) = e -inψ d l mn (ϑ)e imϕ , m, n = -(2l + 1) , ..., 2l + 1

where the indices m, n indicate, respectively, columns and rows, and

d l mn (ϑ) = (-1) l-n [(l + m)!(l -m)!(l + n)!(l -n)!] 1/2 × k (-1) k cos ϑ 2 m+n+2k sin ϑ 2 2l-m-n-2k k!(l -m -k)!(l -n -k)!(m + n + k)! ,
and the sum runs over all k such that the factorials are non-negative; see [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF]Chapter 4] for a huge collection of alternative expressions. Here we simply recall that the elements of D l (ψ, ϑ, ϕ) are related to the spherical harmonics by the relationship

D l 0m (ϕ, ϑ, ψ) = (-1) m 4π 2l + 1 Y l-m (ϑ, ϕ) = 4π 2l + 1 Y * lm (ϑ, ϕ) . (3.9)
In other words, the spherical harmonics correspond (up to a constant) to the elements of the "central" column in the Wigner's D matrix. Such matrices operate irreducibly and equivalently on (2l + 1) spaces (the so-called isotypical spaces), each of them spanned by a different column n of the matrix representation itself. The elements of column n correspond to the so-called spin n spherical harmonics, which enjoy a great importance in particle physics and in harmonic expansions for tensor valued random fields. In this paper, we restrict our attention only to the usual n = 0 spherical harmonics, which correspond to usual scalar functions.

Remark. By exploiting relation (3.9), it is not difficult to show that the usual spectral representation for random fields on the sphere, as given in (2.2), is just the stochastic Peter-Weyl Theorem on the quotient space S 2 = SO(3)/SO [START_REF] Bump | Lie Groups[END_REF]. Indeed, by the stochastic Peter-Weyl Theorem (see e.g. [START_REF] Peccati | Decompositions of stochastic processes based on irreducible group representations[END_REF]) we obtain, for any square integrable, isotropic random field {T (g) :

g ∈ SO(3)} T (g) = T (ϕ, ϑ, ψ) = l m,n a lmn 2l + 1 8π 2 D l mn (ϕ, ϑ, ψ) ,
where dg is the Haar (uniform) measure on SO (3) with total mass 8π 2 . Now if we consider the restriction of T (g) to S 2 = SO(3)/SO(2), denoted by T S 2 (ϕ, ϑ), we deduce that

a lmn = SO(3) T S 2 (g) 2l + 1 8π 2 D l mn (g)dg = S 2 T S 2 (ϕ, ϑ) 2π 0 e inψ dψ 2l + 1 8π 2 d l mn (ϑ)e -imϕ sin ϑdϕdϑ , = S 2 T S 2 (ϕ, ϑ)δ 0 n (2π)
2l + 1 8π 2 d l mn (ϑ)e -imϕ sin ϑdϕdϑ , the second equality following from the fact that T S 2 (g) is constant with respect to ψ. We can thus conclude that

a lmn = 0 for n = 0 √ 2πa lm for n = 0 ,
where the array {a lm } is defined by (2.2).

The Clebsch-Gordan matrices

It follows from standard representation theory that we can exploit the family D l l=0,1,,2,... to build alternative (reducible) representations, either by taking the tensor product family D l1 ⊗ D l2 l1,l2 , or by considering direct sums ⊕ l2+l1 l=|l2-l1| D l l1,l2

. These representations have dimensions

(2l 1 + 1)(2l 2 + 1) × (2l 1 + 1)(2l 2 + 1)
and are unitarily equivalent, whence there exists a unitary matrix C l1l2 such that

D l1 ⊗ D l2 = C l1l2 ⊕ l2+l1 l=|l2-l1| D l C * l1l2 . (3.10)
The matrix C l1l2 is a {(2l 1 + 1)(2l 2 + 1) × (2l 1 + 1)(2l 2 + 1)} block matrix, whose blocks, of dimensions (2l 2 + 1) × (2l + 1), are customarily denoted by C l l1(m1)l2 , m 1 = -l 1 , ..., l 1 ; the elements of such a block are indexed by m 2 (over rows) and m (over columns; note that m = -(2l + 1), ..., 2l + 1). More precisely,

C l1l2 = C l. l1(m1)l2. m1=-l1,...,l1;l=|l2-l1|,...,l2+l1 (3.11) C l. l1(m1)l2.

=

C lm l1m1l2m2 m2=-l2,...,l2;m=-l,...,l .

(3.12)

Remark. The fact that the two matrices D l1 ⊗ D l2 and ⊕ l2+l1 l=|l2-l1| D l have the same dimension follows from the elementary relation (valid for any integers l 1 , l 2 ≥ 0):

l1+l2 l=|l2-l1| (2l + 1) = (2l 1 + 1) (2l 2 + 1) . (3.13)
By induction, one also obtains that, for every n ≥ 3,

l1+l2 λ1=|l2-l1| λ1+l3 λ2=|l3-λ1| • • • λn-2+ln λn-1=|ln-λn-2| (2λ n-1 + 1) = n j=1 (2l j + 1) , (3.14)
for any integers l 1 , ..., l n ≥ 0 (relation (3.14) is needed in Section 5.2).

The Clebsch-Gordan coefficients for SO(3) are then defined as the collection C lm l1m1l2m2 of the the elements of the unitary matrices C l1l2 . These coefficients were introduced in Mathematics in the XIX century, as motivated by the analysis of invariants in Algebraic Geometry; in the 20th century, they have gained an enormous importance in the quantum theory of angular momentum, where C lm l1m1l2m2 represents the probability amplitude that two particles with total angular momentum l 1 , l 2 and momentum projection on the z-axis m 1 and m 2 are coupled to form a system with total angular momentum l and projection m (see e.g. [START_REF] Liboff | Introductory Quantum Mechanics[END_REF]). Their use in the analysis of isotropic random fields is much more recent, see for instance [START_REF] Hu | The Angular trispectrum of the CMB[END_REF] and the references therein.

Remark (More on the structure of the Clebsch-Gordan matrices). To ease the reading of the subsequent discussion, we provide an alternative way of building a Clebsch-Gordan matrix C l1l2 , starting from any enumeration of its entries. Fix integers l 1 , l 2 ≥ 0 such that l 1 ≤ l 2 (this is just for notational convenience), and consider the Clebsch-Gordan coefficients C lm l1m1l2m2 given in (3.11)-(3.12). According to the above discussion, we know that:

(i) -l i ≤ m i ≤ l i for i = 1, 2, (ii) l 2 -l 1 ≤ l ≤ l 1 + l 2 , (iii) -l ≤ m ≤ l,
and (iv) the symbols (l 1 , m 1 , l 2 , m 2 ) label rows, whereas the pairs (l, m) are attached to columns. Now introduce the total order ≺ c on the "column pairs" (l, m), by setting that (l, m) ≺ c (l ′ , m ′ ), whenever either l < l ′ or l = l ′ and m < m ′ . Analogously, introduce a total order ≺ r over the "row symbols" (l 1 , m 1 , l 2 , m 2 ), by

setting that (l 1 , m 1 , l 2 , m 2 ) ≺ r (l ′ 1 , m ′ 1 , l ′ 2 , m ′ 2 ), if either m 1 < m ′ 1 , or m 1 = m ′ 1 and m 2 < m ′ 2 (
recall that l 1 and l 2 are fixed). One can check that the set of column pairs (resp. row symbols) can now be written as a saturated chain 1 with respect to ≺ c (resp. ≺ r ) with a least element given by (l

2 -l 1 , -(l 2 -l 1 )) (resp. (l 1 , -l 1 , l 2 , -l 2 )
) and a maximal element given by (l 2 + l 1 , l 2 + l 1 ) (resp. (l 1 , l 1 , l 2 , l 2 )). Then, (A) dispose the columns from west to east, increasingly according to ≺ c , (B) dispose the rows from north to south, increasingly according to ≺ r . For instance, by setting l 1 = 0 and l 2 ≥ 1, one obtains that C l1l2 is the (2l 2 + 1) × (2l 2 + 1) square matrix C l2m 00l2m2 with column indices m = -(2l 2 + 1), ..., (2l 2 + 1) and row indices m 2 = -(2l 2 + 1), ..., (2l 2 + 1) (from the subsequent discussion, one also deduces that, in general, C lm 00l2m2 = δ l2 l δ m2 m ). By selecting l 1 = l 2 = 1, one sees that C 11 is the 9 × 9 matrix with elements

C lm 1m11m2 (for m 1 , m 2 = -1, 0, 1; l = 0, 1, 2, m = -l, ..., l) arranged as follows:                 C 0,0 1,-1;1,-1 C 1,-1 1,-1;1,-1 C 10 1,-1;1,-1 C 11 1,-1;1,-1 C 2,-2 1,-1;1,-1 C 2,-1 1,-1;1,-1 C 2,0 1,-1;1,-1 C 2,1 1,-1;1,-1 C 2,2 1,-1;1,-1 C 0,0 1,-1;1,0 ... ... ... ... ... ... ... ... C 0,0 1,-1;1,1 ... ... C 11 1,-1;1,1 ... ... ... ... ... C 0,0 1,0;1,-1 ... ... ... ... ... C 2,0 1,0;1,-1 ... ... C 0,0 1,0;1,0 ... ... ... C 2,-2 1,0;1,0 ... ... ... ... C 0,0 1,0;1,1 ... ... ... ... ... ... C 2,1 1,0;1,1 ... C 0,0 1,1;1,-1 ... ... ... ... ... ... ... ... C 0,0 1,1;1,0 ... ... ... ... ... ... ... ... C 0,0 1,1;1,1 C 1,-1 1,1;1,1 C 1,0 1,1;1,1 C 1,1 1,1;1,1 C 2,-2 1,1;1,1 C 2,-1 1,1;1,1 C 2,0 1,1;1,1 C 2,1 1,1;1,1 C 2,2 1,1;1,1                
Explicit expressions for the Clebsch-Gordan coefficients of SO(3) are known, but they are in general hardly manageable. We have for instance (see [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF], expression 8.2.1.5)

C l3-m3 l1m1l2m2 := (-1) l1+l3+m2 2l 3 + 1 (l 1 + l 2 -l 3 )!(l 1 -l 2 + l 3 )!(l 1 -l 2 + l 3 )! (l 1 + l 2 + l 3 + 1)! 1/2 × (l 3 + m 3 )!(l 3 -m 3 )! (l 1 + m 1 )!(l 1 -m 1 )!(l 2 + m 2 )!(l 2 -m 2 )! 1/2 × z (-1) z (l 2 + l 3 + m 1 -z)!(l 1 -m 1 + z)! z!(l 2 + l 3 -l 1 -z)!(l 3 + m 3 -z)!(l 1 -l 2 -m 3 + z)! ,
where the summation runs over all z's such that the factorials are non-negative. This expression becomes much neater for m 1 = m 2 = m 3 = 0, where we have

C l30 l10l20 =      0 , for l 1 + l 2 + l 3 odd (-1) l 1 +l 2 -l 3 2 √ 2l3+1[(l1+l2+l3)/2]! [(l1+l2-l3)/2]![(l1-l2+l3)/2]![(-l1+l2+l3)/2]! (l1+l2-l3)!(l1-l2+l3)!(-l1+l2+l3)! (l1+l2+l3+1)! 1/2 , for l 1 + l 2 + l 3 even .
The coefficients, moreover, enjoy a nice set of symmetry and orthogonality properties, playing a crucial role in our results to follow. From unitary equivalence we have the two relations:

m1,m2 C lm l1m1l2m2 C l ′ m ′ l1m1l2m2 = δ l′ l δ m′ m , (3.15)
l,m

C lm l1m1l2m2 C lm l1m ′ 1 l2m ′ 2 = δ m ′ 1 m1 δ m ′ 2 m2 ; (3.16)
1 Given a finite set A = {a j : j = 1, ..., N } and an order ≺ on A, one says that A is a saturated chain with respect to ≺ if there exists a permutation π of {1, ..., N } such that

a π(1) ≺ a π(2) ≺ • • • ≺ a π(N-1) ≺ a π(N) .
In this case, a π(1) and a π(N) are called, respectively, the least and the maximal elements of the chain (see [28, p. 99]) in particular, (3.15) is a consequence of the orthogonality of row vectors, whereas (3.16) comes from the orthogonality of columns. Other properties are better expressed in terms of the Wigner's coefficients, which are related to the Clebsch-Gordan coefficients by the identities (see [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF], Chapter 8)

l 1 l 2 l 3 m 1 m 2 -m 3 = (-1) l3+m3 1 √ 2l 3 + 1 C l3m3 l1-m1l2-m2
(3.17)

C l3m3 l1m1l2m2 = (-1) l1-l2+m3 2l 3 + 1 l 1 l 2 l 3 m 1 m 2 -m 3 . (3.18)
The Wigner's 3j (and, consequently, the Clebsch-Gordan) coefficients are real-valued, they are different from zero only if m 1 + m 2 + m 3 = 0 and l i ≤ l j + l k for all i, j, k = 1, 2, 3 (triangle conditions), and they satisfy the symmetry conditions

l 1 l 2 l 3 m 1 m 2 m 3 = (-1) l1+l2+l3 l 1 l 2 l 3 -m 1 -m 2 -m 3 , l 1 l 2 l 3 m 1 m 2 m 3 = (-1) sign(π) l π(1) l π(2) l π(3) m 2 m 3 m 1 ,
where π is a permutation of {1, 2, 3}, and sign (π) denotes the sign of π. It follows also that for m 1 = m 2 = m 3 = 0, the coefficients C l30 l10l20 are different from zero only when the sum l 1 + l 2 + l 3 is even. Later in the paper, we shall also need the so-called Wigner's 6j coefficients, which are defined by

a b e c d f := α,β,γ ε,δ,φ (-1) e+f +ε+φ a b e α β ε c d e γ δ -ε a d f α δ -φ c b f γ β φ , (3.19
) see [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF], chapter 9 for analytic expressions and a full set of properties; we simply recall here that the Wigner's 6j coefficients can themselves be given an important interpretation in terms of group representations, namely they relate different coupling schemes in the decomposition of tensor product into direct sum representations, see [START_REF] Biedenharn | The Racah-Wigner algebra in quantum theory[END_REF] for further details.

For future reference, we also recall some further standard properties of Kronecker (tensor) products and direct sums of matrices: we have

⊕ n i=1 (A i B i ) = (⊕ n i=1 A i ) (⊕ n i=1 B i ) , (3.20) (⊕ n i=1 A i ) ⊗ B = ⊕ n i=1 (A i ⊗ B) (3.21)
and, provided all matrix products are well-defined,

(AB ⊗ C) = (A ⊗ I n ) (B ⊗ C) . (3.22)
Here, ⊕ n i=1 A i is defined as the block diagonal matrix diag {A 1 , ..., A n } if A i is a set of square matrices of order r i × r i , whereas it is defined as the stacked column vector of order (

n i=1 r i ) × 1 if the A i are r i × 1 column vectors.

Characterization of polyspectra 4.1 Four general statements

The following result is well-known. As it is crucial in our arguments to follow and we failed to locate any explicit reference, we shall provide a short proof for the sake of completeness. Note that, in the sequel, we use the symbol a l. to indicate the (2l + 1)-dimensional complex-valued random vector defined in (2.4).

Lemma 4 Let T be a strongly isotropic field on S 2 , and let the harmonic coefficients {a lm } be defined according to (2.2). Then, for every l ≥ 0 and every g ∈ SO(3), we have

D l (g)a l. d = a l. , l = 0, 1, 2, ... . (4.23)
The equality (4.23) must be understood in the sense of finite-dimensional distributions for sequences of random vectors, that is, (4.23) takes place if, and only if, for every k ≥ 1 and every 0 ≤

l 1 < l 2 < ••• < l k , D l1 (g)a l1. , ..., D l k (g)a l k . d = {a l1. , ..., a l k . } . (4.24)
Proof. We provide the proof of (4.24) only when k = 1 and l 1 = l ≥ 1. The general case is obtained analogously. By strong isotropy, we have that, for every l ≥ 1, every g ∈ SO (3) and every x 1 , ..., x n ∈ S 2 , the equality (2.1) takes place. Now, (2.1) can be rewritten as follows:

l m a lm Y lm (x 1 ), ..., l m a lm Y lm (x n ) d = l m a lm Y lm (gx 1 ), ..., l m a lm Y lm (gx n ) = l m a lm m ′ D l m ′ m (g)Y lm ′ (x 1 ), ..., l m a lm m ′ D l m ′ m (g)Y lm ′ (x n ) = l m ′ a lm ′ Y lm ′ (x 1 ), ..., l m ′ a lm ′ Y lm ′ (x n ) , (4.25)
where we write

a lm ′ m a lm D l m ′ m (g), (4.26)
and we have used

{Y lm (gx 1 ), ..., Y lm (gx n )} ≡ m ′ D l m ′ m (g)Y lm ′ (x 1 ), ..., m ′ D l m ′ m (g)Y lm ′ (x n ) . (4.27)
which follows from the group representation property and the identity (3.9). To conclude, just observe that (4.25) implies that

a lm ′ = S 2 T (gx) Y lm ′ (x)dx, m ′ = -l, ..., l,
yielding that, due to strong isotropy and with obvious notation, a l. d = a l. . The conclusion follows from the fact that, thanks to (4.26), a l. = D l (g) a l. .

The next theorem connects the invariance properties of the vectors {a l. } to the representations of SO(3). We need first to establish some notation. For every 0 ≤ l 1 , l 2 , ..., l n , we shall write

∆ l1...ln SO(3) D l1 (g) ⊗ D l2 (g) ⊗ ... ⊗ D ln (g) dg , (4.28) ∆ l1...ln (g) D l1 (g) ⊗ D l2 (g) ⊗ ... ⊗ D ln (g), g ∈ SO (3) , (4.29)
and use the symbol S l1...ln (whenever is well-defined), as given in formula (2.5). We stress that ∆ l1...ln and ∆ l1...ln (g) are square matrices with (2l 1 + 1) × ... × (2l n + 1) rows and S l1...ln is a column vector with (2l 1 + 1) × ... × (2l n + 1) elements. The following result applies to an arbitrary n ≥ 2: see [START_REF] Hu | The Angular trispectrum of the CMB[END_REF] for some related results in the case n = 3, 4.

Proposition 5 Let T be a strongly isotropic field with moments of order n ≥ 2. Then, for every 0 ≤ l 1 , l 2 , ..., l n and every fixed g * ∈ SO (3) ∆ l1...ln S l1...ln = S l1...ln (4.30) ∆ l1...ln (g * ) S l1...ln = S l1...ln .

(4.31)

On the other hand, fix n ≥ 2 and assume that T (x) is a not necessarily isotropic random field on the sphere s.t.

sup x (E |T (x)| n ) < ∞.
Then T (.) is P -almost surely Lebesgue square integrable and the nth order spectral moments of T exist and are finite. If moreover (4.30) holds for every 0 ≤ l 1 ≤ • • • ≤ l n , then one has that, for every g ∈ SO (3),

E D l1 (g)a l1. ⊗ • • • ⊗ D ln (g)a ln. = E [a l1. ⊗ • • • ⊗ a ln. ] , (4.32)
and T is n-weakly isotropic.

Proof. By strong isotropy and Lemma 4, one has

E D l1 (g)a l1. ⊗ ... ⊗ D ln (g)a ln. = E {a l1. ⊗ ... ⊗ a ln. } for all g ∈ SO(3) , l 1 , ..., l n ∈ N n .
Now assume that g is sampled randomly (and independently of the {a l. }) according to some probability measure, say P 0 , on SO(3). From the property (3.22) of tensor products and trivial manipulations, we obtain (with obvious notation and by independence)

E D l1 (•)a l1. ⊗ ... ⊗ D ln (•)a ln. = E D l1 (•) ⊗ ... ⊗ D ln (•) [a l1. ⊗ ... ⊗ a ln. ] = E 0 D l1 (•) ⊗ ... ⊗ D ln (•) E {a l1. ⊗ ... ⊗ a ln. } .
Now, if one chooses P 0 to be equal to the Haar (uniform) measure on SO (3), one has that

E 0 D l1 (•) ⊗ ... ⊗ D ln (•) = ∆ l1...ln ,
thus giving (4.30). On the other hand, if one chooses P 0 to be equal to the Dirac mass at some g * ∈ SO (3), one has that

E 0 D l1 (•) ⊗ ... ⊗ D ln (•) = ∆ l1...ln (g * ) ,
which shows that (4.31) is satisfied. Now let T satisfy the assumptions of the second part of the statement for some n ≥ 2. We recall first that the representation (2.2) continues to hold, in a pathwise sense. To see that the nth order joint moments of the harmonic coefficients a lm are finite it is enough to use Jensen's inequality, along with a standard version of the Fubini theorem, to obtain that

E |a lm | n = E S 2 T (x)Y lm (x)dx n ≤ E S 2 |T (x)| n |Y lm (x)| n dx ≤ sup x∈S 2 |Y lm (x)| n sup x∈S 2 E|T (x)| n ≤ 2l + 1 4π n/2 sup x∈S 2 E|T (x)| n < ∞ .
It is then straightforward that, if S l1...ln satisfies (4.30), one also has that for any fixed g ∈ SO(3) Note that relation (4.30) can be rephrased by saying that, for a strongly isotropic field, the joint moment vector E {a l1. ⊗ a l2. ⊗ ... ⊗ a ln. } must be an eigenvector of the matrix (4.28) for every n ≥ 2 and every 0 ≤ l 1 ≤ • • • ≤ l n . A similar characterization holds for cumulants polyspectra. Recall the notation S c l1...ln introduced in (2.6).

E D l1 (g) ⊗ ... ⊗ D ln (g) [a l1. ⊗ ... ⊗ a ln. ] = D l1 (g) ⊗ ... ⊗ D ln (g) E [a l1. ⊗ ... ⊗ a ln. ] = D l1 (
Proposition 6 Let T be a strongly isotropic field with moments of order n ≥ 2. Then, for every 0 ≤ l 1 , l 2 , ..., l n and every fixed g * ∈ SO (3), On the other hand, fix n ≥ 2 and assume that T (x) is a not necessarily isotropic random field on the sphere s.t.

sup x (E |T (x)| n ) < ∞.
Then T (.) is P -almost surely Lebesgue square integrable and the nth order spectral moments of T exist and are finite. If moreover (4.33) holds for every 0 ≤ l 1 ≤ • • • ≤ l n , then one has that, for every g ∈ SO (3), relation (4.32) holds, and T is n-weakly isotropic.

Proof. For every x 1 , ..., x n ∈ S 2 , write Cum {T (x 1 ) , ..., T (x n )} the joint cumulant of the random variables T (x 1 ) , ..., T (x n ). By using isotropy, one has that, for every g ∈ SO (3),

Cum {T (x 1 ) , ..., T (x n )} = Cum {T (gx 1 ) , ..., T (gx n )} . ( 4
)} ≡ m ′ D l1 m ′ m1 (g)Y l1m ′ (x 1 ), ..., m ′ D ln m ′ mn (g)Y lnm ′ (x n ) .
The second part of the statement is proved by arguments analogous to the ones used in the proof of Proposition 5.

We now present an alternative (and more involved) characterization of the cumulant polyspectra associated with an isotropic field. Given n ≥ 2 and a partition π = {b 1 , ..., b k } ∈ Π (n), we build a permutation v π = (v π (1) , ..., v π (n)) ∈ S n as follows: (i) write the partition

π = {b 1 , ..., b k } = i 1 1 , ..., i 1 |b1| , ..., i k 1 , ..., i k |b k | (4.37)
(where |b j | ≥ 1 stands for the size of b j ) by means of the convention outlined in Section 2 (that is, order the blocks and the elements within each block according to the lexicographic order); (ii) define v π = S n by simply removing the brackets in (4.37), that is, set

v π = (v π (1) , ..., v π (n)) = i 1 1 , ..., i 1 |b1| , i 2 1 , ..., i 2 |b2| , ..., i k 1 , ..., i k |b k | .
For instance, if a partition π of {1, ..., 6} is composed of the blocks {1, 3} , {6, 4} and {2, 5}, one first writes π in the form π = {{1, 3} , {2, 5} , {4, 6}}, and then defines v π = (v π (1) , ..., v π (6)) = (1, 3, 2, 5, 4, 6).

Given n ≥ 2, 0 ≤ l 1 ≤ • • • ≤ l n ,
D l v π (1) (g) ⊗ D l v π (2) (g) ⊗ ... ⊗ D l v π (n) (g) dg , (4.38)
obtained from the matrix ∆ l1...ln in (4.28), by permuting the indexes l i according to v π . Plainly, if v π is equal to the identity permutation, then ∆ π l1...ln = ∆ l1...ln . We also set, for every fixed g ∈ SO (3),

∆ π l1...ln (g) D l v π (1) (g) ⊗ D l v π (2) (g) ⊗ ... ⊗ D l v π (n) (g).
Proposition 7 Let T be a strongly isotropic field with finite moments of order n ≥ 2. For 0 ≤ l 1 , l 2 , ..., l n , define S c l1...ln according to (2.6). Then, for every 0 ≤ l 1 , l 2 , ..., l n , and every g ∈ SO (3)

S c l1,...ln = π={b1,...,b k }∈Π(n) (-1) k-1 (k -1)!∆ π l1...ln E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ] (4.39) = π={b1,...,b k }∈Π(n) (-1) k-1 (k -1)!∆ π l1...ln (g) E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ] (4.40)
On the other hand, fix n ≥ 2 and assume that T (x) is a (not necessarily isotropic) random field on the sphere s.t.

sup x (E |T (x)| n ) < ∞.
Then, the nth order spectral moments and cumulants of T exist and are finite. If moreover (4.40) holds for every 0 ≤ l 1 , l 2 , ..., l n and every g ∈ SO (3), then one has that T is n-weakly isotropic.

Proof. Fix π = {b 1 , ..., b k } ∈ Π (n). By strong isotropy and Lemma 4, one has that, for a fixed g * ∈ SO (3), the quantity

E ⊗ i∈b1 D li (g) a li. ⊗ • • • ⊗ E ⊗ i∈b k D li (g) a li. = ∆ π l1...ln (g * ) E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ]
does not depend on g * , so that

E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ] = ∆ π l1...ln (g * ) E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ] = SO(3) ∆ π l1...ln (g) E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ] dg = ∆ π l1...ln E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ] .
To prove the second part of the statement, suppose that T (x) verifies sup x (E |T (x)| n ) < ∞, and that its associated harmonic coefficients verify (4.40). Then, for every fixed rotation g * ∈ SO(3),

π={b1,...,b k }∈Π(n) (-1) k-1 (k -1)!E ⊗ i∈b1 D li (g * ) a li. ⊗ • • • ⊗ E ⊗ i∈b k D li (g * ) a li. = π={b1,...,b k }∈Π(n) (-1) k-1 (k -1)! × ×[D l v π (1) (g * ) ⊗ • • • ⊗ D l v π (n) (g * )]E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ] = π={b1,...,b k }∈Π(n) (-1) k-1 (k -1)! × ∆ π l1...ln (g * ) E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ] = π={b1,...,b k }∈Π(n) (-1) k-1 (k -1)!E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ]
By the definition of cumulants, this last equality gives that

E D l1 (g * )a l1. ⊗ • • • ⊗ D ln (g * )a ln. = E [a l1. ⊗ • • • ⊗ a ln. ] .
Since g * is arbitrary, the n-weak isotropy follows from (2.2).

Remark. By combining (4.33) and (4.39) we obtain for instance that the nth cumulant polyspectrum of an isotropic field verifies the identity

S c l1...ln = ∆ l1...ln S c l1...ln = π={b1,...,b k }∈Π(n) (-1) k-1 (k -1)!∆ π l1...ln E [⊗ i∈b1 a li. ] ⊗ • • • ⊗ E [⊗ i∈b k a li. ] .
5 Angular polyspectra and the structure of ∆ l 1 ...l n

Spectra of strongly isotropic fields

Our aim in this section is to investigate more deeply the structure of the matrix ∆ l1...ln appearing in (4.28), in order to derive an explicit characterization for the angular polyspectra. As a preliminary example, we deal with the case n = 2.

Proposition 8 For integers l 1 , l 2 ≥ 0, one has that

∆ l1l2 = SO(3) D l1 (g) ⊗ D l2 (g) dg = δ l2 l1 C 00 l1.l2. (C 00 l1.l2. ) ′ , (5.41
)

that is: if l 1 = l 2 , then ∆ l1l2 is a (2l 1 + 1) (2l 2 + 1) × (2l 1 + 1) (2l 2 + 1) zero matrix; if l 1 = l 2 , then ∆ l1l2 = ∆ l1l1
is given by C 00 l1.l1. (C 00 l1.l1. ) ′ .

Proof. Using the equivalence of the two representations D l1 (g) ⊗ D l2 (g) and ⊕ l2+l1 λ=|l2-l1| D λ (g), as well as the definition of the Clebsch-Gordan matrices, we obtain that

SO(3) D l1 (g) ⊗ D l2 (g) dg = C l1l2 SO(3) ⊕ l2+l1 λ=|l2-l1| D λ (g) dg C * l1l2 .
(5.42) Now, if l 1 = l 2 , then the RHS of (5.42) is equal to the zero matrix since, as a consequence of the Peter-Weyl theorem and for λ = 0, the entries of D λ (•) are orthogonal to the constants. If l 1 = l 2 , then the integrated matrix on the RHS of (5.42) becomes SO(3) ⊕ 2l1 λ=0 D λ (g) dg, that is, a (2l 1 + 1) 2 × (2l 1 + 1)

2
matrix which is zero everywhere, except for the entry in the top-left corner, which is equal to one (since

SO(3) dg = 1
). The proof is concluded by checking that

C l1l1 SO(3) ⊕ 2l1 λ=0 D λ (g) dg C * l1l1 = C 00 l1.l1. (C 00 l1.l1. ) ′ .
Remark. Recall that C 00 l1.l2. is a column vector of dimension (2l 1 + 1) (2l 2 + 1), corresponding to the first column of the matrix C l1l2 . Also, according e.g. to [31, .

Proposition 8 provides a characterization of the spectrum of a strongly isotropic field.

Corollary 9 Let T be a strongly isotropic field with second moments, and let the vectors of the harmonic coefficients {a l. } be defined according to (2.2). Then, for any integers l 1 , l 2 ≥ 0, one has that

E {a l1. ⊗ a l2. } = (-1) m1 2l 1 + 1 δ l2 l1 δ -m2 m1 C l1 (5.43)
for some C l1 ≥ 0 depending uniquely on l 1 .

Proof. According to (4.30), one has that

E {a l1. ⊗ a l2. } = δ l2 l1 C 00 l1.l2. (C 00 l1.l2. ) ′ E {a l1. ⊗ a l2. } ,
implying that E {a l1. ⊗ a l2. } is (a) equal to the zero vector for l 1 = l 2 , and (b) of the form C l1 × C 00 l1.l2. , for some constant C l1 , when l 1 = l 2 . To see that C l1 cannot be negative, just observe that a l10 is real-valued for every l 1 ≥ 0, so that (5.43) yields that

C l1 = (2l 1 + 1) × E a 2 l10 .
In the subsequent two subsections, we shall obtain, for every n ≥ 3, a characterization of ∆ l1...,ln and E{a l1. ⊗ • • • ⊗ a ln. }, respectively analogous to (5.41) and (5.43).

The structure of ∆ l 1 ...ln

We first need to establish some further notation. Definition B. Fix n ≥ 3. For integers l 1 , ..., l n ≥ 0, we define C l1...ln to be the unitary matrix, of dimension Remarks.

(1) Fix l 1 , ..., l n ≥ 0, as well as g ∈ SO (3). Then, the matrix

⊕ l2+l1 λ1=|l2-l1| ⊕ l3+λ1 λ2=|l3-λ1| ... ⊕ ln+λn-2 λn-1=|ln-λn-2| D λn-1 (g) (5.46)
is a block-diagonal matrix, obtained as follows. (a) Consider vectors of integers (λ 1 , ..., λ n-1 ) satisfying the relations

|l 2 -l 1 | ≤ λ 1 ≤ l 1 + l 2 , and |l k+1 -λ k-1 | ≤ λ k ≤ l k+1 + λ k-1 , for k = 2, ..., n -1. (b)
Introduce a (total) order ≺ 0 on the collection of these vectors by saying that (λ 1 , ..., λ n-1 ) ≺ 0 λ ′ 1 , ..., λ ′ n-1 , (5.47) whenever either λ 1 < λ ′ 1 , or there exists k = 2, ..., n -2 such that λ j = λ ′ j for every j = 1, ..., k, and λ k+1 < λ ′ k+1 . (c) Associate to each vector (λ 1 , ..., λ n-1 ) the matrix D λn-1 (g). (d) Construct a blockdiagonal matrix by disposing the matrices D λn-1 (g) from the top-left corner to the bottom-right corner, in increasing order with respect to ≺ 0 . As an example, consider the case where n = 3 and l 1 = l 2 = l 3 = 1. Here, the vectors (λ 1 , λ 2 ) involved in the direct sum (5.45) are (in increasing order with respect to ≺ 0 ) (0, 1) , (1, 0) , (1, 1) , (1, 2) , (2, 1) , (2, 2) and (2, 3) , and the matrix (5.46) is therefore given by 

         D 1 (
          (5.48)
where the dots indicate zero entries, and we have used the fact that D 0 (g) ≡ 1.

(2) The fact that the representation (5.45) has dimension n j=1 (2l j + 1) is a direct consequence of formula (3.14).

(3) The fact that the two representations (5.44) and ( 5.45) are equivalent can be proved by iteration. Indeed, by standard representation theory, on has that (5.44) is equivalent to

⊕ l2+l1 λ1=|l2-l1| D λ1 (.) ⊗ D l3 (•) ⊗ • • • ⊗ D ln (•) ,
which is in turn equivalent to

⊕ l2+l1 λ1=|l2-l1| ⊕ l3+λ1 λ2=|l3-λ1| D λ2 (.) ⊗ D l4 (•) ⊗ • • • ⊗ D ln (•) .
By iterating the same procedure until all tensor products have disappeared (that is, by successively replacing the tensor product

D λ k (.) ⊗ D l k+2 (•) with ⊕ l k+2 +λ k λ k+1 =|l k+2 -λ k | D λ2 (.) for k = 2, ..., n -1)
, one obtains the desired conclusion.

For every n ≥ 3 and every l 1 , ..., l n ≥ 0, the elements of the matrix C l1...ln , introduced in Definition B, can be written in the form C λ1...λn-1,µn-1 l1m1...lnmn . The indices (m 1 , ..., m n ) are such that -l i ≤ m i ≤ l i (i = 1, ..., n) and label rows; on the other hand, the indices (λ 1 ...λ n-1 , µ n-1 ) label columns, and verify the relations

|l 2 -l 1 | ≤ λ 1 ≤ l 1 + l 2 , |l k+1 -λ k-1 | ≤ λ k ≤ l k+1 + λ k-1 (k = 2, ..., n -1) and -λ n-1 ≤ µ n-1 ≤ λ n-1 .
It is well known (see e.g. [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF]) that the quantity C

= µn-2    µ1...µn-3 C λ1µ1 l1m1l2m2 C λ2µ2 λ1µ1l3m3 ...C λn-2µn-2 λn-3µn-3ln-1mn-1    C λn-1µn-1 λn-2µn-2lnmn = µ1...µn-2 C λ1µ1 l1m1l2m2 C λ2µ2 λ1µ1l3m3 ...C λn-2µn-2 λn-3µn-3ln-1mn-1 C λn-1µn-1 λn-2µn-2lnmn .
Remark. Given an enumeration of the coefficients C λ1...λn-1,µn-1 l1m1...lnmn , the matrix C l1...ln can be built (analogously to the case of the Clebsch-Gordan matrices of Section 3.2) by disposing rows (from top to bottom) and columns (from left to right) increasingly according to two separate total orders. The order ≺ r on the symbols (m 1 , ..., m n ) is obtained by setting that (m 1 , ..., m n ) ≺ r (m ′ 1 , ..., m ′ n ) whenever either m 1 < m ′ 1 , or there exists k = 2, ..., n -1 such that m j = m ′ j for every j = 1, ..., k, and m k+1 < m ′ k+1 . The order ≺ c on the symbols (λ 1 ...λ n-1 , µ n-1 ) is obtained by setting that (λ

1 ...λ n-1 , µ n-1 ) ≺ c λ ′ 1 ...λ ′ n-1 , µ ′ n-1
whenever either (λ 1 , ..., λ n-1 ) ≺ 0 λ ′ 1 , ..., λ ′ n-1 , as defined in (5.47), or λ i = λ ′ i for every i = 1, ..., n -1 and µ n-1 < µ ′ n-1 .

One has also the following (useful) alternative representation of generalized Clebsch-Gordan matrices.

Proposition 10 For every n ≥ 3 and every l 1 , ..., l n ≥ 0, one can represent the matrix C l1...ln , as follows

C l1...ln = C l1l2l3...ln-1 ⊗ I 2ln+1 (⊕ l2+l1 λ1=|l2-l1| ... ⊕ ln+λn-3 λn-2=|ln-λn-3| C λn-2ln ,
where I m indicates a m × m identity matrix. Also, one has that

C l1...ln = (C l1l2 ⊗ I 2l3+1 ⊗ ... ⊗ I 2ln+1 ) × (⊕ l2+l1 λ=|l2-l1| C λl3 ) ⊗ ... ⊗ I 2ln+1 ×... × (⊕ l2+l1 λ1=|l2-l1| ... ⊕ ln+λn-3 λn-2=|ln-λn-3| C λn-2ln ,
where × stands for the usual product between matrices.

Definition C. For every n ≥ 3 and every l 1 , ..., l n ≥ 0, we define E l1...ln to be the Π n j=1 (2l j + 1) × Π n j=1 (2l j + 1) square matrix

E l1...ln := ⊕ l2+l1 λ1=|l2-l1| ... ⊕ ln+λn-2 λn-1=|ln-λn-2| δ 0 λn-1 I 2λn-1+1 .
(5.49)

In other words, E l1...ln is the diagonal matrix built from the matrix (5.46), by replacing every block of the type D λn-1 (g), with λ n-1 > 0, with a (2λ n-1 + 1) × (2λ n-1 + 1) zero matrix, and by letting the 1 × 1 blocks D 0 (g) = 1 unchanged. For instance, by setting n = 3 and l 1 = l 2 = l 3 = 1 (and by using (5.48)) one obtains a 27 × 27 matrix E 111 whose entries are all zero, except for the fourth element (starting from the top-left corner) of the main diagonal.

The following result states that the matrix ∆ l1.

D l1 (g) ⊗ D l2 (g) ⊗ • • • ⊗ D ln (g) dg (5.51) = SO(3) C l1...ln ⊕ l2+l1 λ1=|l2-l1| ⊕ l3+λ1 λ2=|l3-λ1| ... ⊕ ln+λn-2 λn-1=|ln-λn-2| D λn-1 (g)C * l1...ln dg.
By linearity and by the definition of the integral of a matrix-valued function, one has that the last line of (5.51) equals

C l1...ln ⊕ l2+l1 λ1=|l2-l1| ⊕ l3+λ1 λ2=|l3-λ1| ... ⊕ ln+λn-2 λn-1=|ln-λn-2| SO(3) D λn-1 (g)dg C * l1...ln .
Now observe that, if λ n-1 > 0, then SO(3) D λn-1 (g)dg equals a (2λ n-1 + 1) × (2λ n-1 + 1) zero matrix, whereas SO(3) D 0 (g)dg = SO(3) 1dg = 1. The conclusion is obtained by resorting to the definition of E l1...ln given in (5.49).

Existence and characterization of reduced polyspectra of arbitrary orders

Combining the previous Proposition with [START_REF] Baldi | On the characterization of isotropic random fields on homogeneous spaces of compact groups[END_REF], we obtain the main result of this paper.

Theorem 12 If a random field is strongly isotropic with finite moments of order n ≥ 3, then for every l 1 , ..., l n there exists two arrays P l1....ln (λ 1 , ..., λ n-3 ) and P C l1....ln (λ 1 , ..., λ n-3 ), with C λ1µ1 l1m1l2m2 C λ2µ2 λ1µ1l3m3 ...C ln,-mn λn-3µn-3ln-1mn-1 .

|l 2 -l 1 | ≤ λ 1 ≤ l 2 + l 1 , |l 3 -λ 1 | ≤ λ 2 ≤ l 3 + λ 1 , ..., |l n-2 -λ n-4 | ≤ λ n-3 ≤ l n-2 + λ n-4 , such that
(5.54)

Remark. For a fixed n ≥ 2, the real-valued arrays {P l1...ln (•) : l 1 , ..., l n ≥ 0} and P C l1...ln (•) : l 1 , ..., l n ≥ 0 are, respectively, the reduced polyspectrum of order n -1 and the reduced cumulant polyspectrum of order n -1 associated with the underlying strongly isotropic random field.

Proof of Theorem 12. We shall prove only (5.52), since the proof of (5.53) is entirely analogous. By Proposition 5 and Proposition 11, if the random field is isotropic, then It follows that S l1...ln is a solution if and only if the column vector C * l1...ln S l1...ln has zeroes corresponding to the zeroes of E l1...ln , whereas the elements corresponding to unity can be arbitrary. In view of the orthonormality properties of C * l1...ln , this condition is met if, and only if, S l1...ln is a linear combination of the columns in the matrix C * l1...ln corresponding to non-zero elements of the diagonal E l1...ln . These linear combinations can be written explicitly as

l2-l1 λ1=l2-l1 l3+λ1 λ2=l3-λ1 ... ln+λn-2 λn-1=ln-λn-2 C λ1...λn-2lm l1m1....lnmn P l1....ln (λ 1 , ..., λ n-3 , λ n-2 )δ 0 l = l2-l1 λ1=l2-l1 l3+λ1 λ2=l3-λ1 ... ln+λn-2 λn-1=ln-λn-2    µ1...µn-2 C λ1µ1 l1m1l2m2 C λ2µ2 λ1µ1l3m3 ...C lm λn-2µn-2.lnmn δ 0 l    × × P l1....ln (λ 1 , ..., λ n-3 , λ n-2 ) = l2-l1 λ1=l2-l1 l3+λ1 λ2=l3-λ1 ... ln+λn-2 λn-1=ln-λn-2    µ1...µn-2 C λ1µ1 l1m1l2m2 C λ2µ2 λ1µ1l3m3 ...C 00 λn-2µn-2.lnmn    × × P l1....ln (λ 1 , ..., λ n-3 , λ n-2 ).
Recalling again that

C 0m l1m1l2m2 = (-1) m1 2l 1 + 1 δ l2 l1 δ -m2 m1 δ 0 m ,
(see [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF], 8.5.1.1), we obtain that

= l2-l1 λ1=l2-l1 l3+λ1 λ2=l3-λ1 ... ln+λn-2 λn-1=ln-λn-2    µ1...µn-2 C λ1µ1 l1m1l2m2 C λ2µ2 λ1µ1l3m3 ... (-1) mn 2l n + 1 δ ln λn-2 δ -mn µn-2    × × P l1....ln (λ 1 , ..., λ n-3 , λ n-2 ) = l2-l1 λ1=l2-l1 l3+λ1 λ2=l3-λ1 ... ln+λn-2 λn-1=ln-λn-2    µ1...µn-2 C λ1µ1 l1m1l2m2 C λ2µ2 λ1µ1l3m3 ...C ln-mn λn-3µn-3.ln-1mn-1 (-1) mn    × × P l1....ln (λ 1 , ..., λ n-3 ) = l2-l1 λ1=l2-l1 l3+λ1 λ2=l3-λ1
... where we have set

P l1....ln (λ 1 , ..., λ n-3 ) := 1 2l n + 1 P l1....ln (λ 1 , ..., λ n-3 , l n ) .
All there is left to show is that the coefficients of this linear combination are necessarily real. To see this, it is sufficient to specialize the previous discussion to the case where m 1 = m 2 = ... = m n = 0,, and to observe that, in this case Examples. For n = 3, Theorem 12 implies that, under isotropy Ea l1m1 a l2m2 a l3m3 = (-1) m3 C l3-m3 l1m1l2m2 P l1l2l3 .

From this last relation, we can recover the so-called reduced bispectrum, noted b l1l2l3 , defined for instance in [START_REF] Hu | The Angular trispectrum of the CMB[END_REF], [START_REF] Marinucci | High-resolution asymptotics for the angular bispectrum of spherical random fields[END_REF] and [START_REF] Marinucci | A Central Limit Theorem and Higher Order Results for the Angular Bispectrum[END_REF], which satisfies indeed the relationship

P l1l2l3 = b l1l2l3 C l30 l10l20 (2l 1 + 1)(2l 2 + 1) (2l 3 + 1)4π .
For n = 4 (i.e. the trispectrum, [START_REF] Hu | The Angular trispectrum of the CMB[END_REF]) we obtain the expression

Ea l1m1 a l2m2 a l3m3 a l4m4 = (-1) m4 l2+l1 λ=|l2-l1| C λl4-m4 l1m1l2m2l3m3 P l1l2l3l4 (λ) = l2+l1 λ=|l2-l1| λ µ=-λ C λµ l1m1l2m2 C l4-m4 λµl3m3 P l1l2l3l4 (λ) .
The next result gives a further probabilistic characterization of the reduced bispectrum.

Proposition 13 Fix n ≥ 2. A real-valued array {A l1...ln (•) : l 1 , ..., l n ≥ 0} is the reduced polyspectrum of order n -1 (resp. the reduced cumulant polyspectrum of order n -1) of some strongly isotropic random field if, and only if, there exists a sequence {X l : l ≥ 0} of zero-mean real-valued random variables such that l≥0 (2l + 1) E X 2 l < +∞ and, for every l 1 , ..., l n ≥ 0

E (X l1 • • • X ln ) = l2+l1 λ1=l2-l1
... Proof. We shall only prove (5.55). For the necessity it is enough to take X l = a l0 , where a l0 is the harmonic coefficient of index (l, 0) associated with a strongly isotropic field with moments of all orders. For the sufficiency, we consider first the (anisotropic) random field

Z (x) = l≥0 X l Y l0 (x) .
Then, by taking T (x) = Z (gx), where g is sampled randomly with the uniform Haar measure on SO (3), one obtains a random field with the desired characteristics.

There are two very important issues that are left open by Theorem 12. As a first issue, it seems natural to look for characterizations of the reduced polyspectra P l1...ln , at least under natural models of physical interest. As a second point, we note that the explicit expressions provided in Theorem 12 depend on the ordering l 1 , ..., l n we chose for the decomposition of ∆ l1...ln . In the next two sections, we try to address these (and other) points.

Some Explicit Examples

In this section we provide explicit computations for the reduced polyspectra P l1...ln (n ≥ 2), or P C l1...ln , for some models of physical interest. Of course, the Gaussian isotropic fields can be easily dealt with. Indeed, in this case one has that P C l1...ln = 0 for all n ≥ 3. In what follows, we shall therefore be concerned with polyspectra of Gaussian subordinated isotropic fields, that is, random fields that can be written as a deterministic and non-linear function of some collection of Gaussian isotropic fields. In general, this class of random fields allow for a clear-cut mathematical treatment, whilst covering a great array of empirically relevant circumstances.

A simple physical model

The general Gaussian-subordinated model has the form

T = q j=1 f j H j T G / E (T 2 G ) = f 1 T G + f 2 (T 2 G /E T 2 G -1) + ..., (6.57)
where f j is a real constant, H j (.) denotes the jth Hermite polynomial (see e.g. [START_REF] Surgailis | CLTs for polynomials of linear sequences: Diagram formula with illustrations[END_REF]), and T G is a Gaussian, zero-mean isotropic random field. Note that we have implicitly defined the sequence of Hermite polynomials in such a way that H 1 (x) = x, H 2 (x) = x 2 -1, H 3 (x) = x 3 -3x, and so on. In this section, when no further specification is needed, the spectral decomposition of the underlying Gaussian field T G is written

T G = lm a lm Y lm .
We shall sometimes use the following notation

T = lm a lm Y lm = q j=1
f j a lm (j)Y lm , (6.58)

a lm (j) = S 2 H j T G (x) / E (T 2 G ) Y lm (x) dx, (6.59
)

a lm = q j=1
a lm (j) . (6.60)

For instance, models of Cosmic Microwave Background radiation are currently dominated by assumptions such as the Sachs-Wolfe model with the so-called Bardeen's potential (see e.g. [START_REF] Bartolo | Non-Gaussianity from Inflation: Theory and Observations[END_REF] or [START_REF] Dodelson | Modern Cosmology[END_REF]). The latter can be written down explicitly as

T = T G + f N L (T 2 G -ET 2 G ) , (6.61)
where f N L is a nonlinearity parameters which depends upon physical constants in the associated "slowroll" inflationary model (see e.g. [START_REF] Bartolo | Non-Gaussianity from Inflation: Theory and Observations[END_REF]). Note that (6.61) has can be written in the form (6.57), by setting

f 1 = 1, f 2 = f N L × E T 2
G and f j = 0, for j ≥ 3. The value of the constant f N L × E T 2 G is expected to be very small, namely of the order 10 -4 [START_REF] Bartolo | Non-Gaussianity from Inflation: Theory and Observations[END_REF]. To simplify the discussion, we now assume that ET 2 G = 1. In this case, by using (6.58)-( 6.60), one has that

a lm = a lm + f N L a lm (2) , a lm (2) = S 2 T 2 Y lm dx = S 2 ℓ1ℓ2 m1m2 a ℓ1m1 a ℓ2m2 Y ℓ1m1 Y ℓ2m2 Y lm dx = ℓ1ℓ2 m1m2 a ℓ1m1 a ℓ2m2 (2ℓ 1 + 1)(2ℓ 2 + 1) (2l + 1)4π C l0 ℓ10ℓ20 C lm ℓ1m1ℓ2m2 .
It follows that

C l := E| a lm | 2 = C l + 2f 2 N L l1l2 C l1 C l2 (2l 1 + 1)(2l 2 + 1) 4π(2l + 1) C l0 l10l20 2 , so that V ar(T ) = l 2l + 1 4π C l = l 2l + 1 4π C l + 2f 2 N L l1l2 C l1 C l2 (2l 1 + 1)(2l 2 + 1) (4π) 2 l C l0 l10l20 2 = l 2l + 1 4π C l + 2f 2 N L l1 C l1 (2l 1 + 1) 4π 2 = V ar(T G ) + f 2 N L V ar(H 2 (T G )) ,
as expected, due to the orthogonality properties of Hermite polynomials. For the bispectrum, we obtain therefore

E a l1m1 a l2m2 a l3m3 = E {(a l1m1 + f 2 a l1m1 (2))(a l2m2 + f 2 a l2m2 (2))(a l3m3 + f 2 a l3m3 (2))} = f 2 Ea l1m1 (2)a l2m2 a l3m3 + f 2 Ea l1m1 a l2m2 (2)a l3m3 +f 2 Ea l1m1 a l2m2 a l3m3 (2) + f 3 2 Ea l1m1 (2)a l2m2 (2)a l3m3 (2) = (-1) m3 C l3-m3 l1m1l2m2 P l1l2l3 ,
where

P l1l2l3 = 6f 2 (2l 1 + 1)(2l 2 + 1) (2l 3 + 1)4π C l30 l10l20 {C l1 C l2 + C l1 C l3 + C l2 C l3 } (6.62) +f 3 2 ℓ1ℓ2ℓ3 C l10 ℓ10ℓ20 C l20 ℓ10ℓ30 C l30 ℓ20ℓ30 (2ℓ 1 + 1)(2ℓ 2 + 1)(2ℓ 3 + 1) (4π) 3 × (6.63) × 8(-1) l3 √ 2l 3 + 1 ℓ 1 ℓ 2 ℓ 3 l 3 l 2 l 1 {C ℓ1 C ℓ2 C ℓ3 } .
The lack of symmetry with respect to the l 3 term is only apparent and can be easily dispensed with by permuting the multipoles in C l3m3 l1m1l2m2 or using expression (3.18). Formula (6.62) is consistent with the cosmological literature, where (6.63) is considered a higher order term and hence neglected (see again ([15])).

The Connection with Higher Order Moments

We now provide a simple result, connecting the reduced polyspectrum with the higher order moments of the associated spherical random field.

Proposition 14

The following identity holds for every isotropic field with finite moments of order p and with a reduced polyspectrum P l1...lp (•) : l 1 , ..., l p ≥ 0 : for every x ∈ S 2 , ET (x) p ≡ l1...lp Proof. We use the trivial fact that

(2l 1 + 1) • • • (2l p + 1) (4π)
T (x) d = T (0) = l a l0 Y l0 (0) = l a l0 2l + 1 4π ,
where 0 is the North Pole and we used the fact that, for m = 0, Y lm (0) = 0 and Y l0 (0) = 2l+1 4π (see e.g. [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF]Chapter 5]). Hence, Example. Take T = H q (T G ), where H q is the qth Hermite polynomial. Then ET p = c pq ET 2 qp/2 , where c pq ∈ N denotes the number of Gaussian diagrams without flat edges with p rows and q columns (see [START_REF] Surgailis | CLTs for polynomials of linear sequences: Diagram formula with illustrations[END_REF]). Therefore, one has the identity

ET p = l1...lp (2l 1 + 1) • • • (2l p + 1) (4π) p E a l10 ...a lp0 = l1...lp (2l 1 + 1) • • • (2l p + 1) (4π)

The χ 2 ν polyspectrum

Previously in (6.63), we have implicitly derived the "χ 2 1 bispectrum", that is, the bispectrum associated with a field of the type T = H 2 (T G ), where T G is Gaussian, centered, isotropic and with unit variance. More precisely, with the notation (6.58)-( 6.60), one deduces from (6.63) that

Ea l1m1 (2)a l2m2 (2)a l3m3 (2) (6.64) = ℓ1ℓ2ℓ3× ℓ4ℓ5ℓ6 µ1...µ6 C l10 ℓ10ℓ20 C l1m1 ℓ1µ1ℓ2µ2 C l20 ℓ30ℓ40 C l2m2 ℓ3µ3ℓ4µ4 C l30 ℓ50ℓ60 C l3m3 ℓ5µ5ℓ6µ6 × × (2ℓ 1 + 1)(2ℓ 2 + 1) (2l 1 + 1)4π (2ℓ 3 + 1)(2ℓ 4 + 1) (2l 2 + 1)4π (2ℓ 5 + 1)(2ℓ 6 + 1) (2l 3 + 1)4π × ×E {a ℓ1µ1 a ℓ2µ2 a ℓ3µ3 a ℓ4µ4 a ℓ5µ5 a ℓ6µ6 } = 8(-1) l3-m3 ℓ1ℓ2ℓ3 C l10 ℓ10ℓ20 C l20 ℓ10ℓ30 C l30 ℓ20ℓ30 (2ℓ 1 + 1)(2ℓ 2 + 1)(2ℓ 3 + 1) (4π) 3 × × C l3-m3 l1m1l2m2 √ 2l 3 + 1 ℓ 1 ℓ 2 ℓ 3 l 3 l 2 l 1 {C ℓ1 C ℓ2 C ℓ3 } , (6.65)
see [31, p. 260 ; p. 454]. We now wish to extend these results to polyspectra of order p = 4, 5, 6 for random fields of the type T = H 2 (T G ), where (as above) T G is Gaussian, centered, isotropic and with unit variance . As anticipated, here we focus on cumulants instead of moments. We have the following result. where the reduced cumulant polyspectrum P C l1...lp (•) : l 1 , ..., l p ≥ 0 is given by

P C;1 l1l2l3l4 (λ) = 48 (2λ + 1) (4π) 4 (2l 4 + 1) ℓ1...ℓ4 C ℓ1 ...C ℓ4 C l10 ℓ10ℓ20 C l30 ℓ20ℓ30 C l40 ℓ30ℓ40 C l20 ℓ40ℓ10 ×(2ℓ 1 + 1)...(2ℓ 4 + 1)(-1) l1+l2+ℓ2+ℓ4 l 1 l 2 λ ℓ 4 ℓ 2 ℓ 1 λ l 3 l 4 ℓ 3 ℓ 4 ℓ 2 for p = 4 , P C;1 l1...l5 (λ 1 , λ 2 ) = 384 (2λ 1 + 1) (2λ 2 + 1) (4π) 5 (2l 5 + 1) ℓ1...ℓ5 C ℓ1 ...C ℓ5 C l10 ℓ10ℓ20 C l20 ℓ20ℓ30 C l40 ℓ30ℓ40 C l50 ℓ40ℓ50 C l30 ℓ50ℓ10 × ×(2ℓ 1 + 1)...(2ℓ 5 + 1)(-1) ℓ1+ℓ5+l3 l 1 l 2 λ 1 ℓ 3 ℓ 1 ℓ 2 λ 1 l 3 λ 2 ℓ 5 ℓ 3 ℓ 1 λ 2 l 4 l 5 ℓ 4 ℓ 5 ℓ 3
,for p = 5 , and

P C;1 l1...l6 (λ 1 , λ 2 , λ 3 ) = 3840 (2λ 1 + 1) (2λ 2 + 1) (2λ 3 + 1) (4π) 6 (2l 5 + 1) ℓ1...ℓ5 C ℓ1 ...C ℓ6 C l10 ℓ10ℓ20 C l20 ℓ20ℓ30 C l30 ℓ30ℓ40 C l50 ℓ40ℓ50 C l60 ℓ50ℓ60 C l40 ℓ60ℓ10 × ×(2ℓ 1 + 1)...(2ℓ 6 + 1)(-1) λ1+ℓ3+ℓ6+l4 l 1 l 2 λ 1 ℓ 3 ℓ 1 ℓ 2 λ 1 l 5 λ 2 ℓ 5 ℓ 3 ℓ 1 λ 2 l 3 l 4 ℓ 4 ℓ 5 ℓ 3 ,for p = 6 .
Proof. The result can be proved by means of the standard graphical techniques for convolutions of Clebsch-Gordan coefficients, as described in [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF]Chapters 11 and 12]. Here, we only provide the complete proof for the case p = 6. Let {a ℓm } be the random harmonic coefficients associated with the underlying Gaussian field T G . By definition, the field H 2 (T G ) admits the expansion ). The proof now follows directly from graphical techniques. In particular, the previous term can be associated with an hexagon, having in each vertex an outward line corresponding to a "free" (i.e. not summed up) index l i m i , i = 1, ..., 6. An expression for convolutions of Clebsch-Gordan coefficients corresponding to such a configuration can be found in [31, p. 461], eq. 12.1.6.30. From this, standard combinatorial arguments and a convenient relabelling of the indexes, we obtain that

H 2 (T G ) = l≥0 l m=-l a lm (2) Y lm , where a lm (2) = ℓ1m1ℓ2m2 a ℓ1m1 a ℓ2m2 S 2 Y ℓ1m1 (x) Y ℓ2m2 (x) Y lm (x)dx = ℓ1m1ℓ2m2 a ℓ1m1 a ℓ2m2 ℓ 1 ℓ 2 l m 1 m 2 -m × (-1) m × × ℓ 1 ℓ 2 l 0 0 0 (2ℓ 1 + 1) (2ℓ 2 + 1) (2l + 1) 4π = ℓ1m1ℓ2m2 a ℓ1m1 a ℓ2m2 C lm ℓ1m1ℓ2m2 C lm ℓ10ℓ20 ( 2ℓ
P C;1 l1...l6 (λ 1 , λ 2 , λ 3 ) = 3840 3 j=1 (2λ j + 1) (4π) 6 (2l p + 1) × (-1) λ1+ℓ3+ℓ6+l4 × ℓ1...ℓ6 (2ℓ 1 + 1) • • • (2ℓ 6 + 1)C ℓ1 ...C ℓ6 C l10 ℓ10ℓ20 C l20 ℓ20ℓ30 C l30 ℓ30ℓ40 C l50 ℓ40ℓ50 C l60 ℓ50ℓ60 C l40 ℓ60ℓ10 × l 1 l 2 λ 1 ℓ 3 ℓ 1 ℓ 2 λ 1 λ 2 l 3 ℓ 4 ℓ 3 ℓ 1 λ 2 l 4 λ 3 ℓ 6 ℓ 4 ℓ 1 λ 3 l 5 l 6 ℓ 5 ℓ 6 ℓ 4 .
Note that 3840 = 2 p-1 (p -1)! = 2 5 5! is the number of automorphisms between graphs belonging to M (Λ (lm)).

We recall that the Clebsch-Gordan coefficients C c0 a0b0 are identically zero unless a + b + c is even; it is hence easy to see that the previous polyspectra are non-zero only if the sum {l 1 + ... + l p } is even as well.

From the previous Proposition, we can derive the corresponding expressions for the cumulant polyspectra for χ 2 ν random field.

Definition B. We say the random field T χ 2 ν has a chi-square law with ν ≥ 1 degrees of freedom if there exist ν independent and identically distributed Gaussian random fields T i such that

T χ 2 ν law = T 2 1 + ... + T 2 ν .
It is trivial to show that T χ 2 ν is mean-square continuous and isotropic if T i is. We have the following Proposition 16 The cumulant polyspectra of T χ 2 ν (for p ≥ 2) are given by P C;ν l1...lp (λ 1 , ..., λ p-3 ) = νP C;1 l1...lp (λ 1 , ..., λ p-3 ).

Proof. Note that the cumulant polyspectra of order p ≥ 2 of T χ 2 ν coincide with those of the centered field T χ 2 ν -ET χ 2 ν (due to the translation-invariance properties of cumulants). Then, the proof is an immediate consequence of Proposition 15 and the of the standard multinearity properties of cumulants.

Further Issues and Applications

The purpose of this final Section is to introduce what we view as promising directions for further research, where the ideas of this paper may perhaps yield further insights. We shall delay to future work a more thorough investigation of the issues which are left open below.

Representations of the Symmetric Group

As a further link between representation theory and higher order angular power spectra, we mention the following. It is to be stressed that the decomposition of ∆ l1...ln that we achieved in the previous Proposition 11 is by no means unique. In particular, what we did was to choose a particular sequence of "couplings", i.e. we partitioned tensor products of the Wigner's matrices D l in a specific order before decomposing them into direct sums. Alternative partitions yield different eigenvectors and therefore, different expressions for the polyspectra/joint moments . Alternatively, we could maintain the same coupling scheme (for instance, "start always from the first pair on the left", as we did earlier) but acting on (l 1 , ..., l n ) by the symmetric group S n . However, not all coupling schemes can be achieved by simply permuting the elements of (l 1 , l 2 , ..., l n ). This is the well-known problem of parentheses in Mathematical Physics (see for instance [START_REF] Biedenharn | The Racah-Wigner algebra in quantum theory[END_REF]).

We suggest here that one can establish a link between alternate expressions for the angular polyspectra and representations of the symmetric group. More precisely the alternate expressions that we find for the polyspectra P l1....ln (λ 1 , ..., λ n-3 ) of a strongly isotropic field (with n-moments) must be such that, for every permutation π ∈ S n , λ1 ... Similarly as in the previous section, the sum of products of Clebsch-Gordan coefficients on the right hand side can be expressed in terms of higher order Wigner's coefficients. Since this section is just informal, for brevity's sake we do not give explicit expressions (see e.g. [START_REF] Varshalovich | Quantum Theory of Angular Momentum[END_REF]Chapter 10]). The two expressions (7.66) and (7.67) imply that, for every fixed (l 1 , ..., l n ) and every permutation π, there exists a square matrix A ((l 1 , ..., l n ) ; π) such that P l1....ln = A {(l 1 , ..., l n ) ; π} P π(l1)....π(ln) , where P l1...ln is the vector with entries P l1...ln (λ 1 , ..., λ n ). We conjecture that in this way one can build a representation of the symmetric group S n on the vector space generated by admissible polyspectra P l1....ln . If this is indeed the case, some important questions are left open: for instance, whether or not the representation is faithful (see [START_REF] Diaconis | Group Representations in Probability and Statistics[END_REF]), and whether these ideas can lead to algorithms for the numerical simulation of representation matrices, along the lines of what we shall pursue in the next subsection.

Random data compression

In this subsection we shall show how we can exploit the previous results to develop a probabilistic algorithm to compress information on Clebsch-Gordan coefficients. Note first that # C l3m3 l1m1l2m2 : l 1 , l 2 , l 3 ≤ L, C l3m3 l1m1l2m2 = 0 ≈ O(L 6 ) ; it is therefore clear how for most applications the storage of Clebsch-Gordan coefficients for future usage is simply unfeasible, whatever the supercomputing facilities (for instance, for CMB data analysis, L ≈ 3×10 3 is currently required, so that the number of Clebsch-Gordan coefficients to be saved would exceed 10 20 ).

Let us consider again a chi-square field as defined before, i.e. Let us assume we simulate B times T χ 2 (x), which is trivially done by simply squaring a Gaussian field: the latter is obtained by sampling independent complex Gaussian variables with variance C l . We store the triangular arrays a i lm l=1,...,L;m=-l,...,l , i = 1, ..., B; here the dimension is of order B × L 2 . We can then recover any value C l3m3 l1m1l2m2 by means of the Monte Carlo estimate

C l3m3 l1m1l2m2 = h -1 l1l2l3 B i=1 a (i) l1m1 a (i) l2m2 a (i) l3m3
B , which requires B steps and B × L 2 + L 3 storage capacity, as opposed to L 6 storage capacity by the direct method. We leave for further research a more thorough investigation on the convergence properties of this algorithm; we stress, however, that the procedure we advocate is completely general, i.e. it does not depend on peculiar features of the group SO(3) we are currently considering. We believe, hence, that similar ideas can be implemented for the numerical estimation of Clebsch-Gordan coefficients for other compact groups of interest for theoretical physicists. We leave this and the previous issues in this Section as topics for further research.

{T (x 1

 1 ), ...T (x k )} d = {T (gx 1 ), ...T (gx k )} ,(2.1)where d = denotes equality in distribution.

  g) ⊗ ... ⊗ D ln (g) ∆ l1...ln S l1...ln = D l1 (g) ⊗ ... ⊗ D ln (g) SO(3) D l1 (g) ⊗ ... ⊗ D ln (g) dg S l1...ln = SO(3) D l1 (gg) ⊗ D l2 (gg) ⊗ ... ⊗ D ln (gg) dg S l1...ln = ∆ l1...ln S l1...ln = E {a l1. ⊗ ... ⊗ a ln. } , which proves the n-th spectral moment is invariant to rotations. The fact that T is n-weakly isotropic is a consequence of the spectral representation (2.2).

∆

  l1...ln S c l1...ln = S c l1...ln (4.33) ∆ l1...ln (g * ) S c l1...ln = S c l1...ln . (4.34)

  and π ∈ Π (n), we define the matrix ∆ π l1...ln SO(3)

  two equivalent representations of SO (3)D l1 (.) ⊗ D l2 (.) ⊗ • • • ⊗ D ln (.) (5.44) and ⊕ l2+l1 λ1=|l2-l1| ⊕ l3+λ1 λ2=|l3-λ1| ... ⊕ ln+λn-2λn-1=|ln-λn-2| D λn-1 (.).(5.45)

  λ1...λn-1,µn-1 l1m1...lnmn can be represented as a convolution of the Clebsch-Gordan coefficients introduced in Section 3.2, namely: C λ1,...,λn-1,µn-1 l1m1...lnmn = C λ1,...,λn-2,. l1m1...ln-1mn-1 C λn-1µn-1 λn-2lnmn

3 C 3 C

 33 Ea l1m1 ...a lnmn = (-1) λ1...λn-3ln-mn l1m1....ln-1mn-1 P l1....ln (λ 1 , ..., λ n-3 ) (5.52) Cum {a l1m1 , ..., a lnmn } = (-1) λ1...λn-3ln-mn l1m1....ln-1mn-1 P C l1....ln (λ 1 , ..., λ n-3 ) (5.53) C λ1...λn-3;ln-mn l1m1....ln-1mn-1

S

  l1...ln = C l1...ln E l1...ln C * l1...ln S l1...ln , that is, because C l1...ln is unitary C * l1...ln S l1...ln = E l1...ln C * l1...ln S l1...ln .

2 C

 2 λ1...λn-3ln-mn l1m1....ln-1mn-1 P l1....ln (λ 1 , ..., λ n-3 ),

3 C

 3 Ea l10 ...a ln0 = λ1 ... λn-λ1...λn-3ln0 l10....ln-10 P l1....ln (λ 1 , ..., λ n-3 ) is real by definition (note indeed that the columns of C l1...ln are linearly independent). Let us illustrate the previous results by some more examples.

λn- 3 C 3 C

 33 λ1...λn-3ln0 l10....ln-10 A l1....ln (λ 1 , ..., λ n-3 ) (5.55) (resp. Cum {X l1 , • • •, X ln } = λ1...λn-3ln0 l10....ln-10 A l1....ln (λ 1 , ..., λ n-3 ). ) (5.56)

3 P

 3 l1...lp (λ 1 , ..., λ p-3 )C λ1...λp-3lp0 l10...lp-20 .

p λ1...λp- 3 P

 3 l1...lp (λ 1 , ..., λ p-3 )C λ1...λp-3lp0 l10...lp-20 .

3 P

 3 l1...lp (2l 1 + 1)...(2l p + 1) (4π) p λ1...λp-l1...lp (λ 1 , ..., λ p-3 )C λ1...λp-3lp0 l10...lp-20 = c pq l (2l + 1) 4π C l pq/2 .

Proposition 15 3 C

 153 The cumulant χ a l1m1 (2) , ..., a lpmp (2) (p = 4, 5, 6) associated with the harmonic coefficients of an isotropic random field of the type H 2 (T G ) (where T G is Gaussian and isotropic, with angular power spectrum {C l : l ≥ 0}) given by χ a l1m1 (2) , ..., a lpmp (2) = (-1) lp-mp λ1...λp-λ1...λp-3lp-mp l1m1...lp-1mp-1 × P C;1 l1...lp (λ 1 , ..., λ p-3 ) ,

λn- 3 C

 3 λ1...λn-3;ln-mn l1m1....ln-1mn-1 P l1....ln (λ 1 , ..., λ n-3 )

1 . 3 C 1 .λn- 3 P

 1313 ..λ ′ n-3 ;ln-mn π(l1)m1....π(ln-1)mn-1 P π(l1)....π(ln) (λ ′ 1 , ..., λ ′ n-3 ) . Now let us multiply both sides by C λ ′′ 1 ...λ ′′ n-3 ;lnm ′ n l1m1....ln-1mn-1 , where (λ ′′ 1 , ..., λ ′′ n-3 ) is fixed, and sum over (m 1 , ...m n ). In view of the unitary properties of Clebsch-Gordan coefficients we obtain for the left-hand side m1...mn C λ ′′ 1 ...λ ′′ n-3 ;ln-mn l1m1....ln-1mn-1 λ1...λn-3;ln-mn l1m1....ln-1mn-1 P l1....ln (λ 1 , ..., λ n-3 ) ..λ ′′ n-3 ;ln-mn l1m1....ln-1mn-1 C λ1...λn-3;ln-mn l1m1....ln-1mn-1 P l1....ln (λ 1 , ..., λ n-3 ) l1....ln (λ 1 , ..., λ n-3 ) = P l1....ln (λ ′′ 1 , ..., λ ′′ n

  )m1....π(ln-1)mn-1 P π(l1)....π(ln) (λ ′ 1 , ..., λ ′ n-3 )

  )m1....π(ln-1)mn-1 C λ ′ 1 ...λ ′ n-3 ;lnmn π(l1)m1....π(ln-1)mn-1 P π(l1)....π(ln) (λ ′ 1 , ..., λ ′ n-3 ). (7.67)

T χ 2 1 ℓ 1 ℓ 2 ℓ 3 l 1 l 2 l 3

 213 (x) = H 2 (T G (x)) = lm a lm (2)Y lm (x) ;we have proved earlier in (6.65) thatEa l1m1 (2)a l2m2 (2)a l3m3 (2) = (-1) m3 C l3m3 l1m1l2m2 h l1l2l3 where h l1l2l3 {C ℓ1 C ℓ2 C ℓ3 } ,which can be calculated analytically and stored, with storage dimension # h l1l2l3 : l 1 , l 2 , l 3 ≤ L, C l30 l10l20 = 0 ≈ O(L 3 ) .

  For a given lm = (ℓ 11 m 11 , ℓ 12 m 12 ; ...; ℓ 61 m 61 , ℓ 62 m 62 ), the quantity Cum {a ℓ11m11 a ℓ12m12 , ..., a ℓ61m61 a ℓ62m62 } is computed as follows: 11 m 11 ℓ 12 m 12 ℓ 21 m 21 ℓ 22 m 22 ℓ 31 m 31 ℓ 32 m 32 ℓ 41 m 41 ℓ 42 m 42 ℓ 51 m 51 ℓ 52 m 52 ℓ 61 m 61 ℓ 62 m 62 Define the class M (Λ (lm)) of connected, Gaussian non-flat diagrams over Λ, that is, every γ ∈ M (Λ (lm)) is a partition of the entries of Λ (lm), into pairs belonging to different rows; moreover, such a partition has to be connected, in the sense that γ cannot be divided into two separate diagrams. For instance, an element of M (Λ (lm)) is γ = {{ℓ 11 m 11 , ℓ 21 m 21 } {ℓ 22 m 22 , ℓ 32 m 32 } {ℓ 31 m 31 , ℓ 41 m 41 } {ℓ 42 m 42 , ℓ 52 m 52 } {ℓ 51 m 61 , ℓ 61 m 61 } {ℓ 62 m 62 , ℓ 12 m 12 }} Use the standard diagram formula (see again [30]), to obtain that Cum {a ℓ11m11 a ℓ12m12 , ..., a ℓ61m61 a ℓ62m62 } = where the first sum runs over all vectors of the type lm = (ℓ 11 m 11 , ℓ 12 m 12 ; ...; ℓ 61 m 61 , ℓ 62 m 62

					δ (γ) .
					γ∈M(Λ(lm))
	It follows that
	=	Cum {a l1m1 (2) , ..., a l6m6 (2)} lm γ∈M(Λ(lm)) δ (γ) 6 j=1 C lj mj ℓj1mj1ℓj2mj2 C ℓj10ℓj20 lj mj	(2ℓ j1 + 1) (2ℓ j2 + 1) 4π (2l j + 1)	,
					1 + 1) (2ℓ 2 + 1) 4π (2l + 1)	.
		6 j=1	C ℓj1mj1ℓj2mj2 C lj mj ℓj10ℓj20 lj mj	(2ℓ j1 + 1) (2ℓ j2 + 1) 4π (2l j + 1)	.
	• Build the 6 × 2 matrix	Λ (lm) =	       	ℓ        
	• • For every γ ∈ M (Λ (lm)), write
			δ (γ) =

By using once again the multilinearity of cumulants, one obtains that

Cum {a l1m1 (2) , ..., a l6m6 (2)} = ℓ11m11ℓ12m12 • • • ℓ61m61ℓ61m61 Cum {a ℓ11m11 a ℓ12m12 , ..., a ℓ61m61 a ℓ62m62 } × × {ℓ ab m ab ,ℓ cd m cd }∈γ δ ℓ ab ℓ cd δ -m cd m ab (-1) m ab C ℓ ab

(where δ b a is the usual Kronecker symbol)

•