Aliasing eﬀects for random ﬁelds over spheres of arbitrary dimension

: In this paper, aliasing eﬀects are investigated for random ﬁelds deﬁned on the d -dimensional sphere S d and reconstructed from discrete samples. First, we introduce the concept of an aliasing function on S d . The aliasing function allows one to identify explicitly the aliases of a given harmonic coeﬃcient in the Fourier decomposition. Then, we exploit this tool to establish the aliases of the harmonic coeﬃcients approximated by means of the quadrature procedure named spherical uniform sampling. Subsequently, we study the consequences of the aliasing errors in the approximation of the angular power spectrum of an isotropic random ﬁeld, the harmonic de- composition of its covariance function. Finally, we show that band-limited random ﬁelds are aliases-free, under the assumption of a suﬃciently large amount of nodes in the quadrature rule.


Overview
We are concerned with the study of the aliasing effects for the harmonic expansion of a random field defined on the d-dimensional sphere S d . A spherical random field T is a stochastic process defined over the unit sphere S d and thus depending on the location x = (ϑ, ϕ) = ϑ (1) , . . . , ϑ (d−1) , ϕ ∈ S d , where ϑ (i) ∈ [0, π], for i = 1, . . . , d − 1, and ϕ ∈ [ 0, 2π) . Harmonic analysis has been proved to be an insightful tool to study several issues related to random fields on the sphere and the development of spherical random fields in a series of spherical harmonics has many applications in several branches of probability and statistics. We are referring, for example, to the study of the asymptotic behavior of the bispectrum of spherical random fields (see [Mar06]), their Euler-Poincaré characteristic (see [CM18]), the estimation of their spectral parameters ( [DLM14]), and the development of quantitative central limit theorems for nonlinear functional of corresponding random eigenfunctions (see [MR15]). Under some integrability conditions on T (see Section 2.2), the following harmonic expansion holds: where ∈ N and m = (m 1 , . . . , m d−1 ) ∈ N d−2 ⊗ Z are the harmonic (or wave) numbers.
The set of spherical harmonics Y ,m = Y ,m1,...,m d−1 : S d → C provides an orthonormal basis for the space L 2 S d = L 2 S d , dx , where dx is the uniform Lebesgue measure over S d (see Section 2.1). The harmonic coefficients a ,m = a ,m1,...,m d−1 are given by and contain all the stochastic information of T (ϑ, ϕ). Nevertheless, the explicit computation of the integral (1) is an unachievable target in many experimental situations. Indeed, the measurements of T (ϑ, ϕ) can be in practical cases collected only over a finite sample of locations As a consequence, for any choice of and m the integral producing the harmonic coefficient a ,m is approximated by the sum of finitely many elements T (x i ), i = 1 . . . , N, the samples of the random field. As well-known in the literature, an exact reconstruction of the harmonic coefficients by means of finite sums represents a reachable target when considering band-limited random processes. Band-limited random processes are characterized by a bandwidth L 0 , so that all the harmonic coefficients for ≥ L 0 are null. A suitable choice of a sampling theorem and the cardinality of the sampling points yields the exact reconstruction for the non-null coefficients (see also, for example, [Mü07,SB93]). Further details will be discussed in Section 6.
However, if the random field is not band-limited or if the sampling theorem is not properly selected, the approximation of the integral in (1) by a finite sum can produce the so-called aliasing errors, that is, different coefficients become indistinguishable -aliases -of one another (see, for example, [Mü07,SB93]). The set of coefficients, acting as aliases of each other, depends specifically on the chosen sampling procedure.
The concept of aliasing (also known as confounding) comes from signal processing theory and related disciplines. In general, aliasing makes different signals indistinguishable when sampled, and it can be produced when the reconstruction of the signal from samples is different from the original continuous one (see, for example, [PM96, Chapter 1]).
The aliasing phenomenon arising in the harmonic expansion of a 2-dimensional spherical random field has been investigated by [LN97]. On the one hand, band-limited random fields over S 2 , which can be roughly viewed as linear combinations of finitely many spherical harmonics, can be uniquely reconstructed with a sufficiently large sample size. On the other hand, an explicit definition of the aliasing function, a crucial tool to identify the aliases of a given harmonic coefficient, is developed when the sampling is based on the combination of a Gauss-Legendre quadrature formula and a trapezoidal rule (see Section 4 for further details). In many practical applications, this sampling procedure is the most convenient scheme to perform numerical analysis over the sphere (see, for example, [AH12,SB93,Sze75]). Further reasons of interest to study the aliasing effects in S 2 have arisen in the field of optimal design of experiments. In [DMP05], designs over S 2 based on this sampling scheme have been proved to be optimal with respect to the whole set of Kiefer's Φ p -criteria, presented in [Kie74], that is, they are the most efficient among all the approximate designs for regression problems with spherical predictors.
Recently, interest has occurred in regression problems in spherical frameworks of arbitrary dimension and the related discretization problems (see, for example, [LS15]). In particular, in [DKSG18], experimental designs obtained by the discretization of the uniform distribution over S d by means of the combination of the so-called Gegenbauer-Gauss quadrature rules (see Section 3.2 for further details) and a trapezoidal rule, have been proved to be optimal with respect not only to the aforementioned Kiefer's Φ p -criteria, but also to another class of orthogonally invariant information criteria, the Φ Es -criteria. Given the increasing interest for spheres of dimension larger than 2 (see Subsection 1.2 for further details), it is therefore pivotal to carry out further investigations into the aliasing effects for random fields sampled over S d , d > 2. On the one hand, this research improves the understanding of the behavior of the approximated harmonic coefficients when computed over discrete samplings, in particular over a spherical uniform sampling (see Section 3.3). On the other hand, our investigations make extensive use of the properties of the hyperspherical harmonics, thus providing a deeper insight on their structure, carrying on with the results presented in [DKSG18].
In this paper, we work under the following assumption: a spherical random field T is observed over a finite set of locations x i ∈ S d : i = 1, . . . , N , the so-called sampling points, associated to the weights {w i : i = 1, . . . , N}. Thus, for any set of harmonic numbers and m, the approximated -or aliasedharmonic coefficient is given bỹ where τ ( , m; , m ) is the aforementioned aliasing function and is given by Further details can be found in Section 4.1. The coefficient a ,m is said to be an alias of a ,m with intensity |τ ( , m; , m )| if τ ( , m; , m ) = 0.
First, we study the general structure of the aliasing function under the very mild assumption that the sampling scheme is separable with respect to the angular coordinates, that is, the sampling points {x i : i = 1, . . . , N} can be written as follows Heuristically, a sample scheme is separable if a different discretization procedure is developed for each distinct coordinate. Then, we investigate on the explicit structure of this function and, consequently, on the identification of aliases assuming a spherical uniform design as the sampling procedure.
Second, under the assumption of isotropy, we consider the aliasing effects for the angular power spectrum of a random field, which describes the decomposition of the covariance function in terms of the frequency ≥ 0 (see Section 2.2), providing information on the dependence structure of the random field.
Third, we investigate also on the aliasing effects for band-limited random fields. More specifically, we establish suitable conditions on the sample size in order to guarantee the annihilation of the aliasing phenomenon.

