Regularity properties and simulations of Gaussian random fields on the Sphere cross Time

We study the regularity properties of Gaussian fields defined over spheres cross time. In particular, we consider two alternative spectral decompositions for a Gaussian field on $\mathbb{S}^d \times \mathbb{R}$. For each decomposition, we establish regularity properties through Sobolev and interpolation spaces. We then propose a simulation method and study its level of accuracy in the $L^2$ sense. The method turns to be both fast and efficient.


Introduction
Spatio-temporal variability is of major importance in many fields, in particular for anthropogenic and natural processes, such as earthquakes, geographic evolution of diseases, income distributions, mortality fields, atmospheric pollutant concentrations, hydrological basin characterization and precipitation fields, among others. For many natural phenomena involving, for instance, climate change and atmospheric variables, many branches of applied sciences have been increasingly interested in the analysis of data distributed over the whole sphere representing planet Earth and evolving through time. Hence the need for random fields models where the spatial location is continuously indexed through the sphere, and where time can be either continuous or discrete. It is common to consider the observations as a partial realization of a spatio-temporal random field which is usually considered to be Gaussian (see [6,7,10]). Thus, the dependence structure in space-time is governed by the covariance of the spatiotemporal Gaussian field, and we refer the reader to [13,24,27] for significant contributions in this direction.
Specifically, let d be a positive integer, and let S d = {x ∈ R d+1 , x = 1} be the d-dimensional unit sphere in the Euclidean space R d+1 , where · denotes the Euclidean norm of x ∈ R d+1 . We denote Z = {Z(x, t), (x, t) ∈ S d × R} a Gaussian field on S d × R.
The tour de force in [17] characterizes isotropic Gaussian random fields on the sphere S d through Karhunen-Loève expansions with respect to the spherical harmonics functions and the angular power spectrum. They show that the smoothness of the covariance is connected to the decay of the angular power spectrum and then discuss the relation to sample Hölder continuity and sample differentiability of the random fields.
The present paper extends part of the work of [17] to space-time. Such extension is non-trivial and depends on two alternative spectral decompositions of a Gaussian field on spheres cross time. In particular, we propose either Hermite or classical Karhunen-Loève expansions, and show how regularity properties can evolve dynamically over time. The crux of our arguments rely on recent advances on the characterization of covariance functions associated to Gaussian fields on spheres cross time (see [3] and [21]). Notably, the Berg-Porcu representation in terms of Schoenberg functions inspires the proposal of alternative spectral decompositions for the temporal part, which become then crucial to establish the regularity properties of the associated Gaussian field.
The second part of the paper is devoted to simulation methods which should be computationally fast while keeping a reasonable level of accuracy, resulting in a notable step forward. Efficient simulation methods for random fields defined on the sphere cross time are, currently, almost unexplored. Cholesky decomposition is an appealing alternative since it is an exact method. However, the method has an order of computation of O(N 3 ), with N denoting the sample size. This makes its implementation computationally challenging for large scale problems (the so called Big "n"problem). It is therefore mandatory to investigate efficient simulation methods. Here, we propose a simulation method based on a suitable truncation of the proposed double spectral decompositions. We establish its accuracy in the L 2 sense and illustrate how the model keeps a reasonable level of accuracy while being considerably fast, even when the number of spatiotemporal locations is very high.
The remainder of the article is as follows. Section 2 provides the basic material for a self-contained exposition. The expansions for the kernel covariances and the random field are presented in Section 3. Section 4 presents the regularity results for the kernel covariance functions in terms of weighted Sobolev spaces and weighted bi-sequence spaces. In Section 5 our simulation method is developed and its accuracy is studied. Also, we provide some numerical experiments for illustrative purposes. In Appendix 7.1 we also provide a rather general version of the Karhunen-Loève theorem.
The manuscript is intended for complex-valued random fields over S d × R with d ∈ N, except for Section 5 where the simulations considered are for real-valued random fields over S 2 × R.

