Spherically Restricted Random Hyperbolic Diffusion

This paper investigates solutions of hyperbolic diffusion equations in R3 with random initial conditions. The solutions are given as spatial-temporal random fields. Their restrictions to the unit sphere S2 are studied. All assumptions are formulated in terms of the angular power spectrum or the spectral measure of the random initial conditions. Approximations to the exact solutions are given. Upper bounds for the mean-square convergence rates of the approximation fields are obtained. The smoothness properties of the exact solution and its approximation are also investigated. It is demonstrated that the Hölder-type continuity of the solution depends on the decay of the angular power spectrum. Conditions on the spectral measure of initial conditions that guarantee short- or long-range dependence of the solutions are given. Numerical studies are presented to verify the theoretical findings.


Introduction
Numerous environmental, biological and astrophysical applications require the modelling of changes in data on the unit sphere S 2 or in the 3D space R 3 [1][2][3][4][5][6].One of the conventional tools for such modelling is stochastic partial differential equations (SPDEs)-see, for example [1,5,7,8], and the references therein.Random fields that are solutions of such SPDEs often exhibit dynamics dependent on initial conditions.Investigating the properties of these random fields is important for theoretical insight and practical applications.
SPDEs on surfaces and Riemannian manifolds have found numerous applications to problems in cosmology, physics, biology, fluid dynamics and pattern formation on surfaces, just to mention a few.See [7,[9][10][11] and the references therein.
Random fields on a sphere, one of the simplest two-dimensional manifolds, have been used as a standard model in the astrophysical and cosmological literature in the last several decades [3,4,8,12].NASA and ESA space missions [4] obtained very detailed measurements of Cosmic Microwave Background radiation (CMB), which are interpreted as a realisation of a spherical random field superimposed on an underlying signal of large-scale acoustic waves in plasma near the time of recombination.The theory of the standard inflation scenario uses a Gaussian model for the density fluctuation of this field [3,4,6].Several new cosmological models were proposed using non-Gaussian assumptions and employed sophisticated statistical tests to justify possible departures from Gaussianity.
Entropy 2020, 22, 217 2 of 31 The understanding of changes in CMB temperature fluctuations is important to predict future cosmological evolution and accurately reconstruct past states of the Universe.It also can help in the estimation and statistical inference of physical parameters obtained from the CMB data.SPDEs on the sphere can be used to model changes in CMB temperature fluctuations; see [7,8].
The pronounced spectral peaks at very large wavelengths in CMB temperature data are evidence of acoustic waves that were seeded by earlier superluminal inflation, with remnant coherent waves remaining until last scattering of photons and recombination of atoms around 340,000 years after the big bang [6,13].In the plasma Universe there was chaotic mixing, but it is problematic to represent the underlying particle kinematics as standard Brownian motion, which leads to the standard diffusion equation for particle densities.Under standard diffusion, density disturbances have unbounded propagation speeds, which is unacceptable in relativistic cosmological contexts wherein remnant structures that are coherent over space-like domains have not been smeared away by diffusion.Reimberg [14] has directly modelled a sequence of photon-electron collisions backwards in time from the last scattering, with random changes of direction, and with the same distance travelled over equal time steps-unlike in usual random flight theory.Giona has developed a Feynman-Kac stochastic dynamics by which a particle undergoes a succession of collisions with speed-limited jumps.Consequently, the diffusion coefficients decrease as a power of the Lorentz-Fitzgerald contraction factor of an inertial observer.This has the interesting consequence that different observers may disagree on whether a process is deterministic or stochastic.Non-trivial continuous relativistic Markov processes on position space are simply not possible [15].A simpler alternative phenomenological model is effected by replacing the standard diffusion equation by the simplest hyperbolic diffusion equation that has a variable but bounded speed of propagation.Ali and Zhang [16] recast the hyperbolic diffusion equation as a Lorentz-invariant Liouville conservation equation in one time and four space dimensions, before restricting x 4 to be ict.Ali and Zhang then retain the second law of thermodynamics but only as a reaction-diffusion equation in 1+4 dimensions.Section 8 shows how information entropy may decrease by a small amount during the propagation of a point source by hyperbolic diffusion, whereas the overall increase is much larger.
The so-called Cattaneo hyperbolic diffusion equation [17,18] has been used to explain outcomes of heat conduction experiments in liquid He 4 in the super-fluid state [19,20], and solid He 3 and solid He 4 at very low temperatures.In these materials, and in some nanotubes and other graphite structures [21], heat energy propagates as a "second sound" wave mediated by phonons, at a propagation speed of around one-tenth the normal speed of sound.Since the experiments on graphite have been conducted on nano-scale samples, we expect that waves of second sound could likewise be detected in the spherical surface of a C 60 or larger fullerene ball, as formulated in our previous paper on hyperbolic diffusion on a spherical surface [8].
SPDEs on R 3 have been extensively studied.However, SPDEs on manifolds have only recently attracted a great deal of attention [7,8,22,23].The results in these papers demonstrate that the continuity properties of solutions and the convergence rates of approximations to solutions are determined by decay rates of the angular power spectrum of initial random conditions.This article continues studies of solutions of SPDEs on the sphere.However, in contrast to the above publications that directly model spherical random fields using Laplace or Laplace-Beltrami operators on the sphere, we employ another approach.Namely, we consider the restriction of the stochastic hyperbolic diffusion in R 3 to the unit sphere.Compared to the available literature this approach is more consistent with real CMB observations that exist in 3D space but are measured only on S 2 .From a mathematical point of view, additional investigations are required to show that solutions of known models on the sphere admit physically meaningful extensions to R 3 that are consistent with 3D observations.By its construction, our model directly provides this consistency.The proposed model may find new applications for the next generation of CMB experiments, CMB-S 4 , which will be collecting 3D observations.A very detailed discussion of SPDEs on manifolds and their physical and mathematical justification for CMB problems can be found in [8].The hyperbolic diffusion equation prohibits the superluminal Entropy 2020, 22, 217 3 of 31 propagation of density disturbances, which is an unwanted feature of pure diffusion models over super-galactic distances.In addition, the linear hyperbolic diffusion equation, expressed in terms of co-moving material space coordinates and conformal time coordinate, is a good approximation to the field equation of a scalar field minimally coupled to an expanding Robertson-Walker space-time.However, speed-limited diffusion raises some interesting questions about the dynamics of Shannon entropy.For physical concentrations governed by linear or nonlinear heat diffusion equations of parabolic type, Shannon entropy is fully analogous to thermodynamic entropy and it increases monotonically [24].It will be explained that at low wave numbers, the hyperbolic diffusion equation behaves as a dissipative diffusion equation.However, above some cut-off wave number it behaves as a bi-directional wave equation which has increasing entropy when twin pulses separate but has decreasing entropy when pulses approach each other constructively.
This paper is organised as follows.Section 2 presents definitions and results about spatial-temporal random fields in R 3 .It also introduces hyperbolic diffusion equations with random initial conditions and their solutions.Section 3 investigates the spatial-temporal hyperbolic diffusion field on the unit sphere.The Hölder-type continuity of the exact solution of the spatial-temporal hyperbolic diffusion field on the sphere is investigated in Section 4. In Section 5 we study the dependence structures of the spherical hyperbolic diffusion random fields.Section 6 obtains the mean-square convergence rate to the diffusion field in terms of the angular power spectrum.Section 7 provides some numerical results.Finally, Shannon entropy behaviour is discussed in Section 8, followed by some conclusions.
In the following sections we will use the symbol C to denote constants that are not important for our exposition.The same symbol may be used for different constants appearing in the same proof.

Spatial Random Hyperbolic Diffusion
This section reviews the basic theory of random fields in R 3 and introduces a hyperbolic diffusion with random initial conditions.Then, the solution of the diffusion equation is derived and analysed.
We consider the hyperbolic diffusion equation where q(x, t) = q(x, t, ω), ω ∈ Ω, is a random field satisfying the random initial conditions: where ∆ is the Laplacian in R 3 and the random field η(x) = η(x, ω), x ∈ R 3 , ω ∈ Ω, defined on a suitable complete probability space (Ω, F , P), is assumed to be measurable, mean-square continuous, wide-sense homogeneous and isotropic with zero mean and the covariance function The covariance function has the following representation: for some bounded, non-negative measures Entropy 2020, 22, 217 4 of 31 See [25], pp.1-5 and [26], pp.10-15 for more details.Then, there exists a complex-valued orthogonally scattered random measure Z(•) such that for every x ∈ R 3 , the field η(x) itself has the spectral representation Let Y lm (θ, ϕ), θ ∈ [0, π], ϕ ∈ [0, 2π), l = 0, 1, . . ., m = −l, . . ., l, be complex spherical harmonics defined by the relation where P m l (•) are the associated Legendre polynomials with indices l and m.For spherical harmonics it holds that where the symbol * denotes the complex conjugation and δ l ′ l is the Kronecker delta function.The addition formula for spherical harmonics gives The Bessel function J ν (•) of the first kind of order ν is defined by where Γ(•) is the Gamma function.
It admits the following representation by the Poisson integral, see (10.9.4) in [27]: By the addition theorem for Bessel functions(e.g., [26], p. 14), where Z lm (•) is a family of random measures on (R 1 + , B(R 1 + )), such that The stochastic integrals in (3) and ( 4) are viewed as L 2 (Ω) integrals with the structural measures F and G, correspondingly.
Let us consider the initial conditions of the form: where δ(x) is the Dirac delta function.
Let Q(x, t), x ∈ R 3 , t ≥ 0, be the fundamental solution (or the Green's function) of the initial-value problem (1) and (6), and let be its Fourier transform.
The following theorem derives the Fourier transform H(κ, t).Contrary to many models in the literature, for the initial-value problem (1) and (2) it can be explicitly written in terms of elementary functions.Later, this result will be used to obtain the solution q(x, t, ω), x ∈ R 3 , t ≥ 0, ω ∈ Ω, and its covariance function.
Theorem 1.The Fourier transform (7) of the initial-value problem (1) and (6) is given by the formula where I {•} denotes the indicator function.
Proof of Theorem 1.The Fourier transform (7) is the solution of the initial-value problem The characteristic equation for the ordinary differential equation in (11) is with the roots Therefore, the general solution of the ordinary differential equation in (11) has the form where K 1 (κ), K 2 (κ) are some functions that do not depend on t and z 1 (κ), z 2 (κ) are given by (12).
The following theorem provides the solution of the initial-value problem and its covariance function in terms of the Fourier transform H(κ, t).As the explicit expression of H(κ, t) in terms of elementary functions is given in Theorem 1, it can be used to obtain an explicit expression for the solution and then easily investigate various properties of q(x, t).Theorem 2. The solution q(x, t) = q(x, t, ω), x ∈ R 3 , t ≥ 0, ω ∈ Ω, of the initial-value problem (1) and (2) can be written as the convolution The covariance function of the spatio-temporal random field (20) is Proof of Theorem 2. Notice that where H(κ, t) is given by ( 8), assuming that the random initial condition has the spectral measure F, such that Under the condition (22), the stochastic integral (20) exists in the L 2 (Ω)-sense.
By Lemma 1 the function |H(κ, t)| can be bounded by a constant C(t) which depends only on t. 0) we obtain (22).The representation (21) immediately follows from (20) and the orthogonality of Z(•).

Spherical Random Hyperbolic Diffusion
In this section we investigate a restriction of the spatial-temporal hyperbolic diffusion field from Section 2 to the unit sphere.
Consider the sphere S 2 = {x ∈ R 3 : x = 1} in the three-dimensional Euclidean space R 3 with the Lebesgue measure A spatio-temporal spherical random field defined on a probability space (Ω, F , P) is a stochastic function We consider a real-valued spatio-temporal spherical random field T with zero mean and finite second-order moments and being continuous in the mean-square sense (see, e.g., Marinucci and Peccati [3] for definitions and other details).Under these conditions, the zero-mean random field T can be expanded in the mean-square sense as the Laplace series [25]: where the functions Y lm (θ, ϕ) represent the spherical harmonics and the coefficients a lm (t) are given by the formula We assume that the field T is isotropic (in the weak sense), that is, ET 2 (x, t) < ∞, and ET(x, t)T(y, t ′ ) = ET(gx, t)T(gy, t ′ ) for every g ∈ SO(3), the group of rotations in R 3 .This is equivalent to the condition that the covariance function E T(θ, ϕ, t) T(θ ′ , ϕ ′ , t ′ ) depends only on the angular distance γ = γ PQ between the points P = (θ, ϕ) and Q = (θ ′ , ϕ ′ ) on S 2 for every t, t ′ ≥ 0. The field is isotropic if and only if Hence, The functional series {C l (t, t ′ ), l = 0, 1, . . .} is called the angular time-dependent power spectrum of the isotropic random field T(θ, ϕ, t).
We can define a covariance function between two locations with the angular distance γ at times t and t ′ by where Entropy 2020, 22, 217 9 of 31 If T(θ, ϕ, t) is a zero-mean isotropic Gaussian field, then the coefficients a lm (t), m = −l, . . ., l, l ≥ 1, are complex-valued Gaussian stochastic processes with By Remark 1, the random field q(x, t), x ∈ R 3 given by ( 20) is homogeneous and isotropic in x, and hence its covariance function (21) can be written in the form: where (r, θ, ϕ) and (r ′ , θ ′ , ϕ ′ ) are spherical coordinates of x and x ′ respectively.Using the Karhunen theorem we obtain the following spectral representation of the random field: where the random measures Z lm (•) are given in (5).
Similarly to the condition ( 22) the isotropic measure G(•) satisfies the following condition if the field has a finite variance: Subclasses of covariance functions of the isotropic fields on the sphere can be obtained from covariance functions of homogeneous isotropic random fields in Euclidean space, since a restriction of the homogeneous and isotropic random field to the sphere yields an isotropic spherical field (e.g., [25], p. 76).
Consider two locations x and x ′ on the unit sphere S 2 with the angle γ ∈ [0, π] between them.Then, the Euclidean distance between these two points is 2 sin γ 2 , the inner product is x, x ′ = cos γ, which gives a direct correspondence between the covariance function R 0 ( x − x ′ , t, t ′ ) in the Euclidean space and the covariance function R(cos γ, t, t ′ ) = R 0 (2 sin γ 2 , t, t ′ ) on the sphere for every fixed t, t ′ ≥ 0. Thus, the restriction of the homogeneous and isotropic hyperbolic diffusion field (25) to S 2 is an isotropic spherical random field for every fixed t, t ′ ≥ 0. We will call it the spherical hyperbolic diffusion isotropic random field Its covariance function is of the form: By the addition theorem for Bessel functions, the random field T H (x, t) = TH (θ, ϕ, t) has the following spectral representation: where and the random measure Z lm (•) satisfies (5).
Entropy 2020, 22, 217 10 of 31 Thus, the angular spectrum of the isotropic spherical random field T H (x, t) is given by the formula Therefore, we obtained the following result.
Then, the restriction of the spatio-temporal hyperbolic-diffusion random field (25) to the sphere S 2 is an isotropic spatio-temporal spherical random field with the following angular spectrum: The field and its covariance functions are given by ( 27) and (26), respectively.
This result investigates the restriction of the spatio-temporal hyperbolic-diffusion random field to the sphere S 2 .It shows how the angular power spectrum C l (t, t ′ ) of the restriction depends on the Fourier transform H(µ, t).Hence, one can explicitly compute coefficients C l (t, t ′ ) and study the contributions of different spherical harmonics to the spatial-temporal field T H (x, t).
Remark 3. It follows from Lemma 2 and the estimate |P l (cos θ)| ≤ 1 that the solution's covariance function given by ( 24) is finite if the initial condition η(x), x ∈ S 2 , has a finite variance.

Smoothness of Solutions
Numerous problems in mathematical physics and geosciences require studying the regularity properties of solutions of differential equations.Smoothness, boundedness of derivatives or Hölder continuity conditions are often used to describe and investigate local changes and growth rates of solutions.Knowing regularity properties is also essential for an adequate approximation of SPDE solutions.In those cases where solutions are given by infinite series, it is a rather difficult mathematical problem, as the tail terms of such series can accumulate.
In this section, we investigate the Hölder-type continuity of the solution T(θ, ϕ, t) given by ( 27) on the sphere.Estimations of the closeness of T values at spherical points (θ, ϕ) and (θ ′ , ϕ ′ ) are obtained.We demonstrate how they depend on the decay of the angular power spectrum and provide some specifications in terms of the spectral measure G(•).
First, we obtain the continuity of the solution with respect to the geodesic distance on the sphere.To prove this we use the approach from Corollary 5 in [8].
Theorem 4. Let TH (θ, ϕ, t) be the solution of the initial value problem (1) and (2) and the random initial condition η(x), x ∈ S 2 , has the angular power spectrum {C l , l = 0, 1, 2, . . .} satisfying the assumption where γ is the angle between directions (θ, ϕ) and (θ Proof of Theorem 4. (a) It follows from ( 23), ( 24), ( 27) and ( 32) that Applying the property of Legendre polynomials, [23], p. 16, one obtains the statement (a) of the theorem.(b) It follows from the proof of (32) that in the case of G([0, c 2D ]) = 0 it holds The remaining steps are similar to the proof in (a).
The next two results provide conditions on the field's spectrum that guarantee the Hölder-type regularity of TH (θ, ϕ, t).
even for the case of α = 0 in (33).
The next result gives sufficient conditions to guarantee (33).

Short and Long Memory
Investigating statistical dependence between measurements at two points with increasing time or spatial distance between them is an important issue for practical temporal or spatial predictions.The spatial domain of the considered random fields is restricted to the sphere S 2 with the geodesic distance γ.Note that this distance is bounded to the interval [0, π], but time t can unboundedly increase and takes values in [0, +∞).Hence, this section investigates only temporal statistical dependencies, namely, slow or fast decays of covariance functions in time.The corresponding cases represent long or short memory scenarios.
In this section we use the representation (26) of covariance functions to investigate the structure of dependences of T H (x, t) over time.We demonstrate that conditional on the spectral isotropic measure G(•) of the initial random condition η(x), x ∈ R 3 , the random field T H (x, t) can exhibit short-or long-range dependence.
The random field T H (x, t) will be called short-range dependent if for all t ≥ 0 and γ ∈ [0, π].If the integral in (35) is divergent, the field is called long-range dependent.
Results that link behaviours of covariance functions at infinity and spectral measures at the origin are called Abelian-Tauberian theorems.A very detailed overview of such results for random fields can be found in [28].
First we investigate the case of x = x ′ in (26) (i.e., the behaviour of R(1, t + h, t)).
Theorem 7.For x = x ′ the random field T H (x, t) exhibits short-range dependence if and only if µ −2 G(dµ) is integrable in a neighbourhood of zero.
Proof of Theorem 7. It follows from ( 16), ( 17) and ( 26) that Using the upper bound from (19) we get or, by ( 16) and cosh Noting that sin(h) A on [0, A], A > 0, we obtain that (37) is finite if and only if the following integral converges: The last integral is finite only if Now, note that for γ ∈ (0, π) the interval [0, c/2D) can be split into a finite number of subintervals Note that by lim µ→0 sin(µ) µ = 1 this condition is equivalent to the one in Theorem 7.This completes the proof.

Approximations to Solutions
The results in the previous sections were based on the series representation of the random field TH (θ, ϕ, t).In applications and numerical studies, only a finite number of series terms are available.
Hence, one has to investigate the behaviours of finite cumulative sums.This section provides an analysis of truncated series expansions of the solution field TH (θ, ϕ, t) and shows the role of the decay rate of the angular power spectrum.These results can be used to determine the number of terms for a given accuracy of approximate solutions.
This section introduces and studies approximate solutions of the initial value problem (1) and (2).A mean-square convergence rate to the diffusion field in terms of the angular power spectrum C l is obtained.Then, several specifications in terms of the measure G(•) are discussed.
We define the approximation TH,L (θ, ϕ, t) of the truncation degree L ∈ N to the solution TH (θ, ϕ, t) given by (27) as The next result provides the convergence rate of TH,L (θ, ϕ, t) to TH (θ, ϕ, t) when L → ∞.Theorem 9. Let TH (θ, ϕ, t) be the solution to the initial value problem (1) and (2) and TH,L (θ, ϕ, t) the corresponding approximation of truncation degree L ∈ N.Then, Proof of Theorem 9. Note that by properties of a lm (t) we get Then, by ( 23) and ( 27) it follows that Using the addition formula for spherical harmonics one gets Finally, by ( 32) For the SPDE model studied in [8] it was shown that its solution has an exponential decay in t and the corresponding approximation error can be bounded as see (36) in [8].
The following result shows that the considered model is more complex.In the general case of an arbitrary measure G(•) it is impossible to get a bound similar to (39) even for a sufficiently large L. Theorem 10.For any fixed C > 0 and L ∈ N there exist t > 0 and an initial random condition η(x), x ∈ R 3 , such that the norm of the approximation error TH (θ, ϕ, t) − TH,L (θ, ϕ, t) does not satisfy (39).
Then, 1 Hence, by (38) and Theorem 3 for any C, L > 0, there exists t, ε > 0, and the measure G(•) such that for the corresponding TH (θ, ϕ, t) and TH,L (θ, ϕ, t) it holds However, it is possible to obtain a rate of convergence that is exponential in t if the measure G(•) has a bounded support.Theorem 11.Let η(x), x ∈ R 3 , have the measure G(•) such that G([0, δ]) = 0 for some δ ∈ (0, c 2D ).Then, for the solution TH (θ, ϕ, t) of the initial value problem (1)-(2) and its approximation TH,L (θ, ϕ, t) it holds that x is an increasing function on (0, ∞) it follows from (16) that for µ ≥ δ Notice that for x ≥ 0 and a ∈ (0, 1) it holds that 1 + x ≤ 1 a exp(xa).Then, using the definition of H2 (µ, t) in ( 17) we get for t ≥ 0 Applying this bound to (38) we obtain the statement of the theorem.
Remark 4. The rates of convergence in Theorems 9, 11 and Corollary 1 are sharp.Indeed, for t = 0 one obtains The angular power spectrum {C l , l = 0, 1, . . .} of the initial random field η(x) is determined by the measure G(•).The following results provide some insight into the behaviour of ∑ ∞ l=L (2l + 1)C l in terms of the spectral measure G(•).

(a)
Then, it holds that Proof of Theorem 12. (a) It follows from the representation By von Lommel's formula (see (2.60) in [29]), where µ ∈ R and ν > −1, we obtain Now, (40) follows by substituting the last expression in (43).(b) Using the inequality from [30] |J we obtain that for L ≥ 2 which after the substitution in ( 43) gives ( 41).
(c) By the Poisson integral formula and the identity one obtains If [0, δ], δ > 0, is the support of the measure G(•), then it follows from (43)-(45) that which completes the proof.

Numerical Studies
This section presents numerical studies of the solution T H (x, t), its angular spectrum, and covariance functions over time.We also provide some numerical analysis of approximation errors.
All numerical computations and simulations in this paper were performed using the software R version 3.6.1 and Python version 3.7.5.The results were derived using the HEALPix representation of spherical data (see [31] and http://healpix.sourceforge.net).In particular, the R package rcosmo [32,33] was used for computations and visualisations of the obtained results.The Python package heal py was used for fast spherical harmonics generation of spherical maps from Laplace series coefficients.The R and Python code used for numerical examples in Section 7 are freely available in the folder "Research materials" from the website https://sites.google.com/site/olenkoandriy/.
It is important to clarify that the numerical analysis in this paper is rather different from the one in [8] and requires more advanced approximation approaches.Namely, the stochastic model in [8] yielded the representation of the Laplace series coefficients a lm (t) = C[A l (t) + B l (t)]a lm (0) for some functions A l (t) and B l (t) which can be explicitly computed in terms of elementary functions.However, for the model ( 1)-(2) there is no such simple functional relation that links a lm (t) and a lm (0).As a result, there are no explicit elementary functional relations between C l (t, t ′ ), R(cos γ, t, t ′ ) and C l (0, 0), R(cos γ, 0, 0) respectively.To compute spectral and covariance functions of T H (x, t) at time t > 0 one has to use formulae ( 26), ( 28) and (29).These integral representations are given in terms of the spectral measure G(•) and stochastic measures Z lm (•) of the initial random condition field η(x).
By (1.2.5) in [26], it follows from which can be used to compute ( 26), ( 29) and simulate Z lm (•) for computations in (28).However, obtaining a reliable approximation of the integral in (46) and stochastic measures Z lm (•) requires the estimation of the empirical covariance function R(r) on a dense grid.Moreover, for data observed on bounded subsets of R 3 , covariance functions can be estimated only for distances that do not exceed their diameters.Thus, it is important to verify that empirical covariance functions are sufficiently quickly decaying to be assumed negligible for distances greater than these diameters.We postpone the solution of these technical problems and analysis of real data to future publications.
In the following examples, we study the properties of solutions and their approximations using simulated data.The examples were constructed to demonstrate that the model is sufficiently powerful to imitate behaviours of the empirical CMB covariance function and oscillating angular spectrum (see [8,34]).The actual CMB covariance function and angular spectrum are shown in Figure 1a,b.Note that the estimated angular CMB spectrum shown in Figure 1b was obtained by a piecewise fitting of several physical models and interpolation techniques for different intervals of the spectrum [6,13].Some actual spectrum estimates deviate substantially from the fitted curve in Figure 1b [34].Therefore, in predicting CMB and spectrum changes over time, small details can be ignored and one needs to focus on a general pattern.Thus, the following examples with the analysis of simulated data which spectrum is analogous to the real one can offer important insights on the future evolution of CMB and its spectral properties.In the examples, the case of a discrete measure G(•) is considered (i.e., the support of G(•) is a finite set {µ i , i = 1, . . ., I}).We employ real-valued stochastic measures Z lm (•) that are concentrated on this set and satisfy the condition We assume that the random field η(x) is centred Gaussian.Hence, we can choose Z lm (µ i ) ∼ N(0, σ 2 i ) that are independent for different l, m and i.
In these settings, formulae ( 26), ( 28) and ( 29) take the following discrete forms: which are convenient for simulations.This approach can also be used to approximate absolutely continuous spectral measures G(•) by considering a sufficiently large I, Example 1.This example illustrates changes over time of the covariance function R(cos γ, 0, t) defined by (26) and the power spectrum C l (t, t) defined by (29).To produce plots and computations we used the corresponding discrete Equations ( 47) and ( 48) with values σ i = 100 i by i ∈ {1, 2, . . ., 10} and a discrete spectrum concentrated on the interval [1,40].All computations and plots in this example are presented for the values c = 1 and D = 1 of the parameters in Equation (1).
Figure 2a shows the covariance R(cos γ, t, t) at the time lags t = 0, t = 0.1 and t = 0.5 as functions of the angular distance γ.To understand the effect of time and the angular distance γ on the covariance function we provided 3D-plots (see Figure 2b) showing the covariance as a function of the time lag t.The plot in Figure 2b is normalized by dividing each value by max γ∈[0,π] R(cos γ, 0, 0).It is obvious that the covariance decays through time and changes very little except for values of γ, which are close to 0.
To understand the effect of the parameters c and D on the covariance function, we also produced Figure 3.It illustrates changes of the covariance function R(cos γ, t, t) at a specific time t as functions of the angular distances γ and the parameters c or D. To produce this figure we used t = 0.1.with the above values of σ i , i = 1, . . ., 10. From this figure it is clear that the power spectrum C l (t, t) decays very quickly to 0 when l increases.To investigate the effect of the parameter l we provide a plot of the ratio R 0.1,0,l = C l (0.1, 0.1)/C l (0, 0) for the first 70 coefficients C l (see Figure 4b).This figure confirms that the ratio R 0.1,0,l is bounded by 1 and changes very little when l increases.
Figure 5a plots the tail sums ∑ l≥L (2l + 1)C l (0, 0) and ∑ l≥L (2l + 1)C l (0.1, 0.1) as functions of L, while Figure 5b displays the corresponding ratio RR 0.1,0,L = ∑ l≥L (2l+1)C l (0.1,0.1) ∑ l≥L (2l+1)C l (0,0) .From Figure 5a it is clear that  Example 2. In this example we use a discrete spectrum concentrated on the two intervals [0, 20] and [80,90].Thus, the initial condition random field η(x) has low-and high-frequency components.To produce realisations of η(x) and T H (x, t), x ∈ S 2 , which are similar to small real CMB values, we used σ 2 i = 0.00003 and 0.0001 for low-and high-frequency components, respectively.These small values allow us to employ the visualisation tools and colour palettes used for CMB plotting in the R package rcosmo [33] and the Python package heal py.
To produce the plots and computations in this paper we used the first 100 coefficients C l obtained by applying (48) to the above discrete spectrum.They are shown in Figure 6 in red.In this example we use the values c = 1 and D = 2 of the parameters in Equation (1).The coefficients C l (t, t) for t = 0.05 and 0.1 are plotted in blue and green respectively.The graph indicates two regions with relatively large values of C l that correspond to the spectral measure G(•) used for these computations.It can be seen that values C l (t, t) decrease over time.However, the corresponding spherical maps change rather slowly.Therefore, only two maps, for t = 0 and 0.05, are plotted in Figure 7.For the following numerical studies we used simulated data from two windows shown in Figure 7b.The estimated means in Table 1 confirm that T H (θ, x, t) has a zero mean.It can be observed from Figure 7 and the estimated interquartile ranges (IQRs) in Table 1 that the magnitude of T H (x, t) values decrease with time.However, the distribution type of the combined values does not change substantially.Namely, the combined values of T H (x, t) exhibit an approximately bell-shaped behaviour with tails that are heavier than in the Gaussian case (see Figures 8 and 9).Similar results were obtained for various observation windows of S 2 .For example, for the second rectangular window shown in Figure 7b Q-Q plots and histograms of observations in this window are given in Figures 8 and 9, respectively.These results about distributions of combined values were also confirmed by computing the Shannon entropy for the empirical distributions { pi } given by the histograms in Figure 9. Values of Ĥ do not change much over time t (Table 1).They are not substantially different from the entropy upper bound log(16) ≈ 2.77.The q-statistic, see [35], was used to investigate heterogeneity between values of T H (θ, ϕ, t) in windows 1 and 2 from Figure 7b.Table 1 indicates that heterogeneity is absent at time 0 and the evolution due to the model (1) does not introduce heterogeneity, at least for short time periods.

Entropy and Hyperbolic Diffusion
This section discusses the evolution of Shannon entropy for hyperbolic diffusion.Theoretical analysis and several numerical examples are presented.To simplify the exposition and plots, only the case of x ∈ R and various non-random initial conditions are studied.
For diffusive transport that arises from the random motion of particles, the mass distribution may indeed be regarded as a probability distribution, after which the Shannon entropy may be calculated.For a simple thermodynamic system governed purely by linear or nonlinear heat conduction, there is a close analogy between thermodynamic entropy and Shannon entropy (e.g., [24,36]).When the transport mechanism is modified to hyperbolic diffusion, the behaviour of entropy requires more scrutiny.In order to illustrate this, consider one-dimensional solutions q(x, t) on [−ℓ, ℓ] × R + , subject to Neumann boundary conditions This may represent transport in the x-direction through a linear conduit of cross-sectional area A, with the variation of density in each cross section being effectively zero.It will be seen that the total mass M is constant.Therefore, the scaled density q * = qA/M has constant unit integral on [−L, L], from which physically relevant non-negative solutions q * (x, t) may be regarded as distributions.By choosing length scale D/c and time scale D/c 2 , it may be assumed that the coefficients in the hyperbolic diffusion equation are normalised to ±1.
Let t * = tc 2 /D, x * = xc/D and L * = Lc/D.Then, subject to boundary conditions q * x * = 0, x * = ±L * and initial conditions q * (x * , 0) = u 0 (x * ), q * t * (x * , 0) = v 0 (x * ).Defining Shannon entropy density to be s = −q * log q * , the hyperbolic diffusion equation for q * (x, t) implies The case of unbounded speed of propagation is obtained by taking the limit c → ∞, which results in a positive entropy production rate Dq * 2 x /q * .This is familiar from the theory of heat conduction, for which the entropy production rate is L e DT 2 x /T, where T is absolute temperature and L e is the Lewis number, which is the order-1 ratio of thermal diffusivity to mass diffusivity.
For uni-directional waves of velocity ±c, the entropy production rate is zero.For bi-directional waves, the total Shannon entropy is constant when opposite-travelling waves are not superposing, increasing when opposite-travelling superposing waves are separating, and decreasing when they are superposing and approaching.However, non-constant travelling wave solutions of the hyperbolic diffusion equation must have speed less than c and they must have an amplitude that decreases with time.For the remainder of this section, the asterisk superscripts will be conveniently omitted.Some solutions of the hyperbolic diffusion equation may be of dissipative diffusive type, while others may be dissipative bi-directional waves.In order to illustrate this, by the completeness of the Fourier transform, the general even solution by the separation of variables is a n e −0.5t cos(ω n t) cos(k n x), where The first summation covers modes that are purely dissipative in character, just as for the linear heat diffusion equation.However, in this case, the dissipative modes exist only when L ≥ 2π.The second summation covers standing wave modes with decaying amplitude.These may be regarded as a superposition of a decaying left-travelling wave and a decaying right-travelling wave.Note that the dissipative mode with logarithmic decay rate α − 1 decays more slowly than all other modes.The above solution is mass-conserving with mean value a 0 and constant mass integral 2La 0 = 1 by normalisation.For a single decaying standing wave mode of a hyperbolic diffusion equation distribution, for some value of t, Then, the total Shannon entropy is At times t = (2m + 1)π/2ω n , m ∈ Z, the distribution is uniform, which is the state of maximum entropy S = log(2L).Overall, the total entropy oscillates as it approaches the limiting equilibrium state.However the negative excursions of entropy may be quite small since the amplitude of oscillation decreases exponentially.Figure 10 plots the total entropy for a wave with single harmonic, calculated by trapezoidal integration with 400 intervals, versus time.It would be helpful to have a point-source solution for the hyperbolic diffusion equation.As far as we are aware, there is no known simple expression for the point source evolution but it has the standard uniform Fourier spectrum that evolves according to (50).It is plotted in Figure 11 after truncating the Fourier series at 100 terms.As in the d'Alembert wave equation, two separating travelling delta waves emerge but now the amplitudes of the truncated spikes are decreasing and there is an additional central symmetric hump due to the purely diffusive terms.The leading edges of the spikes are travelling at maximum speed c.In two and three dimensions there would be similar solutions with a single travelling cylindrical or spherical shock wave surrounding a central hump.It is instructive also to view the motion of an initial rectangular disturbance of finite amplitude.This is approximated in Figures 12 and 13 by a Fourier series of 200 terms.The truncated Fourier series is an exact solution, but due to the truncation and the boundary conditions, the solution is negative at some values of the domain, so that Shannon entropy cannot be calculated.However, the solution is indicative of the behaviour of a non-negative solution with initial rectangle.As in the bidirectional wave equation, the symmetric solution consists of two superposed rectangles that increase entropy as they begin to separate by travelling in opposite directions.After they have separated, their amplitude decreases, which leads to further entropy increase.The height of the leading edge decreases more rapidly than the trailing edge, so each rectangle evolves to a trapezoid.The leading edge-which is the boundary of the disturbance-continues to move at maximum speed c. Between the trapezoids, there is a central hump that eventually dominates, and resembles a diffusive Gaussian, increasing entropy further.With this kind of peaked initial condition, there is no indication of any significant period of entropy decrease.

Future Research Problems
This paper investigated evolutions of random fields determined by hyperbolic diffusion equations with random initial conditions.Spherical random fields were modelled as restrictions of 3D solution fields to the sphere.Compared to the previous publications, it resulted in a more realistic physical model.However, the solution field for the new model cannot be represented by using Laplace series coefficients a lm (0) of the initial condition directly.The more complicated representation involves spectral measures of initial random conditions.
Detailed studies of the solutions and their approximations were presented.In particular, regularity properties and temporal dependencies of solutions were investigated.Approximations to the SPDE solutions were proposed, and the upper bound analysis of approximation errors was provided.It was demonstrated that the magnitude of approximation errors is determined by the angular power spectrum C l and decreases at the rate of the cumulative tail sums (∑ ∞ l=L (2l + 1)C l ) 1/2 .The numerical studies investigated the dependence of solution fields on parameters of the SPDE model and provided some insight into the evolution of Shannon entropy for hyperbolic diffusion.Some important problems and extensions for future research are: • Investigating the sharpness of the obtained upper bounds on approximation errors (see [8]); • Developing statistical estimators of the equation parameters and studying their asymptotic properties; • Extending the methodology to tangent spherical vector fields (see [37]); • Developing numerical methods for the obtained representations to deal with spectra of initial conditions; • Extending the analysis and numerical studies in Section 8 to other scenarios; • In line with the theme of this Special Issue, in future we intend to study the effect of nonlinear diffusivity in the equation q t + 1 c 2 q tt = ∇ • [D(q)∇q].For example, if q is the electron density in a plasma, D(q) is typically decreasing [38].
Figure 3a displays R(cos γ, 0.1, 0.1) for D = 1 as a function of c and the angular distances γ.
Figure 3b displays R(cos γ, 0.1, 0.1) for c = 1 as a function of D and the angular distances γ.The plots in Figure 3 are normalised by dividing each value by max γ∈[0,π]R(cos γ, 0.1, 0.1).It is clear form Figure3athat the covariance decays through c (also through D; see Figure3b) and changes very little except values of γ which are close to 0. Figure3bdemonstrates that the normalised covariance function exhibits decaying periodic behaviour when D increases.

Figure 4a displays the
Figure4adisplays the power spectrum C l (t, t) as a function of t ≥ 0. To produce this figure we used t ∈ [0, 1] and l = 2, 5 and 10.The first 70 coefficients C l were computed by applying the Equation (48) with the above values of σ i , i = 1, . . ., 10. From this figure it is clear that the power spectrum C l (t, t) decays very quickly to 0 when l increases.To investigate the effect of the parameter l we provide a plot of the ratio R 0.1,0,l = C l (0.1, 0.1)/C l (0, 0) for the first 70 coefficients C l (see Figure4b).This figure confirms that the ratio R 0.1,0,l is bounded by 1 and changes very little when l increases.Figure5aplots the tail sums ∑ l≥L (2l + 1)C l (0, 0) and ∑ l≥L (2l + 1)C l (0.1, 0.1) as functions of L, while
By the Poisson integral representation of the Bessel function it follows that Hence, to study the integrability of the covariance function |R(1, t + h, t)| one has to investigate the integral which completes the proof.Now we extend Theorem 7 to the case of arbitrary x and x ′ from S 2 .
Therefore, similar to the proof of Theorem 7 we obtain the sufficient and necessary condition for the integrability of |R(cos γ, t ′ , t)|

Table 1 .
Statistics for windows 1 and 2.
Total entropy for standing wave with single harmonic.L = 3π, wave number k 2 = 2π/L.