Some applications and further research
An accurate characterization of the aliasing phenomena has great significance from both the points of view of theoretical statistics and its practical applications. More specifically, the analysis of spherical random fields over S d is strongly motivated by a growing set of applications in several scientific disciplines, such as Cosmology and Astrophysics for d = 2 (see, for example, [BM07,MP10] As already mentioned, aliasing phenomena can be detected in all the experimental situations where harmonic coefficients are measured by means of a discretization of the integral given by Equation (1). In this case, the presence of aliases can bring some crucial disadvantages for the experimenter.
In the classical optimal design approach (see for example [DKSG18]), the construction of experiments concerning spherical data is very sensitive to the aliasing effects. The outcomes of these experiments can be indeed affected by the aliasing of some terms belonging to the experimental design with other ones, potentially important but not included in the chosen model (see, for example, [JN11]). According for instance to [Mü07], in the construction of experimental designs for the regression of random fields, the experimenter can exploit a firstorder regression model, where interactions and aliasing are not considered. On the one hand, these designs are optimal to estimate primary effects. On the other hand, they can still present some undesirable aliasing effects producing some alias-depending bias. In this case, the information on the aliasing effects for each term is developed by means of the aforementioned aliasing function. The intensities of the aliases can be then collected in the so-called alias matrix (see, for further details, [JN11]). The alias matrix depends specifically on the experimental design; for further details, the reader is referred, for example, to [GJ11]. The construction of optimal designs minimizing the alias-depending bias subject to constraints on design efficiency, in the sense of the aforementioned optimality criteria (see, again, [DKSG18]), is therefore a topic of extreme interest (see also [JN11]).
Hyperspherical random fields in S d can be furthermore exploited to study random fields defined over the unit ball B d−1 in R d−1 , which currently represent a very challenging topic in data analysis. On the one hand, random fields defined over the unit ball B 3 are a very useful tool, aimed to generate realistic three-dimensional models from observational data in several research branches of Cosmology, and other disciplines, such as, for instance, Medical Brain Imaging, and Seismology (see, respectively, [DFH + 14, BSX + 07, LM12] and the references therein). On the other hand, the construction of a cubature formula on the unit ball is a complicated task. Indeed, even if it is theoretically known that the cubature points must correspond to the zeroes of Bessel functions of increasing degrees, in practice these points are not explicitly calculable (see, for example, [LM12]). As proved in [PX08], under some mild smoothing conditions, this issue can be overcome by linking the construction of frames and the definition of cubature formulas on B d−1 with the ones on S d . More specifically, orthogonal polynomials on the unit sphere an those on unit ball can be related by the following map linking the points in B d−1 with the one in the upper hemisphere of S d (see [PX08][Equation (4.5)]). We can thus define a distance B d−1 which corresponds to the geodesic distance on S d . The map given by (2) provides also a connection between (weighted) L p -spaces on B d−1 and L p S d . This allows one to study random fields over the unit ball by means of objects defined over spheres of higher dimension. The understanding of aliasing effects over S d becomes crucial to produce useful measurements related thus to these random fields.
By the point of view of applications, in Medical Image Analysis the statistical representation of the shape of a brain region is commonly modeled as the realization of a Gaussian random field, defined across the entire surface of the region (see for example [BSX + 07]). Many shape modeling frameworks in computational anatomy apply shape particularization techniques for cortical structures based on the spherical harmonic representation, to encode global shape features into a small number of coefficients (see [HCW + 13]). This data reduction technique, however, can not provide a proper representation with a single parametrization of multiple disconnected sub-cortical structures, specifically the left and right hippocampus and amygdala. The so-called 4Dhyperspherical harmonic representation of surface anatomy aims to solve this issue by means of a stereographic projection of an entire collection of disjoint 3-dimensional objects onto the hypersphere of dimension 4. Indeed, as aforementioned, a stereographic projection embeds a 3-dimensional volume onto the surface of a 4-dimensional hypersphere, avoiding thus, the issues related to flatten 3-dimensional surfaces to the 3-dimensional sphere. Subsequently, any disconnected objects of dimension 3 can be projected onto a connected surface in S 4 , and, thus, represented as the linear combination of hyperspherical harmonics of dimension 4 (see [HCK + 15]). Finally, further investigations can be done to study the aliasing effects arising when alternative sampling schemes to the Gauss-Gegenbauer quadrature are taken into account. For example, we refer to the so-called equiangular sampling schemes, which involve a uniform discretization of all the angular coordinates, introduced by [Sku86], and then developed, among others, by [DH94,MW11]. Another relevant sampling scheme concerns the decomposition of the hypersphere into Voronoi cells (see, for example, [NPW06]). This sampling scheme allows one to build the so-called spherical needlets, a class of spherical wavelets featuring a wide range of applications in Statistics (see, for example, [BKMP09a,DLM13,Dur16]). In view of these applications, the aliasing effects related to this sampling procedure are of vibrant interest.

Organization of the paper
This paper is structured as follows. In Section 2, we introduce some fundamental results on the harmonic analysis over the d-dimensional sphere as well as a short review of spherical random fields. Section 3 includes a short overview on the so-called Gegenbauer-Gauss quadrature formula, crucial to build a spherical uniform sampling, and provides some auxiliary results. In Section 4, we present the main findings of this work. In particular, Theorem 1 describes the construction of the aliasing function τ ( , m; , m ) under the assumption of the separability of the sampling with respect to the angular components, while Theorem 2 identifies the aliases for any harmonic coefficient a ,m when the sampling is uniform. In Section 5, we study the aliasing effects for the angular power spectrum of an isotropic random field (see Theorem 3), while in Section 6 we provide an algorithm to remove the aliasing effects for a band-limited random field sampled over a spherical uniform design (see Theorem 4). Section 7 presents an explanatory example, while Section 8 collects all the proofs.

Preliminaries
This section collects some introductory results, concerning harmonic analysis and its application to spherical random fields. It also includes a quick overview on the Gegenbauer-Gauss formula. The reader is referred to [SW71,AH12,VK91] for further details about the harmonic analysis on the sphere, to [AT07] for a detailed description of random fields and their properties, while [MP11] provides an extended description of spherical random fields over S 2 . Further details concerning the Gegenbauer-Gauss quadrature rule can be found in [AS64, AH12, SB93, Sze75].

Harmonic analysis on the sphere
Let ϑ (i) ∈ [0, π], for i = 1, . . . , d−1, and ϕ ∈ [ 0, 2π) be the spherical polar coordinates over S d . From now on, we will denote by x = (ϑ, ϕ) = ϑ (1) , . . . , ϑ (d−1) , ϕ the generic spherical coordinate, that is, the direction of a point on S d . Let the Thus, the uniform Lebesgue measure dx over S d , namely, the element of the solid angle, is defined by such that the surface area of the hypersphere corresponds to where Γ denotes the Gamma function. Let us denote by H the restriction of the space of harmonic homogeneous polynomials of order to S d . As well-known in the literature (see, for example, [AH12,SW71]), the space of square-integrable functions over S d can be described as the direct sum of the spaces H , that is, H .
For any integer ≥ 0, from now on called frequency, we define the following set where, for x ∈ S d , Y ,m = Y ,m1,...,m d−1 : S d → C denotes the so-called spherical -or hyperspherical -harmonic of degree and order m. In other words, fixed ≥ 0, M appoints the finitely many vectors m which identify the spherical harmonics spanning the space H .
Another common approach to introduce spherical harmonics exploits the socalled d-spherical Laplace-Beltrami operator Δ S d (see, for example, [MP11]). Fixed ≥ 0, the spherical harmonics Y ,m (x) corresponding to any m ∈ M are the eigenfunctions of As proved for example in [AW82], for any ≥ 0, the size of {Y ,m : m ∈ M }, namely, the multiplicity of the set of spherical harmonics with eigenvalue ε ;d , is given by The set {Y ,m (x) : ≥ 0; m ∈ M } provides therefore an orthonormal basis for L 2 S d . For any g ∈ L 2 S d , the following Fourier -or harmonic -expansion holds where {a ,m : ≥ 0; m ∈ M } are the so-called harmonic coefficients, given by the integral From now on, for the sake of notational simplicity, we fix m 0 = . Furthermore, we will use indifferently the two equivalent short and long notations , the hyperspherical harmonics are defined by , is the Gegenbauer (or ultraspherical) polynomial of degree n and parameter α. Following for example [AS64,Sze75], they are orthogonal with respect to the measure Roughly speaking, each hyperspherical harmonic in (6) can be viewed as product of a complex exponential function and a set of Gegenbauer polynomials, whose orders and parameters are properly nested and normalized to guarantee orthonormality, that is,

C. Durastanti and T. Patschkowski
Hyperspherical harmonics feature also the following property, known as addition formula (see, for example, [AW82] where ·, · is the standard inner product in L 2 R d+1 . Note that K can be viewed as the kernel of the projector over the harmonic space H , the restriction to the sphere of the space of homogeneous and harmonic polynomials of order . The projection P of g ∈ L 2 S d onto H is given by It follows that and that any function g ∈ L 2 S d can be rewritten as the sum of projections over the spaces H ,

Spherical random fields
Given a probability space {Ω, F, P}, a spherical random field T ω (x), ω ∈ Ω and x ∈ S d , describes a stochastic process defined the sphere S d . From now on, the dependence on ω ∈ Ω will be omitted and the random field will be denoted by T (x), x ∈ S d , for the sake of the simplicity (see also [AT07]).
If T has a finite second moment, that is, spherical random field can be decomposed in terms of the projections over the space H , ≥ 0, so that where T (x) = P [T ] (x). Each projector onto H can be described as a linear combination of finitely many hyperspherical harmonics, As in the deterministic case described in Section 2.1, for any ≥ 0 and m ∈ M , the random harmonic coefficient is defined by The random harmonic coefficients contain all the stochastic information of the random field T , namely, a ,m = a ,m (ω), for ω ∈ Ω, ≥ 0 and m ∈ M . A random field is said to be band-limited if there exists a bandwidth L 0 ∈ N, so that a ,m = 0 for any > L 0 , whenever m ∈ M . In this case, it holds that By the practical point of view, band-limited random fields provide a useful approximation of fields with harmonic coefficients decaying fast enough as the frequency grows. Let us define the expectation μ (x) = E [T (x)]; the covariance function Υ : where, for z ∈ C,z denotes its complex conjugate. Without losing any generality, assume that T is centered, so that, for x, x ∈ S d , it holds that A spherical random field is said to be isotropic if it is invariant in distribution with respect to rotations of the coordinate system or, more precisely, where d = denotes equality in distribution, and SO (d + 1) is the so-called special group of rotations in R d+1 . Following [BKMP09b,BM07,MP11], if the random field is isotropic, then Υ depends only on γ and its variance σ 2 (x) = Υ (x, x) does not depend on the location x ∈ S d , so that it holds that where σ 2 ∈ R + . The covariance function itself can be therefore rewritten in terms of its dependence on the distance between x and x , so that Let us finally define the correlation function ρ : [−1, 1] → [−1, 1], which is invariant with respect to rotations when the random field is isotropic, that is As far as the random harmonic coefficients {a ,m : ≥ 0, m ∈ M } are concerned, since μ (x) = 0 for x ∈ S d , we have that E [a ,m ] = 0. On the one hand, the Fourier expansion of T can be read as a decomposition of the field into a sequence of uncorrelated random variables, preserving its spectral characteristics, that is, where {C : ≥ 0} is the so-called angular power spectrum of T .
On the other hand, the spectral decomposition of the covariance function is given by where we rewrite the covariance function in terms of the projection kernel corresponding to the frequency level . Combining (9), (14) and (16), the angular power spectrum of a random field can be viewed as the harmonic decomposition of its covariance function and can be rewritten as the average where Ξ d ( ) is given by (5) (see, for example, [Mar06] for d = 2).

The Gauss-Gegenbauer quadrature formula and the spherical uniform design
This section includes a quick overview on the Gegenbauer-Gauss formula. We also introduce the spherical uniform sampling and two related auxiliary results. Further details concerning the Gegenbauer-Gauss quadrature rule can be found in [AS64, AH12, SB93, Sze75], while the spherical uniform sampling is presented by [DKSG18].

Separability of the sampling
We first introduce a very mild condition on the sampling procedure. Generalizing the proposal introduced by [LN97] on S 2 to S d , d > 2, here we consider a discretization scheme produced by the combination of d one-dimensional quadrature rules, with respect to the coordinates ϑ (j) , j = 1, . . . , d − 1, and ϕ. More specifically, we introduce the following condition on the sampling points and weights.
Condition 1 (Separability of the sampling scheme).
so that The sampling points {x i : i = 1, . . . , N} are component-wise given by Roughly speaking, each sequence in (18) corresponds to the set of weights for a quadrature formula with respect to the j-th angular component of the angle Each value of the index i * ∈ {1, . . . , N} corresponds uniquely to a suitable choice of values k * 0 , . . . , k * d−1 , while the related weight w i * is given by

The Gauss-Gegenbauer quadrature formula
In general, a quadrature rule denotes an approximation of a definite integral of a function by means of a weighted sum of function values, estimated at specified points within the domain of integration (see, for example, [SB93]). In particular, a r-point Gaussian quadrature rule is a formula specifically built to yield an exact result for polynomials of degree smaller or equal to 2r − 1, after a suitable choice of the points and weights {t k , ω k : k = 0, . . . , r − 1}. For this reason, it is also called quadrature formula of degree 2r − 1. The domain of integration is conventionally taken as [−1, 1], and the choice of points and weights usually depends on the so-called weight function a, whereas the integral can be written in the form Here p (t) is approximately polynomial, and a (t) ∈ L 1 ([−1, 1]) is a well-known function. In this case, a proper selection of From now on, while the letter ω will concern weights related to quadrature formulas for coordinates on the interval [−1, 1], the letter w will denote weights related to quadrature formulas for angular coordinates.
Following for example [SB93], it can be shown that the quadrature points can be chosen as the roots of some polynomial belonging to some suitable class of orthogonal polynomials, depending on the function a.
When a (t) = 1 for all t ∈ [−1, 1], the associated polynomials are the Legendre polynomials. In this case, the method is then known as Gauss-Legendre quadrature (see [AS64,Formula 25.4.29]). Such a method is widely used in the 2-dimensional spherical framework (see, for example, [AH12]), and the aliases produced by this formula were largely investigated in [LN97]).
More in general, as stated in [AS64,Formula 25.4.33], when a (t) = a α, , the method is known as the Gauss-Jacobi quadrature formula, since it makes use of the Jacobi polynomials (see also [Sze75,p.47]). Since it is well-known that Jacobi polynomials reduce to Gegenbauer polynomials when α = β (see, for example, [Sze75, Formula 4.1.5]), we refer to the quadrature rule denoted by a weight function ν α (t) (equal to a α,β (t) for α = β) as the Gauss-Gegenbauer quadrature (see, for example, [ESM14]).
Subsequently, the discrete uniform sampling over the sphere is obtained by combining a trapezoidal rule for the angle ϕ and (d − 1) Gauss-Gegenbauer quadrature rules for the coordinates ϑ (j) , for j = 1, This method has been described in details by [DKSG18,Lemma 3.1] in the framework of optimal design for regression problems with spherical predictors. Indeed, by the theoretical point of view, the (continuous) uniform distribution on the sphere provides an optimal design for experiments on the unit sphere, but this distribution is not implementable as a design in real experiments (for more details, see [DKSG18, Theorem 3.1]). Thus, a set of equivalent discrete designs is established by means of the combination of the following quadrature formulas over the sphere, written as in [DKSG18, Lemma 3.1], to which we refer to for a proof.
if and only if the following conditions are satisfied: 1. The polynomial is orthogonal to all polynomials of degree smaller or equal to z − r with respect to a (t), 2. the weights ω k are given by where λ k (t) is the k-th Lagrange interpolation formula with nodes t 0 , . . . , t r−1 , given by

The spherical uniform sampling
Assume now z = 2Q 0 in Definition 1. Following [Sze75, Formula 4.7.15] (see also (8)), the Gegenbauer polynomials C In Lemma 1 above, we have recalled a set of quadrature formulas for the interval [−1, 1], each of those associated to the corresponding weight function ν α(j) , for j = 1, . . . , d − 1. The following Condition exploits properly these quadrature formulas for ϑ, combined with a trapezoidal rule for ϕ, to establish a well-defined uniform distribution over the sphere of arbitrary dimension d (see also, for example, [AH12,DKSG18]). Observe that this choice yields a suitable quadrature formula for each angular component in S d .
As already discussed in [AH12,DKSG18], the Gauss-Gegenbauer quadrature in Lemma 1 is characterized by a unitary sum of the weights for each component, while Condition 2 guarantees orthonormality for spherical harmonics Y ,m and Y ,m so that + ≤ 2Q 0 , that is, We present now two auxiliary results crucial to prove Theorem 2, referring to the aliasing effects under Condition 2. Their proofs can be found in Section 8.2.
The first Lemma establishes the parity properties of the cubature points and weights for each angular component ϑ (j) with respect to ϑ (j) = π/2, for j = 1, . . . , d − 1. Indeed, due to the parity formula C The next result exploits Lemma 2 to develop parity properties on the Gauss-Gegenbauer quadrature formula.
Then it holds that

Aliasing effects on the sphere
This section presents our main results concerning the aliasing phenomenon for d-dimensional spherical random fields. First, we define the aliasing function, the key tool to explicitly determine the aliases for any given harmonic coefficient. Then, we study the aliasing function and the set of harmonic numbers identifying the aliases for any given coefficient a ,m in two different cases. The proof of the theorems presented in this section are collected in Section 8.1. As a first step, we just assume that the aliasing function is separable with respect to the angular components. This assumption is very mild, as it reflects both the separability of the spherical harmonics and the practical convenience of choosing separable sampling points, with respect to the angular coordinates.
As a second step, we study the aliasing effects under the assumption that the sample comes from a spherical uniform design.

The aliasing function
In practical applications, the measurements of the random fields can be sampled only over a finite number of locations on S d . As a straightforward consequence, the integral (12) can not be explicitly computed, but it has to be replaced by a sum of finitely many samples of T .
where f (ϑ) is given by (3). Combining (10) and (11) with (28) yields where τ ( , m; , m ) is given by From now on, we will refer to τ ( , m; , m )  As stated by [LN97] for the case d = 2, on the one hand, the following equality is a necessary and sufficient condition to identify a ,m andã ,m . This equality does not hold in general (see Section 6). On the other hand, fixed , , m and m , if τ ( , m; , m ) = 0, that is, a ,m is an alias of a ,m , its intensity, denoting how large is the contribution of this alias, is given by |τ ( , m; , m )|. The total amount of aliases in (29) and the corresponding intensity depends specifically on the choice of the sampling points {x i : i = 1, . . . , N} over S d , which characterizes entirely the subsequent structure of (30). In other words, every setting chosen for the sampling points leads to a specific set of aliases, described by the corresponding aliasing function.
Here we study the aliasing function τ ( , m; , m ) first in a more general framework, under the assumption of a separable sampling with respect to the angular coordinates in Section 4.2, and then for a discrete version of the spherical uniform distribution in Section 4.3.

The separability of the aliasing function
Let us assume now that the assumptions of Condition 1 hold. Thus, given Q 0 , Q 1 , . . . , Q d−1 ∈ N, so that N = d−1 j=0 Q j , for j = 1, . . . , d − 1, the corresponding set of quadrature points and weights is given by As a straightforward consequence, we obtain the following result.

Aliasing and spherical uniform designs
As already mentioned in Section 1.1, the motivations behind the study of this particular setting come from two different sources. On the one hand, the uniform design is largely used in the framework on numerical analysis over the sphere (see [AH12,SB93,Sze75]). On the other hand, in the field of mathematical statistics, the spherical uniform sampling has be proved to be the the most efficient design with respect to a large set of optimality criteria such as the Kiefer's Φ p -as well as the Φ Es -criteria, in the framework of optimal designs of experiments (see [DKSG18]). Furthermore, in Remark 3, we show that our findings align with the results established [LN97]) for the two-dimensional case. The example described in Section 7 establishes explicitly the set of aliases of a given harmonic coefficient.
The main results of this section, stated in the forthcoming Theorem 2, require some further notation, produced in Remark 2.
Let us now define the following sets, and, for j = 1, . . . , d − 2, Observe that the definition of A j and B j is formally correct to take into account all the possible combinations of s j−1 and Q j . It is straightforward to observe that

Define now the following sets
while In other words, when s j ∈ Δ j , it can take any value in H

will be labeled as secondary locations. According to Corollary 1, a proper choice of the sampling points can annihilate the aliases having secondary locations. The same does not hold for the ones in the primary locations. It is indeed impossible to remove all the aliases in primary locations just by choosing the sampling points and parameters. As we will discuss in Section 6, these aliases can be completely erased, after a proper selection of sampling points, only if the random field is band-limited.
Finally, note that under the assumptions of Corollary 1, it holds that

Aliasing for angular power spectrum
In this section, our purpose is to investigate on the aliasing effects as far as the spectral approximation of an isotropic random field is concerned. More specifically, we establish a method to identify the aliases of each element of the power spectrum {C : ≥ 0}. Assume to have an isotropic random field on S d , so that (15) and (16) hold. When the integral (12) is replaced with the sum (29) under the Condition 2, we want to study how the aliasing errors arising in (29), affect the estimation of C = Var (a ,m ) (see (16)). In particular we are interested in developing the presence of aliases when C is approximated by the averagẽ where Ξ d ( ) is given by (5) (cf, for example, (17)). Let us recall that D 0 ( ) is given by (34), and let V Q , m ( ) be defined by Our findings, which extend to the d-dimensional sphere the outcomes of [LN97, Theorem 3.1] (cf. Remark 3), are produced in the following theorem.
Theorem 3. Let T be an isotropic random field on S d with angular power spectrum given by (16). Under the assumption given in Condition 2, it holds thatC = s0∈D0( ) The proof of Theorem 3 can be found in Section 8.1.

Band-limited random fields
In this section, we establish the condition on the sample size, leading to an exact reconstruction of the harmonic coefficients a ,m for band-limited random fields, in the paradigm of the spherical uniform design. In other words, for band-limited random fields and for a suitable choice of Q, the approximation of the integral (12) by the sum (28) is exact and, then, there are no aliases, analogously to the findings described in [LN97, Section 4] for d = 2. The reader is referred to Section 8.1 for the proofs of the theorems collected in this section. If the number of sampling points is sufficiently large with respect to the bandwidth characterizing the random field, we obtain two crucial results, stated in the next theorem. On the one hand, the band-limited random fields are aliasfree inã ,m and, on the other, they are exactly reconstructed by means of the Gaussian quadrature procedure described above.
Theorem 4. Assume that T (x) is band-limited with bandwidth L 0 , that is, the harmonic expansion given by (13) holds. If also Condition 2 holds, with Furthermore, for any L ∈ N satisfying Q ≥ L ≥ L 0 , the following reconstruction holds exactly: where x k = ϑ k0,...,k d−2 , ϕ k d−1 and K is given by (9).
Remark 6. In view of the results presented in Theorem 4, the sample size N has to satisfy the following condition, in order to avoid aliasing effects for band-limited random processes with bandwidth L 0 .
Remark 7. If the random field is band-limited, the only possible aliases belong to secondary locations (see Remark 5). Thus, a suitable choice of the parameters Q 0 , . . . , Q d−2 , M annihilates all the potential aliases.
A random field has a band-limited power spectrum with bandwidth P L if C = 0 for any > P L . The following theorem shows that these random fields are aliases-free inC , employing a Gauss sampling under Condition 2 and given a suitable sample size.

C. Durastanti and T. Patschkowski
s0,s1,rM I Q 0,0 (2s 0 , 2s 1 ) I Q 0,0 (2s 1 , 2rM ) Observe that the first line in (54) describes the aliases obtained for s 0 ∈ A 0 , while the other two lines contain the aliases corresponding to s 0 ∈ B 0 . Notice that if s 0 ∈ A 0 , then B 1 = ∅. As a consequence, it follows that both the indexes s 1 and s 2 can not take the null-value. When s 0 ∈ B 0 , we have that A 1 = {0, . . . , Q − 1} and B 1 = {Q, . . . , s 0 }. Hence, we obtain the second and the third sums in (54). We want to establish here the locations of the aliases that affect a 0,0,0 for some choices of Q and M . Let us take Q = 2, 4 and M = Q/2, Q. Here, for the sake of the computational simplicity, we will take into account only s 0 = 1, . . . , 4. All the aliases of a 0,0,0 for the considered range of s 0 are collected in Table 1 and in Table 2. For any choice of Q and M , each column contains aliases belonging to the sets {s 0 ∈ A 0 , s 1 ∈ A 1 }, {s 0 ∈ B 0 , s 1 ∈ A 1 }, and {s 0 ∈ B 0 , s 1 ∈ B 1 } respectively. The locations of the aliases are also shown in Figure 1, for Q = 2, and Figure 2, for Q = 4.
According to the results here produced, we can notice that • the minimum distance of the aliases increases when Q grows, following Remark 4 and Equation (50) in Remark 5; • all the aliases with secondary locations (see Remark 5), belonging thus to the subsets {s 0 ∈ A 0 , s 1 ∈ A 1 } and {s 0 ∈ B 0 , s 1 ∈ A 1 }, vanish for M = Q, as stated in Corollary 1; • the coefficient a 0,0,0 is not affected by aliasing if it is the harmonic coefficient of a band-limited function with band width L 0 < Q, as stated in Theorem 4.

Proofs
In this section, we provide proofs for the main and auxiliary results.

Proofs of the main results
Proof of Theorem 1. Using (3) . . .
as claimed. For both cases, we follow a backward induction step, studying first the aliasing effects due to the trapezoidal sampling for coordinate j = d, using the results holding for the j-th component to prove the statement for the j − 1-th component, until we reach j = 1.
Part 1 -Here our purpose it to exploit either properties due to the uniform sampling and the ones related to the harmonic numbers of spherical harmonics, to establish lower and, where possible, upper bounds for the indices s 0 , . . . , s d−2 , r. These indices identify the aliases of the harmonic coefficient a ,m , given in the form a +2s0,m+2s . Let us consider initially j = d and apply to the coordinate ϕ the standard trapezoidal rule. As well as in [LN97] (see also [DKSG18]), using (22) and (23) in (32) yields where r ∈ Z is such that Consider now j = d − 1. The component ϑ (d−1) is subject to the aforementioned Gauss-Legendre quadrature formula (cf. the case d = 2 in [LN97]). Indeed, by using (55) jointly with the definition of the sampling points and weights given by (24) and (25) respectively with j = d − 1, the (d − 1)-th aliasing factor is given by (56) Observe now that the Legendre polynomials can be expressed in terms of a Gegenbauer polynomial by means of the formula , see for example [Sze75,Formula 4.7.35]. Hence, we obtain that where In analogy to [LN97, Theorem 2.1], using (26), given in Lemma 3, for j = d − 1, in (57) leads to In other words, the d − 1-th aliasing factor is not null only for even values of where s d−2 ∈ D m d−2 , given by which guarantees that m d−2 ≥ 0 and, thus, a well-defined aliasing factor in (56). On the one hand, using m d−2 = m d−2 + 2s d−2 in the set concerning the d-th aliasing factor, we have that r ∈ R M m d−1 (m d−2 + 2s d−2 ), as given by (36). On the other hand, following (4) and (6), it holds that Therefore we obtain that Consider now 2 ≤ j ≤ d − 2. For each component, we use a suitable Gauss-Gegenbauer quadrature rule described above (see also [DKSG18, On the one hand, Formula (27) in Lemma 3 with m j = m j + 2s j yields I Qj−1 mj−1,mj m j−1 , m j + 2s j = 0 only for m j−1 = m j−1 + 2s j−1 , so that the aliases with respect to the j-th component are identified by the function It is straightforward to set s j−1 ∈ D mj−1 , where so that the polynomials in I On the other hand, taking into account (4) and (6), it follows that m j−1 = m j−1 + 2s j−1 ≤ m j−2 . Thus we obtain that s j−1 ∈ R mj−1 m j−2 , where with m j−2 = m j−2 + 2s j−2 . Combining these two results and recalling (35), for j = 2, . . . , d − 1, it holds that Furthermore, the following step of the backward procedure yields m j−2 = m j−2 + 2s j−2 , so that for j = 2, . . . , d− 1. Consider, finally, the case j = 1. This aliasing factor is given by I Q0 ,m1 ( , m 1 + 2s 1 ) for s 1 ∈ H (1) m1 ( ) .
Here we can thus select = + 2s 0 , s 0 ∈ D 0 ( ), where D 0 ( ) is given by (34). Note that s 0 is the only index that is not selected from a set of finitely many elements.
Part 2 -Here our aim is to use the order of the used quadrature formula to convert, when possible, the sums of I Qj−1 mj−1,mj m j−1 , m j to integrals. Then, we exploit the orthogonality of the Gegenbauer polynomials (see Section 2) to establish further combinations of indices s 0 , . . . , s d−1 , r which lead to a null aliasing function.
Recombining all these results for j = 1, . . . , d yields the fact that the aliases a +2s0,m+2s exist for s ∈ Z Q ,m , where Z Q ,m is defined by (47), as well as for s 0 ∈ D 0 ( ) (cf. Part 1), as claimed.
Proof of Theorem 3. Let us fix ≥ 0 and m ∈ M , and recall furthermore that the random variables a +2s0,m+s , s 0 ∈ D 0 ( ) , s ∈ Z Q ,m are uncorrelated with variance C +2s0 . The variance ofã ,m is, thus, given by Var (ã ,m ) = Using now (60) leads to s d−2 = 0. Reiterating this backward procedure for the other harmonic numbers m j , j = d − 3, . . . , 1 and yields (52).
To prove (53), it suffices to use the band-width in the expansion (11), that is, Using now in the equation above (29), (49), and (52) yields the claimed result.

Proofs of the auxiliary results
Proof of Lemma 2. The symmetry of the sampling angles follow the symmetry of the roots of the Gegenbauer polynomials. Furthermore, note that sin ϑ Then, we have that Qj−1−kj−1−1 , as claimed. Proof of Lemma 3. First of all, note that this result for d = 2, involving thus Legendre polynomials, has been already claimed in [LN97, Theorem 2.1].
As far as d > 2 is concerned, let us preliminarily recall that, for t ∈ [−1, 1], C as claimed.
In order to prove (27), consider initially only even values of Q. Hence, by means of Lemma 2, we have that Moreover, if Q is odd, since sampling points have to be symmetric with respect to π/2, the additional point with respect to the previous case has to coincide with π/2. Thus G (π/2) = 0 and (27) holds, as claimed.
Proof of Corollary 1. This proof follows directly from the proof of Theorem 2-Part 2. Indeed, if M ≥ Q, it follows that r = 0. Then, combining (58), (59) and (60) yields the claimed result.