Preliminaries
This section is largely expository and devoted to the illustration of the framework and notations that will be of major use throughout the manuscript. All the tools presented in this section are valid in S d , for any d ∈ N. Some particular references to the case d = 2 are exposed, as they will be of use in Section 5.

Spherical Harmonics Functions and Gegenbauer Polynomials
Spherical harmonics are restrictions to the unit sphere S d of real harmonics polynomials in R d+1 . They are also the eigenfunctions of the Laplace-Beltrami operator on S d . A deeper overview on spherical harmonics along with the properties listed in this subsection can be found in [9].
where σ d denotes the surface area measure, and σ d denotes the surface area of S d , with Γ being the Gamma function.
For j = 0, 1, . . . , let H d j denote the linear space of spherical harmonics of degree j over S d . For different degrees, spherical harmonics are orthogonal with respect to the inner product of where δ j,j is the Kronecker delta function, being identically equal to one if j = j , and zero otherwise.
constitutes an orthonormal basis of L 2 (S d ) and Theorem 2.42 of [19] shows that Besides, the addition formula for spherical harmonics states where ·, · R d+1 denotes the inner product in R d+1 . Here C r j (·) are the Gegenbauer (or ultraspherical) polynomials of degree j and order r, defined by where P (α,β) j (·) denotes the Jacobi polynomial of parameters α, β > −1 and degree j, and (·) j denotes the Pochhammer symbol (rising factorial) defined by where (r) j = Γ(r + j)/Γ(r), provided r is not a negative integer.
Gegenbauer polynomials constitute a basis of L 2 (−1, 1) and satisfy the orthogonality relation (section 3.15 in By Stirling's inequalities, for n ∈ N fixed, there exists constants c 1 (n) and c 2 (n) such that Hence, assuming 2r ∈ N, relation (2.3) becomes Besides, for n ≤ j we observe that (see section 10.9 in [2]) In what follows c j (d, x) denotes the standardized Gegenbauer polynomial, being identically equal to one for x = 1 and r = (d − 1)/2, that is It is straightforward to see that Remark 2.1. When d = 2, i.e., the S 2 case, C 1/2 j (·) = P j (·), where P j (·) is the Legendre polynomial of degree j [9].
See [9] and [19] for a deeper overview of spherical harmonics, and [2] for a detailed description of Gegenbauer, Jacobi and Legendre polynomials.

Isotropic Stationary Gaussian Random Fields on the Sphere cross Time
Let (Ω, F, P) be a complete probability space. Consider S d and S d × R as manifolds contained in R d+1 and in R d+2 respectively. Definition 2.1. The geodesic distance θ on S d ("great circle" or "spherical" distance) is defined by Throughout, unless it is explicitly presented in a different way, we write S d instead of (S d , σ d ), σ d being the surface area measure, which is equivalent to other uniformly distributed measures on S d , such as the Haar measure or the Lebesgue spherical measure [8]. Analogously, we write R instead of (R, µ), for µ being the Lebesgue measure.
is multivariate Gaussian distributed. Here, denotes the transpose operator.
for all finite systems of pairwise distinct points {(x k , t k )} n k=1 ⊂ S d × R and constants c 1 , . . . , c n ∈ C. A positive definite function h is strictly positive definite if the inequality (2.9) is strict, unless c 1 = · · · = c n = 0.
We call a function h : S d × R × S d × R −→ C spatially isotropic and temporally stationary if there exists a function ψ : (2.10) Hence, a spatially isotropic and temporally stationary function depends on its arguments via the great circle distance θ(x, y) and the time lag, or equivalently, via the inner product x, y R d+1 and the time lag. Definition 2.3. We call a random field Z = Z(x, t) : is constant for all (x, t) ∈ S d × R, and if the covariance h is a spatially isotropic and temporally stationary function on (S d × R) 2 . The associated function ψ in (2.10) is called the covariance kernel or simply kernel.
Remark 2.2. A Gaussian random field (GRF) which is 2-weakly isotropic stationary on S d ×R, is in fact isotropic in the spatial variable and stationary in the time variable (see [18]), hence, it has an invariant distribution under rotations on the spatial variable, and under translations on the temporal variable.
Throughout the manuscript we work with zero-mean random fields, with no loss of generality.

Kernel Covariance Functions on the Sphere cross Time
In his seminal paper, [23] In particular, the case G = R offers a characterization of spatio-temporal covariance functions of centred 2-weakly isotropic stationary random fields over the sphere cross time.
Let P(R) denote the set of continuous and positive definite functions on R.
For d = 1, 2, . . . , we consider the class P(S d , R) of continuous functions ψ : [−1, 1]×R → C such that the associated spatially isotropic temporally stationary function h in (2.10) is positive definite. The following result, rephrased from [3], allows to identify the class P(S d , R) with the covariances functions of 2-weakly isotropic stationary random fields on S d × R.
• We note that, in comparison with the representation of covariance functions for 2-weakly isotropic random fields on the sphere, representation (3.1) does not consider Schoenberg coefficients ϕ j,d but Schoenberg functions ϕ j,d (·), which will play a fundamental role subsequently.

Expansions for Isotropic Stationary GRFs and Kernel Covariance Functions
According to Theorem 2.1, the kernel ψ of an isotropic stationary GRF Z on S d × R admits the following representation: where {ϕ j,d } j∈N is a sequence of functions in P(R), such that the series is uniformly convergent.
Expression (3.1) allows to consider different expansions for the kernel. Before introducing these representations, we present the expansion for the random field which motivates the simulation methodology.

Karhunen-Loève Expansions for Isotropic Stationary GRFs on the Sphere cross Time
In [15], the following Karhunen-Loève representation for an isotropic stationary GRF Z over S 2 × R is proposed, where {X j,m (t)} j,m is a sequence of one-dimensional complex-valued mutually independent stochastic processes. The set of all X j,m (t) forms a denumerable infinite dimensional stochastic process which completely defines the process on the sphere, and Y j,m (x) are the elements of an orthonormal basis of H 2 j . Representations (3.1) for the spatio-temporal covariance and (3.2) for the random field, allow us to introduce the following family of GRFs on S d × R for any d ∈ N. Definition 3.1. Let Z be a random field on S d × R defined, in the mean square sense, as Here, for each j, m and d, {X j,m,d (t), t ∈ R} is a complex-valued zero-mean stationary Gaussian process such that where {ϕ j,d } j∈N represents the Schoenberg's functions associated to the mapping ψ in Equation (3.1).
Remark 3.1. Examples of random field satisfying Definition (3.1) are found in the Appendix of [21]. Proposition 3.1. Let Z be a random field as in Definition 3.1. Then, Z is an isotropic stationary GRF, with zero-mean and covariance function given by ψ as in (3.1).
Proof. The proof follows straight by using the properties of the process X j,m,d , and the addition formula in Equation (2.2).
For each j, m, d, the process {X j,m,d (t), t ∈ R} has the following Karhunen-Loève expansion (see Appendix 7.1), Remark 3.2. By (3.4), an alternative way to write (3.3) is: Expressions (3.3) or (3.5) represent a way to construct isotropic stationary GRFs on the sphere cross time, and suggest a spectral simulation method. However, it is not yet proved that any isotropic stationary GRF on the sphere cross time can be written in this way.

Double Karhunen-Loève Expansion of Kernel Covariance Functions
By stationarity of the process X and Karhunen-Loève theorem we have, Hence, a slight abuse of notation allows to reformulate the last expression to Therefore, the kernel covariance function ψ in (3.1) also admits the expansion Following [17], we call the series {a j,k,d } j,k∈N ⊂ R the spatio-temporal angular power spectrum.
Let ν be the standard Gaussian measure. As the Schoenberg's functions {ϕ j,d } j∈N associated to ψ in Equation (3.1) belong to the class P(R) ⊂ L 2 (R, ν), they can be expanded in terms of normalized Hermite polynomials.
The series converges in L 2 (R, ν), and each H k is a normalized Hermite polynomial of degree k given by Consequently, the kernel ψ in (3.1) can be reformulated as We call the series {α k,j,d } j,k∈N ⊂ C the spatio-temporal Hermite power spectrum.

Regularity Properties
This section is devoted to study the behaviour of the kernel covariance functions associated to an isotropic stationary GRF Z. It will be shown that the regularity of such kernels is closely related to the decay of the Hermite power spectrum or the angular power spectrum. Moreover, the latter characterizes also the (K −J)term truncation of a GRF Z as in Equation (3.5). We recall that we have introduced two expansions for the kernel covariance function of an isotropic stationary GRF Z: i.-The Spatio-temporal power spectrum, by using a double Kaurhunen-Loève expansion according to formula (3.6), valid for the kernel covariance function of any isotropic stationary GRF as in (3.5) over S d × R, σ d ⊗ µ with µ the Lebesgue measure.
ii.-The Hermite power spectrum, by using Hermite polynomials according to formula (3.7), valid for the kernel covariance function of any isotropic stationary GRF Z over S d × R, σ d ⊗ ν with ν the standard Gaussian measure.
Remark 4.1. Considering the relations between Gegenbauer and Legendre polynomials, in the case d = 2, the kernel (3.1) turns out to be In what follows we present the regularity analysis of the kernel in terms of the behaviour of the two proposed expansions (3.6) and (3.7). We address first the relation with the Hermite power spectrum (3.7).

Regularity analysis for the Hermite expansion
In this part of the manuscript we consider the measure space (R, ν) with ν the standard Gaussian measure.
For n, m ∈ N let the spaces H n (−1, 1) ⊂ L 2 (−1, 1) and H m (R) ⊂ L 2 (R) be the standard Sobolev spaces. We extend the proposal in [17] and consider the function spaces W n,m := W n,m ((−1, 1)×R) as the closure of H n (−1, 1)×H m (R) with respect to the weighted norm · W n,m ((−1,1)×R) given by Note that (W n,n , n ∈ N) is a decreasing scale of separable Hilbert spaces, i.e. By Theorem 5.2 in [1], the norm of W n,n is equivalent to the first and the last element of the sum, i.e.
We now derive another equivalent norm of W n,n in terms of summability of the spectrum. We first observe that, as the normalized Hermite polynomials {H k } k∈N constitute a orthonormal basis of L 2 (R), and that for any fixed r > −1/2 the Gengenbauer polynomials {C r j } j∈N a basis for L 2 (−1, 1), it is apparent that {H k ·C r j } j,k∈N is a basis for L 2 ((−1, 1)×R) ∼ = L 2 (−1, 1)⊗L 2 (R). Therefore, any f ∈ L 2 can be expanded in the series that is, f can be written as a covariance kernel ψ of the type (3.7).
This allows to tackle the problem in a different way. Instead of using spectral techniques, regularity of kernels might be shown through an isomorphism between the spaces W n,m ((−1, 1) × R) and the weighted square summable bisequence spaces (kj) m,n : where k m/2 · j (d−1)/2+n ; k, j ∈ N denotes the sequence of weights.
From now on, for the sake of simplicity, we only consider weighted Sobolev spaces W n := W n,n , and the weighted square summable bi-sequence spaces (kj) n := (kj) n,n , obtained as a special case of Equation (4.2) when n = m.
In order to extend the isomorphism to spaces W η with η being not an integer, following [26] we now introduce the interpolation spaces W η := W η ((−1, 1)×R) for n < η < n + 1, defined through: equipped with the norm f W η given by where the functional J is defined by The definition of the interpolation spaces (kj) η for η non-integer is carried out in analogous way.
The interpolation property (see section 2.4.1 in [25]), implies that, if the spaces W n and (kj) n are isomorphic for all n ∈ N, then they are isomorphic for all η ∈ R + . Theorem 4.1. Let f ∈ L 2 having expansion (4.1) and η ∈ R + be given. Then, is an equivalent norm in W η .
For ψ ∈ W n , n ∈ N, the equivalence is reduced to: α k,j k n/2 j (d−1)/2+n , k, j ∈ N is in 2 Proof. Assume first that the claim is already proved for η ∈ N, i.e., W n is isomorphic to the weighted bi-sequence space (kj) n for all n ∈ N.
Fix n ∈ N, let n < η < n + 1 and set κ := η − n. By the interpolation theorem of Stein-Weiss (see Theorem 5.4.1 in [4]), the weights of (kj) η are given by Now, we prove the isomorphism between W n and (kj) n for n ∈ N, which is equivalent to prove the second formulation of the theorem. We have that, where and Standard properties of normalized Hermite polynomials show that I k,k = 0 if k < n, and

On one hand, Stirling inequality
On the other hand, Hence, for k ≥ n, there exists constants c 1 (n) and c 2 (n) such that Therefore, I k,k k n δ k,k , k ≥ n. Therefore, from (4.3), (4.4) and (4.5) we deduce that which concludes the proof.

Regularity analysis for the double Karhunen-Loève expansion
With the exception of minor details on the definitions of the spaces, this part of the manuscript follows similarly to Section 4.1. On the other hand, we now consider the measure space (R, µ) with µ the Lebesgue measure.
For n ∈ N consider the function spaces V n T := V n T ((−1, 1) × R) as the closure of H n (−1, 1) × L 2 (R) with respect to the weighted norms · V n T given by (V n T ) n∈N is a decreasing sequence of separable Hilbert spaces, and Now we look at the weighted square summable bi-sequence spaces (kj) T,n := 2 j (d−1)/2+n ; k, j ∈ N .
As in Section 4.1, we consider the interpolation spaces V η T for n < η < n+1. The proof of next result follows exactly the same lines as the Theorem 4.1, hence it is omitted. Theorem 4.2. Let f ∈ L 2 and η ∈ R + be given. i.e.
Remark 4.2. By taking into account the normalizing constants, all the previous results encompasses the results in Section 3 of [17].

Spatio-Temporal Spectral Simulation
We now study a spectral simulation method for random fields on S 2 × [0, T ], where T denotes the time horizon. For a neater exposition, along this section, we omit the subscripts associated to the spatial dimension d = 2.

(5.3)
We consider Re(X j,m (t)) and Im(X j,m (t)) as independent processes with the following Fourier expansions Here, {A q j,k,m } and {B q j,k,m }, for q = 1, 2, are sequences of independent centred real-valued Gaussian random variables, such that var(A 1 j,k,0 ) = var(B 1 j,k,0 ) = a j,k and var(A q j,k,m ) = var(B q j,k,m ) = a j,k /2, for m = 0, where {a j,k } k,j∈N is a summable bi-sequence of non-negative coefficients. A direct calculation shows that the covariance function of Z(x, t) is spatially isotropic and temporally stationary. More precisely, we have that Finally, given two positive integers J and K, we truncate expression (5.3) in the index j and k, respectively. Thus, we simulate a space-time Gaussian random field on S 2 × [0, T ] using the explicit approximation where P j,m (·) = (2j + 1)(j − m)!/(4π(j + m))!P j,m (·).
We assess the L 2 (Ω × S 2 × [0, T ]) error associated to the truncated expansion Z(x, t) given in Equation (5.4), in terms of the positive integers J and K. We follow the scheme used by [17] in the spatial context and extend their result to the space-time case. Next, we state the main result of this section. Theorem 5.1. Let ν 1 , ν 2 ≥ 2. Suppose that there exist positive constants C i , for i = 1, 2, 3, and positive integers j 0 and k 0 , such that a j,0 ≤ C 1 j −ν1 , a 0,k ≤ C 2 k −ν2 and a j,k ≤ C 3 j −ν1 k −ν2 , for all j ≥ j 0 and k ≥ k 0 . Then, the following inequality holds

5)
for some positive constants C 1 and C 2 .
Proof. We decompose Z(x, t) − Z(x, t) in Equation (5.4) into two independent terms, where and with ∆ j,k (x, t) defined as For the second term, we have used the identity which is satisfied for any summable bi-sequence {ξ j,k } j,k∈N .
Simple examples can be generated from the following space-time angular power spectrum with ν i > 2, for i = 1, 2. We illustrate space-time realizations on S 2 × {1, 2}, over 24000 spatial locations, with coefficients (5.6), in two cases: (a) ν 1 = 3 and ν 2 = 5, and (b) ν 1 = ν 2 = 5. Figures 1 and 2 show the corresponding realizations for Scenarios (a) and (b), respectively. For each case, we truncate the series using K = J = 50. Note that the parameter ν 1 is the responsible of the spatial scale and smoothness of the realization. In [17], some realizations are illustrated using a similar spectrum, in a merely spatial context.
We now compare the empirical and theoretical convergence rates for the cases (a) and (b) described above. In our experiment, we consider J = K and study the (Log) error in terms of (Log) J, taking as the exact solution J = 50. Note that, under this choice, the bound (5.5) implies that the order of convergence is min{(ν 1 − 2)/2, (ν 2 − 2)/2}. Following [17], instead of calculating the L 2 -error, we take the maximum error over all the points on the space-time grid. The empirical errors are calculated on the basis of 100 independent samples. Our studies reflect the theoretical results (see Figure 3).

Conclusions and discussion
The present work has provided a deep look at the regularity properties of Gaussian fields evolving temporally over spheres. We hope this effort will put the basis for facing important challenges related to space-time processes. There are in fact many open problems related to mathematical modeling as well as to statistical inference and to optimal prediction. A list of open problems is included in the recent survey [20]. Amongst them, our paper is certainly related to Problem 1, that is to the construction of non-stationary processes on spheres cross time. Our work could also put the basis to solve Problem 2, related to the construction of multivariate space-time processes. It might be interesting to extend the study of regularity properties to the vector valued case. This would imply the use of a pretty different machinery. Problem 10 is closely related to our approach, because regularity properties are crucial to study Gaussian fields under infill asymptotics.
On the other hand, a question arise naturally: is it possible to make inference with a representation like (3.5) ? Or with its respective spectral decomposition?
The answer is, a priori, no. Establishing a clear relation between the parameters of any random field and its spectrum is not an obvious task, in fact, up today the only familiar stochastic process with known spectrum is Brownian motion, and his closer generalization, the fractional Brownian motion, doesn't have yet a known spectrum. However, under relatively weak hypotheses, the covariance kernel of a GRF turns out to be a Mercer kernel. This opens an alternative to the negative answer previously mentioned, by considering the eigenvalue problem associated to the integral operator induced by the kernel.

Karhunen-Loève Theorem
Recent results on functional analysis (see [11] and [12]) allow to construct Mercer's kernels in more general contexts. A proper interpretation of these results allows to generalize the classic Karhunen-Loève theorem in a very neat way. We first introduce the framework and basics notations from [11].
Let S be a nonempty set and K a positive definite kernel on S, i.e., a function K : S × S −→ C satisfying the inequality whenever n ≥ 1, {x 1 , x 2 , . . . , x n } is a subset of S and {c 1 , c 2 , . . . , c n } is a subset of C. The set of all positive definite kernels over S is denoted by PD(S).
If S is endowed with a measure υ, denote by L 2 PD(S, υ) the class of kernels such that the associated integral operator is positive, that is, when the following conditions holds i.e., K(f ), f L 2 (S,υ) ≥ 0, f ∈ L 2 (S, υ).
Finally we define what is a Mercer's kernel according to [12]: A continuous kernel K on S is a Mercer's kernel when it possesses a series representation of the form where {ζ j } j∈N is an L 2 (S, υ)-orthonormal basis of continuous functions on S, {a j (K)} j∈N decreases to 0 and the series converges uniformly and absolutely on compact subsets of S × S.
For the rest of this manuscript we will consider S to be a topological space endowed with a strictly positive measure υ, that is, a complete Borel measure on S for which two properties hold: every open nonempty subset of S has positive measure and every x ∈ S belongs to an open subset of S having finite measure. Besides, as in Section 2, (Ω, F, P) denotes a complete probability space. Then, the kernel K associated with the covariance function of X is a Mercer's Kernel. Therefore, X admits a Karhunen-Loève expansion where {ζ j } j∈N is an orthonormal basis of L 2 (S, υ), with E[λ j ] = 0, and there exists a sequence {a j } j∈N of non-negative real numbers such that E[λ j λ k ] = δ jk a j and Var[λ j ] = a j .
The series expansion (7.2) converges in L 2 (Ω × S; C), i.e., The series expansion (7.2) converges in L 2 (Ω; C) for all s ∈ S, i.e., for all s ∈ S The convergence of the series expansion (7.2) is absolute and uniform on compact subsets of S in the mean-square sense.
Proof. Let K be the kernel associated to the covariance of the stochastic process X, i.e.
and let K be it's associated integral operator. From hypothesis (7.1) it is direct to see that the mapping κ, such that s ∈ S → C κ(s) := K(s, s), belongs to L 1 (S, υ). Since K is positive definite in the usual sense, K(s, t) = K(t, s) and the matrix K(s, s) K(s, t) K(s, t) K(t, t) is positive definite, hence its determinant K(s, s)K(t, t) − K(s, t)K(s, t) is nonnegative. Thus Now, consider f, g ∈ L 2 (S, υ) and the classic tensor product of functions, i.e., with inner product given by .
In conclusion, the kernel K is continuous and L 2 (S, υ)-positive definite on S, and the mapping κ belongs to L 1 (S, υ). Then, by theorem 3.1 in [12] K is a Mercer's kernel.
The rest of the proof concerns the Karhunen-Loève expansion of the process X and it follows well-known arguments that we reproduce for the convenience of the reader. We have that K has a L 2 (S × S, υ ⊗ υ)-convergent series representation in the form K(t, s) = +∞ j=1 a j (K)ζ j (t)ζ j (s), t, s ∈ S, {a j (K)} j∈N decreases to 0, {ζ j } j∈N is a L 2 (S, υ)-orthonormal basis. The convergence of the series is absolute and uniform on compact subsets of S × S.
Note that E(λ 2 j ) = K(ζ j ), ζ j 2 , hence Cauchy-Scharwz inequality guarantees that λ j ∈ L 2 (Ω) for all j. Also, Fubini's theorem allows to see that E[λ j ] = 0, E[λ j λ k ] = δ jk a j and Var[λ j ] = a j . Now, for any fixed ω ∈ Ω 0 it is clear that By condition (7.1) is clear that E(|f 0 |) < ∞, hence, the dominated convergence theorem allows us to conclude that Now, fix s ∈ S. Again by Fubini's theorem we observe that,  Thus, the proof is concluded.
Remark 7.1. Karhunen-Loève expansion (or Karhunen-Loève theorem), usually require extra hypothesis, like compactness of the associated space S or some kind of invariance of the field. In that line the Stochastic Peter-Weyl theorem (theorem 5.5 introduced in [18]) may be understood as a Karhunen-Loève theorem for 2-weakly isotropic fields over G a topological compact group with associated Haar measure of unit mass. Theorem 7.1 only require condition (7.1) and the continuity of the covariance function.